text:
"""
    calculate_distance (lon1, lat1, lon2, lat2)

Returns the distance between coordinates (lon1, lat1) and (lon2, lat2) in
units of [km].

# Arguments
- 'lon1::Number': Longitude of point 1
- 'lat1::Number': Latitude of point 1
- 'lon2::Number': Longitude of point 2
- 'lat2::Number': Latitude of point 2

# Returns
- 'Number': Distance between (lon1, lat1) and (lon2, lat2) in kilometers.
"""
function calculate_distance(lon1, lat1, lon2, lat2)

    # Perform in-range checks
    if (lon1 < -180) | (lon1 > 180)
        throw(DomainError(lon1, "lon1 must be [-180 180]"))
    end

    if (lon2 < -180) | (lon2 > 180)
        throw(DomainError(lon2, "lon2 must be [-180 180]"))
    end

    if (lat1 < -90) | (lat1 > 90)
        throw(DomainError(lat1, "lat1 must be [-90 90]"))
    end

    if (lat2 < -90) | (lat2 > 90)
        throw(DomainError(lat2, "lat2 must be [-90 90]"))
    end

    deltalon = deg2rad(lon2 - lon1)

    lat1_rad = deg2rad(lat1)
    lat2_rad = deg2rad(lat2)
    deltalat = lat2_rad - lat1_rad

    a = sin(deltalat / 2.0)^2 + cos(lat1_rad) * cos(lat2_rad) * sin(deltalon / 2.0)^2

    return 2.0 * 6378.0 * asin(sqrt(a))

end

"""
    mask_locations_within_radius (lon, lat, radius, locations)

Takes an array of locations {Geolocation} along with a center lon and
center lat coordinate, a radius (in [km]) and returns a Bool array which
represents whether a location is within the radius or not.

# Arguments
- 'lon::Number': Longitude of central coordinate
- 'lat::Number': Latitude of central coordinate
- 'radius::Number': Radius [km] within which points will be kept
- 'locations::Array{Geolocation, 1}': Array of Geolocations to be masked

# Returns
- 'Array{Bool, 1}': Boolean mask array

"""
function check_location_within_radius(lon, lat, radius, location::Geolocation)

    dist = calculate_distance(
            lon, lat,
            location.lon,
            location.lat
        )

    return dist <= radius

end

function check_location_within_radius(lon, lat, radius, loctime::GeolocationTime)

    dist = calculate_distance(
        lon, lat,
        loctime.loc.lon,
        loctime.loc.lat
    )

    return dist <= radius

end

"""

"""
function calculate_solar_angles(loctime::GeolocationTime)

    # Convert to Julian date
    jdate = jdcnv(loctime.time)
    # Obtain solar position
    sun_radec = sunpos(jdate)
    # Get the local position
    loc_radec = zenpos(jdate, loctime.loc.lat, loctime.loc.lon)
    # Compute apparent ALT and AZ of the sun from (lon, lat) position
    sun_alt, sun_az, sun_hour = eq2hor(sun_radec[1], sun_radec[2], jdate, loctime.loc.lat, loctime.loc.lon)

    if sun_alt < 0
        return NaN, sun_az
    else
        # Return SZA and SAA in degrees
        return 90.0 - sun_alt, sun_az
    end

end

"""
    Function to return bin indices in one dimension, taken from:
https://stackoverflow.com/questions/54879412/get-the-mapping-from-each-element-of-input-to-the-bin-of-the-histogram-in-julia

"""
binindices(edges, data) = searchsortedlast.(Ref(edges), data)

"""

Function to produce a 2D histogram which also returns the indices
of the 2D grid in which certain coordinates fall into (stdlib Julia functions
don't return the indices sadly, so we built our own)

"""
function calculate_regular_2d_histogram(lons::Array{<:Number, 1},
                                        lats::Array{<:Number, 1},
                                        delta_lon,
                                        delta_lat)

    # Construct the lon and lat boundaries
    # ------------------------------------

    # The grids are constructed such that the space between first
    # bin edge and the smallest data point is the same as the
    # last bin edge and the largest data point.

    Nlon = convert(Int, ceil((maximum(lons) - minimum(lons)) / delta_lon))
    Nlat = convert(Int, ceil((maximum(lats) - minimum(lats)) / delta_lat))

    lon_bin0 = (minimum(lons) + maximum(lons) - (Nlon) * delta_lon) / 2
    lat_bin0 = (minimum(lats) + maximum(lats) - (Nlat) * delta_lat) / 2

    lon_grid = collect(range(lon_bin0, step=delta_lon, length=Nlon+1))
    lat_grid = collect(range(lat_bin0, step=delta_lat, length=Nlat+1))

    # Assign to bin indices
    lon_idx = binindices(lon_grid, lons)
    lat_idx = binindices(lat_grid, lats)

    cell_idx = collect(zip(lon_idx, lat_idx))

    return (lon_grid, lat_grid), cell_idx

end

function subset_scene_time(IS::InstrumentSampling,
                           date_min::DateTime,
                           date_max::DateTime)

    T = typeof(IS)

    # Understand which scene time fall into the time window
    idx_times = findall(
        (get_time(IS) .>= date_min) .&
        (get_time(IS) .<= date_max)
    )

    # Which of those resulting locations drop out of the locations list?
    newloc = intersect(
        get_loc(IS),
        (p -> p.loctime.loc).(IS.scenes[idx_times])
    )

    return T(
        IS.info,
        IS.instrument,
        newloc,
        IS.scenes[idx_times]
    )

end

ext: jl


text:
# This file was generated by JuDoc, do not modify it. # hide
evaluate(pipe, X, height, resampling=Holdout(),
         measure=rms) |> pprint

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
# using MPI
import MPI

MPI.Init()
rank = MPI.Comm_rank(MPI.COMM_WORLD)
size = MPI.Comm_size(MPI.COMM_WORLD)

manager = MPI.start_main_loop(MPI.MPI_TRANSPORT_ALL)

println("Hello from rank $rank")
println("Size $size")

@everywhere function print_cpus()
    run(pipeline(`ps -eo pid,psr,pcpu,cmd`, `grep julia`))
    run(`hostname`)
    # run(`echo $LD_LIBRARY_PATH`)
end

println("Before map")
print_cpus()

@everywhere function pcpus(mat)
    println("In map")
    print_cpus()
    svd(mat)
end

M = [rand(100,100) for i=1:size]

pmap(pcpus, M)

println("After map.")
print_cpus()

MPI.stop_main_loop(manager)


text:
using SpikingNetworks

parseParameters"""
N=100
τ = 1
W=0.1*randn(N,N)"""

parseEquations"""
dr/dt = (Iₑ(t)-r)/τ+max(W⊙r,0)
Iₑ(t)=0.1"""

net=init"""
N=N
r=zeros(N)"""

stateMonitor"[1:30;]"

@time run!(net,20/net.dt)
using MatlabPlot
mplot(net.rRecord')


text:
"""
    frst(signal, α, p)

Computing the α order fractional sine transform of the input **signal**.

# Example

```julia-repl
julia> frst([1,2,3], 0.5, 2)
3-element Vector{ComplexF64}:
1.707106781186548 + 1.207106781186547im
1.9999999999999998 - 1.7071067811865481im
-1.1213203435596437 - 1.2071067811865468im
```

### References

```tex
@article{article,
author = {Pei, Soo-Chang and Yeh, Min-Hung},
year = {2001},
month = {07},
pages = {1198 - 1207},
title = {The discrete fractional cosine and sine transforms},
volume = {49},
journal = {Signal Processing, IEEE Transactions on},
doi = {10.1109/78.923302}
}
```
"""
function frst(signal, α, p)
    N = length(signal)
    @views signal = signal[:]
    p = min(max(2, p), N-1)
    E = dFRST(N,p)
    result = E *(exp.(-im*pi*α*collect(0:N-1)) .*(E' *signal))
    return result
end

function dFRST(N, p)
    N1 = 2*N+2
    d2 = [1, -2, 1]

    d_p = 1
    s = 0
    st = zeros(1, N1)

    for k = 1:floor(Int, p/2)
        if isa(d_p, Number)
            d_p = @. d2*d_p
        else
            d_p = conv(d2, d_p)
        end
        st[vcat(collect(N1-k+1:N1), collect(1:k+1))] = d_p
        st[1]=0
        temp=vcat(union(1, collect(1:k-1)), union(2,collect(1:k-1)))
        temp = temp[:] ./collect(1:2*k)
        s=s.+(-1)^(k-1)*prod(temp)*2*st
    end

    H = Toeplitz(s[:], s[:]) +diagm(real.(fft(s[:])))

    V = hcat(zeros(N), reverse(zeros(N, N)+I, dims=1), zeros(N), -(zeros(N, N)+I)) ./sqrt(2)

    Od = V*H*V'

    _, vo = eigen(Od)

    return reverse(reverse(vo, dims=2), dims=1)
end


text:
"""
    TruncSVD

ITensor factorization type for a truncated singular-value
decomposition, returned by `svd`.
"""
struct TruncSVD
  U::ITensor
  S::ITensor
  V::ITensor
  spec::Spectrum
  u::Index
  v::Index
end

# iteration for destructuring into components `U,S,V,spec,u,v = S`
iterate(S::TruncSVD) = (S.U, Val(:S))
iterate(S::TruncSVD, ::Val{:S}) = (S.S, Val(:V))
iterate(S::TruncSVD, ::Val{:V}) = (S.V, Val(:spec))
iterate(S::TruncSVD, ::Val{:spec}) = (S.spec, Val(:u))
iterate(S::TruncSVD, ::Val{:u}) = (S.u, Val(:v))
iterate(S::TruncSVD, ::Val{:v}) = (S.v, Val(:done))
iterate(S::TruncSVD, ::Val{:done}) = nothing

@doc """
    svd(A::ITensor, inds::Index...; <keyword arguments>)

Singular value decomposition (SVD) of an ITensor `A`, computed
by treating the "left indices" provided collectively
as a row index, and the remaining "right indices" as a
column index (matricization of a tensor).

The first three return arguments are `U`, `S`, and `V`, such that
`A ≈ U * S * V`.

Whether or not the SVD performs a trunction depends on the keyword
arguments provided.

# Examples

```julia
i = Index(2)
j = Index(5)
k = Index(2)

A = randomITensor(i, j, k)
U, S, V = svd(A, i, k);
@show norm(A - U * S * V) <= 10 * eps() * norm(A)

# This will truncate the last 2 singular values.
# The norm of the difference with the original tensor
# will be the sqrt root of the sum of the squares of the
# singular values that get truncated.
Utrunc, Strunc, Vtrunc = svd(A, i, k; maxdim=2);
@show norm(A - Utrunc * Strunc * Vtrunc) ≈ sqrt(S[3, 3]^2 + S[4, 4]^2)

# Alternatively we can specify that we want to truncate
# the weights of the singular values up to a certain cutoff,
# so the error will be no larger than the cutoff.
Utrunc2, Strunc2, Vtrunc2 = svd(A, i, k; cutoff=1e-10);
@show norm(A - Utrunc2 * Strunc2 * Vtrunc2) <= 1e-10
```

# Keywords
- `maxdim::Int`: the maximum number of singular values to keep.
- `mindim::Int`: the minimum number of singular values to keep.
- `cutoff::Float64`: set the desired truncation error of the SVD, by default defined as the sum of the squares of the smallest singular values.
- `lefttags::String = "Link,u"`: set the tags of the Index shared by `U` and `S`.
- `righttags::String = "Link,v"`: set the tags of the Index shared by `S` and `V`.
- `alg::String = "divide_and_conquer"`. Options:
  - `"divide_and_conquer"` - A divide-and-conquer algorithm. LAPACK's gesdd. Fast, but may lead to some innacurate singular values for very ill-conditioned matrices. Also may sometimes fail to converge, leading to errors (in which case "qr_iteration" or "recursive" can be tried).
  - `"qr_iteration"` - Typically slower but more accurate for very ill-conditioned matrices compared to `"divide_and_conquer"`. LAPACK's gesvd.
  - `"recursive"` - ITensor's custom svd. Very reliable, but may be slow if high precision is needed. To get an `svd` of a matrix `A`, an eigendecomposition of ``A^{\\dagger} A`` is used to compute `U` and then a `qr` of ``A^{\\dagger} U`` is used to compute `V`. This is performed recursively to compute small singular values.
- `use_absolute_cutoff::Bool = false`: set if all probability weights below the `cutoff` value should be discarded, rather than the sum of discarded weights.
- `use_relative_cutoff::Bool = true`: set if the singular values should be normalized for the sake of truncation.

See also: [`factorize`](@ref), [`eigen`](@ref)
"""
function svd(A::ITensor, Linds...; kwargs...)
  utags::TagSet = get(kwargs, :lefttags, get(kwargs, :utags, "Link,u"))
  vtags::TagSet = get(kwargs, :righttags, get(kwargs, :vtags, "Link,v"))

  # Keyword argument deprecations
  #if haskey(kwargs, :utags) || haskey(kwargs, :vtags)
  #  @warn "Keyword arguments `utags` and `vtags` are deprecated in favor of `leftags` and `righttags`."
  #end

  Lis = commoninds(A, indices(Linds))
  Ris = uniqueinds(A, Lis)

  if length(Lis) == 0 || length(Ris) == 0
    error(
      "In `svd`, the left or right indices are empty (the indices of `A` are ($(inds(A))), but the input indices are ($Lis)). For now, this is not supported. You may have accidentally input the wrong indices.",
    )
  end

  CL = combiner(Lis...)
  CR = combiner(Ris...)

  AC = A * CR * CL

  cL = combinedind(CL)
  cR = combinedind(CR)
  if inds(AC) != (cL, cR)
    AC = permute(AC, cL, cR)
  end

  USVT = svd(tensor(AC); kwargs...)
  if isnothing(USVT)
    return nothing
  end
  UT, ST, VT, spec = USVT
  UC, S, VC = itensor(UT), itensor(ST), itensor(VT)

  u = commonind(S, UC)
  v = commonind(S, VC)

  if hasqns(A)
    # Fix the flux of UC,S,VC
    # such that flux(UC) == flux(VC) == QN()
    # and flux(S) == flux(A)
    for b in nzblocks(UC)
      i1 = inds(UC)[1]
      i2 = inds(UC)[2]
      newqn = -dir(i2) * flux(i1 => Block(b[1]))
      setblockqn!(i2, newqn, b[2])
      setblockqn!(u, newqn, b[2])
    end

    for b in nzblocks(VC)
      i1 = inds(VC)[1]
      i2 = inds(VC)[2]
      newqn = -dir(i2) * flux(i1 => Block(b[1]))
      setblockqn!(i2, newqn, b[2])
      setblockqn!(v, newqn, b[2])
    end
  end

  U = UC * dag(CL)
  V = VC * dag(CR)

  settags!(U, utags, u)
  settags!(S, utags, u)
  settags!(S, vtags, v)
  settags!(V, vtags, v)

  u = settags(u, utags)
  v = settags(v, vtags)

  return TruncSVD(U, S, V, spec, u, v)
end

"""
    TruncEigen

ITensor factorization type for a truncated eigenvalue
decomposition, returned by `eigen`.
"""
struct TruncEigen
  D::ITensor
  V::ITensor
  Vt::ITensor
  spec::Spectrum
  l::Index
  r::Index
end

# iteration for destructuring into components `D, V, spec, l, r = E`
iterate(E::TruncEigen) = (E.D, Val(:V))
iterate(E::TruncEigen, ::Val{:V}) = (E.V, Val(:spec))
iterate(E::TruncEigen, ::Val{:spec}) = (E.spec, Val(:l))
iterate(E::TruncEigen, ::Val{:l}) = (E.l, Val(:r))
iterate(E::TruncEigen, ::Val{:r}) = (E.r, Val(:done))
iterate(E::TruncEigen, ::Val{:done}) = nothing

"""
    eigen(A::ITensor[, Linds, Rinds]; <keyword arguments>)

Eigendecomposition of an ITensor `A`, computed
by treating the "left indices" `Linds` provided collectively
as a row index, and remaining "right indices" `Rinds` as a
column index (matricization of a tensor).

If no indices are provided, pairs of primed and unprimed indices are
searched for, with `Linds` taken to be the primed indices and
`Rinds` taken to be the unprimed indices.

The return arguments are the eigenvalues `D` and eigenvectors `U`
as tensors, such that `A * U ∼ U * D` (more precisely they are approximately
equal up to proper replacements of indices, see the example for details).

Whether or not `eigen` performs a trunction depends on the keyword
arguments provided. Note that truncation is only well defined for
positive semidefinite matrices.

# Arguments
- `maxdim::Int`: the maximum number of singular values to keep.
- `mindim::Int`: the minimum number of singular values to keep.
- `cutoff::Float64`: set the desired truncation error of the eigenvalues, by default defined as the sum of the squares of the smallest eigenvalues. For now truncation is only well defined for positive semi-definite eigenspectra.
- `ishermitian::Bool = false`: specify if the matrix is Hermitian, in which case a specialized diagonalization routine will be used and it is guaranteed that real eigenvalues will be returned.
- `plev::Int = 0`: set the prime level of the Indices of `D`. Default prime levels are subject to change.
- `leftplev::Int = plev`: set the prime level of the Index unique to `D`. Default prime levels are subject to change.
- `rightplev::Int = leftplev+1`: set the prime level of the Index shared by `D` and `U`. Default tags are subject to change.
- `tags::String = "Link,eigen"`: set the tags of the Indices of `D`. Default tags are subject to change.
- `lefttags::String = tags`: set the tags of the Index unique to `D`. Default tags are subject to change.
- `righttags::String = tags`: set the tags of the Index shared by `D` and `U`. Default tags are subject to change.
- `use_absolute_cutoff::Bool = false`: set if all probability weights below the `cutoff` value should be discarded, rather than the sum of discarded weights.
- `use_relative_cutoff::Bool = true`: set if the singular values should be normalized for the sake of truncation.

# Examples
```julia
i, j, k, l = Index(2, "i"), Index(2, "j"), Index(2, "k"), Index(2, "l")
A = randomITensor(i, j, k, l)
Linds = (i, k)
Rinds = (j, l)
D, U = eigen(A, Linds, Rinds)
dl, dr = uniqueind(D, U), commonind(D, U)
Ul = replaceinds(U, (Rinds..., dr) => (Linds..., dl))
A * U ≈ Ul * D # true
```

See also: [`svd`](@ref), [`factorize`](@ref)
"""
function eigen(A::ITensor, Linds, Rinds; kwargs...)
  @debug_check begin
    if hasqns(A)
      @assert flux(A) == QN()
    end
  end

  N = ndims(A)
  NL = length(Linds)
  NR = length(Rinds)
  NL != NR && error("Must have equal number of left and right indices")
  N != NL + NR &&
    error("Number of left and right indices must add up to total number of indices")

  ishermitian::Bool = get(kwargs, :ishermitian, false)

  tags::TagSet = get(kwargs, :tags, "Link,eigen")
  lefttags::TagSet = get(kwargs, :lefttags, tags)
  righttags::TagSet = get(kwargs, :righttags, tags)

  plev::Int = get(kwargs, :plev, 0)
  leftplev::Int = get(kwargs, :leftplev, plev)
  rightplev::Int = get(kwargs, :rightplev, plev)

  if lefttags == righttags && leftplev == rightplev
    leftplev = rightplev + 1
  end

  # Linds, Rinds may not have the correct directions
  Lis = indices(Linds)
  Ris = indices(Rinds)

  # Ensure the indices have the correct directions,
  # QNs, etc.
  # First grab the indices in A, then permute them
  # correctly.
  Lis = permute(commoninds(A, Lis), Lis)
  Ris = permute(commoninds(A, Ris), Ris)

  for (l, r) in zip(Lis, Ris)
    if space(l) != space(r)
      error("In eigen, indices must come in pairs with equal spaces.")
    end
    if hasqns(A)
      if dir(l) == dir(r)
        error("In eigen, indices must come in pairs with opposite directions")
      end
    end
  end

  CL = combiner(Lis...; dir=Out, tags="CMB,left")
  CR = combiner(dag(Ris)...; dir=Out, tags="CMB,right")

  AC = A * dag(CR) * CL

  cL = combinedind(CL)
  cR = dag(combinedind(CR))
  if inds(AC) != (cL, cR)
    AC = permute(AC, cL, cR)
  end

  AT = ishermitian ? Hermitian(tensor(AC)) : tensor(AC)

  DT, VT, spec = eigen(AT; kwargs...)
  D, VC = itensor(DT), itensor(VT)

  V = VC * dag(CR)

  # Set right index tags
  l = uniqueind(D, V)
  r = commonind(D, V)
  l̃ = setprime(settags(l, lefttags), leftplev)
  r̃ = setprime(settags(l̃, righttags), rightplev)

  replaceinds!(D, (l, r), (l̃, r̃))
  replaceind!(V, r, r̃)

  l, r = l̃, r̃

  # The right eigenvectors, after being applied to A
  Vt = replaceinds(V, (Ris..., r), (Lis..., l))

  @debug_check begin
    if hasqns(A)
      @assert flux(D) == QN()
      @assert flux(V) == QN()
      @assert flux(Vt) == QN()
    end
  end

  return TruncEigen(D, V, Vt, spec, l, r)
end

function eigen(A::ITensor; kwargs...)
  Ris = filterinds(A; plev=0)
  Lis = Ris'
  return eigen(A, Lis, Ris; kwargs...)
end

function noinds_error_message(decomp::String)
  return "$decomp without any input indices is currently not defined. In the future it may be defined as performing a $decomp decomposition treating the ITensor as a matrix from the primed to the unprimed indices."
end

qr(A::ITensor; kwargs...) = error(noinds_error_message("qr"))

# TODO: write this in terms of combiners and then
# call qr on the order-2 tensors directly
function qr(A::ITensor, Linds...; kwargs...)
  tags::TagSet = get(kwargs, :tags, "Link,qr")
  Lis = commoninds(A, indices(Linds))
  Ris = uniqueinds(A, Lis)
  Lpos, Rpos = NDTensors.getperms(inds(A), Lis, Ris)
  QT, RT = qr(tensor(A), Lpos, Rpos; kwargs...)
  Q, R = itensor(QT), itensor(RT)
  q = commonind(Q, R)
  settags!(Q, tags, q)
  settags!(R, tags, q)
  q = settags(q, tags)
  return Q, R, q
end

polar(A::ITensor; kwargs...) = error(noinds_error_message("polar"))

# TODO: allow custom tags in internal indices?
# TODO: return the new common indices?
function polar(A::ITensor, Linds...; kwargs...)
  U, S, V = svd(A, Linds...; kwargs...)
  u = commoninds(S, U)
  v = commoninds(S, V)
  δᵤᵥ′ = δ(u..., v'...)
  Q = U * δᵤᵥ′ * V'
  P = dag(V') * dag(δᵤᵥ′) * S * V
  return Q, P, commoninds(Q, P)
end

function factorize_qr(A::ITensor, Linds...; kwargs...)
  ortho::String = get(kwargs, :ortho, "left")
  if ortho == "left"
    L, R, q = qr(A, Linds...; kwargs...)
  elseif ortho == "right"
    Lis = uniqueinds(A, indices(Linds))
    R, L, q = qr(A, Lis...; kwargs...)
  else
    error(
      "In factorize using qr decomposition, ortho keyword $ortho not supported. Supported options are left or right.",
    )
  end
  return L, R
end

function factorize_svd(A::ITensor, Linds...; kwargs...)
  ortho::String = get(kwargs, :ortho, "left")
  alg::String = get(kwargs, :svd_alg, "divide_and_conquer")
  USV = svd(A, Linds...; kwargs..., alg=alg)
  if isnothing(USV)
    return nothing
  end
  U, S, V, spec, u, v = USV
  if ortho == "left"
    L, R = U, S * V
  elseif ortho == "right"
    L, R = U * S, V
  elseif ortho == "none"
    sqrtS = S
    sqrtS .= sqrt.(S)
    L, R = U * sqrtS, sqrtS * V
    replaceind!(L, v, u)
  else
    error(
      "In factorize using svd decomposition, ortho keyword $ortho not supported. Supported options are left, right, or none.",
    )
  end
  return L, R, spec
end

function factorize_eigen(A::ITensor, Linds...; kwargs...)
  ortho::String = get(kwargs, :ortho, "left")
  delta_A2 = get(kwargs, :eigen_perturbation, nothing)
  if ortho == "left"
    Lis = commoninds(A, indices(Linds))
  elseif ortho == "right"
    Lis = uniqueinds(A, indices(Linds))
  else
    error(
      "In factorize using eigen decomposition, ortho keyword $ortho not supported. Supported options are left or right.",
    )
  end
  simLis = sim(Lis)
  A2 = A * replaceinds(dag(A), Lis, simLis)
  if !isnothing(delta_A2)
    # This assumes delta_A2 has indices:
    # (Lis..., prime(Lis)...)
    delta_A2 = replaceinds(delta_A2, Lis, dag(simLis))
    noprime!(delta_A2)
    A2 += delta_A2
  end
  F = eigen(A2, Lis, simLis; ishermitian=true, kwargs...)
  D, _, spec = F
  L = F.Vt
  R = dag(L) * A
  if ortho == "right"
    L, R = R, L
  end
  return L, R, spec
end

factorize(A::ITensor; kwargs...) = error(noinds_error_message("factorize"))

"""
    factorize(A::ITensor, Linds::Index...; <keyword arguments>)

Perform a factorization of `A` into ITensors `L` and `R` such that `A ≈ L * R`.

# Arguments
- `ortho::String = "left"`: Choose orthogonality properties of the factorization.
  - `"left"`: the left factor `L` is an orthogonal basis such that `L * dag(prime(L, commonind(L,R))) ≈ I`.
  - `"right"`: the right factor `R` forms an orthogonal basis.
  - `"none"`, neither of the factors form an orthogonal basis, and in general are made as symmetrically as possible (depending on the decomposition used).
- `which_decomp::Union{String, Nothing} = nothing`: choose what kind of decomposition is used.
  - `nothing`: choose the decomposition automatically based on the other arguments. For example, when `nothing` is chosen and `ortho = "left"` or `"right"`, and a cutoff is provided, `svd` or `eigen` is used depending on the provided cutoff (`eigen` is only used when the cutoff is greater than `1e-12`, since it has a lower precision). When no truncation is requested `qr` is used for dense ITensors and `svd` for block-sparse ITensors (in the future `qr` will be used also for block-sparse ITensors in this case).
  - `"svd"`: `L = U` and `R = S * V` for `ortho = "left"`, `L = U * S` and `R = V` for `ortho = "right"`, and `L = U * sqrt.(S)` and `R = sqrt.(S) * V` for `ortho = "none"`. To control which `svd` algorithm is choose, use the `svd_alg` keyword argument. See the documentation for `svd` for the supported algorithms, which are the same as those accepted by the `alg` keyword argument.
  - `"eigen"`: `L = U` and ``R = U^{\\dagger} A`` where `U` is determined from the eigendecompositon ``A A^{\\dagger} = U D U^{\\dagger}`` for `ortho = "left"` (and vice versa for `ortho = "right"`). `"eigen"` is not supported for `ortho = "none"`.
  - `"qr"`: `L=Q` and `R` an upper-triangular matrix when `ortho = "left"`, and `R = Q` and `L` a lower-triangular matrix when `ortho = "right"` (currently supported for dense ITensors only).
In the future, other decompositions like QR (for block-sparse ITensors), polar, cholesky, LU, etc. are expected to be supported.

For truncation arguments, see: [`svd`](@ref)
"""
function factorize(A::ITensor, Linds...; kwargs...)
  ortho::String = get(kwargs, :ortho, "left")
  tags::TagSet = get(kwargs, :tags, "Link,fact")
  plev::Int = get(kwargs, :plev, 0)
  which_decomp::Union{String,Nothing} = get(kwargs, :which_decomp, nothing)
  cutoff = get(kwargs, :cutoff, nothing)
  eigen_perturbation = get(kwargs, :eigen_perturbation, nothing)
  if !isnothing(eigen_perturbation)
    if !(isnothing(which_decomp) || which_decomp == "eigen")
      error("""when passing a non-trivial eigen_perturbation to `factorize`,
               the which_decomp keyword argument must be either "automatic" or
               "eigen" """)
    end
    which_decomp = "eigen"
  end

  # Deprecated keywords
  if haskey(kwargs, :dir)
    error(
      """dir keyword in factorize has been replace by ortho.
Note that the default is now `left`, meaning for the results L,R = factorize(A), L forms an orthogonal basis.""",
    )
  end

  if haskey(kwargs, :which_factorization)
    error("""which_factorization keyword in factorize has been replace by which_decomp.""")
  end

  # Determines when to use eigen vs. svd (eigen is less precise,
  # so eigen should only be used if a larger cutoff is requested)
  automatic_cutoff = 1e-12
  dL, dR = dim(indices(Linds)), dim(indices(setdiff(inds(A), Linds)))
  maxdim = get(kwargs, :maxdim, min(dL, dR))
  might_truncate = !isnothing(cutoff) || maxdim < min(dL, dR)

  if isnothing(which_decomp)
    if !might_truncate && !hasqns(A) && ortho != "none"
      which_decomp = "qr"
    elseif isnothing(cutoff) || cutoff ≤ automatic_cutoff
      which_decomp = "svd"
    elseif cutoff > automatic_cutoff
      which_decomp = "eigen"
    end
  end

  if which_decomp == "svd"
    LR = factorize_svd(A, Linds...; kwargs...)
    if isnothing(LR)
      return nothing
    end
    L, R, spec = LR
  elseif which_decomp == "eigen"
    L, R, spec = factorize_eigen(A, Linds...; kwargs...)
  elseif which_decomp == "qr"
    hasqns(A) && error("QR factorization of an ITensor with QNs is not yet supported.")
    L, R = factorize_qr(A, Linds...; kwargs...)
    spec = Spectrum(nothing, 0.0)
  else
    throw(
      ArgumentError(
        """In factorize, factorization $which_decomp is not currently supported. Use `"svd"`, `"eigen"`, `"qr"` or `nothing`.""",
      ),
    )
  end

  # Set the tags and prime level
  l = commonind(L, R)
  l̃ = setprime(settags(l, tags), plev)
  replaceind!(L, l, l̃)
  replaceind!(R, l, l̃)
  l = l̃

  return L, R, spec, l
end


text:
using InteractionNetworkModels, Distributions, Distances, StructuredDistances
using BenchmarkTools, StatsBase, Plots, StatsPlots

V = 50
Sᵐ = [rand(1:V,rand(2:5)) for i in 1:10]
d = FastEditDistance(FastLCS(100), 100)
γ = 3.0
K_inner, K_outer = (DimensionRange(1,10), DimensionRange(1,25))
# K_inner, K_outer = (Inf,Inf)
model = SIS(
    Sᵐ, γ,
    d,
    1:V,
    K_inner, K_outer
)

aux_mcmc = SisMcmcInsertDelete(
    ν_ed=1, ν_td=1, β=0.7, len_dist=TrGeometric(0.9,1,K_inner.u),
    lag=1, burn_in=0,
    K=200
)

S_curr = deepcopy(Sᵐ)
S_prop = S_curr[1:(end-1)]

aux_eval = SisAuxTermEvaluator(aux_mcmc, d, 1:V, K_inner, K_outer)
initialise!(aux_eval, S_prop)

n = 100
aux_data = [[Int[]] for i in 1:n]

@time aux_term_current(
    S_curr, S_prop,
    γ, d, 1:V, K_inner, K_outer,
    aux_mcmc, aux_data
)

@time aux_term_new(
    S_curr, S_prop,
    γ,
    n,
    aux_eval
)

plot(map(x->d(x, S_prop), y), label="New")
plot!(map(x->d(x, S_prop), aux_data), label="Old")

aux_eval.S_curr

aux_mcmc.dist_curr

@btime aux_term_current(
    $S_curr, $S_prop,
    $γ, $d, 1:$V, $K_inner, $K_outer,
    $aux_mcmc, $aux_data
)

@btime aux_term_new(
    $S_curr, $S_prop,
    $γ,
    $n,
    $aux_eval
)

plot(map(x->d(x, S_prop), y), label="New")
plot!(map(x->d(x, S_prop), aux_data), label="Old")

m = 100
current = [aux_term_current(
                S_curr, S_prop,
                γ, d, 1:V, K_inner, K_outer,
                aux_mcmc, aux_data
            ) for i in 1:m
];
new = [aux_term_new(
            S_curr, S_prop,
            γ,
            n,
            aux_eval
        ) for i in 1:m
];

dotplot([current, new])


text:
module SmallSets

using StaticArrays

import Base: show_vector, join_eltype, issubset,
             _default_eltype, grow_to!, eltype,
             similar, show, isempty, length, in,
             push!, pop!,
             delete!, sizehint!, empty!, rehash!,
             start, done, next,
             union, union!, intersect, setdiff!,
             ==, <, <=

export SmallSet, issmall

mutable struct SmallSet{N, T} <: Base.AbstractSet{T} where N <: Integer
  vsize::Int
  vector::MVector{N, T} # TODO Union Vector/MVector slows down everything...
  set::Union{Set{T}, Void}

  SmallSet{N, T}() where {N, T} = begin
    if N > 32
      warn("Consider using Set when using N > 32.")
    end
    new(0, MVector{N, T}(), nothing)
  end
end

SmallSet{N}() where N <: Integer = SmallSet{N, Any}()
SmallSet{N}(itr) where {N} = SmallSet{N, eltype(itr)}(itr)
SmallSet{N, T}(itr) where {N, T} = union!(SmallSet{N, T}(), itr)
function SmallSet{N}(g::Base.Generator) where {N}
  T = _default_eltype(typeof(g))
  (isleaftype(T) || T === Union{}) || return grow_to!(SmallSet{N, T}(), g)
  return SmallSet{N, T}(g)
end

# private
join_sizes() = 0
join_sizes(v1::SmallSet{N, T}, vs::SmallSet...) where {N, T} = max(N, join_sizes(vs...))

# public
issmall(s::SmallSet) = typeof(s.set) === Void
eltype(::Type{SmallSet{N, T}}) where {N, T} = T
similar(s::SmallSet{N, T}) where {N, T} = SmallSet{N, T}()
similar(s::SmallSet{N, T}, NT::Type) where {N, T} = SmallSet{N, NT}()

function show(io::IO, s::SmallSet{N, T}) where {N, T}
  print(io, "SmallSet")
  if isempty(s)
    print(io, "{", N, ",", T, "}()")
    return
  end
  print(io,"(")
  show_vector(io, issmall(s) ? s.vector[1 : s.vsize] : s.set, "[","]") # FIXME (currently creates copy)
  print(io,")")
end

isempty(s::SmallSet) = s.vsize === 0
length(s::SmallSet{N, T}) where {N, T} = issmall(s) ? s.vsize : length(s.set::Set{T})

function in(x, s::SmallSet{N, T}) where {N, T}
  if issmall(s)
    for i = 1 : s.vsize
      if x == s.vector[i]
        return true
      end
    end
    return false
  else
    return x in s.set::Set{T}
  end
end

function push!(s::SmallSet{N, T}, x) where {N, T}
  if !(x in s)
    if s.vsize < N
      s.vsize += 1
      s.vector[s.vsize] = x
      return s
    end

    if issmall(s)
      s.set = Set{T}(s.vector)
    end

    push!(s.set::Set{T}, x)
  end
  return s
end

#TODO: Switch back to small? When?
function pop!(s::SmallSet{N, T}, x) where {N, T}
  if issmall(s)
    i = findfirst(s.vector, x)
    if i === 0
      throw(KeyError(x))
    else
      s.vector[i] = s.vector[s.vsize]
      s.vsize -= 1
    end
  else
    pop!(s.set::Set{T}, x)
  end
  return x
end
pop!(s::SmallSet, x, deflt) = x in s ? pop!(s, x) : deflt
pop!(s::SmallSet{N, T}) where {N, T} = issmall(s) ? pop!(s, s.vector[1]) : pop!(s, s.set::Set{T})
delete!(s::SmallSet, x) = (pop!(s, x, 0); s)

copy(s::SmallSet) = union!(similar(s), s)

sizehint!(s::SmallSet, newsz) = issmall(s) ? s : sizehint!(s.set) # FIXME?
empty!(s::SmallSet) = (s.set = nothing; s.vsize = 0; s)
rehash!(s::SmallSet) = issmall(s) ? s : rehash!(s.set)

start(s::SmallSet)       = issmall(s) ? start(s.vector) : start(s.set)
done(s::SmallSet, state) = issmall(s) ? done(s.vector, state) : done(s.set, state)
next(s::SmallSet, i)     = issmall(s) ? next(s.vector, i) : next(s.set, i)

union(s::SmallSet) = copy(s)
function union(s::SmallSet, sets::SmallSet...)
  u = SmallSet{join_sizes(s, sets...), join_eltype(s, sets...)}()
  union!(u, s)
  for t in sets
    union!(u, t)
  end
  return u
end
const ∪ = union
union!(s::SmallSet, xs) = (for x = xs; push!(s, x); end; s)
union!(s::SmallSet, xs::AbstractArray) = (sizehint!(s, length(xs)); for x = xs; push!(s, x); end; s) # FIXME

intersect(s::SmallSet) = copy(s)
#function intersect(s::Set, sets::Set...)
const ∩ = intersect
#function setdiff(a::Set, b::Set)
setdiff!(s::SmallSet, xs) = (for x=xs; delete!(s,x); end; s)

==(l::SmallSet, r::SmallSet) = (length(l) == length(r)) && (l <= r)
<( l::SmallSet, r::SmallSet) = (length(l) < length(r)) && (l <= r)
<=(l::SmallSet, r::SmallSet) = issubset(l, r)

const ⊆ = issubset
⊊(l::SmallSet, r::SmallSet) = <(l, r)
⊈(l::SmallSet, r::SmallSet) = !⊆(l, r)
⊇(l, r) = issubset(r, l)
⊉(l::SmallSet, r::SmallSet) = !⊇(l, r)
⊋(l::SmallSet, r::SmallSet) = <(r, l)

end


text:
mutable struct Node
    name::String
    parent::Union{Node, Nothing}
    children::Array{Node,1}
end
Base.show(io::IO, z::Node) = print(io, z.name)

function find_or_create(name::AbstractString, arr::Array{Node,1})::Node
    i = findfirst(x -> x.name == name, arr)
    if i === nothing
        n = Node(name, nothing, [])
        push!(arr, n)
        return n
    else
        return arr[i]
    end
end

function climb(node::Node, arr = Node[])::Array{Node,1}
    if node.parent === nothing
        return arr
    else
        climb(node.parent, push!(arr, node.parent))
    end
end

function climb_to(node::Node, target::Node, arr = Node[])::Array{Node,1}
    if node.name == target.name
        return arr
    elseif node.parent === nothing
        throw(ErrorException("target node not found"))
    else
        climb_to(node.parent, target, push!(arr, node.parent))
    end
end

function build_orbitmap()::Array{Node,1}
    data = open("orbitmap.txt") do f
        split.(readlines(f), ")")
    end
    objects = Node[]
    for o in data
        parent_name, child_name = o
        parent = find_or_create(parent_name, objects)
        child = find_or_create(child_name, objects)
        push!(parent.children, child)
        child.parent = parent
    end

    return objects
end

function orbit_count(objects::Array{Node,1})::Int
    orbit_count = 0
    for o in objects
        orbit_count = orbit_count + length(climb(o))
    end

    return orbit_count
end

function orbit_distance(from::Node, to::Node)::Int
    root = intersect(climb(from), climb(to))[1]

    return length(climb_to(you, root)) + length(climb_to(santa, root)) - 2
end

orbitmap = build_orbitmap()

println(orbit_count(orbitmap))

you = find_or_create("YOU", orbitmap)
santa = find_or_create("SAN", orbitmap)
println(orbit_distance(you, santa))


text:
# no module defined, so testsets added to parentmodule
using ReTest

@testset "FakePackage: no module" begin
    @test true
end

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
export basis, dbasis, dbasisu

"""
    basis(npts, ord, t, x[])

Generate a B-spline basis functions `N[]` for a given open knot vector `x[]`.

A B-Spline basis is a collection of functions of a parameter `t` wich form
 a basis for the vectorial space of functions. The conformation of this set
 higly depends on the choosing of the knots `x[]` the curve is bound to.
The basis is computed with _Cox-de Boor_ recursive function applied to the
 basis dependency tree in order to optimise the computation.

---

# Arguments
- `npts::Int64`: the number of the points of the controll polygon.
- `ord::Int64`: the order of the B-Spline (`deg = ord-1`).
- `t::Float64`: the parameter value of the parametric curve.
- `x::Array{Float64}`: the knot vector.

---

# Examples

```jldoctest
julia> basis(5,3,0.5,[0. 0 0 1 2 3 3 3])
5-element Array{Float64,1}:
 0.25
 0.625
 0.125
 0.0
 0.0
```

---

_By Elia Onofri, Giuseppe Santorelli_
"""

function basis(npts::Int64, ord::Int64, t::Float64, x::Array{Float64})::Array{Float64}

    local m::Int64 = npts + ord	    # Dimension of the knot vector

    local N::Array{Float64} = zeros(m)   # Basis progressive vector

    # Local check of the knot vector correctness
    @assert length(x) == m ("ERROR: incompatibile knot vector with given parameters n+1 = $(npts), k = $(ord)")

    # Eval N_{i,1} for i = 1:max_B
    for i = 1 : m - 1
        if (t >= x[i]) && (t < x[i+1])
            N[i] = 1
        else
            N[i] = 0
        end
    end

    # Eval higher basis N_{i,deg} for deg = 2:ord and i = 1:max_B-deg
    for deg = 2 : ord
        for i = 1 : m - deg
            # Eval of the direct dependency
            if N[i] == 0
                ddep = 0.0
            else
                ddep = ((t - x[i]) * N[i]) / (x[i+deg-1] - x[i])
            end
            # Eval of the forward dependency
            if N[i+1] == 0
                fdep = 0.0
            else
                fdep = ((x[i+deg] - t) * N[i+1]) / (x[i+deg] - x[i+1])
            end
            # Collection of the dependencies
            N[i] = ddep + fdep
        end
    end

    # Otherwise last point is zero
    if t == x[m]
        N[npts] = 1
    end

    return N[1:npts]
end

#--------------------------------------------------------------------------------------------------------------------------------

"""
    dbasis(npts, ord, t, x[])

Generate a B-spline basis functions `N[]` and its derivatives `D1[], D2[]` for a given open knot vector `x[]`.

A B-Spline basis is a collection of functions of a parameter `t` wich form
 a basis for the vectorial space of functions. The conformation of this set
 higly depends on the choosing of the knots `x[]` the curve is bound to.
The basis is computed with _Cox-de Boor_ recursive function applied to the
 basis dependency tree in order to optimise the computation.

---

# Arguments
- `npts::Int64`: the number of the points of the controll polygon.
- `ord::Int64`: the order of the B-Spline (`deg = ord-1`).
- `t::Float64`: the parameter value of the parametric curve.
- `x::Array{Float64}`: the knot vector.

---

# Examples

```jldoctest
julia> dbasis(5,3,0.5,[0. 0 0 1 2 3 3 3])[1]
5-element Array{Float64,1}:
 0.25
 0.625
 0.125
 0.0
 0.0

julia> dbasis(5,3,0.5,[0. 0 0 1 2 3 3 3])[2]
5-element Array{Float64,1}:
 -1.0
  0.5
  0.5
  0.0
  0.0

julia> dbasis(5,3,0.5,[0. 0 0 1 2 3 3 3])[3]
5-element Array{Float64,1}:
  2.0
 -3.0
  1.0
  0.0
  0.0
```

---

_By Elia Onofri, Giuseppe Santorelli_
"""

function dbasis(npts::Int64, ord::Int64, t::Float64, x::Array{Float64})::Tuple{Array{Float64}, Array{Float64}, Array{Float64}}

    local m::Int64 = npts + ord			# Dimension of the knot vector

    local N::Array{Float64} = zeros(m)   # Basis progressive vector
    local D1::Array{Float64} = zeros(m)  # First Derivative progressive vector
    local D2::Array{Float64} = zeros(m)  # Second Derivative progressive vector

    # Local check of the knot vector correctness
    @assert length(x) == m ("ERROR: incompatibile knot vector with given parameters n+1 = $(npts), k = $(ord)")

    # Eval N_{i,1} for i = 1:max_B
    for i = 1 : m - 1
        if (t >= x[i]) && (t < x[i+1])
            N[i] = 1
        else
            N[i] = 0
        end
    end

    if t == x[m]    # last(x)
        N[npts] = 1
    end

    # Eval higher basis N_{i,deg} for deg = 2:ord and i = 1:max_B-deg
    for deg = 2 : ord
        for i = 1 : m - deg

            # ----Eval of the direct dependency----
            if N[i] == 0
                ddep = 0.0
                f1 = 0.0
            else
                ddep = ((t - x[i]) * N[i]) / (x[i+deg-1] - x[i])
                f1 = N[i] / (x[i+deg-1] - x[i])
            end

            # ----Eval of the forward dependency----
            if N[i+1] == 0
                fdep = 0.0
                f2 = 0.0
            else
                fdep = ((x[i+deg] - t) * N[i+1]) / (x[i+deg] - x[i+1])
                f2 = -N[i+1] / (x[i+deg] - x[i+1])
            end

            #----Eval of the direct first derivative dependency----
            if D1[i] != 0    # if the lower order basis function is zero skip the calculation
                f3 = ((t - x[i]) * D1[i]) / (x[i+deg-1] - x[i])
                s1 = (2 * D1[i]) / (x[i+deg-1] - x[i])
            else
                f3 = 0.0
                s1 = 0.0
            end

            #----Eval of the forward first derivative dependency----
            if D1[i+1] != 0    # if the lower order basis function is zero skip the calculation
                f4 = ((x[i+deg] - t) * D1[i+1]) / (x[i+deg] - x[i+1])
                s2 = ((-2) * D1[i+1]) / (x[i+deg] - x[i+1])
            else
                f4 = 0.0
                s2 = 0.0
            end

            #----Eval of the direct second derivative dependency----
            if D2[i] != 0    # if the lower order basis function is zero skip the calculation
                s3 = ((t - x[i]) * D2[i] ) / (x[i+deg-1]-x[i])
            else
                s3 = 0.0
            end

            #----Eval of the forward second derivative dependency----
            if D2[i+1] != 0    # if the lower order basis function is zero skip the calculation
                s4 = ((x[i+deg] - t) * D2[i+1]) / (x[i+deg] - x[i+1])
            else
                s4 = 0.0
            end

            # Collection of the dependencies
            N[i] = ddep + fdep
            D1[i] = f1 + f2 + f3 + f4
            D2[i] = s1 + s2 + s3 + s4
        end
    end

    return (N[1:npts], D1[1:npts], D2[1:npts])
end

#--------------------------------------------------------------------------------------------------------------------------------

"""
    dbasisu(npts, ord, t, x[])

Generate a B-spline basis functions `N[]` and its derivatives `D1[], D2[]` for a given uniform periodic knot vector `x[]`.

A B-Spline basis is a collection of functions of a parameter `t` wich form
 a basis for the vectorial space of functions. The conformation of this set
 higly depends on the choosing of the knots `x[]` the curve is bound to.
The basis is computed with _Cox-de Boor_ recursive function applied to the
 basis dependency tree in order to optimise the computation.

---

# Arguments
- `npts::Int64`: the number of the points of the controll polygon.
- `ord::Int64`: the order of the B-Spline (`deg = ord-1`).
- `t::Float64`: the parameter value of the parametric curve.
- `x::Array{Float64}`: the knot vector.

---

# Examples

```jldoctest
julia> dbasisu(5,3,0.5,[0. 1 2 3 4 5 6 7])[1]
5-element Array{Float64,1}:
 0.125
 0.0
 0.0
 0.0
 0.0

julia> dbasisu(5,3,0.5,[0. 1 2 3 4 5 6 7])[2]
5-element Array{Float64,1}:
 0.5
 0.0
 0.0
 0.0
 0.0

julia> dbasisu(5,3,0.5,[0. 1 2 3 4 5 6 7])[3]
5-element Array{Float64,1}:
 1.0
 0.0
 0.0
 0.0
 0.0
```

---

_By Elia Onofri, Giuseppe Santorelli_
"""

function dbasisu(npts::Int64, ord::Int64, t::Float64, x::Array{Float64})::Tuple{Array{Float64}, Array{Float64}, Array{Float64}}
    local m::Int64 = npts+ord			# Dimension of the knot vector
    local N::Array{Float64} = zeros(m)   # Basis progressive vector
    local D1::Array{Float64} = zeros(m)  # First Derivative progressive vector
    local D2::Array{Float64} = zeros(m)  # Second Derivative progressive vector

    # Local check of the knot vector correctness
    @assert length(x) == m ("ERROR: incompatibile knot vector with given parameters n+1 = $(npts), k = $(ord)")

    # Eval N_{i,1} for i = 1:max_B
    for i = 1 : m-1
        if (t >= x[i]) && (t < x[i+1])
            N[i] = 1
        else
            N[i] = 0
        end
    end

    if t == x[npts+1]    # handle the end specially
        N[npts] = 1   	 # resetting the first order basis functions.
        N[npts+1] = 0
    end

    # Eval higher basis N_{i,deg} for deg = 2:ord and i = 1:max_B-deg
    for deg = 2 : ord
        for i = 1 : m-deg

            # ----Eval of the direct dependency----
            if N[i] == 0
                ddep = 0.0
                f1 = 0.0
            else
                ddep = ((t - x[i]) * N[i]) / (x[i+deg-1] - x[i])
                f1 = N[i] / (x[i+deg-1] - x[i])
            end

            # ----Eval of the forward dependency----
            if N[i+1] == 0
                fdep = 0.0
                f2 = 0.0
            else
                fdep = ((x[i+deg] - t) * N[i+1]) / (x[i+deg] - x[i+1])
                f2 = -N[i+1] / (x[i+deg] - x[i+1])
            end

            #----Eval of the direct first derivative dependency----
            if D1[i] != 0    # if the lower order basis function is zero skip the calculation
                f3 = ((t - x[i]) * D1[i]) / (x[i+deg-1] - x[i])
                s1 = (2 * D1[i]) / (x[i+deg-1] - x[i])
            else
                f3 = 0.0
                s1 = 0.0
            end

            #----Eval of the forward first derivative dependency----
            if D1[i+1] != 0    # if the lower order basis function is zero skip the calculation
                f4 = ((x[i+deg] - t) * D1[i+1]) / (x[i+deg] - x[i+1])
                s2 = ((-2) * D1[i+1]) / (x[i+deg] - x[i+1])
            else
                f4 = 0.0
                s2 = 0.0
            end

            #----Eval of the direct second derivative dependency----
            if D2[i] != 0    # if the lower order basis function is zero skip the calculation
                s3 = ((t - x[i]) * D2[i] ) / (x[i+deg-1]-x[i])
            else
                s3 = 0.0
            end

            #----Eval of the forward second derivative dependency----
            if D2[i+1] != 0    # if the lower order basis function is zero skip the calculation
                s4 = ((x[i+deg] - t) * D2[i+1]) / (x[i+deg] - x[i+1])
            else
                s4 = 0.0
            end

            # Collection of the dependencies
            N[i] = ddep + fdep
            D1[i] = f1 + f2 + f3 + f4
            D2[i] = s1 + s2 + s3 + s4
        end
    end

    return (N[1:npts], D1[1:npts], D2[1:npts])
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -3.517322036224115 <= q <= 6.144134546328734)

                     add_NL_constraint(m, :(1/(1 + exp(-(0.9006950404367178 + -0.8785024386485549*1/(1 + exp(-(0.872987896263127 + 0.13249617278479375*$(x[1]) + 0.7878182386475117*$(x[2]) + -0.4203533662190604*$(x[3]) + -0.7649366056085984*$(x[4]) + -0.57871509222098*$(x[5])))) + 0.18543992065441683*1/(1 + exp(-(0.8256824299447256 + -0.6033988061517293*$(x[1]) + -0.15324994390362567*$(x[2]) + 0.6933635130851963*$(x[3]) + 0.1238690975426957*$(x[4]) + 0.2325901400223671*$(x[5]))))))) + 1/(1 + exp(-(0.10821024944347046 + 0.16585705930469619*1/(1 + exp(-(0.872987896263127 + 0.13249617278479375*$(x[1]) + 0.7878182386475117*$(x[2]) + -0.4203533662190604*$(x[3]) + -0.7649366056085984*$(x[4]) + -0.57871509222098*$(x[5])))) + 0.9368217650268278*1/(1 + exp(-(0.8256824299447256 + -0.6033988061517293*$(x[1]) + -0.15324994390362567*$(x[2]) + 0.6933635130851963*$(x[3]) + 0.1238690975426957*$(x[4]) + 0.2325901400223671*$(x[5]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function roccurve(scores, truelabels)
    # truelabels l ϵ {0, 1}
    if all(unique(truelabels) .== [1, 2])
        truelabels .-= 1
    end

    descendingidx = sortperm(scores, rev = true)
    scores = scores[descendingidx]
    truelabels = truelabels[descendingidx]

    distincvalueidx = findall(diff(scores) .!= 0)
    thresholdidx = vcat(distincvalueidx, length(truelabels))

    tps = cumsum(truelabels)[thresholdidx]
    fps = thresholdidx .- tps

    if length(tps) == 0 || fps[1] != 0 || tps[1] != 0
        # Add an extra threshold position if necessary
        # to make sure that the curve starts at (0, 0)
        tps = vcat(0, tps)
        fps = vcat(0, fps)
    end

    if fps[end] <= 0
        warn("No negative samples in y_true, false positive value should be meaningless")
        fpr = nothing
    else
        fpr = fps ./ fps[end]
    end

    if tps[end] <= 0
        warn("No positive samples in y_true, true positive value should be meaningless")
        tpr = nothing
    else
        tpr = tps ./ tps[end]
    end

    return fpr, tpr
end

auc(ŷ::Vector{T}, y::Vector{Int}) where {T<: Real} = auc(roccurve(ŷ, y)...)

function auc(x::Vector{T}, y::Vector{T}) where {T<:Real}
    if all(unique(y) == [0,1]) || all(unique(y) == [1, 2])
        @warn "it seems like you are passing a function labels, invoking roccurve"
        x, y = roccurve(x, y)
    end

    dx = diff(x)
    dy = y[2:end] + y[1:end - 1]
    return dx' * dy / 2
end


text:
"""
    lift(func, arg)
Calls `func(arg)`, propagating `missing` values
"""
lift(func, ::Missing)=missing
lift(func, arg) = func(arg)

miss_null(::Nothing) = missing
miss_null(x) = x

"""
    getfirst(dict, keys...)

Returns the element coresponding to the first key that is found.
Returns `missing` if no key is found.
"""
function getfirst(json, key, otherkeys...)
    get(json,  key) do
        getfirst(json, otherkeys...)
    end
end
getfirst(json) = missing

###############################
"""
    getpage_raw(url)

Downloads the page from the URL,
returning the raw (unparsed) text of the body.
"""
getpage_raw(url) = String(HTTP.request("GET", url).body)

"""
    getpage(url)

downloads and parses the page from the URL
"""
const getpage = parsehtml ∘ getpage_raw

###

"""
    text_only(doc)

Extracts just the unformatted text (no attributes etc),
from a HTML document or fragment(/s)
"""
text_only(doc::HTMLDocument) = text_only(doc.root)
text_only(frag) = join([replace(text(leaf), "\r" => "") for leaf in Leaves(frag) if leaf isa HTMLText], " ")
text_only(frags::Vector) = join(text_only.(frags), " ")

"""
    filter_html(text)

Strips any HTML tags out of the `text`.
If that is required.
"""
function filter_html(text)
    # Check if the text is a HTML or not
    # Note we are not parsing it, just checking if we should parse it
    if occursin(r"<(\"[^\"]*\"|'[^']*'|[^'\">])*>", text)
        # It seems like it may be HTML, so now parse it.
        text_only(parsehtml(text))
	else
    	text
    end
end

filter_html(::Missing) = missing

###########

"
    indent(str)

Indents each line in a string
"
indent(str, indentwith="\t") = join(indentwith.*strip.(split(str, "\n")), "\n")

"""
    escape_multiline_string

like Escape string, but does not escape newlines
"""
function escape_multiline_string(s::AbstractString)
    escaped = s
    escaped = replace(escaped, '\$' => raw"\$")
    escaped = replace(escaped, '\\' =>  raw"\\")
    escaped
end

function remove_cite_version(code)
    replace(code, r"\(Version \d+\) " => "")
end

"""
    leaf_subtypes(T)
Returns all the nonabstract types decedent from `T`.
"""
function leaf_subtypes(T)
       if isconcretetype(T)
           T
       else
           vcat(leaf_subtypes.(subtypes(T))...)
       end
end

##############################################

"""
	match_doi(uri::String
"""
function match_doi(uri::String)
    identifier_match = match(r"\b(10[.][0-9]{4,}(?:[.][0-9]+)*\/(?:(?![\"&\'<>])\S)+)\b", uri)
    identifier_match === nothing ? nothing : identifier_match.match
end


text:
# Autogenerated wrapper script for libsigcpp_jll for x86_64-apple-darwin
export libsigc

JLLWrappers.@generate_wrapper_header("libsigcpp")
JLLWrappers.@declare_library_product(libsigc, "@rpath/libsigc-3.0.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libsigc,
        "lib/libsigc-3.0.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using IterTools: product
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(path::String, n::Int)::Int
    text = parse.(Int, split(replace(read(path, String), "\"" => ""), ","))
    keys = Dict(i => Int[] for i ∈ 1:n)
    letters = 97:122
    for i ∈ 1:n
        for j ∈ letters
            for k ∈ i:n:length(text)
                keys[i] = ∪(keys[i], j)
                if !(32 ≤ text[k] ⊻ j ≤ 122)
                    keys[i] = setdiff(keys[i], j)
                    break
                end
            end
        end
    end
    for key ∈ product(map(x -> x[2], sort(collect(keys)))...)
        decrypted_text = ""
        result = 0
        for (i, j) ∈ enumerate(text)
            xor = j ⊻ key[(i - 1) % n + 1]
            decrypted_text = string(decrypted_text, Char(xor))
            result += xor
        end
        if occursin(" the ", decrypted_text)
            return result
        end
    end
end

compute("problems/0059/p059_cipher.txt", 3)

@benchmark compute("problems/0059/p059_cipher.txt", 3)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using Base.Test
using QuantumOptics

ωc = 1.2
ωa = 0.9
g = 1.0
γ = 0.5
κ = 1.1

T = Float64[0.,1.]

fockbasis = FockBasis(7)
spinbasis = SpinBasis(1//2)
basis = tensor(spinbasis, fockbasis)

sx = sigmax(spinbasis)
sy = sigmay(spinbasis)
sz = sigmaz(spinbasis)
sp = sigmap(spinbasis)
sm = sigmam(spinbasis)

Ha = embed(basis, 1, 0.5*ωa*sz)
Hc = embed(basis, 2, ωc*number(fockbasis))
Hint = sm ⊗ create(fockbasis) + sp ⊗ destroy(fockbasis)
H = Ha + Hc + Hint

Ja = embed(basis, 1, sqrt(γ)*sm)
Jc = embed(basis, 2, sqrt(κ)*destroy(fockbasis))
J = [Ja, Jc]

Ψ₀ = spinup(spinbasis) ⊗ fockstate(fockbasis, 5)
ρ₀ = Ψ₀⊗dagger(Ψ₀)

op1 = DenseOperator(spinbasis, [1.2+0.3im 0.7+1.2im;0.3+0.1im 0.8+3.2im])
op2 = DenseOperator(spinbasis, [0.2+0.1im 0.1+2.3im; 0.8+4.0im 0.3+1.4im])
@test tracedistance(spre(op1)*op2, op1*op2) < 1e-12
@test tracedistance(spost(op1)*op2, op2*op1) < 1e-12

@test spre(SparseOperator(op1))*op2 == op1*op2
@test spost(SparseOperator(op1))*op2 == op2*op1

L = liouvillian(H, J)
ρ = -1im*(H*ρ₀ - ρ₀*H)
for j=J
    ρ += j*ρ₀*dagger(j) - 0.5*dagger(j)*j*ρ₀ - 0.5*ρ₀*dagger(j)*j
end
@test tracedistance(L*ρ₀, ρ)<1e-10

tout, ρt = timeevolution.master([0.,1.], ρ₀, H, J; reltol=1e-7)
@test tracedistance(expm(full(L))*ρ₀, ρt[end])<1e-6


text:
module DivConformingFESpacesTests

using Test
using Gridap.Geometry
using Gridap.ReferenceFEs
using Gridap.FESpaces
using Gridap.CellData
using Gridap.TensorValues
using Gridap.Fields
using Gridap.Io

domain =(0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

order = 1

u(x) = x

reffe = ReferenceFE(raviart_thomas,order)

V = TestFESpace(model,reffe,dirichlet_tags = [1,6])
test_single_field_fe_space(V)

U = TrialFESpace(V,u)

uh = interpolate(u,U)

e = u - uh

Ω = Triangulation(model)
dΩ = Measure(Ω,2*order)

el2 = sqrt(sum( ∫( e⋅e )*dΩ ))
@test el2 < 1.0e-10

#using Gridap.Visualization
#
#writevtk(Ω,"trian",nsubcells=10,cellfields=["uh"=>uh])

order = 1

reffe = ReferenceFE(TET,raviart_thomas,order)

domain =(0,1,0,1,0,1)
partition = (3,3,3)
model = simplexify(CartesianDiscreteModel(domain,partition))

labels = get_face_labeling(model)
dir_tags = Array{Integer}(undef,0)

V = FESpace(model,reffe,conformity=DivConformity())

v(x) = VectorValue(-0.5*x[1]+1.0,-0.5*x[2],-0.5*x[3])
vh = interpolate(v,V)

e = v - vh

Ω = Triangulation(model)
dΩ = Measure(Ω,2*order)

el2 = sqrt(sum( ∫( e⋅e )*dΩ ))
@test el2 < 1.0e-10

#using Gridap.Visualization
#
#writevtk(trian,"test",order=3,cellfields=["vh"=>vh])

end # module


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

@doc raw"""Status is a return value for calls that don't return other objects.

    IoK8sApimachineryPkgApisMetaV1Status(;
        apiVersion=nothing,
        code=nothing,
        details=nothing,
        kind=nothing,
        message=nothing,
        metadata=nothing,
        reason=nothing,
        status=nothing,
    )

    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - code::Int32 : Suggested HTTP return code for this status, 0 if not set.
    - details::IoK8sApimachineryPkgApisMetaV1StatusDetails : Extended data associated with the reason.  Each reason may define its own extended details. This field is optional and the data returned is not guaranteed to conform to any schema except that defined by the reason type.
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - message::String : A human-readable description of the status of this operation.
    - metadata::IoK8sApimachineryPkgApisMetaV1ListMeta : Standard list metadata. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - reason::String : A machine-readable description of why this operation is in the \"Failure\" status. If this value is empty there is no information available. A Reason clarifies an HTTP status code but does not override it.
    - status::String : Status of the operation. One of: \"Success\" or \"Failure\". More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#spec-and-status
"""
mutable struct IoK8sApimachineryPkgApisMetaV1Status <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    code::Any # spec type: Union{ Nothing, Int32 } # spec name: code
    details::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1StatusDetails } # spec name: details
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ListMeta } # spec name: metadata
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason
    status::Any # spec type: Union{ Nothing, String } # spec name: status

    function IoK8sApimachineryPkgApisMetaV1Status(;apiVersion=nothing, code=nothing, details=nothing, kind=nothing, message=nothing, metadata=nothing, reason=nothing, status=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1Status, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApimachineryPkgApisMetaV1Status, Symbol("code"), code)
        setfield!(o, Symbol("code"), code)
        validate_property(IoK8sApimachineryPkgApisMetaV1Status, Symbol("details"), details)
        setfield!(o, Symbol("details"), details)
        validate_property(IoK8sApimachineryPkgApisMetaV1Status, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApimachineryPkgApisMetaV1Status, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(IoK8sApimachineryPkgApisMetaV1Status, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApimachineryPkgApisMetaV1Status, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        validate_property(IoK8sApimachineryPkgApisMetaV1Status, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1Status

const _property_map_IoK8sApimachineryPkgApisMetaV1Status = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("code")=>Symbol("code"), Symbol("details")=>Symbol("details"), Symbol("kind")=>Symbol("kind"), Symbol("message")=>Symbol("message"), Symbol("metadata")=>Symbol("metadata"), Symbol("reason")=>Symbol("reason"), Symbol("status")=>Symbol("status"))
const _property_types_IoK8sApimachineryPkgApisMetaV1Status = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("code")=>"Int32", Symbol("details")=>"IoK8sApimachineryPkgApisMetaV1StatusDetails", Symbol("kind")=>"String", Symbol("message")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ListMeta", Symbol("reason")=>"String", Symbol("status")=>"String")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1Status }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1Status))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1Status }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1Status[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1Status }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1Status[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1Status)
    true
end

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


text:
# Autogenerated wrapper script for SPIRV_LLVM_Translator_jll for powerpc64le-linux-gnu-cxx11
export libLLVMSPIRV, llvm_spirv

JLLWrappers.@generate_wrapper_header("SPIRV_LLVM_Translator")
JLLWrappers.@declare_library_product(libLLVMSPIRV, "libLLVMSPIRVLib.so.9jl")
JLLWrappers.@declare_executable_product(llvm_spirv)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libLLVMSPIRV,
        "lib/libLLVMSPIRVLib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        llvm_spirv,
        "bin/llvm-spirv",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function jamsd_reg_eqns(ctx, m::JAMSDMathProgModel, equil=false)
    if has_objective(m)
        obj_offset = 1
    else
        obj_offset = 0
    end

    jamsd_declare_eqns(ctx, m, equil)

    if has_objective(m)
        jamsd_add_obj_nl(ctx, m)
        jamsd_add_obj_lin(ctx, m)
        jamsd_add_obj_quad(ctx, m)
    end

    if m.ncon > 0 || length(m.quad_equs) > 0
#        write_nl_k_block(f, m)
        if length(m.nonquad_idx) != m.ncon
            error("The number of constraint and the equation index do not match")
        end
        jamsd_add_cons_lin(ctx, m, obj_offset)
        jamsd_add_cons_nl(ctx, m, obj_offset)
        jamsd_add_cons_quad(ctx, m, obj_offset)
        jamsd_add_contraint_sense(ctx, m, obj_offset)
    end

    # This is most likely broken, we need another way of getting this data if
    # it is important
#    return m.ncon+obj_offset+length(m.quad_equs)
end

function create_jamsd_ctx(m::JAMSDMathProgModel)
    # Not all model have an objective function
    if has_objective(m)
        obj_offset = 1
    else
        obj_offset = 0
    end

    ctx = ctx_create(m.nvar, length(m.quad_equs)+m.ncon+obj_offset)

    jamsd_declare_vars(ctx, m)
    jamsd_add_var_guess(ctx, m)

    jamsd_reg_eqns(ctx, m)

    return ctx
end

function has_objective(m::JAMSDMathProgModel)
    CONFIG[:debug] && println("DEBUG: has_objective $(m.lin_obj) $(m.obj) $(m.quad_obj)")
    res = !isempty(m.lin_obj) || (isa(m.obj, Expr)) || length(m.quad_obj) == 3
    if !res return res end
    res = any(v != 0. for v in values(m.lin_obj)) || isa(m.obj, Expr) && (m.obj != 0.) || length(m.quad_obj) == 3
    return res
end

# Nonlinear constraint trees
function jamsd_add_cons_nl(ctx, m::JAMSDMathProgModel, offset)
    for idx in 1:m.ncon
        if isa(m.constrs[idx], Expr) && (m.constrs[idx] != 0.)
            eidx = m.nonquad_idx[idx] + offset - 1 + m.offset
            tree, node = jamsd_get_treedata(ctx, eidx)
            CONFIG[:debug] && println("jamsd_add_cons_nl: storing expression $(eidx)")
            jamsd_add_nlexpr(ctx, tree, node, m, m.constrs[idx])
        end
    end
end

function jamsd_quad(ctx, m::JAMSDMathProgModel, idx, equ, offset, isObj::Bool=false)
    ##########################################################################
    #
    ##########################################################################
    if isObj
        eidx = m.offset
    else
        eidx = m.quad_idx[idx] - 1 + m.offset + offset
    end

    CONFIG[:debug] && println("jamsd_quad: $equ")
    lidx, lval, rowidx, colidx, qval = equ[1:5]
    CONFIG[:debug] && println("jamsd_quad: $lidx $lval $rowidx $colidx $qval")

    ##########################################################################
    # lin_dict collects all the linear terms
    #
    # we don't construct the dict directly, since it may have duplicate values
    ##########################################################################

    lin_dict = Dict(zip(lidx, Iterators.repeated(0.)))
    for i in 1:length(lidx)
        lin_dict[lidx[i]] += lval[i]
    end

    ###########################################################################
    # MPB is a bit insane here: via setquadobj! it provides a triangular matrix
    # representing .5 x^TQx ... But this makes no sense since we assume that
    # the multiplication is commutative here. Therefore the values of the off-
    # diagonal element have to be multiplied by 2
    ###########################################################################

    quad_dict = Dict(zip(zip(rowidx, colidx), Iterators.repeated(0.)))

    if isObj
        for i in 1:length(rowidx)
            if rowidx[i] != colidx[i]
                quad_dict[(rowidx[i], colidx[i])] += 2*qval[i]
            else
                quad_dict[(rowidx[i], colidx[i])] += qval[i]
            end
        end
    else
        for i in 1:length(rowidx)
            quad_dict[(rowidx[i], colidx[i])] += qval[i]
        end
    end

    qidxC = keys(quad_dict)
    rowidxU = [elt[1]-1 for elt in qidxC]
    colidxU = [elt[2]-1 for elt in qidxC]
    qvalU = collect(values(quad_dict))

    ##########################################################################
    # Add the linear part <c,x> for the purely linear variables
    ##########################################################################

    lidxS = Set(lidx)
    lin_vars = setdiff(lidxS, union(Set(rowidx), Set(colidx)))
    qidxS = setdiff(lidxS, lin_vars)
    for lindx in lin_vars
        # add as linear variable
        v = lin_dict[lindx]
        if abs(v) > 0.
            CONFIG[:debug] && println("DEBUG: in jamsd_quad, lin var $lindx with value $v")
            ctx_add_lin_var(ctx, eidx, m.v_index_map[lindx], v)
        end
    end

    ##########################################################################
    # Add the term <c,x> for the quadratic variables
    ##########################################################################

    if length(lin_dict) > 0 && length(qidxS) > 0
        # TODO(xhub) filter for qvals[i] = 0.
        qvals = collect(lin_dict[qidx] for qidx in qidxS)
        CONFIG[:debug] && println("DEBUG: in jamsd_quad, qidxS = $qidxS; qval= $qvals")
        abs_var = jamsd_avar(length(qidxS), collect(qidxS) .- 1)
        CONFIG[:debug] && println("DEBUG: in jamsd_quad, quad var $(collect(qidxS)) with value $qvals")
        jamsd_equ_add_lin_tree(ctx, eidx, qvals, abs_var, 1.)
        jamsd_avar_free(abs_var)
    end

    ##########################################################################
    # Add the quadratic terms <x, Mx>
    ##########################################################################

    if length(qvalU) > 0
        CONFIG[:debug] && println("DEBUG: in jamsd_quad, quad term $rowidxU $colidxU $qvalU")
        @assert length(qvalU) == length(rowidxU)
        @assert length(qvalU) == length(colidxU)
        mat = jamsd_mat_coo(rowidxU, colidxU, qvalU)
        midxS = union(BitSet(rowidx), BitSet(colidx))
        avar = jamsd_avar(length(midxS), collect(i-1 for i in midxS))
        if isObj
            c = 1.
        else
            c = 2.
        end
        # the coeff is 2 since this function adds .5 x^TMx
        jamsd_equ_add_quadratic(ctx, eidx, mat, avar, c)
        jamsd_avar_free(avar)
        jamsd_mat_free(mat)
    end
end

function jamsd_add_obj_quad(ctx, m::JAMSDMathProgModel)
    CONFIG[:debug] && println("DEBUG: quad_obj = $(m.quad_obj)")
    if length(m.quad_obj) == 3
        rowidx = m.quad_obj[1]
        colidx = m.quad_obj[2]
        quadval = m.quad_obj[3]
        jamsd_quad(ctx, m, 0, (Int[], Float64[], rowidx, colidx, quadval), 0, true)
    elseif length(m.quad_obj) == 0
        #doing noting here
    else
        error("jamsd_add_obj_quad :: invalid quad_obj object $(m.quad_obj)")
    end
end

function jamsd_add_cons_quad(ctx, m::JAMSDMathProgModel, offset)
    for (idx, equ) in enumerate(m.quad_equs)
        jamsd_quad(ctx, m, idx, equ, offset)
    end
end

# Nonlinear objective tree
function jamsd_add_obj_nl(ctx, m::JAMSDMathProgModel)
    # Get tree, node, ...
    if m.obj != 0.
        tree, node = jamsd_get_treedata(ctx, m.offset)
        CONFIG[:debug] && println("jamsd_add_obj_nl: storing expression $(m.offset)")
        jamsd_add_nlexpr(ctx, tree, node, m, m.obj)
    end
end

# Initial primal guesses
function jamsd_add_var_guess(ctx, m::JAMSDMathProgModel)
    for idx in 0:(m.nvar - 1)
        i = m.v_index_map_rev[idx]
        CONFIG[:debug] && println("Setting level value variable $i to $(m.x_0[i])")
        ctx_setvarlone(ctx, idx, m.x_0[i])
    end
end

# Constraint bounds
function jamsd_add_contraint_sense(ctx, m::JAMSDMathProgModel, offset)
    for idx in 1:m.ncon
        lower = m.g_l[idx]
        upper = m.g_u[idx]
        rel = m.r_codes[idx]
        # TODO(xhub) document
        if rel == -1
            error("Doubly constraint equation is not yet supported")
        elseif rel == 0         # ==
            value = lower
        elseif rel == 1         # >=
            value = lower
        elseif rel == 2         # <=
            value = upper
        elseif rel == 4         # ???
            value == upper
        else
            error("unsupported rel = $(rel) for equation index $i")
        end
        eidx = m.nonquad_idx[idx] + offset - 1 + m.offset
        CONFIG[:debug] && println("Setting sense and rhs for equation $eidx: $rel $value")
        jamsd_set_rhs(ctx, eidx, value)
        jamsd_set_equtype(ctx, eidx, rel)
    end

    for (idx, equ) in enumerate(m.quad_equs)
        rel, value = equ[end-1:end]
        eidx = m.quad_idx[idx] + offset - 1 + m.offset
        jamsd_set_rhs(ctx, eidx, value)
        jamsd_set_equtype(ctx, eidx, relation_to_jamsd[rel])
    end
end

function jamsd_declare_var(ctx, vtype, lower, upper)
    if lower == -Inf
        if upper == Inf
            jamsd_add_free_var(ctx, 1)
        elseif upper == 0
            jamsd_add_neg_var(ctx, 1)
        else
            jamsd_add_box_var(ctx, lower, upper)
        end
    else
        if lower == upper
            jamsd_add_box_var(ctx, lower, upper)
        elseif upper == Inf
            if lower == 0
                jamsd_add_pos_var(ctx, 1)
            else
                jamsd_add_box_var(ctx, lower, upper)
            end
        else
            jamsd_add_box_var(ctx, lower, upper)
        end
    end
    idx = hack_last_vidx(ctx)
    CONFIG[:debug] && println("DEBUG: var $(idx) has type $(vtype)")
    # We have no consistency check
    # TODO SOS, semicont
    if vtype == :Bin
        jamsd_set_vartype(ctx, idx, 1)
    elseif vtype == :Int
        jamsd_set_vartype(ctx, idx, 2)
    end

end

# Variable bounds
function jamsd_declare_vars(ctx, m::JAMSDMathProgModel)
    for idx in 0:(m.nvar - 1)
        i = m.v_index_map_rev[idx]
        lower = m.x_l[i]
        upper = m.x_u[i]
        jamsd_declare_var(ctx, m.vartypes[i], lower, upper)
    end

    # Set variable names
    if m.d.hasvalue
        ctx_setvarnames(ctx, m.d.value.m.colNames)
    end
end

# Jacobian counts
# TODO(xhub) resuse this to prealloc the model_repr?
#function write_nl_k_block(f, m::JAMSDMathProgModel)
#    println(f, "k$(m.nvar - 1)")
#    total = 0
#    for index = 0:(m.nvar - 2)
#        i = m.v_index_map_rev[index]
#        total += m.j_counts[i]
#        println(f, total)
#    end
#end

function jamsd_declare_eqns(ctx, m::JAMSDMathProgModel, equil=false)
    if has_objective(m)
        cor = 0
    else
        cor = -1
    end

    for idx in 0:m.ncon+cor+length(m.quad_equs)
        jamsd_decl_eqn(ctx, idx+m.offset)
    end
    if equil || !has_objective(m)
        jamsd_set_objeqn(ctx, -1)
    else
        jamsd_set_objeqn(ctx, 0)
    end
#    return m.ncon+cor+length(m.quad_equs)
end

# Linear constraint expressions
function jamsd_add_cons_lin(ctx, m::JAMSDMathProgModel, offset)
    for idx in 1:m.ncon
        num_vars = length(m.lin_constrs[idx])
        if num_vars > 0
#TODO(Xhub)            jamsd_reserve_lequ(ctx, idx, num_vars)
            eidx = m.nonquad_idx[idx] + offset - 1 + m.offset
            for (k,v) in m.lin_constrs[idx]
                if abs(v) > 0.
                    CONFIG[:debug] && println("DEBUG: add var $k with value $v in eqn $(eidx)")
                    ctx_add_lin_var(ctx, eidx, m.v_index_map[k], v)
                end
            end
        end
    end
end

# Linear objective expression
function jamsd_add_obj_lin(ctx, m::JAMSDMathProgModel)
    for (k,v) in m.lin_obj
        if abs(v) > 0.
            ctx_add_lin_var(ctx, m.offset, m.v_index_map[k], v)
        end
    end
end

function jamsd_add_quad_lin(ctx, m::JAMSDMathProgModel)
end

##############################################################################
# Low level functions
##############################################################################

function write_arithm_op(ctx, tree, node, m, fn, args)
    # TODO(xhub) implement better management of variable/constant base
    len = length(args)
    equtree_arithm(tree, node, arithm_ops[fn], len)
    pnode = equnode_deref(node)
    for i in 0:(len-1)
        child = equnode_get_child_addr(pnode, i)
        jamsd_add_nlexpr(ctx, tree, child, m, args[len-i])
    end
end

# Convert an expression tree
jamsd_add_nlexpr(ctx, tree, node, m, c) = error("Unrecognized expression $c")
# Handle numerical constants e.g. pi
jamsd_add_nlexpr(ctx, tree, node, m, c::Symbol) =  jamsd_add_nlexpr(ctx, tree, node, m, float(eval(c)))

# write down constant
function jamsd_add_nlexpr(ctx, tree, node, m, c::Real)
#    if abs(c) > 0.
        equtree_cst(ctx, tree, node, c)
#    end
end

jamsd_add_nlexpr(ctx, tree, node, m, c::LinearityExpr) = jamsd_add_nlexpr(ctx, tree, node, m, c.c)
function jamsd_add_nlexpr(ctx, tree, node, m, c::Expr)
    CONFIG[:debug] && println("DEBUG: encoding expr $c")
    if c.head == :ref
        CONFIG[:debug] && println("DEBUG: variable case: $c")
        # This is a variable
        if c.args[1] == :x
            @assert isa(c.args[2], Int)
            equtree_var(ctx, tree, node, m.v_index_map[c.args[2]], 1.)
        else
            error("Unrecognized reference expression $c")
        end
    elseif c.head == :call
        if c.args[1] in keys(arithm_ops)
            write_arithm_op(ctx, tree, node, m, c.args[1], c.args[2:end])
        elseif c.args[1] == :neg
            equtree_umin(ctx, tree, node)
            jamsd_add_nlexpr(ctx, tree, node, m, c.args[2])
        # TODO(xhub) support nary_functions
        elseif c.args[1] in nary_functions
            error("Unsupported n-ary function $(c.args[1])")
        else
            equtree_call(ctx, tree, node, func_to_jamsd[c.args[1]])
            len = length(c.args)
            pnode = equnode_deref(node)

            for i in 0:(len - 2)
                child = equnode_get_child_addr(pnode, i)
                jamsd_add_nlexpr(ctx, tree, child, m, c.args[i+2])
            end
        end

    elseif c.head == :comparison
        # .nl only handles binary comparison
        @assert length(c.args) == 3
        # Output comparison type first, followed by args
        CONFIG[:debug] && println(f, nl_operator(c.args[2]))
        map(arg -> jamsd_add_nlexpr(ctx, tree, node, m, arg), c.args[1:2:end])

    elseif c.head in [:&&, :||]
        error("unsupported binary condition")
    else
        error("Unrecognized expression $c")
    end
end


text:
#############################################################################
#############################################################################
#
# This file implements polynomial addition
#
#############################################################################
#############################################################################

"""
Add two polynomials.
"""
function +(p1::Polynomial, p2::Polynomial)::Polynomial #same as with polynomials not mod p
    p1, p2 = deepcopy(p1), deepcopy(p2)
    p3 = Term[]
    while !iszero(p1) && !iszero(p2)
        t1, t2 = leading(p1), leading(p2)
        if t1.degree == t2.degree
            push!(p3, (prepop!(p1)+prepop!(p2)))
        elseif t1.degree < t2.degree
            push!(p3,pop!(p2))
        else
            push!(p3,pop!(p1))
        end
    end
    while !iszero(p1)
        push!(p3,pop!(p1))
    end
    while !iszero(p2)
        push!(p3,pop!(p2))
    end

    return Polynomial(p3)
end

function plus(p1::Polynomial, p2::Polynomial)::Polynomial
    iszero(p1) && return p2
    iszero(p2) && return p1
    (iszero(p1) && iszero(p1)) && return p2

    degree(p1) > degree(p2) ? n = degree(p1)+1 : n = degree(p2)+1

    p1 = coeffvector(p1,n)
    p2 = coeffvector(p2,n)

    output = [Term(0,0) for _ in 1:n]
    for i in 1:n
        output[i] = Term(p1[i]+p2[i],n-i)
    end

    return Polynomial(output,true)

end

function +(p1::PolynomialModP, p2::PolynomialModP)::PolynomialModP #New, fast addition method
    @assert p1.mod == p2.mod
    p = p1.mod
    iszero(p1) && return p2
    iszero(p2) && return p1
    (iszero(p1) && iszero(p1)) && return p2

    #Same basic structure as old method accept that the polynomil are turned into vectors

    degree(p1) > degree(p2) ? n = degree(p1)+1 : n = degree(p2)+1

    p1 = coeffvector(p1.polynomial,n)
    p2 = coeffvector(p2.polynomial,n)

    output = [Term(0,0) for _ in 1:n]
    i = 1
    for j in 1:n
        c = mod(p1[j]+p2[j],p)
        if c == 0
            popat!(output,i)
        else
            output[i] = Term(c,n-j)
            i = i+1
        end

    end
    poly = Polynomial(output,true)
    return PolynomialModP(poly,p,true)

end

function plus(p1::PolynomialModP, p2::PolynomialModP)::PolynomialModP #old, slow addition method
    @assert p1.mod == p2.mod
    p = p1.mod
    p1, p2 = deepcopy(p1), deepcopy(p2)
    p3 = Term[]
    while !iszero(p1) && !iszero(p2)
        t1, t2 = leading(p1), leading(p2)
        if t1.degree == t2.degree
            push!(p3, mod((prepop!(p1)+prepop!(p2)),p))
        elseif t1.degree < t2.degree
            push!(p3,pop!(p2))
        else
            push!(p3,pop!(p1))
        end
    end
    while !iszero(p1)
        push!(p3,pop!(p1))
    end
    while !iszero(p2)
        push!(p3,pop!(p2))
    end

    return PolynomialModP(p3,p)
end

"""
Add a polynomial and a term.
"""
+(p::Polynomial, t::Term) = p + Polynomial(t)
+(t::Term, p::Polynomial) = p + t

+(p::PolynomialModP, t::Term) = p + PolynomialModP(t,p.mod)
+(t::Term, p::PolynomialModP) = p.polynomial + t

"""
Add a polynomial and an integer.
"""
+(p::Polynomial, n::Int) = p + Term(n,0)
+(n::Int, p::Polynomial) = p + Term(n,0)

+(p::PolynomialModP, n::Int) = p + Term(n,0)
+(n::Int, p::PolynomialModP) = p + Term(n,0)


text:
using RLEVectors

function disjoinA(x::RLEVector, y::RLEVector)
    length(x) != length(y) && error("RLEVectors of unequal length.")
    runends = disjoin(x.runends, y.runends)
    runvalues_x = x[runends]
    runvalues_y = y[runends]
    return( (runends, runvalues_x, runvalues_y ) )
end

function disjoinB(x::RLEVector, y::RLEVector)
    length(x) != length(y) && error("RLEVectors of unequal length.")
    runends = disjoin(x.runends, y.runends)
    #nrun = disjoin_length(x.runends, y.runends) # Alternative
    #nx = nrun(x)
    #ny = nrun(y)
    xv = x.runvalues
    yv = y.runvalues
    xe = x.runends
    ye = y.runends
    runends = disjoin(xe, ye)
    nrun = length(runends)

    runvalues_x = Array(eltype(x), nrun)
    runvalues_y = Array(eltype(x), nrun)
    i = j = 1
    while true
        if xe[i] > ye[j]
            println("i bigger")
            current_end = runends[runind] = xe[i]
#            runvalues_x[runind] = xv[i]
#            runvalues_y[runind] = yv[j]
#            i = i - 1
        elseif xe[i] < ye[j]
            println("j bigger")
            runends[runind] = ye[j]
#            runvalues_x[runind] = xv[i]
#            runvalues_y[runind] = yv[j]
#            j = j - 1
        else
            println("tie")
            runends[runind] = xe[i]
#            runvalues_x[runind] = xv[i]
#            runvalues_y[runind] = yv[j]
#            i = i - 1
#            j = j - 1
        end

        runind = runind - 1
        println(runvalues_x)
        println(runvalues_y)
    end
    return( (runends, runvalues_x, runvalues_y ) )
end

function disjoinC(x::RLEVector, y::RLEVector)
    length(x) != length(y) && error("RLEVectors of unequal length.")
    runends = disjoin(x.runends, y.runends)
    runvalues_x = x[runends]
    runvalues_y = y[runends]
#    #nrun = disjoin_length(x.runends, y.runends) # Alternative
#    #nx = nrun(x)
#    #ny = nrun(y)
#    xv = x.runvalues
#    yv = y.runvalues
#    xe = x.runends
#    ye = y.runends
#    runends = disjoin(xe, ye)
#    nrun = length(runends)
#
#    runvalues_x = Array(eltype(x), nrun)
#    runvalues_y = Array(eltype(x), nrun)
#    runvalues_x[:] = 0 # debug
#    runvalues_y[:] = 0 # debug
#    i = j = 1
#    while true
#        if xe[i] > ye[j]
#            println("i bigger")
#            current_end = runends[runind] = xe[i]
##            runvalues_x[runind] = xv[i]
##            runvalues_y[runind] = yv[j]
##            i = i - 1
#        elseif xe[i] < ye[j]
#            println("j bigger")
#            runends[runind] = ye[j]
##            runvalues_x[runind] = xv[i]
##            runvalues_y[runind] = yv[j]
##            j = j - 1
#        else
#            println("tie")
#            runends[runind] = xe[i]
##            runvalues_x[runind] = xv[i]
##            runvalues_y[runind] = yv[j]
##            i = i - 1
##            j = j - 1
#        end
#
#        runind = runind - 1
#        println(runvalues_x)
#        println(runvalues_y)
#    end
    return( (runends, runvalues_x, runvalues_y ) )
end

foo = IntegerRle(Vector{Int32}(collect(1:1000)), Vector{Int32}(collect(5:5:5000)))


text:
include("dependencies_for_runtests.jl")

using MPI

# # Distributed model tests
#
# These tests are meant to be run on 4 ranks. This script may be run
# stand-alone (outside the test environment) via
#
# mpiexec -n 4 julia --project test_distributed_models.jl
#
# provided that a few packages (like TimesDates.jl) are in your global environment.
#
# Another possibility is to use tmpi ():
#
# tmpi 4 julia --project
#
# then later:
#
# julia> include("test_distributed_models.jl")
#
# When running the tests this way, uncomment the following line

MPI.Init()

# to initialize MPI.

using Oceananigans.Distributed: reconstruct_global_grid
using Oceananigans.Distributed: ZXYPermutation, ZYXPermutation

@kernel function set_distributed_solver_input!(permuted_ϕ, ϕ, ::ZYXPermutation)
    i, j, k = @index(Global, NTuple)
    @inbounds permuted_ϕ[k, j, i] = ϕ[i, j, k]
end

@kernel function set_distributed_solver_input!(permuted_ϕ, ϕ, ::ZXYPermutation)
    i, j, k = @index(Global, NTuple)
    @inbounds permuted_ϕ[k, i, j] = ϕ[i, j, k]
end

function random_divergent_source_term(grid)
    # Generate right hand side from a random (divergent) velocity field.
    Ru = XFaceField(grid)
    Rv = YFaceField(grid)
    Rw = ZFaceField(grid)
    U = (u=Ru, v=Rv, w=Rw)

    Nx, Ny, Nz = size(grid)
    set!(Ru, (x, y, z) -> rand())
    set!(Rv, (x, y, z) -> rand())
    set!(Rw, (x, y, z) -> rand())

    arch = architecture(grid)
    fill_halo_regions!(Ru)
    fill_halo_regions!(Rv)
    fill_halo_regions!(Rw)

    # Compute the right hand side R = ∇⋅U
    ArrayType = array_type(arch)
    R = zeros(Nx, Ny, Nz) |> ArrayType
    event = launch!(arch, grid, :xyz, divergence!, grid, U.u.data, U.v.data, U.w.data, R,
                    dependencies=Event(device(arch)))
    wait(device(arch), event)

    return R
end

function divergence_free_poisson_solution_triply_periodic(grid_points, ranks)
    topo = (Periodic, Periodic, Periodic)
    arch = MultiArch(CPU(), ranks=ranks, topology = topo)
    local_grid = RectilinearGrid(arch, topology=topo, size=grid_points, extent=(1, 2, 3))

    bcs = FieldBoundaryConditions(local_grid, (Center, Center, Center))
    bcs = inject_halo_communication_boundary_conditions(bcs, arch.local_rank, arch.connectivity)

    # The test will solve for ϕ, then compare R to ∇²ϕ.
    ϕ   = CenterField(local_grid, boundary_conditions=bcs)
    ∇²ϕ = CenterField(local_grid, boundary_conditions=bcs)
    R   = random_divergent_source_term(local_grid)

    global_grid = reconstruct_global_grid(local_grid)
    solver = DistributedFFTBasedPoissonSolver(global_grid, local_grid)

    # Solve it
    ϕc = first(solver.storage)

    event = launch!(arch, local_grid, :xyz,
                    set_distributed_solver_input!, ϕc, R, solver.input_permutation,
                    dependencies = device_event(arch))

    wait(device(arch), event)

    solve!(ϕ, solver)

    # "Recompute" ∇²ϕ
    compute_∇²!(∇²ϕ, ϕ, arch, local_grid)

    return R ≈ interior(∇²ϕ)
end

@testset "Distributed FFT-based Poisson solver" begin
    @info "  Testing 3D distributed FFT-based Poisson solver..."
    @test divergence_free_poisson_solution_triply_periodic((44, 44, 8), (1, 4, 1))
    @test divergence_free_poisson_solution_triply_periodic((44, 16, 8), (1, 4, 1))
    @test divergence_free_poisson_solution_triply_periodic((16, 44, 8), (1, 4, 1))
    @test divergence_free_poisson_solution_triply_periodic((44, 16, 8), (2, 2, 1))
    @test divergence_free_poisson_solution_triply_periodic((16, 44, 8), (2, 2, 1))

    @info "  Testing 2D distributed FFT-based Poisson solver..."
    @test divergence_free_poisson_solution_triply_periodic((44, 16, 1), (1, 4, 1))
    @test divergence_free_poisson_solution_triply_periodic((44, 16, 1), (4, 1, 1))
    @test divergence_free_poisson_solution_triply_periodic((16, 44, 1), (1, 4, 1))
    @test divergence_free_poisson_solution_triply_periodic((16, 44, 1), (4, 1, 1))

    @test_throws ArgumentError divergence_free_poisson_solution_triply_periodic((16, 44, 1), (2, 2, 1))
    @test_throws ArgumentError divergence_free_poisson_solution_triply_periodic((44, 16, 1), (2, 2, 1))
end


text:
function summarize_regression_test(fields, correct_fields)
    for (field_name, φ, φ_c) in zip(keys(fields), fields, correct_fields)
        Δ = φ .- φ_c

        Δ_min      = minimum(Δ)
        Δ_max      = maximum(Δ)
        Δ_mean     = mean(Δ)
        Δ_abs_mean = mean(abs, Δ)
        Δ_std      = std(Δ)

        matching    = sum(φ .≈ φ_c)
        grid_points = length(φ_c)

        @info @sprintf("Δ%s: min=%+.6e, max=%+.6e, mean=%+.6e, absmean=%+.6e, std=%+.6e (%d/%d matching grid points)",
                       field_name, Δ_min, Δ_max, Δ_mean, Δ_abs_mean, Δ_std, matching, grid_points)
    end
end

function get_fields_from_checkpoint(filename)
    file = jldopen(filename)

    tracers = keys(file["tracers"])
    tracers = Tuple(Symbol(c) for c in tracers)

    velocity_fields = (u = file["velocities/u/data"],
                       v = file["velocities/v/data"],
                       w = file["velocities/w/data"])

    tracer_fields =
        NamedTuple{tracers}(Tuple(file["tracers/$c/data"] for c in tracers))

    current_tendency_velocity_fields = (u = file["timestepper/Gⁿ/u/data"],
                                        v = file["timestepper/Gⁿ/v/data"],
                                        w = file["timestepper/Gⁿ/w/data"])

    current_tendency_tracer_fields =
        NamedTuple{tracers}(Tuple(file["timestepper/Gⁿ/$c/data"] for c in tracers))

    previous_tendency_velocity_fields = (u = file["timestepper/G⁻/u/data"],
                                         v = file["timestepper/G⁻/v/data"],
                                         w = file["timestepper/G⁻/w/data"])

    previous_tendency_tracer_fields =
        NamedTuple{tracers}(Tuple(file["timestepper/G⁻/$c/data"] for c in tracers))

    close(file)

    solution = merge(velocity_fields, tracer_fields)
    Gⁿ = merge(current_tendency_velocity_fields, current_tendency_tracer_fields)
    G⁻ = merge(previous_tendency_velocity_fields, previous_tendency_tracer_fields)

    return solution, Gⁿ, G⁻
end

include("regression_tests/thermal_bubble_regression_test.jl")
include("regression_tests/rayleigh_benard_regression_test.jl")
include("regression_tests/ocean_large_eddy_simulation_regression_test.jl")

@testset "Regression" begin
    @info "Running regression tests..."

    for arch in archs
        @testset "Thermal bubble [$(typeof(arch))]" begin
            @info "  Testing thermal bubble regression [$(typeof(arch))]"
            run_thermal_bubble_regression_test(arch)
        end

        @testset "Rayleigh–Bénard tracer [$(typeof(arch))]" begin
            @info "  Testing Rayleigh–Bénard tracer regression [$(typeof(arch))]"
            run_rayleigh_benard_regression_test(arch)
        end

        @testset "Ocean large eddy simulation [$(typeof(arch))]" begin
            for closure in (AnisotropicMinimumDissipation(), ConstantSmagorinsky())
                closurename = string(typeof(closure).name.wrapper)
                @info "  Testing oceanic large eddy simulation regression [$closurename, $(typeof(arch))]"
                run_ocean_large_eddy_simulation_regression_test(arch, closure)
            end
        end
    end
end


text:
module ControllerHelper

using Helpers, Renderer

end

max_stars_repo_name: anirudhmurali/hello-julia-genie

max_issues_repo_name: anirudhmurali/hello-julia-genie

max_forks_repo_name: anirudhmurali/hello-julia-genie


text:
function block_sizes(::Type{T}, _α, _β, R₁, R₂) where {T}
    W = pick_vector_width(T)
    α = _α * W
    β = _β * W
    L₁ₑ = first_cache_size(T) * R₁
    L₂ₑ = second_cache_size(T) * R₂
    block_sizes(W, α, β, L₁ₑ, L₂ₑ)
end
function block_sizes(W, α, β, L₁ₑ, L₂ₑ)
    mᵣ, nᵣ = matmul_params()
    MᵣW = mᵣ * W

    Mc = floortostaticint(√(L₁ₑ)*√(L₁ₑ*β + L₂ₑ*α)/√(L₂ₑ) / MᵣW) * MᵣW
    Kc = roundtostaticint(√(L₁ₑ)*√(L₂ₑ)/√(L₁ₑ*β + L₂ₑ*α))
    Nc = floortostaticint(√(L₂ₑ)*√(L₁ₑ*β + L₂ₑ*α)/√(L₁ₑ) / nᵣ) * nᵣ

    Mc, Kc, Nc
end
function block_sizes(::Type{T}) where {T}
    block_sizes(T, W₁Default(), W₂Default(), R₁Default(), R₂Default())
end

"""
    split_m(M, Miters_base, W)

Splits `M` into at most `Miters_base` iterations.
For example, if we wish to divide `517` iterations into roughly 7 blocks using multiples of `8`:

```julia
julia> split_m(517, 7, 8)
(72, 2, 69, 7)
```

This suggests we have base block sizes of size `72`, with two iterations requiring an extra remainder of `8 ( = W)`,
and a final block of `69` to handle the remainder. It also tells us that there are `7` total iterations, as requested.
```julia
julia> 80*2 + 72*(7-2-1) + 69
517
```
This is meant to specify roughly the requested amount of blocks, and return relatively even sizes.

This method is used fairly generally.
"""
@inline function split_m(M, _Mblocks, W)
    Miters = cld_fast(M, W)
    Mblocks = min(_Mblocks, Miters)
    Miter_per_block, Mrem = divrem_fast(Miters, Mblocks)
    Mbsize = Miter_per_block * W
    Mremfinal = M - Mbsize*(Mblocks-1) - Mrem * W
    Mbsize, Mrem, Mremfinal, Mblocks
end

"""
  solve_block_sizes(::Type{T}, M, K, N, α, β, R₂, R₃)

This function returns iteration/blocking descriptions `Mc`, `Kc`, and `Nc` for use when packing both `A` and `B`.

It tries to roughly minimize the cost
```julia
MKN/(Kc*W) + α * MKN/Mc + β * MKN/Nc
```
subject to constraints
```julia
Mc - M ≤ 0
Kc - K ≤ 0
Nc - N ≤ 0
Mc*Kc - L₁ₑ ≤ 0
Kc*Nc - L₂ₑ ≤ 0
```
That is, our constraints say that our block sizes shouldn't be bigger than the actual dimensions, and also that
our packed `A` (`Mc × Kc`) should fit into the first packing cache (generally, actually the `L₂`, and our packed
`B` (`Kc × Nc`) should fit into the second packing cache (generally the `L₃`).

Our cost model consists of three components:
1. Cost of moving data in and out of registers. This is done `(M/Mᵣ * K/Kc * N/Nᵣ)` times and the cost per is `(Mᵣ/W * Nᵣ)`.
2. Cost of moving strips from `B` pack from the low cache levels to the highest cache levels when multiplying `Aₚ * Bₚ`.
   This is done `(M / Mc * K / Kc * N / Nc)` times, and the cost per is proportional to `(Kc * Nᵣ)`.
   `α` is the proportionality-constant parameter.
3. Cost of packing `A`. This is done `(M / Mc * K / Kc * N / Nc)` times, and the cost per is proportional to
   `(Mc * Kc)`. `β` is the proportionality-constant parameter.

As `W` is a constant, we multiply the cost by `W` and absorb it into `α` and `β`. We drop it from the description
from  here on out.

In the full problem, we would have Lagrangian, with μ < 0:
f((Mc,Kc,Nc),(μ₁,μ₂,μ₃,μ₄,μ₅))
MKN/Kc + α * MKN/Mc + β * MKN/Nc - μ₁(Mc - M) - μ₂(Kc - K) - μ₃(Nc - N) - μ₄(Mc*Kc - L2) - μ₅(Kc*Nc - L3)
```julia
0 = ∂L/∂Mc = - α * MKN / Mc² - μ₁ - μ₄*Kc
0 = ∂L/∂Kc = - MKN / Kc² - μ₂ - μ₄*Mc - μ₅*Nc
0 = ∂L/∂Nc = - β * MKN / Nc² - μ₃ - μ₅*Kc
0 = ∂L/∂μ₁ = M - Mc
0 = ∂L/∂μ₂ = K - Kc
0 = ∂L/∂μ₃ = N - Nc
0 = ∂L/∂μ₄ = L₁ₑ - Mc*Kc
0 = ∂L/∂μ₅ = L₂ₑ - Kc*Nc
```
The first 3 constraints complicate things, because they're trivially solved by setting `M = Mc`, `K = Kc`, and `N = Nc`.
But this will violate the last two constraints in general; normally we will be on the interior of the inequalities,
meaning we'd be dropping those constraints. Doing so, this leaves us with:

First, lets just solve the cost w/o constraints 1-3
```julia
0 = ∂L/∂Mc = - α * MKN / Mc² - μ₄*Kc
0 = ∂L/∂Kc = - MKN / Kc² - μ₄*Mc - μ₅*Nc
0 = ∂L/∂Nc = - β * MKN / Nc² - μ₅*Kc
0 = ∂L/∂μ₄ = L₁ₑ - Mc*Kc
0 = ∂L/∂μ₅ = L₂ₑ - Kc*Nc
```
Solving:
```julia
Mc = √(L₁ₑ)*√(L₁ₑ*β + L₂ₑ*α)/√(L₂ₑ)
Kc = √(L₁ₑ)*√(L₂ₑ)/√(L₁ₑ*β + L₂ₑ*α)
Nc = √(L₂ₑ)*√(L₁ₑ*β + L₂ₑ*α)/√(L₁ₑ)
μ₄ = -K*√(L₂ₑ)*M*N*α/(L₁ₑ^(3/2)*√(L₁ₑ*β + L₂ₑ*α))
μ₅ = -K*√(L₁ₑ)*M*N*β/(L₂ₑ^(3/2)*√(L₁ₑ*β + L₂ₑ*α))
```
These solutions are indepedent of matrix size.
The approach we'll take here is solving for `Nc`, `Kc`, and then finally `Mc` one after the other, incorporating sizes.

Starting with `N`, we check how many iterations would be implied by `Nc`, and then choose the smallest value that would
yield that number of iterations. This also ensures that `Nc ≤ N`.
```julia
Niter = cld(N, √(L₂ₑ)*√(L₁ₑ*β + L₂ₑ*α)/√(L₁ₑ))
Nblock, Nrem = divrem(N, Niter)
Nblock_Nrem = Nblock + (Nrem > 0)
```
We have `Nrem` iterations of size `Nblock_Nrem`, and `Niter - Nrem` iterations of size `Nblock`.

We can now make `Nc = Nblock_Nrem` a constant, and solve the remaining three equations again:
```julia
0 = ∂L/∂Mc = - α * MKN / Mc² - μ₄*Kc
0 = ∂L/∂Kc = - MKN / Kc² - μ₄*Mc - μ₅*Ncm
0 = ∂L/∂μ₄ = L₂ₑ - Mc*Kc
```
yielding
```julia
Mc = √(L₁ₑ)*√(α)
Kc = √(L₁ₑ)/√(α)
μ₄ = -K*M*N*√(α)/L₁ₑ^(3/2)
```
We proceed in the same fashion as for `Nc`, being sure to reapply the `Kc * Nc ≤ L₂ₑ` constraint:
```julia
Kiter = cld(K, min(√(L₁ₑ)/√(α), L₂ₑ/Nc))
Kblock, Krem = divrem(K, Ki)
Kblock_Krem = Kblock + (Krem > 0)
```
This leaves `Mc` partitioning, for which, for which we use the constraint `Mc * Kc ≤ L₁ₑ` to set
the initial number of proposed iterations as `cld(M, L₁ₑ / Kcm)` for calling `split_m`.
```julia
Mbsize, Mrem, Mremfinal, Mblocks = split_m(M, cld(M, L₁ₑ / Kcm), StaticInt{W}())
```

Note that for synchronization on `B`, all threads must have the same values for `Kc` and `Nc`.
`K` and `N` will be equal between threads, but `M` may differ. By calculating `Kc` and `Nc`
independently of `M`, this algorithm guarantees all threads are on the same page.
"""
@inline function solve_block_sizes(::Type{T}, M, K, N, _α, _β, R₂, R₃, Wfactor) where {T}
    W = pick_vector_width(T)
    α = _α * W
    β = _β * W
    L₁ₑ =  first_cache_size(T) * R₂
    L₂ₑ = second_cache_size(T) * R₃

    # Nc_init = round(Int, √(L₂ₑ)*√(α * L₂ₑ + β * L₁ₑ)/√(L₁ₑ))
    Nc_init⁻¹ = √(L₁ₑ) / (√(L₂ₑ)*√(α * L₂ₑ + β * L₁ₑ))

    Niter = cldapproxi(N, Nc_init⁻¹) # approximate `ceil`
    Nblock, Nrem = divrem_fast(N, Niter)
    Nblock_Nrem = Nblock + One()#(Nrem > 0)

    ((Mblock, Mblock_Mrem, Mremfinal, Mrem, Miter), (Kblock, Kblock_Krem, Krem, Kiter)) = solve_McKc(T, M, K, Nblock_Nrem, _α, _β, R₂, R₃, Wfactor)

    (Mblock, Mblock_Mrem, Mremfinal, Mrem, Miter), (Kblock, Kblock_Krem, Krem, Kiter), promote(Nblock, Nblock_Nrem, Nrem, Niter)
end
# Takes Nc, calcs Mc and Kc
@inline function solve_McKc(::Type{T}, M, K, Nc, _α, _β, R₂, R₃, Wfactor) where {T}
    W = pick_vector_width(T)
    α = _α * W
    β = _β * W
    L₁ₑ =  first_cache_size(T) * R₂
    L₂ₑ = second_cache_size(T) * R₃

    Kc_init⁻¹ = Base.FastMath.max_fast(√(α/L₁ₑ), Nc*inv(L₂ₑ))
    Kiter = cldapproxi(K, Kc_init⁻¹) # approximate `ceil`
    Kblock, Krem = divrem_fast(K, Kiter)
    Kblock_Krem = Kblock + One()

    Miter_init = cldapproxi(M * inv(L₁ₑ), Kblock_Krem) # Miter = M * Kc / L₁ₑ
    Mbsize, Mrem, Mremfinal, Miter = split_m(M, Miter_init, W * Wfactor)
    Mblock_Mrem = Mbsize + W * Wfactor

    promote(Mbsize, Mblock_Mrem, Mremfinal, Mrem, Miter), promote(Kblock, Kblock_Krem, Krem, Kiter)
end

@inline cldapproxi(n, d⁻¹) = Base.fptosi(Int, Base.FastMath.add_fast(Base.FastMath.mul_fast(n, d⁻¹), 0.9999999999999432)) # approximate `ceil`

"""
  find_first_acceptable(M, W)

Finds first combination of `Miter` and `Niter` that doesn't make `M` too small while producing `Miter * Niter = num_cores()`.
This would be awkard if there are computers with prime numbers of cores. I should probably consider that possibility at some point.
"""
@inline function find_first_acceptable(M, W)
    Mᵣ, Nᵣ = matmul_params()
    factors = calc_factors()
    for (miter, niter) ∈ factors
        if miter * (StaticInt(2) * Mᵣ * W) ≤ M + (W + W)
            return miter, niter
        end
    end
    last(factors)
end
"""
  divide_blocks(M, Ntotal, _nspawn, W)

Splits both `M` and `N` into blocks when trying to spawn a large number of threads relative to the size of the matrices.
"""
@inline function divide_blocks(M, Ntotal, _nspawn, W)
    _nspawn == num_cores() && return find_first_acceptable(M, W)
    mᵣ, nᵣ = matmul_params()
    Miter = clamp(div_fast(M, W*mᵣ * MᵣW_mul_factor()), 1, _nspawn)
    nspawn = div_fast(_nspawn, Miter)
    if (nspawn ≤ 1) & (Miter < _nspawn)
        # rebalance Miter
        Miter = cld_fast(_nspawn, cld_fast(_nspawn, Miter))
        nspawn = div_fast(_nspawn, Miter)
    end
    Miter, cld_fast(Ntotal, max(2, cld_fast(Ntotal, nspawn)))
end