text:
export dvv_lstsq

"""
    dvv_lstsq(time_axis, dt, w)

Linear regression of phase shifts to time to give dt/t=-dv/v

# Arguments
- `time_axis::AbstractArray`: Array of lag times
- `dt::AbstractArray`: Time shifts for each lag time
- `w::AbstractArray`: Weights for the linear regression

# Returns
- `m::Float64`: dt/t for current correlation
- `em::Float64`: Error for calculation of `m`
- `a::Float64`: Intercept for regression calculation
- `ea::Float64`: Error on intercept
- `m0::Float64`: dt/t for current correlation with no intercept
- `em0::Float64`: Error for calculation of `m0`
"""
function dvv_lstsq(time_axis::AbstractArray, dt::AbstractArray; w::AbstractArray=ones(length(time_axis)))
    # regress data using least squares
    # force line through origin
    model0 = glm(@formula(Y ~0 + X),DataFrame(X=time_axis,Y=dt),Normal(),
                IdentityLink(),wts=w)
    # allow nonzero intercept
    model = glm(@formula(Y ~ X),DataFrame(X=time_axis,Y=dt),Normal(),
                IdentityLink(),wts=w)

    # parameters for the regression with nonzero intercept
    a,m = -coef(model).*100
    ea, em = stderror(model).*100
    # parameters for the regression through the origin
    m0 = -coef(model0)[1]*100
    em0 = stderror(model0)[1]*100

    return m, em, a, ea, m0, em0
end


text:
# Units
const ENERGY_UNIT = "eV"
const LENGTH_UNIT = "Å"
const TIME_UNIT = "ps"

const PARAMDIR = abspath(@__DIR__, "..", "params")
#const GLOBAL_SKT = loaddir(joinpath(PARAMDIR, "ob2-1-1", "base"))
const GLOBAL_SKT = loaddir(joinpath(PARAMDIR, "mio-1-1"))
const SLATERKOSTER_HAMILTONIAN = SlaterKoster.hamiltonian(GLOBAL_SKT)
const SLATERKOSTER_OVERLAP = SlaterKoster.overlap(GLOBAL_SKT)

const CC_SINGLE_BOND_LENGTH = 1.54
const CC_DOUBLE_BOND_LENGTH = 1.34
const CC_TRIPLE_BOND_LENGTH = 1.20
const CH_BOND_LENGTH = 1.09

const DEFAULT_α = 1.0
# In units of `electronvolt`
const DEFAULT_λ = 6e-3

# In units of `electronvolt * picosecond`
const HBAR = 6.582119569e-4

ext: jl


text:
struct UnidirectionalFlux{D, F}
    flux::F
end
UnidirectionalFlux{D}(f) where D = UnidirectionalFlux{D, eltype(f)}(f)

function directional_splitting(f::LinearAdvectionFlux{T}) where T<:SVector{2}
    z = zero(eltype(f.velocity))
    (UnidirectionalFlux{1}(LinearAdvectionFlux(SVector(f.velocity[1], z))),
     UnidirectionalFlux{2}(LinearAdvectionFlux(SVector(z, f.velocity[2]))))
end

function directional_splitting(f::LinearAdvectionFlux{T}) where T<:SVector{3}
    z = zero(eltype(f.velocity))
    (UnidirectionalFlux{1}(LinearAdvectionFlux(SVector(f.velocity[1], z, z))),
     UnidirectionalFlux{2}(LinearAdvectionFlux(SVector(z, f.velocity[2], z))),
     UnidirectionalFlux{3}(LinearAdvectionFlux(SVector(z, z, f.velocity[3]))))
end

function directional_splitting(f::AdvectionFlux{2})
    (
    UnidirectionalFlux{1}(AdvectionFlux{2}((m, i_face) -> SVector(f.velocity_at_face(m, i_face)[1], 0.0))),
    UnidirectionalFlux{2}(AdvectionFlux{2}((m, i_face) -> SVector(0.0, f.velocity_at_face(m, i_face)[2])))
    )
end

function numerical_flux(wrapped::UnidirectionalFlux{D}, mesh::CartesianMeshes.AbstractCartesianMesh, w, scheme, i_face, dt) where D
    if CartesianMeshes._direction(i_face) == D
        numerical_flux(wrapped.flux, mesh, w, scheme, i_face, dt)
    else
        zero(eltype(w))
    end
end

# ANOTHER OPTION: OVERLOADING DIV TO LOOP ONLY ON SOME FACES
# function faces_in_direction(::Val{D}, mesh::CartesianMeshes.AbstractCartesianMesh) where D
#     ...
# end
# function div!(Δw, wrapped::UnidirectionalFlux{D}, mesh::CartesianMeshes.AbstractCartesianMesh, w, scheme, i_face, dt) where D
#     @inbounds for i_face in faces_in_direction(Val(D), mesh)
#         ϕ = numerical_flux(flux, mesh, w, scheme, i_face, dt)
#         i_cell_1, i_cell_2 = cells_next_to_inner_face(mesh, i_face)
#         Δw[i_cell_1] += ϕ * face_area(mesh, i_face) / cell_volume(mesh, i_cell_1)
#         Δw[i_cell_2] -= ϕ * face_area(mesh, i_face) / cell_volume(mesh, i_cell_2)
#     end
# end

function courant(Δt, wrapped::UnidirectionalFlux{D}, mesh, w, i_face) where D
    if CartesianMeshes._direction(i_face) == D
        courant(Δt, wrapped.flux, mesh, w, i_face)
    else
        zero(Δt)
    end
end


text:
@reexport module MountainCar

export MountainCarEnv

using ReinforcementLearningBase

const ACTIONS = [-1, 0, 1]
const POSITION_MIN = -1.2
const POSITION_MAX = 0.5
const VELOCITY_MIN = -0.07
const VELOCITY_MAX = 0.07

Base.@kwdef mutable struct MountainCarEnv <: AbstractEnv
    p::Float64 = -0.6 + rand() / 5
    v::Float64 = 0.0
end

RLBase.get_observation_space(env::MountainCarEnv) = MultiContinuousSpace(
    [POSITION_MIN, VELOCITY_MIN],
    [POSITION_MAX, VELOCITY_MAX]
)

RLBase.get_action_space(env::MountainCarEnv) = DiscreteSpace(length(ACTIONS))

function (env::MountainCarEnv)(a)
    action = ACTIONS[a]
    v′ = min(
        max(VELOCITY_MIN, env.v + 0.001 * action - 0.0025 * cos(3 * env.p)),
        VELOCITY_MAX,
    )
    p′ = min(max(POSITION_MIN, env.p + v′), POSITION_MAX)
    if p′ == POSITION_MIN
        v′ = 0.0
    end
    env.p = p′
    env.v = v′
    nothing
end

RLBase.observe(env::MountainCarEnv) =
    (
        reward = env.p == POSITION_MAX ? 0.0 : -1.0,
        terminal = env.p == POSITION_MAX,
        state = [env.p, env.v],
    )

function RLBase.reset!(env::MountainCarEnv)
    env.p = -0.6 + rand() / 5
    env.v = 0.0
    nothing
end

end


text:
module SpinGlassTensors
    using LinearAlgebra
    using TensorOperations, TensorCast
    #using LowRankApprox
    using Memoize

    using DocStringExtensions

    include("base.jl")
    include("linear_algebra_ext.jl")
    include("compressions.jl")
    include("identities.jl")
    include("contractions.jl")
end # module


text:
function BeidouAlmanac(arg0::jint, arg1::jint, arg10::jdouble, arg11::jdouble, arg12::jdouble, arg13::jint, arg14::AbsoluteDate, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return BeidouAlmanac((jint, jint, jdouble, jdouble, jdouble, jint, AbsoluteDate, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg13, arg14, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function BeidouAlmanac(arg0::jint, arg1::jint, arg10::jdouble, arg11::jdouble, arg12::jdouble, arg13::jint, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return BeidouAlmanac((jint, jint, jdouble, jdouble, jdouble, jint, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg13, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_af0(obj::BeidouAlmanac)
    return jcall(obj, "getAf0", jdouble, ())
end

function get_af1(obj::BeidouAlmanac)
    return jcall(obj, "getAf1", jdouble, ())
end

function get_af2(obj::GNSSOrbitalElements)
    return jcall(obj, "getAf2", jdouble, ())
end

function get_aodc(obj::BeidouOrbitalElements)
    return jcall(obj, "getAODC", jint, ())
end

function get_aode(obj::BeidouOrbitalElements)
    return jcall(obj, "getAODE", jint, ())
end

function get_cic(obj::BeidouAlmanac)
    return jcall(obj, "getCic", jdouble, ())
end

function get_cis(obj::BeidouAlmanac)
    return jcall(obj, "getCis", jdouble, ())
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_crc(obj::BeidouAlmanac)
    return jcall(obj, "getCrc", jdouble, ())
end

function get_crs(obj::BeidouAlmanac)
    return jcall(obj, "getCrs", jdouble, ())
end

function get_cuc(obj::BeidouAlmanac)
    return jcall(obj, "getCuc", jdouble, ())
end

function get_cus(obj::BeidouAlmanac)
    return jcall(obj, "getCus", jdouble, ())
end

function get_date(obj::BeidouAlmanac)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_e(obj::BeidouAlmanac)
    return jcall(obj, "getE", jdouble, ())
end

function get_health(obj::BeidouAlmanac)
    return jcall(obj, "getHealth", jint, ())
end

function get_i0(obj::BeidouAlmanac)
    return jcall(obj, "getI0", jdouble, ())
end

function get_i_dot(obj::BeidouAlmanac)
    return jcall(obj, "getIDot", jdouble, ())
end

function get_iod(obj::BeidouOrbitalElements)
    return jcall(obj, "getIOD", jint, ())
end

function get_m0(obj::BeidouAlmanac)
    return jcall(obj, "getM0", jdouble, ())
end

function get_mean_motion(obj::BeidouAlmanac)
    return jcall(obj, "getMeanMotion", jdouble, ())
end

function get_omega0(obj::BeidouAlmanac)
    return jcall(obj, "getOmega0", jdouble, ())
end

function get_omega_dot(obj::BeidouAlmanac)
    return jcall(obj, "getOmegaDot", jdouble, ())
end

function get_pa(obj::BeidouAlmanac)
    return jcall(obj, "getPa", jdouble, ())
end

function get_prn(obj::BeidouAlmanac)
    return jcall(obj, "getPRN", jint, ())
end

function get_sma(obj::BeidouAlmanac)
    return jcall(obj, "getSma", jdouble, ())
end

function get_tgd1(obj::BeidouOrbitalElements)
    return jcall(obj, "getTGD1", jdouble, ())
end

function get_tgd2(obj::BeidouOrbitalElements)
    return jcall(obj, "getTGD2", jdouble, ())
end

function get_time(obj::BeidouAlmanac)
    return jcall(obj, "getTime", jdouble, ())
end

function get_toc(obj::GNSSOrbitalElements)
    return jcall(obj, "getToc", jdouble, ())
end

function get_week(obj::BeidouAlmanac)
    return jcall(obj, "getWeek", jint, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
using CustomUnitRanges
include(CustomUnitRanges.filename_for_urange)  # defines URange

type BidirectionalVector{T} <: AbstractVector{T}
    data::Vector{T}
    offset::Int
end
BidirectionalVector{T}(v::AbstractVector{T}, inds::AbstractUnitRange) =
    BidirectionalVector(copyelts(v), first(inds)-1)
BidirectionalVector(v::AbstractVector) = BidirectionalVector(v, Base.indices1(v))

# copies but doesn't preserve the indices
function copyelts{T}(v::AbstractVector{T})
    inds = Base.indices1(v)
    n = length(inds)
    dest = Array{T}(n)
    for (vel, j) in zip(v, 1:n)
        dest[j] = vel
    end
    dest
end

# Don't implement size or length
Base.indices1(v::BidirectionalVector) = URange(1+v.offset, length(v.data)+v.offset)
Base.indices( v::BidirectionalVector) = (Base.indices1(v),)

function Base.similar(v::AbstractArray, T::Type, inds::Tuple{URange})
    inds1 = inds[1]
    n = length(inds1)
    BidirectionalVector(Array{T}(n), first(inds1)-1)
end

function Base.similar(f::Union{Function,DataType}, inds::Tuple{URange})
    inds1 = inds[1]
    n = length(inds1)
    BidirectionalVector(f(Base.OneTo(n)), first(inds1)-1)
end

@inline function Base.getindex(v::BidirectionalVector, i::Int)
    @boundscheck checkbounds(v, i)
    @inbounds ret = v.data[i-v.offset]
    ret
end

@inline function Base.setindex!(v::BidirectionalVector, val, i::Int)
    @boundscheck checkbounds(v, i)
    @inbounds v.data[i-v.offset] = val
    val
end

Base.push!(v::BidirectionalVector, x) = (push!(v.data, x); v)
Base.pop!(v::BidirectionalVector) = pop!(v.data)
Base.append!(v::BidirectionalVector, collection2) = (append!(v.data, collection2); v)
function Base.prepend!(v::BidirectionalVector, collection2)
    v.offset -= _length(collection2)
    prepend!(v.data, collection2)
    v
end
function Base.shift!(v::BidirectionalVector)
    v.offset += 1
    shift!(v.data)
end
function Base.unshift!(v::BidirectionalVector, x)
    v.offset -= 1
    unshift!(v.data, x)
    v
end
@inline function Base.unshift!(v::BidirectionalVector, y...)
    v.offset -= length(y)
    unshift!(v.data, y...)
    v
end

function deletetail!(v::BidirectionalVector, n::Integer)
    deleteat!(v.data, length(v.data)-n+1:length(v.data))
    v
end
function deletehead!(v::BidirectionalVector, n::Integer)
    v.offset += n
    deleteat!(v.data, 1:n)
    v
end

_length(a::AbstractArray) = length(linearindices(a))


text:
"""
bayessupport \\

For calculating the support (number of non-zero probabilities) from a Bayesian point of view  \\
Arguments:  probs = naive probability distribution \\
            n_trials = number of trials \\
Returns:    a bayesian estimate for the support of the probability distribution  \\
"""
function bayessupport(probs::Array{Float64}, n_trials::Int)
  dim = length(probs)
  probs_nonzero = probs[probs .> 0]
  support_naive = length(probs_nonzero)
  support = support_naive
  if support_naive < dim
    support_expected = support_naive - sum((1 - probs_nonzero) .^ n_trials)
    prev_support_delta = dim
    support_delta = abs(support_naive - support_expected)
    xtr = 0.0
    while (support_delta < prev_support_delta) && ((support_naive+xtr)<dim)
      xtr = xtr + 1.0
      support_expected = 0.0
      # occupied bins
      gamma = xtr*(1.0 - (n_trials/(n_trials+support_naive)^(1/n_trials)))
      support_bayes = ((1.0-gamma)/(n_trials+support_naive))*(probs_nonzero*n_trials+1.0)
      suppport_expected = sum(1.0-(1.0-support_bayes)^n_trials)
      # non-occupied bins
      support_bayes = gamma / xtr
      support_expected = support_expected + xtr*(1.0 - (1.0 - support_bayes)^n_trials)
      prev_support_delta = support_delta
      support_delta = abs(support_naive - support_expected)
    end
    support_naive = support_naive + xtr - 1.0
    if support_delta < prev_support_delta
      support_naive += 1.0
    end
  end
return support_naive
end


text: is_unix() && run(`make`)


text:
# Key managment for braiding from PeaceVote.

mutable struct Key
    id::ID
    key
    generated::UInt64
    created::Union{Int, Nothing}
    expired::Union{Int, Nothing}
end

function newkey(notary::Notary)
    signer = newsigner(notary)

    thisid = id(signer)
    key = signer.key
    t = time_ns()

    return Key(thisid, key, t, nothing, nothing)
end

mutable struct KeyChain
    keys::Vector{Key}
    id::ID
    notary::Notary
    index::Int ### Keeps track on the position where the braidchain had been validated
end

KeyChain(keys::Vector{Key}, id::ID, notary::Notary) = KeyChain(keys, id, notary, 0)

PeaceCypher.id(kc::KeyChain) = kc.id

function generate_key!(kc::KeyChain)
    key = newkey(kc.notary)
    push!(kc.keys, key)

    return Signer(key.key, kc.notary)
    #return key.id
end

function KeyChain(notary::Notary)
    key = newkey(notary)
    kc = KeyChain(Key[key], key.id, notary)
    return kc
end

function set_created!(keys::Vector{Key}, id::ID, n::Int)

    for key in keys
        if key.id == id
            #@assert isnothing(key.created)
            key.created = n
        end
    end

    return nothing
end

function set_expired!(keys::Vector{Key}, id::ID, n::Int)

    for key in keys
        if key.id == id
            key.expired = n
        end
    end

    return nothing
end

function Base.push!(kc::KeyChain, x::Transaction)
    kc.index += 1
end

function Base.push!(kc::KeyChain, t::Transaction{Member})
    kc.index += 1

    cid = t.document.id
    set_created!(kc.keys, cid, kc.index)

    return nothing
end

function Base.push!(kc::KeyChain, t::Transaction{T}) where T <: Braid
    kc.index += 1

    for cid in pseudonyms(t.document)
        set_created!(kc.keys, cid, kc.index)
    end

    for eid in id(t.signature, kc.notary)
        set_expired!(kc.keys, eid, kc.index)
    end

    return nothing
end

function get_signer(kc::KeyChain, n::Int)
    for ki in kc.keys
        if ki.created <= n && ( isnothing(ki.expired) || ki.expired > n)
            return Signer(ki.key, kc.notary)
        end
    end
end

#get_signer(kc::KeyChain) = get_signer(kc, kc.index)

Base.lastindex(kc::KeyChain) = kc.index
Base.getindex(kc::KeyChain, n::Int) = get_signer(kc, n)

function update!(kc::KeyChain, bc::BraidChain)
    i0 = kc.index + 1
    for i in i0:length(bc.transactions)
        push!(kc, bc.transactions[i])
    end
end

function addtrigger!(kc::KeyChain, bc::BraidChain)

    transactions = bc.transactions
    trigger = Trigger(transactions, t -> update!(kc, bc))
    bc.transactions = trigger

    return nothing
end

function braid!(kc::KeyChain, bo::BraidOfficer)
    oldsigner = kc[end]
    newsigner = generate_key!(kc)

    braid(oldsigner, newsigner, bo)
end

# The API I would like to have

# vote(proposal, keychain)
# keychain[n] would return a valid signer whoose id is in the pseudonyms.
# @test id(keychain[n]) in pseudonyms(bc, n)
# braid!(keychain, braidchain). The braid command proceeds with last confirmed key.
# The keychain could already store a last index so it would not need to go through all transactions again.


text:
using LinearAlgebra: size
using Test
using SemiLagrangian

# include("mpoly/runtests.jl")
include("test_interpolation.jl")
include("test_lagrange.jl")
include("test_hermite.jl")
include("test_splinelu.jl")
include("test_bspline.jl")
include("testfftbig.jl")
include("test_util.jl")
include("test_mesh.jl")
include("test_advection.jl")
# include("test_poisson1d.jl")
include("test_poisson.jl")
include("test_poisson2d.jl")
# include("test_rotation1d.jl")
include("test_rotation.jl")
# include("test_translation1d.jl")
include("test_translation.jl")
include("test_swirling.jl")
include("test_quasigeostrophic.jl")


text:
module TwentyFour

import Combinatorics.multiset_permutations
import Base: length, show, getindex

export XXIV, solve

"""
`solve(nums...)` solves a Twenty-Four puzzle for the given numbers.
Returns a `String` giving the solution.
Example:
```
julia> solve(4,2,5,1)
"(4+2)*(5-1)"
```
"""
function solve(nums...)::String
  @assert length(nums)>1 "Please give two or more integers or rationals"
  X = XXIV(nums...)
  solver(X)
  if X.solved
    return X.solution[2:end-1]
  else
    return "No solution"
  end
end

"""
The default goal of a *Twenty Four* problem is 24.
"""
global DEFAULT_GOAL = 24

global const OPS = [ (+), (-), (*), (/)]
global const OP_SYMS = [ "+", "-", "*", "/"]

"""
`XXIV` is an instance of a Twenty Four problem.
"""
mutable struct XXIV
  nums::Vector{Rational{Int}}
  numstrs::Vector{String}
  goal::Int
  solved::Bool
  solution::String

  function XXIV(num_list::Vector{Rational{Int}})
    new(deepcopy(num_list),map(str,num_list),DEFAULT_GOAL,false,"")
  end
end

function XXIV(vals...)
  num_list = [ Rational(x) for x in vals ]
  return XXIV(num_list)
end

XXIV() = XXIV(Rational{Int}[])

getindex(X::XXIV,i::Int) = X.nums[i]

function show(io::IO,X::XXIV)
  n = length(X)
  print(io, "[ ")
  for a in X.numstrs
    print(io,a*" ")
  end
  print(io,"] with goal $(X.goal)")
  if X.solved
    print(io," and solution $(X.solution)")
  end
  nothing
end

length(X::XXIV) = length(X.nums)

"""
`get_goal(X::XXIV)` returns the targe goal of this instance
of a *Twenty Four* problem. This is typically 24.
"""
get_goal(X::XXIV) = X.goal

"""
`set_goal(X::XXIV, g)` sets the goal of this *Twenty Four* problem
to `g`. Default `g` is 24.
"""
function set_goal(X::XXIV, g::Int=DEFAULT_GOAL)
  if X.goal == g
    return nothing
  end
  X.goal = g
  X.solved=false
  X.solution=""
  nothing
end

"""
`str(x::Rational{Int})` returns a nice string form of `x`.
"""
function str(x::Rational{Int})
  if denominator(x)==1
    return "$(numerator(x))"
  end
  return "$(numerator(x))/$(denominator(x))"
end

function solver(X::XXIV)
  n = length(X)
  # basis cases
  if n==0
    X.solved = false
    X.solution = ""
    return 0
  end

  if n == 1
    a = X[1]
    if a == X.goal
      X.solved = true
      X.solution = X.numstrs[1]
      return true
    else
      X.solved = false
      X.solution = ""
      return false
    end
  end

  for i=1:n
    for j=1:n
      if i!=j
        ii = i
        jj = j
        if i>j
          ii,jj = jj,ii
        end

        a = X.nums[i]
        sa = X.numstrs[i]
        b = X.nums[j]
        sb = X.numstrs[j]

        for t = 1:4
          XX = deepcopy(X)

          f = OPS[t]
          sf = OP_SYMS[t]

          ab = f(a,b)
          sab = "(" * sa * sf * sb * ")"
          deleteat!(XX.nums,(ii,jj))
          deleteat!(XX.numstrs,(ii,jj))
          append!(XX.nums,[ab])
          append!(XX.numstrs,[sab])
          try
            if solver(XX)
              X.solution = XX.solution
              X.solved = true
              return true
            end
            catch
          end

        end # next t

      end  # endif
    end # next j
  end # next i
  return false

end  # end of function

end # end of module


text:
using Manifolds, Manopt, Test, ManifoldsBase

@testset "Differentials (on Sphere(2))" begin
    # The Adjoint Differentials test using the same variables as the differentials
    # test
    x = [1.0, 0.0, 0.0]
    y = [0.0, 1.0, 0.0]
    M = Sphere(2)
    ξ = log(M, x, y)
    # Text differentials (1) Dx of Log_xy
    @test adjoint_differential_log_basepoint(M, x, x, ξ) == -ξ
    @test adjoint_differential_log_argument(M, x, x, ξ) == ξ
    @test adjoint_differential_exp_basepoint(M, x, zero_tangent_vector(M, x), ξ) == ξ
    @test adjoint_differential_exp_argument(M, x, zero_tangent_vector(M, x), ξ) == ξ
    for t in [0, 0.15, 0.33, 0.66, 0.9]
        @test adjoint_differential_geodesic_startpoint(M, x, x, t, ξ) == (1 - t) * ξ
        @test norm(M, x, adjoint_differential_geodesic_endpoint(M, x, x, t, ξ) - t * ξ) ≈ 0 atol =
            10.0^(-16)
    end
    Mp = PowerManifold(M, NestedPowerRepresentation(), 3)
    xP = [x, y, x]
    yP = [x, x, y]
    ξP = [ξ, zero_tangent_vector(M, x), -ξ]
    @test norm(
        Mp,
        xP,
        adjoint_differential_forward_logs(Mp, xP, ξP) -
        [-ξ, zero_tangent_vector(M, x), zero_tangent_vector(M, x)],
    ) ≈ 0 atol = 4 * 10.0^(-16)
end


text:
module QuantumStateDB

using JSON3
using LibPQ
using DataFrames
using Tables
using InteractiveUtils: subtypes

include("utils.jl")
include("db.jl")
include("table.jl")
include("api.jl")

end # module


text:
abstract type AbstractSparseIdentificationResult end;

mutable struct SparseIdentificationResult <: AbstractSparseIdentificationResult
    coeff::AbstractArray
    parameters::AbstractArray
    equations::Union{Function, Basis}

    opt::Union{Optimise.AbstractOptimiser, Optimise.AbstractSubspaceOptimiser}
    iterations::Int64
    converged::Bool

    error::AbstractArray
    aicc::AbstractArray
    sparsity::AbstractArray
end

Base.show(io::IO, x::SparseIdentificationResult) = print(io, "Sparse Identification Result with $(sum(x.sparsity)) active terms.")

@inline function Base.print(io::IO, x::SparseIdentificationResult)
    println("Sparse Identification Result")
    println("No. of Parameters : $(length(x.parameters))")
    println("Active terms : $(sum(x.sparsity))")
    for (i, si) in enumerate(x.sparsity)
        println("   Equation $i : $si")
    end
    println("Overall error (L2-Norm) : $(sum(x.error))")
    for (i, ei) in enumerate(x.error)
        println("   Equation $i : $ei")
    end
    println("AICC :")
    for (i, ai) in enumerate(x.aicc)
        println("   Equation $i : $ai")
    end

    print("\n$(x.opt)")
    if x.converged
        println(" converged after $(x.iterations) iterations.")
    else
        println(" did not converge after $(x.iterations) iterations.")
    end
end

"""
    print_equations(res; show_parameter)

Print the equations stored inside the `SparseIdentificationResult` `res`. If `show_parameter` is set
to true, the numerical values will be used. Otherwise the symbolic form will appear.
"""
function print_equations(r::SparseIdentificationResult; show_parameter::Bool = false)
    if show_parameter
        eqs = r.equations(variables(r.equations), parameters(r), independent_variable(r.equations))
        for (i, eq) in enumerate(eqs)
            println("f_$i = ", eq)
        end
    else
        print(r.equations)
    end
end

"""
    SparseIdentificationResult()

Contains the result of a sparse identification. Contains the coefficient matrix `Ξ`,
the equations of motion and its associated parameters. It stores also the optimiser, iteration counter and convergence status
is stored.

Additionally, the model is evaluated over the training data and the ``L_2``-error, Akaikes Information Criterion and the ``L_0``-Norm of the coefficients
is stored.
"""
function SparseIdentificationResult(coeff::AbstractArray, equations::Basis, iters::Int64, opt::T , convergence::Bool, Y::AbstractVecOrMat, X::AbstractVecOrMat; p::AbstractArray = [], t::AbstractVector = []) where T <: Union{Optimise.AbstractOptimiser, Optimise.AbstractSubspaceOptimiser}
    Ŷ = coeff'*equations(X, p, t)
    training_error = norm.(eachrow(Y-Ŷ), 2)
    sparsity = Int64.(norm.(eachcol(coeff), 0))

    aicc = similar(training_error)
    for i in 1:length(aicc)
        aicc[i] = AICC(sparsity[i], view(Ŷ, i, :) , view(Y, i, :))
    end
    b_, p_ = derive_parameterized_eqs(coeff, equations, sum(sparsity))
    return SparseIdentificationResult(coeff, [p...;p_...], b_ , opt, iters, convergence,  training_error, aicc,  sparsity)
end

function derive_parameterized_eqs(Ξ::AbstractArray{T, 2}, b::Basis, sparsity::Int64) where T <: Real
    @parameters p[1:sparsity]
    p_ = zeros(eltype(Ξ), sparsity)
    cnt = 1
    b_ = Basis(Operation[], variables(b), parameters = [parameters(b)...; p...], iv = independent_variable(b))

    for i=1:size(Ξ, 2)
        eq = nothing
        for j = 1:size(Ξ, 1)
            if !iszero(Ξ[j,i])
                if eq === nothing
                    eq = p[cnt]*b[j]
                else
                    eq += p[cnt]*b[j]
                end
                p_[cnt] = Ξ[j,i]
                cnt += 1
            end
        end
        if !isnothing(eq)
            push!(b_, eq)
        end
    end
    b_, p_
end

Base.size(r::SparseIdentificationResult) = size(r.sparsity)
Base.length(r::SparseIdentificationResult) = length(r.sparsity)

ModelingToolkit.parameters(r::SparseIdentificationResult) = r.parameters

dynamics(b::SparseIdentificationResult) = dynamics(b.equations)

"""
    get_sparsity(res)

Return the ``L_0``-Norm of the `SparseIdentificationResult`
"""
get_sparsity(r::SparseIdentificationResult) = r.sparsity

"""
    get_error(res)

Return the ``L_2``-Error of the `SparseIdentificationResult` over the training data.
"""
get_error(r::SparseIdentificationResult) = r.error

"""
    get_aicc(res)

Return Akaikakes Information Criterion  of the `SparseIdentificationResult` over the training data.
"""
get_aicc(r::SparseIdentificationResult) = r.aicc

"""
    get_coefficients(res)

Return the coefficient matrix `Ξ` of the `SparseIdentificationResult`.
"""
get_coefficients(r::SparseIdentificationResult) = r.coeff

function ModelingToolkit.ODESystem(b::SparseIdentificationResult)
    return ODESystem(b.equations)
end


text:
module SimpleReassortmentGraph

using Parameters
using Random
using TreeKnit
using TreeTools

import Base: ==
import Base: isequal, push!, write

const SRG = SimpleReassortmentGraph
export SRG

include("objects.jl")
include("misc.jl")
include("prunegraft.jl")
include("construct.jl")
include("trees.jl")
include("IO.jl")

end # module

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text: include("grid.jl")


text:
# Copyright (C) 2020 Argonne National Laboratory
# Written by Alinson Santos Xavier <axavier@anl.gov>

using RELOG

@testset "parse" begin
    basedir = dirname(@__FILE__)
    instance = RELOG.parsefile("$basedir/../../instances/s1.json")

    centers = instance.collection_centers
    plants = instance.plants
    products = instance.products
    location_name_to_plant = Dict(p.location_name => p for p in plants)
    product_name_to_product = Dict(p.name => p for p in products)

    @test length(centers) == 10
    @test centers[1].name == "C1"
    @test centers[1].latitude == 7
    @test centers[1].latitude == 7
    @test centers[1].longitude == 7
    @test centers[1].amount == [934.56, 934.56]
    @test centers[1].product.name == "P1"

    @test length(plants) == 6

    plant = location_name_to_plant["L1"]
    @test plant.plant_name == "F1"
    @test plant.location_name == "L1"
    @test plant.input.name == "P1"
    @test plant.latitude == 0
    @test plant.longitude == 0

    @test length(plant.sizes) == 2
    @test plant.sizes[1].capacity == 250
    @test plant.sizes[1].opening_cost == [500, 500]
    @test plant.sizes[1].fixed_operating_cost == [30, 30]
    @test plant.sizes[1].variable_operating_cost == [30, 30]
    @test plant.sizes[2].capacity == 1000
    @test plant.sizes[2].opening_cost == [1250, 1250]
    @test plant.sizes[2].fixed_operating_cost == [30, 30]
    @test plant.sizes[2].variable_operating_cost == [30, 30]

    p2 = product_name_to_product["P2"]
    p3 = product_name_to_product["P3"]
    @test length(plant.output) == 2
    @test plant.output[p2] == 0.2
    @test plant.output[p3] == 0.5
    @test plant.disposal_limit[p2] == [1, 1]
    @test plant.disposal_limit[p3] == [1, 1]
    @test plant.disposal_cost[p2] == [-10, -10]
    @test plant.disposal_cost[p3] == [-10, -10]

    plant = location_name_to_plant["L3"]
    @test plant.location_name == "L3"
    @test plant.input.name == "P2"
    @test plant.latitude == 25
    @test plant.longitude == 65

    @test length(plant.sizes) == 2
    @test plant.sizes[1].capacity == 1000.0
    @test plant.sizes[1].opening_cost == [3000, 3000]
    @test plant.sizes[1].fixed_operating_cost == [50, 50]
    @test plant.sizes[1].variable_operating_cost == [50, 50]
    @test plant.sizes[1] == plant.sizes[2]

    p4 = product_name_to_product["P4"]
    @test plant.output[p3] == 0.05
    @test plant.output[p4] == 0.8
    @test plant.disposal_limit[p3] == [1e8, 1e8]
    @test plant.disposal_limit[p4] == [0, 0]
end

@testset "parse (geodb)" begin
    basedir = dirname(@__FILE__)
    instance = RELOG.parsefile("$basedir/../../instances/s2.json")

    centers = instance.collection_centers
    @test centers[1].name == "C1"
    @test centers[1].latitude == 41.83956
    @test centers[1].longitude == -88.08857
end

# @testset "parse (invalid)" begin
#     basedir = dirname(@__FILE__)
#     @test_throws ErrorException RELOG.parsefile("$basedir/../fixtures/s1-wrong-length.json")
# end

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
struct LazyMulAdd{M,O,T<:NativeTypesV} <: Number
    data::T
    @inline LazyMulAdd{M,O,T}(data::T) where {M,O,T<:NativeTypesV} = new{M,O,T}(data)
end
# O for offset is kind of hard to read next to default of 0?
@inline LazyMulAdd{M,O}(data::T) where {M,O,T<:Union{Base.HWReal,AbstractSIMD}} = LazyMulAdd{M,O,T}(data)
@inline LazyMulAdd{M}(data::T) where {M,T<:NativeTypesV} = LazyMulAdd{M,0,T}(data)
@inline LazyMulAdd{0,O}(data::T) where {O,T<:Union{Base.HWReal,AbstractSIMD}} = @assert false#StaticInt{O}()
@inline LazyMulAdd{0}(data::T) where {T<:NativeTypesV} = @assert false#StaticInt{0}()
@inline LazyMulAdd(data::T, ::StaticInt{M}) where {M,T} = LazyMulAdd{M,0,T}(data)
@inline LazyMulAdd(data::T, ::StaticInt{M}, ::StaticInt{O}) where {M,O,T} = LazyMulAdd{M,O,T}(data)

@inline data(lm::LazyMulAdd) = data(getfield(lm, :data)) # calls data on inner for use with indexing (normally `data` only goes through one layer)

@inline _materialize(a::LazyMulAdd{M,O,I}) where {M,O,I} = vadd_fast(vmul_fast(StaticInt{M}(), getfield(a, :data)), StaticInt{O}())
@inline _materialize(x) = x
@inline Base.convert(::Type{T}, a::LazyMulAdd{M,O,I}) where {M,O,I,T<:Number} = convert(T, _materialize(a))
@inline Base.convert(::Type{LazyMulAdd{M,O,I}}, a::LazyMulAdd{M,O,I}) where {M,O,I} = a
# @inline Base.convert(::Type{LazyMulAdd{M,O,I}}, a::LazyMulAdd{M}) where {M,O,I} = a
# @inline Base.convert(::Type{LazyMulAdd{M,T,I}}, a::LazyMulAdd{M,StaticInt{O},I}) where {M,O,I,T} = a

Base.promote_rule(::Type{LazyMulAdd{M,O,I}}, ::Type{T}) where {M,O,I<:Number,T} = promote_type(I, T)
Base.promote_rule(::Type{LazyMulAdd{M,O,MM{W,X,I}}}, ::Type{T}) where {M,O,W,X,I,T} = promote_type(MM{W,X,I}, T)
Base.promote_rule(::Type{LazyMulAdd{M,O,Vec{W,I}}}, ::Type{T}) where {M,O,W,I,T} = promote_type(Vec{W,I}, T)

@inline lazymul(a, b) = vmul_fast(a, b)
@inline lazymul(::StaticInt{M}, b) where {M} = LazyMulAdd{M}(b)
@inline lazymul(::StaticInt{1}, b) = b
@inline lazymul(::StaticInt{0}, b) = StaticInt{0}()
@inline lazymul(a, ::StaticInt{M}) where {M} = LazyMulAdd{M}(a)
@inline lazymul(a, ::StaticInt{1}) = a
@inline lazymul(a, ::StaticInt{0}) = StaticInt{0}()
@inline lazymul(a::LazyMulAdd, ::StaticInt{1}) = a
@inline lazymul(::StaticInt{1}, a::LazyMulAdd) = a
@inline lazymul(a::LazyMulAdd, ::StaticInt{0}) = StaticInt{0}()
@inline lazymul(::StaticInt{0}, a::LazyMulAdd) = StaticInt{0}()
# @inline lazymul(a::LazyMulAdd, ::StaticInt{0}) = StaticInt{0}()
# @inline lazymul(::StaticInt{M}, b::MM{W,X}) where {W,M,X} = LazyMulAdd{M}(MM{W}(getfield(b, :data), StaticInt{M}()*StaticInt{X}()))
# @inline lazymul(a::MM{W,X}, ::StaticInt{M}) where {W,M,X} = LazyMulAdd{M}(MM{W}(getfield(a, :data), StaticInt{M}()*StaticInt{X}()))
@inline lazymul(a::MM{W,X}, ::StaticInt{M}) where {W,M,X} = MM{W}(vmul_fast(StaticInt{M}(), data(a)), StaticInt{X}() * StaticInt{M}())
@inline lazymul(::StaticInt{M}, b::MM{W,X}) where {W,M,X} = MM{W}(vmul_fast(StaticInt{M}(), data(b)), StaticInt{X}() * StaticInt{M}())
@inline lazymul(a::MM{W,X}, ::StaticInt{1}) where {W,X} = a
@inline lazymul(::StaticInt{1}, a::MM{W,X}) where {W,X} = a
@inline lazymul(a::MM{W,X}, ::StaticInt{0}) where {W,X} = StaticInt{0}()
@inline lazymul(::StaticInt{0}, a::MM{W,X}) where {W,X} = StaticInt{0}()

@inline lazymul(::StaticInt{M}, ::StaticInt{N}) where {M, N} = StaticInt{M}()*StaticInt{N}()

@inline lazymul(::StaticInt{M}, ::StaticInt{1}) where {M} = StaticInt{M}()
@inline lazymul(::StaticInt,    ::StaticInt{0}) = StaticInt{0}()

@inline lazymul(::StaticInt{1}, ::StaticInt{M}) where {M} = StaticInt{M}()
@inline lazymul(::StaticInt{0}, ::StaticInt   ) = StaticInt{0}()

@inline lazymul(::StaticInt{0}, ::StaticInt{0}) = StaticInt{0}()
@inline lazymul(::StaticInt{0}, ::StaticInt{1}) = StaticInt{0}()
@inline lazymul(::StaticInt{1}, ::StaticInt{0}) = StaticInt{0}()
@inline lazymul(::StaticInt{1}, ::StaticInt{1}) = StaticInt{1}()

@inline function lazymul(a::LazyMulAdd{M,O}, ::StaticInt{N}) where {M,O,N}
    LazyMulAdd(getfield(a, :data), StaticInt{M}() * StaticInt{N}(), StaticInt{O}() * StaticInt{N}())
end
@inline function lazymul(::StaticInt{M}, b::LazyMulAdd{N,O}) where {M,O,N}
    LazyMulAdd(getfield(b, :data), StaticInt{M}() * StaticInt{N}(), StaticInt{O}() * StaticInt{M}())
end

@inline function lazymul(a::LazyMulAdd{M,<:MM{W,X}}, ::StaticInt{N}) where {M,N,W,X}
    LazyMulAdd(MM{W}(getfield(a, :data), StaticInt{M}()*StaticInt{N}()*StaticInt{X}()), StaticInt{M}()*StaticInt{N}())
end
@inline function lazymul(::StaticInt{M}, b::LazyMulAdd{N,<:MM{W,X}}) where {M,N,W,X}
    LazyMulAdd(MM{W}(getfield(b, :data), StaticInt{M}()*StaticInt{N}()*StaticInt{X}()), StaticInt{M}()*StaticInt{N}())
end

@inline lazymul(a::LazyMulAdd{M,<:MM{W,X}}, ::StaticInt{0}) where {M,W,X} = StaticInt{0}()
@inline lazymul(::StaticInt{0}, b::LazyMulAdd{N,<:MM{W,X}}) where {N,W,X} = StaticInt{0}()
@inline lazymul(a::LazyMulAdd{M,<:MM{W,X}}, ::StaticInt{1}) where {M,W,X} = a
@inline lazymul(::StaticInt{1}, b::LazyMulAdd{N,<:MM{W,X}}) where {N,W,X} = b
@inline function lazymul(a::LazyMulAdd{M}, b::LazyMulAdd{N}) where {M,N}
    LazyMulAdd(vmul_fast(getfield(a, :data), getfield(b, :data)), StaticInt{M}()*StaticInt{N}())
end

vmul_fast(::LazyMulAdd{M,O}, ::StaticInt{0}) where {M,O} = Zero()
vmul_fast(::StaticInt{0}, ::LazyMulAdd{M,O}) where {M,O} = Zero()
vmul_fast(a::LazyMulAdd{M,O}, ::StaticInt{1}) where {M,O} = a
vmul_fast(::StaticInt{1}, a::LazyMulAdd{M,O}) where {M,O} = a
@inline function vmul_fast(a::LazyMulAdd{M,O}, ::StaticInt{N}) where {M,N,O}
    LazyMulAdd(getfield(a, :data), StaticInt{M}() * StaticInt{N}(), StaticInt{O}() * StaticInt{N}())
end
@inline function vmul_fast(::StaticInt{N}, a::LazyMulAdd{M,O}) where {M,N,O}
    LazyMulAdd(getfield(a, :data), StaticInt{M}() * StaticInt{N}(), StaticInt{O}() * StaticInt{N}())
end

@inline function Base.:(>>>)(a::LazyMulAdd{M,O}, ::StaticInt{N}) where {M,O,N}
    LazyMulAdd(getfield(a, :data), StaticInt{M}() >>> StaticInt{N}(), StaticInt{O}() >>> StaticInt{N}())
end

# The approach with `add_indices` is that we try and make `vadd_fast` behave well
# but for `i` and `j` type combinations where it's difficult,
# we can add specific `add_indices` methods that increment the pointer.
@inline function add_indices(p::Ptr, i, j) # generic fallback
    p, vadd_fast(i, j)
end
@inline vadd_fast(i::LazyMulAdd, ::Zero) = i
@inline vadd_fast(::Zero, i::LazyMulAdd) = i

# These following two definitions normally shouldn't be hit
@inline function vadd_fast(a::LazyMulAdd{M,O,MM{W,X,I}}, b::Integer) where {M,O,W,X,I}
    MM{W}(vadd_fast(vmul_fast(StaticInt{M}(), data(a)), b), StaticInt{X}() * StaticInt{M}())
end
@inline function vadd_fast(b::Integer, a::LazyMulAdd{M,O,MM{W,X,I}}) where {M,O,W,X,I}
    MM{W}(vadd_fast(vmul_fast(StaticInt{M}(), data(a)), b), StaticInt{X}() * StaticInt{M}())
end
@inline vadd_fast(a::LazyMulAdd{M,O,MM{W,X,I}}, ::Zero) where {M,O,W,X,I} = a
@inline vadd_fast(::Zero, a::LazyMulAdd{M,O,MM{W,X,I}}) where {M,O,W,X,I} = a

# because we should hit this instead:
@inline add_indices(p::Ptr, b::Integer, a::LazyMulAdd{M,O}) where {M,O} = (p + b, a)
@inline add_indices(p::Ptr, a::LazyMulAdd{M,O}, b::Integer) where {M,O} = (p + b, a)
# but in the case of `VecUnroll`s, which skip the `add_indices`, it's useful to still have the former two definitions.
# However, this also forces us to write:
@inline add_indices(p::Ptr, ::StaticInt{N}, a::LazyMulAdd{M,O}) where {M,O,N} = (p, vadd_fast(a, StaticInt{N}()))
@inline add_indices(p::Ptr, a::LazyMulAdd{M,O}, ::StaticInt{N}) where {M,O,N} = (p, vadd_fast(a, StaticInt{N}()))

@inline function vadd_fast(::StaticInt{N}, a::LazyMulAdd{M,O}) where {N,M,O}
    LazyMulAdd(getfield(a, :data), StaticInt{M}(), StaticInt{O}()+StaticInt{N}())
end
@inline function vadd_fast(a::LazyMulAdd{M,O}, ::StaticInt{N}) where {N,M,O}
    LazyMulAdd(getfield(a, :data), StaticInt{M}(), StaticInt{O}()+StaticInt{N}())
end
@inline function vadd_fast(::StaticInt{N}, a::LazyMulAdd{M,O,MM{W,X,I}}) where {N,M,O,W,X,I}
    LazyMulAdd(getfield(a, :data), StaticInt{M}(), StaticInt{O}()+StaticInt{N}())
end
@inline function vadd_fast(a::LazyMulAdd{M,O,MM{W,X,I}}, ::StaticInt{N}) where {N,M,O,W,X,I}
    LazyMulAdd(getfield(a, :data), StaticInt{M}(), StaticInt{O}()+StaticInt{N}())
end

@inline function vadd_fast(a::LazyMulAdd{M,O}, b::LazyMulAdd{M,A}) where {M,O,A}
    LazyMulAdd(vadd_fast(getfield(a, :data), getfield(b, :data)), StaticInt{M}(), StaticInt{O}()+StaticInt{A}())
end

@inline add_indices(p::Ptr, a::LazyMulAdd{M,O,V}, b::LazyMulAdd{N,P,J}) where {M,O,V<:AbstractSIMDVector,N,P,J<:IntegerTypes} = (gep(p, b), a)
@inline add_indices(p::Ptr, b::LazyMulAdd{N,P,J}, a::LazyMulAdd{M,O,V}) where {M,O,V<:AbstractSIMDVector,N,P,J<:IntegerTypes} = (gep(p, b), a)
@inline add_indices(p::Ptr, a::LazyMulAdd{M,O,V}, b::LazyMulAdd{M,P,J}) where {M,O,V<:AbstractSIMDVector,P,J<:IntegerTypes} = (p, vadd_fast(a, b))
@inline add_indices(p::Ptr, b::LazyMulAdd{M,P,J}, a::LazyMulAdd{M,O,V}) where {M,O,V<:AbstractSIMDVector,P,J<:IntegerTypes} = (p, vadd_fast(a, b))

@inline add_indices(p::Ptr, a::AbstractSIMDVector, b::LazyMulAdd{M,O,I}) where {M,O,I<:IntegerTypes} = (gep(p, b), a)
@inline add_indices(p::Ptr, b::LazyMulAdd{M,O,I}, a::AbstractSIMDVector) where {M,O,I<:IntegerTypes} = (gep(p, b), a)
@inline function add_indices(p::Ptr, ::MM{W,X,StaticInt{A}}, a::LazyMulAdd{M,O,T}) where {M,O,T<:IntegerTypes,A,W,X}
    gep(p, a), MM{W,X}(StaticInt{A}())
end
@inline function add_indices(p::Ptr, a::LazyMulAdd{M,O,T}, ::MM{W,X,StaticInt{A}}) where {M,O,T<:IntegerTypes,A,W,X}
    gep(p, a), MM{W,X}(StaticInt{A}())
end

@generated function add_indices(p::Ptr, a::LazyMulAdd{M,O,MM{W,X,StaticInt{I}}}, b::LazyMulAdd{N,P,J}) where {M,O,W,X,I,N,P,J<:IntegerTypes}
    d, r = divrem(M, N)
    if iszero(r)
        quote
            $(Expr(:meta,:inline))
            p, VectorizationBase.LazyMulAdd{$N,$(I*M)}(MM{$W,$d}(getfield(b, :data)))
        end
    else
        quote
            $(Expr(:meta,:inline))
            gep(p, b), a
        end
    end
end
@inline add_indices(p::Ptr, b::LazyMulAdd{N,P,J}, a::LazyMulAdd{M,O,MM{W,X,StaticInt{I}}}) where {M,O,W,X,I,N,P,J<:IntegerTypes} = add_indices(p, a, b)
@generated function vadd_fast(a::LazyMulAdd{M,O,MM{W,X,StaticInt{I}}}, b::LazyMulAdd{N,P,J}) where {M,O,W,X,I,N,P,J<:IntegerTypes}
    d, r = divrem(M, N)
    if iszero(r)
        quote
            $(Expr(:meta,:inline))
            VectorizationBase.LazyMulAdd{$N,$(I*M)}(MM{$W,$d}(getfield(b, :data)))
        end
    else
        quote
            $(Expr(:meta,:inline))
            vadd_fast(a, _materialize(b))
        end
    end
end
@inline vadd_fast(b::LazyMulAdd{N,P,J}, a::LazyMulAdd{M,O,MM{W,X,StaticInt{I}}}) where {M,O,W,X,I,N,P,J<:IntegerTypes} = vadd_fast(a, b)

@inline vadd_fast(a::VecUnroll, b::LazyMulAdd) = VecUnroll(fmap(vadd_fast, getfield(a, :data), b))
@inline vadd_fast(b::LazyMulAdd, a::VecUnroll) = VecUnroll(fmap(vadd_fast, b, getfield(a, :data)))


text:
function packmultihypo(fnc::CommonConvWrapper{T}) where {T<:FunctorInferenceType}
  @warn "packmultihypo is deprecated in favor of Vector only operations"
  fnc.hypotheses !== nothing ? string(fnc.hypotheses) : ""
end
function parsemultihypostr(str::AS) where {AS <: AbstractString}
  @warn "parsemultihypostr is deprecated in favor of Vector only operations"
  mhcat=nothing
  if length(str) > 0
    mhcat = convert(SamplableBelief, str)
  end
  return mhcat
end

## packing converters-----------------------------------------------------------
# heavy use of multiple dispatch for converting between packed and original data types during DB usage

function convert(::Type{PackedFunctionNodeData{P}}, d::FunctionNodeData{T}) where {P <: PackedInferenceType, T <: FactorOperationalMemory}
  # mhstr = packmultihypo(d.fnc)  # this is where certainhypo error occurs
  return PackedFunctionNodeData(d.eliminated, d.potentialused, d.edgeIDs,
          convert(P, _getCCW(d).usrfnc!),
          d.multihypo, _getCCW(d).certainhypo, d.nullhypo, d.solveInProgress)  # extract two values from ccw for storage -- ccw thrown away
end

## unpack converters------------------------------------------------------------

function convert(
            ::Type{IncrementalInference.GenericFunctionNodeData{IncrementalInference.CommonConvWrapper{F}}},
            packed::IncrementalInference.GenericFunctionNodeData{P} ) where {F <: FunctorInferenceType, P <: PackedInferenceType}
  #
  # TODO store threadmodel=MutliThreaded,SingleThreaded in persistence layer
  usrfnc = convert(F, packed.fnc)
  mhcat, nh = parseusermultihypo(packed.multihypo, packed.nullhypo)

  # TODO -- improve prepgenericconvolution for hypotheses and certainhypo field recovery when deserializing
  # reconstitute from stored data
  # FIXME, add threadmodel=threadmodel
  ccw = prepgenericconvolution(DFG.DFGVariable[], usrfnc, multihypo=mhcat, nullhypo=nh)
  ccw.certainhypo = packed.certainhypo

  ret = FunctionNodeData{CommonConvWrapper{typeof(usrfnc)}}(packed.eliminated, packed.potentialused, packed.edgeIDs, ccw,
                                                            packed.multihypo, packed.certainhypo, packed.nullhypo, packed.solveInProgress)
  #
  return ret
end

## Variables

"""
  $(SIGNATURES)
After deserializing a factor using decodePackedType, use this to
completely rebuild the factor's CCW and user data.
"""
function rebuildFactorMetadata!(dfg::AbstractDFG{SolverParams}, factor::DFGFactor)
  # Set up the neighbor data
  neighbors = map(vId->getVariable(dfg, vId), getNeighbors(dfg, factor))
  neighborUserData = map(v->getVariableType(v), neighbors)

  # Rebuilding the CCW
  fsd = getSolverData(factor)
  ccw_new = getDefaultFactorData(dfg, neighbors, getFactorType(factor),
                                 multihypo=fsd.multihypo,
                                 nullhypo=fsd.nullhypo,
                                 eliminated=fsd.eliminated,
                                 potentialused=fsd.potentialused,
                                 edgeIDs=fsd.edgeIDs,
                                 solveInProgress=fsd.solveInProgress)
  setSolverData!(factor, ccw_new)

  #... Copying neighbor data into the factor?
  # JT TODO it looks like this is already updated in getDefaultFactorData -> prepgenericconvolution
  # factormetadata.variableuserdata is deprecated, remove when removing deprecation
  # for i in 1:Threads.nthreads()
  #   ccw_new.fnc.cpt[i].factormetadata.variableuserdata = deepcopy(neighborUserData)
  # end

  return factor
end

# import IncrementalInference: decodefg, loadjld

veeCategorical(val::Categorical) = val.p
veeCategorical(val::Union{Nothing, Vector{Float64}}) = val

function convert(::Type{Tuple{BallTreeDensity,Float64}},
                 p::TreeBelief )
  # @show size(p.val), size(p.bw), p.manifolds
  # (AMP.manikde!(p.val, p.bw[:,1], p.manifolds), p.inferdim)
  (convert(BallTreeDensity, p), p.inferdim)
end

#


text:
"""
    subtract(x, y)

Perform subtraction between two floating point numbers.

# Arguments
- `a::Float64`: Number 1.
- `b::Float64`: Number 2.

# Returns
- `a - b::Float64`

# Examples
```julia-repl
julia> subtract(1.0, 2.0)
-1.0
```
"""
subtract(a::Float64, b::Float64) = a - b


text:
# imports
include("../iron_line_profile.jl")
using Measures
using DelimitedFiles
gr()

function iron_line_heatmap(;
                            tolerance=1e-12,
                            size_multiplier=1,
                            fov=12,
                            dtmax=10000,
                            nbins=100,
                            vars=LinRange(5.0, 85.0, 20)
                            )
    # nbins=600
    # vars = LinRange(5.0, 85.0, 400)
    # vars = 5.0:5.0:85.0

    # values to be plotted
    bins_vals = Array{Float64}[]
    x_vals_vals = []

    # generating data for different spins
    for var in vars
        fvar = Float64(var)
        x_vals, bins = iron_line_profile(
                                        tolerance=tolerance,
                                        size_multiplier=size_multiplier,
                                        fov=fov,
                                        dtmax=dtmax,
                                        obs_angle=fvar,
                                        output="data",
                                        normalised=true,
                                        nbins=nbins
                                        )

        push!(bins_vals, bins)
        push!(x_vals_vals, x_vals)
    end

    new = mapreduce(permutedims, vcat, bins_vals)
    hmp = heatmap(
                new,
                xticks=(LinRange(0,nbins,5),string.(LinRange(0,10,5))),
                # yticks=(LinRange(0,length(vars),5),string.(round.(vars))),
                yticks=(LinRange(0, length(vars), 5), string.(LinRange(vars[1], vars[end], 5))),
                xlabel="Energy (keV)",
                ylabel="Observation Angle",
                cb_title="Flux (Arbitrary Units)"
                )
end

# hmp = iron_line_heatmap(tolerance=1e-8)
# writedlm("heatmap_data.txt", new)
# display(hmp)
# png(hmp, "iron_heatmap.png")


text:
## Load the packages.
using LoopVectorization
using Interpolations
using Statistics
using BenchmarkTools
using PyPlot
using Printf

## Define the upsample functions.
function make_inner_loop!(ex_list, ii, jj, kk, ifac, jfac, kfac, ioff, joff, koff, is_top)
    push!(ex_list, :(i_hi = $ifac*i+$ii+is_hi), :(j_hi = $jfac*j+$jj+js_hi))
    if is_top
        push!(ex_list, :(k_hi = ke_hi))
    else
        push!(ex_list, :(k_hi = $kfac*k+$kk+ks_hi))
    end

    i_pos = -1/2 + ioff + (1/2 - ioff + ii)/ifac
    j_pos = -1/2 + joff + (1/2 - joff + jj)/jfac
    k_pos = -1/2 + koff + (1/2 - koff + kk)/kfac

    i_lo_off = floor(Int, i_pos)
    j_lo_off = floor(Int, j_pos)
    k_lo_off = floor(Int, k_pos)

    if is_top
        push!(ex_list, :(i_lo = i + $i_lo_off + is_lo), :(j_lo = j + $j_lo_off + js_lo), :(k_lo = ke_lo))
    else
        push!(ex_list, :(i_lo = i + $i_lo_off + is_lo), :(j_lo = j + $j_lo_off + js_lo), :(k_lo = k + $k_lo_off + ks_lo))
    end

    fi = mod(i_pos, 1)
    fj = mod(j_pos, 1)
    fk = mod(k_pos, 1)

    coef = zeros(2, 2, 2)
    coef[1, 1, 1] = (((1-fi)*fj+fi-1)*fk+(fi-1)*fj-fi+1)
    coef[2, 1, 1] = ((fi*fj-fi)*fk-fi*fj+fi)
    coef[1, 2, 1] = ((fi-1)*fj*fk+(1-fi)*fj)
    coef[2, 2, 1] = (fi*fj-fi*fj*fk)
    coef[1, 1, 2] = ((fi-1)*fj-fi+1)*fk
    coef[2, 1, 2] = (fi-fi*fj)*fk
    coef[1, 2, 2] = (1-fi)*fj*fk
    coef[2, 2, 2] = fi*fj*fk

    kc_end = is_top ? 1 : 2

    rhs_list = []
    for kc in 1:kc_end, jc in 1:2, ic in 1:2
        if !(coef[ic, jc, kc] ≈ 0)
            push!(rhs_list, :( $(coef[ic, jc, kc])*lo[i_lo+$(ic-1), j_lo+$(jc-1), k_lo+$(kc-1)] ))
        end
    end

    rhs = Expr(:call, :+, rhs_list...)

    push!(ex_list, :(hi[i_hi, j_hi, k_hi] = $rhs))
end

macro upsample_lin_fast(suffix, ifac, jfac, kfac, ioff, joff, koff, add_top)
    ex_inner = []
    for kk in 0:kfac-1, jj in 0:jfac-1, ii in 0:ifac-1
        make_inner_loop!(ex_inner, ii, jj, kk, ifac, jfac, kfac, ioff, joff, koff, false)
    end
    ex_inner_block = quote
        Threads.@threads for k in 0:ktot-1
            for j in 0:jtot-1
                @inbounds @simd for i in 0:itot-1
                    $(Expr(:block, ex_inner...))
                end
            end
        end
    end

    if add_top
        ex_top = []
        kk = 0
        for jj in 0:jfac-1, ii in 0:ifac-1
            make_inner_loop!(ex_top, ii, jj, kk, ifac, jfac, kfac, ioff, joff, koff, true)
        end
        ex_top_block = quote
            Threads.@threads for j in 0:jtot-1
                @inbounds @simd for i in 0:itot-1
                    $(Expr(:block, ex_top...))
                end
            end
        end
    else
        ex_top_block = :()
    end

    name = Symbol(@sprintf "upsample_lin_%s!" suffix)
    ex = quote
        function $name(hi, lo, itot, jtot, ktot, is_hi, js_hi, ks_hi, ke_hi, is_lo, js_lo, ks_lo, ke_lo)
            $ex_inner_block
            $ex_top_block
        end
    end
    println(ex)
    return esc(ex)
end

function upsample_lin_ref!(hi, lo, x_hi, y_hi, z_hi, x_lo, y_lo, z_lo, is_hi, js_hi, ks_hi, ie_hi, je_hi, ke_hi)
    interp = interpolate((x_lo, y_lo, z_lo), lo, (Gridded(Linear()), Gridded(Linear()), Gridded(Linear())))
    hi[is_hi:ie_hi, js_hi:je_hi, ks_hi:ke_hi] .= interp(x_hi[is_hi:ie_hi], y_hi[js_hi:je_hi], z_hi[ks_hi:ke_hi])
end

## Set up the grids.
itot_lo = 256; jtot_lo = 192; ktot_lo = 128
itot_hi = 512; jtot_hi = 384; ktot_hi = 256
igc_lo = 1; jgc_lo = 1; kgc_lo = 1
igc_hi = 1; jgc_hi = 1; kgc_hi = 1

# Derive dimensions
ifac = itot_hi÷itot_lo; jfac = jtot_hi÷jtot_lo; kfac = ktot_hi÷ktot_lo;
icells_lo = itot_lo + 2igc_lo; jcells_lo = jtot_lo + 2jgc_lo; kcells_lo = ktot_lo + 2kgc_lo
icells_hi = itot_hi + 2igc_hi; jcells_hi = jtot_hi + 2jgc_hi; kcells_hi = ktot_hi + 2kgc_hi
is_lo = igc_lo+1; js_lo = jgc_lo+1; ks_lo = kgc_lo+1
is_hi = igc_hi+1; js_hi = jgc_hi+1; ks_hi = kgc_hi+1
ie_lo = igc_lo+itot_lo; je_lo = jgc_lo+jtot_lo; ke_lo = kgc_lo+ktot_lo+1
ie_hi = igc_hi+itot_hi; je_hi = jgc_hi+jtot_hi; ke_hi = kgc_hi+ktot_hi+1

w_lo = rand(icells_lo, jcells_lo, kcells_lo)
w_hi_ref = zeros(icells_hi, jcells_hi, kcells_hi)
w_hi_fast = zeros(icells_hi, jcells_hi, kcells_hi)

dx_lo = 1/itot_lo; dy_lo = 1/jtot_lo; dz_lo = 1/ktot_lo
x_lo = collect(range(0.5-igc_lo, length=icells_lo)) .* dx_lo
y_lo = collect(range(0.5-jgc_lo, length=jcells_lo)) .* dy_lo
z_lo = collect(range(0.5-kgc_lo, length=kcells_lo)) .* dz_lo
xh_lo = collect(range(-igc_lo, length=icells_lo)) .* dx_lo
yh_lo = collect(range(-jgc_lo, length=jcells_lo)) .* dy_lo
zh_lo = collect(range(-kgc_lo, length=kcells_lo)) .* dz_lo

dx_hi = 1/itot_hi; dy_hi = 1/jtot_hi; dz_hi = 1/ktot_hi
x_hi = collect(range(0.5-igc_hi, length=icells_hi)) .* dx_hi
y_hi = collect(range(0.5-jgc_hi, length=jcells_hi)) .* dy_hi
z_hi = collect(range(0.5-kgc_hi, length=kcells_hi)) .* dz_hi
xh_hi = collect(range(-igc_hi, length=icells_hi)) .* dx_hi
yh_hi = collect(range(-jgc_hi, length=jcells_hi)) .* dy_hi
zh_hi = collect(range(-kgc_hi, length=kcells_hi)) .* dz_hi

## Compute a reference using Interpolations.jl
@upsample_lin_fast("222", 2, 2, 2, 0, 0, 0.5, true)
@btime upsample_lin_ref!(
    w_hi_ref, w_lo, x_hi, y_hi, zh_hi, x_lo, y_lo, zh_lo,
    is_hi, js_hi, ks_hi, ie_hi, je_hi, ke_hi)
@btime upsample_lin_222!(
    w_hi_fast, w_lo, itot_lo, jtot_lo, ktot_lo,
    is_hi, js_hi, ks_hi, ke_hi, is_lo, js_lo, ks_lo, ke_lo)

w_lo_int = @view w_lo[is_lo:ie_lo, js_lo:je_lo, ks_lo:ke_lo]
w_hi_ref_int = @view w_hi_ref[is_hi:ie_hi, js_hi:je_hi, ks_hi:ke_hi]
w_hi_fast_int = @view w_hi_fast[is_hi:ie_hi, js_hi:je_hi, ks_hi:ke_hi]

println("Fast equal to ref: ", w_hi_fast_int ≈ w_hi_ref_int)

## Plot the output.
zh_lo_int = @view zh_lo[ks_lo:ke_lo]
zh_hi_int = @view zh_hi[ks_hi:ke_hi]

figure()
plot(zh_hi_int, w_hi_ref_int[1, 1, :], "C1-^")
plot(zh_hi_int, w_hi_fast_int[1, 1, :], "C2-*")
plot(zh_lo_int, w_lo_int[1, 1, :], "k:+")
tight_layout()
display(gcf())
show()


text:
env_config = {
    "dev" => {
        "protocol" => "http",
        "host" => "dev.staffjoy.com",
        "api_key" => "staffjoydev", # This should be from a sudo user
        "sleep" => 10, # Seconds between fetching
        "port" => 8080,
    },
    "stage" => {
        "protocol" => "https",
        "host" => "stage.staffjoy.com",
        "api_key" => get(ENV, "API_KEY", ""),
        "sleep" => get(ENV, "SLEEP", 60),
        "port" => 80,
    },
    "prod" => {
        "protocol" => "https",
        "host" => "www.staffjoy.com",
        "api_key" => get(ENV, "API_KEY", ""),
        "sleep" => get(ENV, "SLEEP", 60),
        "port" => 80,
    },
}

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
module Adapt

  include("stats.jl")
  export l1_norm, l2_norm, lp_norm

end


text:
using FuncPipelines: FixRest

# match utils

struct MatchSplitIterator
    t::Regex
    s::Union{String, SubString}
end
Base.eltype(::Type{MatchSplitIterator}) = Tuple{Bool, SubString}
Base.IteratorSize(::Type{MatchSplitIterator}) = Base.SizeUnknown()

function Base.iterate(itr::MatchSplitIterator, (r, i, e) = (nothing, firstindex(itr.s), lastindex(itr.s)))
    i > e && return nothing
    t, s = itr.t, itr.s
    if !isnothing(r)
        ri, re = first(r), last(r)
        j = isempty(r) ? first(r) : last(r)
        v = (true, SubString(s, ri, re))
        return v, j > e ? (nothing, i, -1) : (nothing, @inbounds(nextind(s, j)), e)
    end

    r = findnext(itr.t, itr.s, i)
    if isnothing(r)
        return (false, SubString(s, i, e)), (nothing, i, -1)
    end

    ri, re = first(r), last(r)
    if i != ri
        return (false, SubString(s, i, @inbounds(prevind(s, ri)))), (r, i, e)
    else
        j = isempty(r) ? first(r) : last(r)
        v = (true, SubString(s, ri, re))
        return v, j > e ? (nothing, i, -1) : (nothing, @inbounds(nextind(s, j)), e)
    end
    nothing
end

matchsplit(t, s) = matchsplit!(Tuple{Bool, SubString}[], t, s)
function matchsplit!(found, t, s)
    i, e = firstindex(s), lastindex(s)

    while true
        r = findnext(t, s, i)
        if isnothing(r)
            push!(found, (false, SubString(s, i, e)))
            break
        end

        ri, re = first(r), last(r)
        i != ri && push!(found, (false, @inbounds SubString(s, i, prevind(s, ri))))
        push!(found, (true, SubString(s, ri, re)))

        j = isempty(r) ? first(r) : last(r)
        j > e && break
        @inbounds i = nextind(s, j)
        i > e && break
    end
    return found
end

function MatchSplits(patterns, x)
    m, ms = first(patterns), @view patterns[2:end]
    sp = MatchSplitIterator(m, x)

    for m in ms
        iters = Iterators.map(sp) do (istoken, s)
            istoken ? ((istoken, s) for _ = 1:1) : MatchSplitIterator(m, s)
        end
        sp = Iterators.Flatten(iters)
    end
    return sp
end

@noinline matchsplits(patterns, x) = collect(Tuple{Bool, SubString}, MatchSplits(patterns, x))

# misc

nestedcall(f) = Base.Fix1(nestedcall, f)
nestedcall(f, x::AbstractArray) = map(nestedcall(f), x)
nestedcall(f, x) = f(x)

function _nestedcall_f!(f, ys, xs)
    @inbounds for i in eachindex(xs, ys)
        ys[i] = nestedcall(f, xs[i])
    end
    return ys
end

function _nestedcall_f_fallback!(f, ys, xs)
    S = Union{}
    @inbounds for i in eachindex(xs, ys)
        ys[i] = nestedcall(f, xs[i])
        S = promote_type(S, typeof(ys[i]))
    end
    return S, ys
end

function nestedcall(f, xs::Array)
    R = Core.Compiler.return_type(nestedcall, Tuple{typeof(f), eltype(xs)})
    if Base.isconcretetype(R)
        return _nestedcall_f!(f, similar(xs, R), xs)
    else
        S, ys = _nestedcall_f_fallback!(f, similar(xs, R), xs)
        if S != R
            zs = similar(xs, S)
            copyto!(zs, ys)
            return zs
        end
        return ys
    end
end

# encode utils

allany(f, x) = mapfoldl(f, _allany, x; init=(true, false))
_allany(a, b) = a[1] & b, a[2] | !b

"""
    with_head_tail(x, head, tail)

Return `[head; x; tail]`. Ignored if `head` or `tail` is `nothing`. `x` can be nested arrays.

# Example

```julia
julia> TextEncodeBase.with_head_tail(1:5, -1, -2)
7-element Vector{Int64}:
 -1
  1
  2
  3
  4
  5
 -2

julia> TextEncodeBase.with_head_tail([1:5, 2:3], -1, -2)
2-element Vector{Vector{Int64}}:
 [-1, 1, 2, 3, 4, 5, -2]
 [-1, 2, 3, -2]

```
"""
with_head_tail(x::AbstractArray, head, tail) = _with_head_tail(x, head, tail)
with_head_tail(x::AbstractArray{<:AbstractVector}, head, tail) = map(with_head_tail(head, tail), x)
function with_head_tail(x::AbstractArray{>:AbstractArray}, head, tail)
    aoa, aov = allany(Base.Fix2(isa, AbstractArray), x)
    if aoa
        map(with_head_tail(head, tail), x)
    elseif aov
        T = mapreduce(typeof, promote_type, x)
        _with_head_tail(T, x, head, tail)
    else
        error("Input array is mixing array and non-array elements")
    end
end

with_head_tail(head, tail) = FixRest(with_head_tail, head, tail)
with_head_tail(x; head=nothing, tail=nothing) = with_head_tail(x, head, tail)
with_head_tail(; head=nothing, tail=nothing) = with_head_tail(head, tail)

@inline function _with_head_tail(::Type{T}, x, head, tail) where T
    S = T
    n = length(x)
    !isnothing(head) && ((n, S) = (n+1, promote_type(S, typeof(head))))
    !isnothing(tail) && ((n, S) = (n+1, promote_type(S, typeof(tail))))
    vec = Vector{S}(undef, n); empty!(vec)
    !isnothing(head) && push!(vec, head)
    append!(vec, x)
    !isnothing(tail) && push!(vec, tail)
    return vec
end
_with_head_tail(x, head, tail) = _with_head_tail(eltype(x), x, head, tail)

"""
    trunc_or_pad(x, n, pad)

Truncate `x` to length `n`, or add `pad` at the end of x until length equal `n`.
 `x` can be either nested or single array. if `n` is `nothing`, the largest length of
 the inner-most array will be used.

    trunc_or_pad(n, pad)

Create a function that will return new array with truncated or padded value of the input.

see also: [trunc_and_pad](@ref)

# Example

```julia
julia> TextEncodeBase.trunc_or_pad(1:5, 7, -1)
7-element Vector{Int64}:
  1
  2
  3
  4
  5
 -1
 -1

julia> TextEncodeBase.trunc_or_pad([1:5, 2:7], 10, -1)
2-element Vector{Vector{Int64}}:
 [1, 2, 3, 4, 5, -1, -1, -1, -1, -1]
 [2, 3, 4, 5, 6, 7, -1, -1, -1, -1]

julia> TextEncodeBase.trunc_or_pad([1:5, [2:7, [1:2]]], nothing, -1)
2-element Vector{Vector}:
 [1, 2, 3, 4, 5, -1]
 Vector[[2, 3, 4, 5, 6, 7], [[1, 2, -1, -1, -1, -1]]]

```
"""
trunc_or_pad(x::AbstractArray, n::Integer, pad) = trunc_or_pad!(similar(x, n), x, n, pad)
trunc_or_pad(x::AbstractArray{<:AbstractArray}, n::Integer, pad) = map(trunc_or_pad(n, pad), x)
function trunc_or_pad(x::AbstractArray{>:AbstractArray}, n::Integer, pad)
    aoa, aov = allany(Base.Fix2(isa, AbstractArray), x)
    if aoa
        map(trunc_or_pad(n, pad), x)
    elseif aov
        trunc_or_pad!(similar(x, n), x, n, pad)
    else
        error("Input array is mixing array and non-array elements")
    end
end
trunc_or_pad(x, ::Nothing, pad) = trunc_or_pad(x, nestedmaxlength(x), pad)

trunc_or_pad(n, pad) = FixRest(trunc_or_pad, n, pad)
trunc_or_pad(x; n=nothing, pad) = trunc_or_pad(x, n, pad)
trunc_or_pad(; n=nothing, pad) = trunc_or_pad(n, pad)

function trunc_or_pad!(vec, x, n, pad)
    if length(x) <= n
        copyto!(vec, x)
        vec[length(x)+1:n] .= pad
    else
        copyto!(vec, 1, x, 1, n)
    end
    return vec
end

"""
    trunc_and_pad(x, maxn, pad)

Truncate `x` if length exceed `maxn`, and add `pad` at the end of x until all length are the same.
 `x` can be either nested or single array. If `maxn` is `nothing`, the largest length of
 the inner-most array will be used, then the behavior equals to `trunc_or_pad` with `nothing`.

    trunc_and_pad(maxn, pad)

Create a function that truncate input to be length <= `maxn`, and add `pad` until all input has equal length.

see also: [trunc_or_pad](@ref)

# Example

```julia
julia> TextEncodeBase.trunc_and_pad(1:5, 7, -1)
5-element Vector{Int64}:
 1
 2
 3
 4
 5

julia> TextEncodeBase.trunc_and_pad([1:5, 2:7], 10, -1)
2-element Vector{Vector{Int64}}:
 [1, 2, 3, 4, 5, -1]
 [2, 3, 4, 5, 6, 7]

julia> TextEncodeBase.trunc_and_pad([1:5, [2:7, [1:2]]], nothing, -1)
2-element Vector{Vector}:
 [1, 2, 3, 4, 5, -1]
 Vector[[2, 3, 4, 5, 6, 7], [[1, 2, -1, -1, -1, -1]]]

```
"""
trunc_and_pad(x, maxn, pad) = (n = nestedmaxlength(x); _trunc_and_pad(x, n, isnothing(maxn) ? n : maxn, pad))

trunc_and_pad(maxn, pad) = FixRest(trunc_and_pad, maxn, pad)
trunc_and_pad(x; maxn=nothing, pad) = trunc_and_pad(x, maxn, pad)
trunc_and_pad(; maxn=nothing, pad) = trunc_and_pad(maxn, pad)

_trunc_and_pad(x, n, maxn, pad) = trunc_or_pad(x, min(n, maxn), pad)

nestedmaxlength(x::AbstractArray{<:AbstractArray}) = mapfoldl(nestedmaxlength, max, x)
nestedmaxlength(x::AbstractArray) = length(x)
function nestedmaxlength(x::AbstractArray{>:AbstractArray})
    aoa, aov = allany(Base.Fix2(isa, AbstractArray), x)
    if aoa
        mapfoldl(nestedmaxlength, max, x)
    elseif aov
        length(x)
    else
        error("Input array is mixing array and non-array elements")
    end
end

_checkeqsize(x, y) = x == y ? x : throw(DimensionMismatch("nested size not the same: $x != $y"))

nestedsize(x::AbstractArray{<:AbstractArray}) = (mapfoldl(nestedsize, _checkeqsize, x)..., size(x)...)
nestedsize(x::AbstractArray) = size(x)
function nestedsize(x::AbstractArray{>:AbstractArray})
    aoa, aov = allany(Base.Fix2(isa, AbstractArray), x)
    if aoa
        (mapfoldl(nestedsize, _checkeqsize, x)..., size(x)...)
    elseif aov
        size(x)
    else
        error("Input array is mixing array and non-array elements")
    end
end

nestedtype(x::AbstractArray{<:AbstractArray}) = mapreduce(nestedtype, promote_type, x)
nestedtype(x::AbstractArray) = mapreduce(typeof, promote_type, x)
function nestedtype(x::AbstractArray{>:AbstractArray})
    aoa, aov = allany(Base.Fix2(isa, AbstractArray), x)
    if aoa
        mapreduce(nestedtype, promote_type, x)
    elseif aov
        mapreduce(typeof, promote_type, x)
    else
        error("Input array is mixing array and non-array elements")
    end
end

"""
    nested2batch(x)

convert nested array into single array

# Example

```julia
julia> TextEncodeBase.nested2batch([[[1 2],[3 4]]])
1×2×2×1 Array{Int64, 4}:
[:, :, 1, 1] =
 1  2

[:, :, 2, 1] =
 3  4

```
"""
function nested2batch(x)
    arr = Array{nestedtype(x)}(undef, nestedsize(x))
    _nested2batch!(arr, 1, x)
    return arr
end

_reduce_nested(dst_offset, xi) = dst_offset[1], _nested2batch!(dst_offset..., xi)[2]
_nested2batch!(arr, offset, x::AbstractArray{<:AbstractArray}) = foldl(_reduce_nested, x; init=(arr, offset))
_nested2batch!(arr, offset, x::AbstractArray) = (copyto!(arr, offset, x, 1, length(x)); (arr, offset+length(x)))
function _nested2batch!(arr, offset, x::AbstractArray{>:AbstractArray})
    aoa, aov = allany(Base.Fix2(isa, AbstractArray), x)
    if aoa
        foldl(_reduce_nested, x; init=(arr, offset))
    elseif aov
        copyto!(arr, offset, x, 1, length(x))
        (arr, offset+length(x))
    else
        error("Input array is mixing array and non-array elements")
    end
end


text:
using Test

using ArnoldiMethod: partialschur, vtype, eigenvalues, SR
using LinearAlgebra

@testset "Zero eigenvalues & low-rank matrices" begin
    # Construct a rank-3 matrix, then with 100% change it'll find the first
    # three eigenvectors in 3 iterations. From then on, all random vectors
    # form a basis for the multiple eigenvalue -- for now we just test that
    # the stopping criterion works for these zero eigenvalues, we don't take
    # into account that it finds the complete eigenspace of the 0 eigenvalue.
    # (we have too few basis vectors for it anyways)
    for T in (Float64, ComplexF64, BigFloat, Complex{BigFloat})
        A = rand(T,10, 3)

        # Rank 3 matrix.
        B = A * A'

        schur, history = partialschur(B, nev = 5, mindim = 5, maxdim = 7, tol = eps())

        @test history.converged
        @test history.mvproducts == 7
        @test norm(schur.Q'schur.Q - I) < 100eps(real(T))
        @test norm(B * schur.Q - schur.Q * schur.R) < 100eps(real(T))
        @test norm(diag(schur.R)[4:5]) < 100eps(real(T))
    end
end

@testset "Right number type" begin
    A = [rand(Bool) ? 1 : 0 for i=1:10, j=1:10]
    @inferred partialschur(A, nev = 2, mindim = 3, maxdim = 8)
    @test vtype(A) == Float64
end

@testset "Find all eigenvalues of a small matrix" begin
    A = rand(3, 3)
    schur, history = partialschur(A)
    @test history.converged
    @test history.mvproducts == 3
end

@testset "Incorrect input" begin
    A = rand(6, 6)

    @test_throws DimensionMismatch partialschur(rand(4, 3))
    @test_throws ArgumentError partialschur(A, mindim = 5, maxdim = 3)
    @test_throws ArgumentError partialschur(A, nev = 5, mindim = 3)
    @test_throws ArgumentError partialschur(A, nev = 5, maxdim = 3)
    @test_throws ArgumentError partialschur(A, nev = 10)
end

@testset "Target non-dominant eigenvalues" begin
    # Dominant eigenvalues 50, 51, 52, 53, but we target the smallest real part
    A = Diagonal([1:0.1:10; 50:53])
    S, hist = partialschur(A, which = SR())
    @test all(x -> real(x) ≤ 10, eigenvalues(S.R))
end


text:
using Base.Test

include("phone-number.jl")

# Tests adapted from `problem-specifications//canonical-data.json` @ v1.2.0
# Returns the cleaned phone number as a digit string if given number is valid,
# else returns `nothing`.

const expected_number = "2234567890"
const valid_10digit_num = (
        "(223) 456-7890",
        "223.456.7890",
        "223 456   7890   ",
)
const valid_11digit_num = (
        "12234567890",
        "  1 223 456 7890 ",
        "+1 (223) 456-7890",
)
const invalid_num = (
        "1223456789",
        "22234567890",
        "321234567890",
        "223-abc-7890",
        "223-@:!-7890",
        "(023) 456-7890",
        "(123) 456-7890",
        "(223) 056-7890",
        "(223) 156-7890",
)

@testset "clean 10-digit number" begin
    @testset "$number" for number in valid_10digit_num
        @test clean(number) == expected_number
    end
end

@testset "clean 11-digit number starting with 1" begin
    @testset "$number" for number in valid_11digit_num
        @test clean(number) == expected_number
    end
end

@testset "detect invalid number" begin
    @testset "$number" for number in invalid_num
        @test clean(number) == nothing
    end
end


text:
export UniquePoints, multiplicities, is_contained, add!, simple_add!, empty!, points, unique_points

const DEFAULT_CAPACITY = Ref(7) # Determined by testing a couple of different values
const NOT_FOUND = -1
const NOT_FOUND_AND_REAL = -2

#############
# SearchBlock
#############
struct SearchBlock{T}
    elements::Vector{Int32}
    children::Vector{Union{Nothing, SearchBlock{T}}}
    capacity::Int
    distances_cache::Vector{Vector{Tuple{T, Int}}} # one per thread
end

function SearchBlock(::Type{T}; capacity = DEFAULT_CAPACITY[]) where T
    children = Vector{Union{Nothing, SearchBlock{T}}}(nothing, capacity)
    distances_cache = [distance_cache(T, capacity) for _=1:Threads.nthreads()]
    SearchBlock(Int32[], children, capacity, distances_cache)
end

function distance_cache(::Type{T}, capacity) where T
    [(typemax(T), i) for i=1:capacity]
end

function SearchBlock(::Type{T}, index::Int; kwargs...) where T
    block = SearchBlock(T; kwargs...)
    push!(block.elements, index)
    block
end

function is_contained(block::SearchBlock{T}, x::NTuple{N,S}, tol::Real, points::Vector, distance::F, threadid=Threads.threadid()) where {T, N, S, F<:Function}
    is_contained(block, SVector{N,S}(x), tol, points, distance)
end

function is_contained(block::SearchBlock{T}, x::AbstractVector{S}, tol::Real, points::Vector, distance::F, threadid=Threads.threadid()) where {T, S, F<:Function}
    if isempty(block.elements)
        return NOT_FOUND
    end

    n = length(block.elements)

    # We compute now the distance to every other element in the block
    # If the distance to one element is smaller than tol, we are done.
    # Otherwise, we need to look into those child whose distance
    # is closest (let's say with distance `d`) and all those
    # children for whose distance `dᵢ` we have |d-dᵢ| < 2*tol
    # This fact can be shown by using the triangle inequality.
    # We therefore simply compute the distances to every point and sort
    # afterwards

    # We go through the elements and compute the distances,
    # keeping track of the three smallest elements
    m₁ = m₂ = m₃ = (typemax(T), 1)
    # we have a distances cache per thread
    distances = block.distances_cache[threadid]
    for i ∈ 1:n
        dᵢ = distance(points[block.elements[i]], x)
        # early exit
        if dᵢ < tol
             # we rely on the distances for look up, so place at the first place the smallest element
            distances[1] = (dᵢ, i)
            return block.elements[i]
        end

        distances[i] = (dᵢ, i)
        # check three smallest elements and update if necessary
        if dᵢ < m₁[1]
            m₃ = m₂
            m₂ = m₁
            m₁ = (dᵢ, i)
        elseif dᵢ < m₂[1]
            m₃ = m₂
            m₂ = (dᵢ, i)
        elseif dᵢ < m₃[1]
            m₃ = (dᵢ, i)
        end
    end

    # Now we computed all distances and also kept track of the two smallest elements
    # Now there are three cases
    # 1) m₁[1] + 2tol < m₂[1] -- > The element can only be in one subtree
    # 2) m₁[1] + 2tol < m₃[1] -- > The element can only be in the first or second subtree
    # 3) else -> The element can also be in more subtrees,
    #            we need to sort the distances vector and look through everything

    # Check smallest element first
    retidx = is_contained(block.children[m₁[2]], x, tol, points, distance, threadid)
    if retidx ≠ NOT_FOUND
        distances[1] = m₁ # we rely on the distances for look up, so place at the first place the smallest element
        return retidx
    end
    # Case 1)
    if m₂[1] - m₁[1] > 2tol
        distances[1] = m₁
        return NOT_FOUND # already checked first tree
    end
    # Case 2) We know m₂[1] - m₁[1] ≤ 2tol
    retidx = is_contained(block.children[m₂[2]], x, tol, points, distance, threadid)
    if retidx ≠ NOT_FOUND
        distances[1] = m₁ # we rely on the distances for look up, so place at the first place the smallest element
        return retidx
    end

    if m₃[1] - m₁[1] > 2tol
        distances[1] = m₁
        return NOT_FOUND # Checked first and second case
    end

    # Since we know als the third element, let's check it
    retidx = is_contained(block.children[m₃[2]], x, tol, points, distance, threadid)
    if retidx ≠ NOT_FOUND
        distances[1] = m₁ # we rely on the distances for look up, so place at the first place the smallest element
        return retidx
    end

    # Case 3)
    # We need to sort distances
    sort!(distances, Base.Sort.InsertionSort, Base.Sort.By(first))

    # We can start at 4 since we already checked the smallest 3
    d = m₃[1]
    for k ∈ 4:n
        dᵢ, i = distances[k]
        if dᵢ - m₁[1] < 2tol
            retidx = is_contained(block.children[i], x, tol, points, distance, threadid)
            if retidx ≠ NOT_FOUND
                return retidx
            end
        else
            break
        end
    end

    return NOT_FOUND
end
is_contained(::Nothing, x::AbstractVector, tol::Real, points::Vector, distance, threadid) = NOT_FOUND

# This assumes that distances_cache is filled
function _insert!(block::SearchBlock{T}, index::Integer, threadid=Threads.threadid()) where {T, V}
    if isempty(block.elements)
        push!(block.elements, index)
        return
    end

    dᵢ, minᵢ = block.distances_cache[threadid][1]
    # if not filled so far, just add it to the current block
    if length(block.elements) < block.capacity
        push!(block.elements, index)
    # we have no children so far, so create a new one
    elseif block.children[minᵢ] === nothing
        block.children[minᵢ] = SearchBlock(T, index; capacity=block.capacity)
    else # a block already exists, so recurse
        _insert!(block.children[minᵢ], index, threadid)
    end
    nothing
end

function Base.empty!(block::SearchBlock)
    empty!(block.elements)
    block.children .= nothing

    nothing
end

#############
# UniquePoints
#############

"""
    UniquePoints{V<:AbstractVector, T, F<:Function}

A data structure which holds points of type `V` where `T=real(eltype(V))`. This data structure
provides an efficient (poly)logarithmic check whether a point already exists where
two points `u,v` are considered equal if `F(u,v)<tol`, where `tol` is a tolerance provided through the [`add!`](@ref) function.

    UniquePoints(v::AbstractVector{<:Number}, distance::F)

Initialize the data structure with just one data point `v`.

    UniquePoints(V::Vector{<:AbstractVector{<:Number}}, distance::F; tol=1e-5)

Initialize the data structure with all points in `v`. These are added in order
by [`add!`](@ref) with the given tolerance `tol`. In particular, 'UniquePoints' structure will contain only points for which the pairwise distance given by `F` is less than `tol`.

    UniquePoints(v; kwargs...) = UniquePoints(v, euclidean_distance; kwargs...)

If `F` is not specialized, [`euclidean_distance`](@ref) is used.

Optional keywords:

* `check_real=true` adds real from points from group orbits (if they exist). The default is `check_real=true`.
* The user can use `group_action=foo` or, if there is more than one group acting, `group_actions=[foo, bar]`. Then, points that are in the same group orbit are considered equal. See [`GroupActions`](@ref) for details regarding the application rules.

## Examples
```julia-repl
julia> points(UniquePoints([[1.0,0.5], [1.0,0.5], [0.5,1.0]]))
2-element Array{Array{Float64,1},1}:
 [1.0, 0.5]
 [0.5, 1.0]

julia> points(UniquePoints([[1.0,0.5], [1.0,0.5], [0.5,1.0]], group_action = x -> [x[2],x[1]]))
1-element Array{Array{Float64,1},1}:
 [1.0, 0.5]
```
"""
struct UniquePoints{V<:AbstractVector, T, F<:Function, MaybeGA<:Union{Nothing, GroupActions}}
    root::SearchBlock{T}
    points::Vector{V}
    distance_function::F
    group_actions::MaybeGA
    check_real::Bool
end

function UniquePoints(::Type{V}, distance::F;
    group_action=nothing,
    group_actions=group_action === nothing ? nothing : GroupActions(group_action),
    check_real::Bool=true) where {T<:Number, V<:AbstractVector{T}, F<:Function}

    if group_actions isa GroupActions
       actions = group_actions
   elseif group_actions === nothing
       actions = nothing
    else
       actions = GroupActions(group_actions)
    end

    root = SearchBlock(real(T))
    points = Vector{V}()
    UniquePoints(root, points, distance, actions, check_real)
end

function UniquePoints(v::AbstractVector{<:Number}, distance::Function; kwargs...)
    data = UniquePoints(typeof(v), distance; kwargs...)
    add!(data, v)
    data
end

function UniquePoints(v::AbstractVector{<:AbstractVector}, distance::F; tol::Real=1e-5, kwargs...) where {F<:Function}
    data = UniquePoints(eltype(v), distance; kwargs...)
    add!(data, v; tol=tol)
end
UniquePoints(v::Type{<:UniquePoints{V}}, distance::F; kwargs...) where {V, F<:Function} = UniquePoints(V, distance; kwargs...)
UniquePoints(v; distance=euclidean_distance, kwargs...) = UniquePoints(v, distance; kwargs...)

function Base.show(io::IO, data::UniquePoints)
    print(io, typeof(data), " with ", length(points(data)), " points")
end
Base.show(io::IO, ::MIME"application/prs.juno.inline", x::UniquePoints) = x

"""
    points(data::UniquePoints)

Return the points stored in `data`.
"""
points(data::UniquePoints) = data.points

Base.getindex(data::UniquePoints, i::Integer) = data.points[i]

"""
    is_contained(data::UniquePoints{V}, x::V; tol=1e-5)::Bool

Check whether `x` is contained in the `data` by using the tolerance `tol` to decide for duplicates.

    is_contained(data::UniquePoints{V}, x::V, Val{true}(); tol=1e-5)::Int

If `x` is contained in `data` by using the tolerance `tol` return the index
of the data point which already exists. If the data point is not existing `-1`
is returned. If `data` has the option `check_real` enabled, a `-2` will be returned once a real vector was added.
"""
function is_contained(data::UniquePoints, x, val=Val{false}(); tol::Real=1e-5)
    is_contained(data, x, val, tol)
end
function is_contained(data::UniquePoints, x::NTuple{N,T}, val, tol::Real) where {N, T}
    is_contained(data, SVector{N,T}(x), val, tol)
end
function is_contained(data::UniquePoints{T}, x::AbstractVector{S}, ::Val{Index}, tol::Real) where {T,S,Index}
    index = is_contained(data.root, x, tol, data.points, data.distance_function)
    if index == NOT_FOUND
        if data.group_actions !== nothing # extra if statement since inference cannot look through &&
            apply_actions(data.group_actions, x) do y
                k = is_contained(data.root, y, tol, data.points, data.distance_function)
                if k ≠ NOT_FOUND
                    index = k
                    return true
                end
                false
            end
        end
    end

    if Index
        return index
    else
        return index ≠ NOT_FOUND
    end
end

"""
    add!(data::UniquePoints{V}, x::V; tol=1e-5)::Bool

Add `x` to `data` if it doesn't already exists by using the tolerance `tol` to decide for duplicates.

    add!(data::UniquePoints{V}, x::V, Val(true); tol=1e-5)::Int

If `x` is contained in `data` by using the tolerance `tol` to decide for duplicates return the index
of the data point which already exists. If the data point is not existing add it to `data` and
return `-1`. If `data` has the option `check_real` enabled, a `-2` will be returned once a real vector was added. The element will be the last element of `points(data)`.
"""
function add!(data::UniquePoints, x::AbstractVector{<:Number}, ::Val{true}; tol::Real=1e-5)
    idx = is_contained(data, x, Val(true), tol)
    idx == NOT_FOUND || return idx

    if !data.check_real
        simple_add!(data, x, tol)
        return NOT_FOUND
    else
        if is_real_vector(x)
            simple_add!(data, x, tol)
            return NOT_FOUND_AND_REAL
        elseif data.group_actions !== nothing
            not_found_and_real = false
            apply_actions(data.group_actions, x) do y
                if is_real_vector(y)
                    simple_add!(data, y, tol)
                    not_found_and_real = true
                    return true
                end
                false
            end
            not_found_and_real && return NOT_FOUND_AND_REAL
        end
        simple_add!(data, x, tol)
        return NOT_FOUND
    end
end
function add!(data::UniquePoints, x::AbstractVector{<:Number}, ::Val{false}=Val(false); tol::Real=1e-5)
    idx = add!(data, x, Val(true); tol=tol)
    idx == NOT_FOUND || idx == NOT_FOUND_AND_REAL
end
function add!(data::UniquePoints, v::AbstractVector{<:AbstractVector}, val::Val=Val(false); tol::Real=1e-5)
    for vᵢ in v
        add!(data, vᵢ; tol=tol)
    end
    data
end

"""
    simple_add!(data::UniquePoints{V}, x::V, tol::Real)::Bool

Similarly to [`add!`](@ref) but does not apply any group actions.
If the data point is not existing add it to `data` and return `-1`. Otherwise
the index of `x` in `data.points` is returned.
"""
function simple_add!(data::UniquePoints, x::AbstractVector, tol::Real)
    idx = is_contained(data.root, x, tol, data.points, data.distance_function)
    if idx ≠ NOT_FOUND
        return idx
    end
    push!(data.points, x)
    _insert!(data.root, length(data.points))
    NOT_FOUND
end

Base.length(data::UniquePoints) = length(data.points)

Base.iterate(data::UniquePoints) = iterate(data.points)
Base.iterate(data::UniquePoints, state) = iterate(data.points, state)
Base.eltype(::Type{<:UniquePoints{V}}) where {V} = V
Base.size(data::UniquePoints) = size(data.points)

"""
    empty!(data::UniquePoints)

Remove all points from `data`.
"""
function Base.empty!(data::UniquePoints)
    empty!(data.root)
    empty!(data.points)
end

"""
    multiplicities(vectors; distance=euclidean_distance, tol::Real = 1e-5, kwargs...)

Returns an array of arrays of integers. Each vector `w` in 'v' contains all indices `i,j` such that `w[i]` and `w[j]` have `distance` at most tol.

Optional keywords:

* `check_real=true` adds real from points from group orbits (if they exist) to the [`UniquePoints`](@ref) data structure used internally. The default is `check_real=false`.
* The user can use `group_action=foo` or, if there is more than one group acting, `group_actions=[foo, bar]`. Then, points that are in the same group orbit are considered equal. See [`GroupActions`](@ref) for details regarding the application rules.

```julia-repl
julia> multiplicities([[1,0.5]; [1,0.5]; [1,1]])
[[1,2]]
```
This is the same as
```julia
multiplicities([[1,0.5]; [1,0.5]; [1,1]]; distance=(x,y) -> LinearAlgebra.norm(x-y))
```
Here is an example for using group actions.
```julia-repl
julia> X = [[1, 2, 3, 4], [2,1,3,4], [1,2,4,3], [2,1,4,3]]
julia> permutation(x) = [x[2], x[1], x[3], x[4]]
julia> m = multiplicities(X, group_action = permutation)
[[1,2], [3,4]]
```
"""
multiplicities(v; kwargs...) = multiplicities(identity, v; kwargs...)
function multiplicities(f, v; kwargs...) where {F<:Function}
    isempty(v) && return Vector{Vector{Int}}()
    _multiplicities(f, v; kwargs...)
end
default_distance(f, v) = isa(f(first(v)), PVector) ? fubini_study : euclidean_distance
function _multiplicities(f, v; distance=default_distance(f,v), tol::Float64=1e-5, check_real=false, kwargs...) where {F<:Function}
    mults = Dict{Int, Vector{Int}}()
    positions = Vector{Int32}()
    x₀ = f(first(v))
    T = typeof(similar(x₀, promote_type(eltype(x₀), Float64)))
    data = UniquePoints(T, distance; check_real=check_real, kwargs...)
    for (i, vᵢ) in enumerate(v)
        k = add!(data, f(vᵢ), Val(true); tol = tol)
        if k > 0
            if haskey(mults, positions[k])
                push!(mults[positions[k]], i)
            else
                mults[positions[k]] = [positions[k], i]
            end
        else
            push!(positions, i)
        end
    end
    collect(values(mults))
end

"""
    unique_points(points::AbstractVector{<:AbstractVector}; options...)

Compute all unique points with respect to the given options. See [`UniquePoints`](@ref)
for possible options. In particular, it is possible to pass group actions.

## Example

```julia
julia> pts = [[1.0,0.5], [1.0,0.5], [0.5,1.0]];
julia> unique_points(pts)
2-element Array{Array{Float64,1},1}:
 [1.0, 0.5]
 [0.5, 1.0]

julia> unique_points(pts; group_action = x -> [x[2],x[1]])
1-element Array{Array{Float64,1},1}:
 [1.0, 0.5]
```
"""
unique_points(v::AbstractVector{<:AbstractVector}; kwargs...) = points(UniquePoints(v; kwargs...))


text:
using LinearAlgebra
using Plots
using Revise

includet("IterativeMethods.jl")
includet("InternalArray.jl")

# Definitions -----------------------------------------------------------------------------

# We take the reference system to be the position of the spring in rest.
# That is, for the one dimensional case: (l, 0).

"""
    SpringProblem

Contains the parameters associated with a typical spring problem.

# Fields
- `t0`: Initial time of the problem (defaults to 0)
- `p0`: Initial position
- `x0`: Initial speed

- `m`: Mass of the object
- `k`: Spring's constant
- `b`: Friction constant
- `A`: External force amplitude
- `ω`: External force angular speed
"""
struct SpringProblem
    t0
    p0
    v0

    m
    k
    b
    A
    ω

    SpringProblem(; t0 = 0, p0, v0, m, k, b = 0, A = 0, ω = 0) = new(t0, p0, v0, m, k, b, A, ω)
end

function eulermethod(prob::SpringProblem, step, nit)
    # We ought to transform the equation to a  differential equation of order one,
    # for which we write x = [position; speed].
    # That is, the variable x will represent the solution rather than the position.
    # And the solution contains the position and the speed.
    n = length(prob.p0)
    function f(t, x)
        p = x[1:n]
        v = x[n+1:2*n]
        return [v; (-prob.k*p -prob.b*v + [prob.A*sin(prob.ω*t)])/prob.m]
    end

    return eulermethod(f, prob.t0, [prob.p0; prob.v0], h = step, iter = nit)
end

function eulermethod_by_approximation(prob::SpringProblem, timelen; error = 1e-3)
    n = length(prob.p0)
    function f(t, x)
        p = x[1:n]
        v = x[n+1:2*n]
        return [v; (-prob.k*p -prob.b*v + [prob.A*sin(prob.ω*t)])/prob.m]
    end

    return eulermethod_by_approximation(f, prob.t0, timelen, [prob.p0; prob.v0]; error = error)
end

max_stars_repo_name: useredsa/numeric-differential-equations

max_issues_repo_name: useredsa/numeric-differential-equations

max_forks_repo_name: useredsa/numeric-differential-equations


text:
using ITensors
using ITensorPySCF
using ITensorGaussianMPS
using LinearAlgebra

mol = gto.M()
N = 10
mol.nelectron = N

U = 4.0
t = 1.0

mf = scf.RHF(mol)
h1 = zeros(N, N)
for i in 1:(N - 1)
  h1[i, i + 1] = h1[i + 1, i] = -t
end
#h1[N, 1] = h1[1, N] = -t  # PBC
eri = zeros(N, N, N, N)
for i in 1:N
  eri[i, i, i, i] = U
end

mf.get_hcore = (_...) -> h1
mf.get_ovlp = (_...) -> Matrix(I(N))
# ao2mo.restore(8, eri, N) to get 8-fold permutation symmetry of the integrals
# ._eri only supports the two-electron integrals in 4-fold or 8-fold symmetry.
mf._eri = ao2mo.restore(8, eri, N)

mf.kernel()

println("Molecular orbital coefficients:")
display(mf.mo_coeff)

Φ = mf.mo_coeff[1:N, 1:(N ÷ 2)]

s = siteinds("Electron", N; conserve_qns=true)
println("Making free fermion starting MPS")
@time ψ0 = slater_determinant_to_mps(
  s, Φ, Φ; eigval_cutoff=1e-5, cutoff=1e-8, maxdim=200
)
@show maxlinkdim(ψ0)

os = OpSum()
for n in 1:(N - 1)
  os .+= -t, "Cdagup", n, "Cup", n + 1
  os .+= -t, "Cdagdn", n, "Cdn", n + 1
  os .+= -t, "Cdagup", n + 1, "Cup", n
  os .+= -t, "Cdagdn", n + 1, "Cdn", n
end
for n in 1:N
  os .+= U, "Nupdn", n
end
H = MPO(os, s)

println("\nFree fermion starting energy")
@show inner(ψ0, H, ψ0)

println("\nRun dmrg with free starting state")
sweeps = Sweeps(10)
maxdim!(sweeps,10,20,40,60)
cutoff!(sweeps,1E-12)
@time dmrg(H, ψ0, sweeps)


text:
"""
    save_tetgen(mesh_all, filename)

Save mesh in the Tetgen output format.
"""
function save_tetgen(mesh_all, filename)

    # Extract mesh
    @unpack points, facets, elements, facetmarkers, elementmarkers = mesh_all

    println("Saving mesh in Tetgen output format at " * filename)

    # Save points
    npoint = size(points, 2)
    open(filename * ".node", "w") do io
        write(io, "$npoint 3 0 0\n")
        for ipoint = 1:npoint
            write(io, @sprintf "%d %f %f %f\n" ipoint points[:, ipoint]...)
        end
    end

    # Read facets and their associated boundaries
    nfacet = size(facets, 2)
    open(filename * ".face", "w") do io
        write(io, "$nfacet 1\n")
        for i = 1:nfacet
            write(io, @sprintf "%d %d %d %d %d\n" i facets[:, i]... facetmarkers[i])
        end
    end

    # Read elements and their associated compartments
    nelement = size(elements, 2)
    open(filename * ".ele", "w") do io
        write(io, "$nelement 4 1\n")
        for i = 1:nelement
            write(io, @sprintf "%d %d %d %d %d %d\n" i elements[:, i]... elementmarkers[i])
        end
    end

end


text:
using OrdinaryDiffEq, Test

f(u,p,t) = 0.98u
u0 = 1.0
tspan = (0.0, 1.0)
prob=ODEProblem(f,u0,tspan)
sol = solve(prob,Tsit5())

# Make sure various differentiation forms work on scalars

sol1 = solve(prob,Rosenbrock23(),abstol=1e-12,reltol=1e-12)
sol2 = solve(prob,Rosenbrock23(autodiff=false),abstol=1e-12,reltol=1e-12)
sol3 = solve(prob,Rosenbrock23(autodiff=false,diff_type=Val{:central}),abstol=1e-12,reltol=1e-12)
sol4 = solve(prob,Rosenbrock23(autodiff=false,diff_type=Val{:complex}),abstol=1e-12,reltol=1e-12)
sol5 = solve(prob,KenCarp4(),abstol=1e-12,reltol=1e-12)
sol6 = solve(prob,KenCarp4(autodiff=false),abstol=1e-12,reltol=1e-12)
sol7 = solve(prob,KenCarp4(autodiff=false,diff_type=Val{:central}),abstol=1e-12,reltol=1e-12)
sol8 = solve(prob,KenCarp4(autodiff=false,diff_type=Val{:complex}),abstol=1e-12,reltol=1e-12)

ts = 0.0:0.1:1.0
@test sol1(ts) ≈ sol2(ts)
@test sol1(ts) ≈ sol3(ts)
@test sol1(ts) ≈ sol4(ts)
@test sol1(ts) ≈ sol5(ts)
@test sol5(ts) ≈ sol6(ts)
@test sol5(ts) ≈ sol7(ts)
@test sol5(ts) ≈ sol8(ts)

# Test that in-place and out-of-place gets the same results

function f_ip(du, u, p, t)
    du[1] = - u[1]
    nothing
end
f_scalar(u, p, t) = -u

prob_ip = ODEProblem(f_ip, [1.0], (0.0, 10.0))
prob_scalar = ODEProblem(f_scalar, 1.0, (0.0, 10.0))
ts = 0:0.1:10.0

rk_algs = [Euler(),Midpoint(),Heun(),Ralston(),RK4(),SSPRK104(),SSPRK22(),SSPRK33(),
        SSPRK432(),BS3(),BS5(),DP5(),DP8(),Feagin10(),Feagin12(),
        Feagin14(),TanYam7(),Tsit5(),TsitPap8(),Vern6(),Vern7(),Vern8(),Vern9()]

@testset "Algorithm $(nameof(typeof(alg)))" for alg in rk_algs
  println(nameof(typeof(alg)))
  sol_ip = solve(prob_ip, alg, dt = 0.1)
  sol_scalar = solve(prob_scalar, alg, dt = 0.1)

  @test sol_ip(ts, idxs=1) ≈ sol_scalar(ts)
  @test sol_ip.t ≈ sol_scalar.t && sol_ip[1, :] ≈ sol_scalar.u
end

working_sdirk_algs = [ImplicitMidpoint(),
                      ImplicitEuler(),
                      SSPSDIRK2()]

sdirk_algs = [Trapezoid(),
              TRBDF2(), SDIRK2(),
              Kvaerno3(), KenCarp3(),
              Cash4(), Hairer4(), Hairer42(), Kvaerno4(), KenCarp4(),
              Kvaerno5(), KenCarp5()]

@testset "Algorithm $(nameof(typeof(alg)))" for alg in working_sdirk_algs
  println(nameof(typeof(alg)))
  sol_ip = solve(prob_ip, alg, dt=0.0125)
  sol_scalar = solve(prob_scalar, alg, dt=0.0125)

  @test sol_ip(ts, idxs=1) ≈ sol_scalar(ts)
  @test sol_ip.t ≈ sol_scalar.t && sol_ip[1, :] ≈ sol_scalar.u
end

@testset "Algorithm $(nameof(typeof(alg)))" for alg in sdirk_algs
  println(nameof(typeof(alg)))
  sol_ip = solve(prob_ip, alg, dt=0.0125)
  sol_scalar = solve(prob_scalar, alg, dt=0.0125)

  @test_broken sol_ip(ts, idxs=1) ≈ sol_scalar(ts)
  @test_broken sol_ip.t ≈ sol_scalar.t && sol_ip[1, :] ≈ sol_scalar.u
end

working_rosenbrock_algs = [Rosenbrock23(),ROS3P(),Rodas3(),
                           RosShamp4(), Veldd4(), Velds4(), GRK4T(), GRK4A(),
                           Ros4LStab(), Rodas4(), Rodas42(), Rodas4P(), Rodas5()]

rosenbrock_algs = [Rosenbrock32(),
              ]

@testset "Algorithm $(nameof(typeof(alg)))" for alg in working_rosenbrock_algs
  println(nameof(typeof(alg)))
  sol_ip = solve(prob_ip, alg, dt=0.0125)
  sol_scalar = solve(prob_scalar, alg, dt=0.0125)

  @test sol_ip(ts, idxs=1) ≈ sol_scalar(ts)
  @test sol_ip.t ≈ sol_scalar.t && sol_ip[1, :] ≈ sol_scalar.u
end

@testset "Algorithm $(nameof(typeof(alg)))" for alg in rosenbrock_algs
  println(nameof(typeof(alg)))
  sol_ip = solve(prob_ip, alg, dt=0.0125)
  sol_scalar = solve(prob_scalar, alg, dt=0.0125)

  @test_broken sol_ip(ts, idxs=1) ≈ sol_scalar(ts)
  @test_broken sol_ip.t ≈ sol_scalar.t && sol_ip[1, :] ≈ sol_scalar.u
end

rkc_algs = [RKC(),ROCK2(),ROCK4(),SERK2()]

@testset "Algorithm $(nameof(typeof(alg)))" for alg in rkc_algs
  println(nameof(typeof(alg)))
  sol_ip = solve(prob_ip, alg)
  sol_scalar = solve(prob_scalar, alg)

  @test sol_ip(ts, idxs=1) ≈ sol_scalar(ts)
  @test sol_ip.t ≈ sol_scalar.t && sol_ip[1, :] ≈ sol_scalar.u
end


text:
using BinaryBuilder, Pkg
using Base.BinaryPlatforms
const YGGDRASIL_DIR = "../.."
include(joinpath(YGGDRASIL_DIR, "platforms", "mpi.jl"))

name = "MPICH"
version_str = "4.0.2"
version = VersionNumber(version_str)

sources = [
    ArchiveSource("https://www.mpich.org/static/downloads/$(version_str)/mpich-$(version_str).tar.gz",
                  "5a42f1a889d4a2d996c26e48cbf9c595cbf4316c6814f7c181e3320d21dedd42"),
]

script = raw"""
################################################################################
# Install MPICH
################################################################################

# Enter the funzone
cd ${WORKSPACE}/srcdir/mpich*

EXTRA_FLAGS=()
# Define some obscure undocumented variables needed for cross compilation of
# the Fortran bindings.  See for example
# * https://stackoverflow.com/q/56759636/2442087
# * https://github.com/pmodels/mpich/blob/d10400d7a8238dc3c8464184238202ecacfb53c7/doc/installguide/cfile
export CROSS_F77_SIZEOF_INTEGER=4
export CROSS_F77_SIZEOF_REAL=4
export CROSS_F77_SIZEOF_DOUBLE_PRECISION=8
export CROSS_F77_TRUE_VALUE=1
export CROSS_F77_FALSE_VALUE=0

if [[ ${nbits} == 32 ]]; then
    export CROSS_F90_ADDRESS_KIND=4
else
    export CROSS_F90_ADDRESS_KIND=8
fi
export CROSS_F90_OFFSET_KIND=8
export CROSS_F90_INTEGER_KIND=4
export CROSS_F90_INTEGER_MODEL=9
export CROSS_F90_REAL_MODEL=6,37
export CROSS_F90_DOUBLE_MODEL=15,307
export CROSS_F90_ALL_INTEGER_MODELS=2,1,4,2,9,4,18,8,
export CROSS_F90_INTEGER_MODEL_MAP={2,1,1},{4,2,2},{9,4,4},{18,8,8},

if [[ "${target}" == i686-linux-musl ]]; then
    # Our `i686-linux-musl` platform is a bit rotten: it can run C programs,
    # but not C++ or Fortran.  `configure` runs a C program to determine
    # whether it's cross-compiling or not, but when it comes to running
    # Fortran programs, it fails.  In addition, `configure` ignores the
    # above exported variables if it believes it's doing a native build.
    # Small hack: edit `configure` script to force `cross_compiling` to be
    # always "yes".
    sed -i 's/cross_compiling=no/cross_compiling=yes/g' configure
    EXTRA_FLAGS+=(ac_cv_sizeof_bool="1")
fi

if [[ "${target}" == aarch64-apple-* ]]; then
    EXTRA_FLAGS+=(
        FFLAGS=-fallow-argument-mismatch
        FCFLAGS=-fallow-argument-mismatch
    )
fi

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --enable-shared=yes --enable-static=no \
    --with-device=ch3 --disable-dependency-tracking \
    --enable-fast=all,O3 \
    --docdir=/tmp \
    --disable-opencl \
    "${EXTRA_FLAGS[@]}"

# Remove empty `-l` flags from libtool
# (Why are they there? They should not be.)
# Run the command several times to handle multiple (overlapping) occurrences.
sed -i 's/"-l /"/g;s/ -l / /g;s/-l"/"/g' libtool
sed -i 's/"-l /"/g;s/ -l / /g;s/-l"/"/g' libtool
sed -i 's/"-l /"/g;s/ -l / /g;s/-l"/"/g' libtool

# Build the library
make -j${nproc}

# Install the library
make install

################################################################################
# Install licenses
################################################################################

install_license $WORKSPACE/srcdir/mpich*/COPYRIGHT
"""

augment_platform_block = """
    using Base.BinaryPlatforms
    $(MPI.augment)
    augment_platform!(platform::Platform) = augment_mpi!(platform)
"""

platforms = expand_gfortran_versions(filter!(!Sys.iswindows, supported_platforms(; experimental=true)))

# Add `mpi+mpich` platform tag
foreach(p -> (p["mpi"] = "MPICH"), platforms)

products = [
    # MPICH
    LibraryProduct("libmpicxx", :libmpicxx),
    LibraryProduct("libmpifort", :libmpifort),
    LibraryProduct("libmpi", :libmpi),
    ExecutableProduct("mpiexec", :mpiexec),
]

dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"), v"0.5.2"),
    Dependency(PackageSpec(name="MPIPreferences", uuid="3da0fdf6-3ccc-4f1b-acd9-58baa6c99267"); compat="0.1"),
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               augment_platform_block, julia_compat="1.6")

max_stars_repo_name: staticfloat/BinaryBuilderRepos

max_issues_repo_name: staticfloat/BinaryBuilderRepos

max_forks_repo_name: staticfloat/BinaryBuilderRepos


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

using Gadfly
using Color

L(R) = x -> R * x * (1 - x)

k = 60
Rs = linspace(0.1, 4.0, 10) # [0.1, 0.533, 0.967, 1.4, ...]
xs = map(R -> iterate(L(R), 0.2, k), Rs)

fps = map(x -> fixedpoint(x, ϵ = 1e-3), xs) # [4, 8, 26, 10, ...]

# Values of R ∈ [3.0, 4.0] have very interesting properties.
# See Wikipedia or explore them yourself.
# Doing a Fourier Transform of the sequence will help illuminate these interesting
# properties


text:
using MATLAB, Base.Test

restart_default_msession()

a = [1. 2. 3.; 4. 5. 6.]
b = [2. 3. 4.; 8. 7. 6.]

@mput a b
@matlab begin
    r1 = a .* b
    r2 = a + b
end
@mget r1 r2

@test isequal(r1, a .* b)
@test isequal(r2, a + b)

@mget r1::Vector
@test isequal(r1, vec(a .* b))

s = sparse([1. 0. 0.; 2. 3. 0.; 0. 4. 5.])
put_variable(:s, s)
s2 = get_variable(:s)
@test isequal(s, s2)

# mxcall

r = mxcall(:plus, 1, a, b)
@test isequal(r, a + b)

(xx, yy) = mxcall(:meshgrid, 2, [1., 2.], [3., 4.])
@test isequal(xx, [1. 2.; 1. 2.])
@test isequal(yy, [3. 3.; 4. 4.])

close_default_msession()