text:
# add method to permute column order (like sample order) -- make 14 histograms of how high the bar was when permuted. if the histogram is only one peak, no correlation.
#   display what correlation is
#   which ordering accentuates correlation for a particular dataset to show someone an ordering that
#   GOAL: "correlation coefficient needs to be 0.8 to see difference in ordering"
# bring two waterfall plots with different orders and show different bar heights

# shuffle columns to show correlation: traditional waterfall
#   average height among all waterfall orderings? if plotting in a particular ordering, have a silouhette of average height.
#       IF plotted height = silouhette, not correlated
#                           < silouhette, something up/
# turn the knob until it's obvious in the parallel coordinates plot.
# add some outliers
# binomial distribution would show two schools of thought -- or sum normals with bimodal distribution

"""
"""
function random_samples(value; kwargs...)
    dim = length(value)
    translation = random_translation(dim; kwargs...)
    # correlation = random_correlation(dim; kwargs...)
    return translation .+ value
end

"""

"""
function random_translation(dim;
    nsample = DEFAULT_NSAMPLE,
    distribution = DEFAULT_DISTRIBUTION,
    fuzziness = DEFAULT_FUZZINESS,
    seed = SEED,
    kwargs...,
)
    Random.seed!(SEED)
    offset = nsample==1 ? fill(0.0,dim) : random_uniform(fuzziness..., dim)
    return _random_translation(nsample, distribution, offset)
end

function _random_translation(N, distribution::Symbol, offset::Vector)
    return _random_translation(N, fill(distribution,length(offset)), offset)
end

function _random_translation(N, distribution::Vector, offset::Vector; seed=SEED)
    # If given multiple distributions to select from, do so at random.
    if length(unique(sort!(distribution)))>1
        Random.seed!(seed)
        distribution .= distribution[rand(1:length(distribution), length(offset))]
    end

    lst = [_random_translation(N, dist, off) for (dist,off) in zip(distribution,offset)]
    return hcat(lst...)'
end

function _random_translation(N, distribution::Symbol, offset::Float64)
    offset!==1.0 && Random.seed!(make_seed(offset))
    return if distribution==:uniform; random_uniform(offset, N)
    elseif distribution==:normal;     random_normal(offset, N)
    else;                             _distribution_error()
    end
end

"""
    random_uniform(x1, x2, N)
    random_uniform(dx, N)
This function returns an `N`-element Array of random values on the interval
`[x1,x2]` or `[-dx,+dx]`
"""
random_uniform(x1::Float64, x2::Float64, N) = rand(Distributions.Uniform(x1, x2), N)
random_uniform(dx::Float64, N) = random_uniform(-dx, dx, N)

"""
    random_normal(var, N)
This function returns an `N`-element Array of random values from a normal distribution
centered at `0` with a variance `var`
"""
random_normal(var::Float64, N) = rand(Distributions.Normal(0, abs(var)), N)

"""
    random_correlation(dim; kwargs...)
    random_correlation(dim, N; kwargs...)
This function returns a scaling matrix ``C \\in \\mathbb{R}^{dim\\times dim}``,
defined such that ``C C^T = R``, where ``R`` is the covariance matrix.

```math
R' = R^T R
\\
C = \\text{Cholesky}(R' / \\max(R'))
```

``R`` is the rotation matrix returned by ``random_rotation``, squared to ensure it is
positive definite, and ``C`` is computed as the lower triangular component of the Cholesky
factorization of ``R'``.

https://cran.r-project.org/web/packages/GauPro/vignettes/IntroductionToGPs.html#:~:text=Gaussian%20correlation,correlation%20function%20is%20the%20Gausian.&text=The%20parameters%20θ%3D(θ1,Gaussian%20process%20model%20to%20data
https://scipy-cookbook.readthedocs.io/items/CorrelatedRandomSamples.html
"""
function random_correlation(dim;
    nrandom = DEFAULT_NCOR,
    method = :eig,
    kwargs...,
)
    rot = random_rotation(dim;
        symmetric=true,
        type=LinearAlgebra.UnitLowerTriangular,
        kwargs...,
    )

    if method==:eig
        eigval = LinearAlgebra.eigvals(rot)
        eigvec = LinearAlgebra.eigvecs(rot)
        cor = eigvec * LinearAlgebra.diagm(sqrt.(eigval))
    else
        cor = LinearAlgebra.cholesky(rot).L
    end

    return cor
end

"""
    random_rotation(dim; kwargs...)
    random_rotation(dim, N; kwargs...)
This function returns a rotation matrix ``R \\in \\mathbb{R}^{dim\\times dim}``,
defined as ``R = I``, with ``2N`` elements filled such that:

```math
R_{i,j} = R_{j,i} = x
,\\;\\text{where}\\;
\\left\\{
    x \\;\\big\\vert\\; |x| \\in \\left[ x_{min},x_{max} \\right],\\, x\\in\\mathbb{R}^{N}
\\right\\}
```

# Keyword Arguments
- `seed=SEED` for random number generator set before picking the values with which to fill
    the rotation matrix.
- `interactivity=(0.0,0.1)` the range of minimum and maximum allowed rotation

# Returns
- `mat::Matrix` rotation matrix

# Example

```jldoctest
julia> random_rotation(4, 3)
4×4 LinearAlgebra.UnitLowerTriangular{Float64,Array{Float64,2}}:
 1.0         ⋅         ⋅    ⋅
 0.766797   1.0        ⋅    ⋅
 0.590845  -0.433763  1.0   ⋅
 0.0        0.0       0.0  1.0
```
"""
function random_rotation(dim::Integer;
    maxdim = missing,
    permutation = missing,
    interactivity = (0.0,0.1),
    nrandom = true,
    seed = SEED,
    type = I,
    kwargs...,
)
    maxdim = coalesce(maxdim, dim)
    permutation = sort(coalesce(permutation, 1:dim))

    rot = matrix(type, maxdim)
    idx = random_index(rot, nrandom)
    N = length(idx)

    # Generate random values.
    Random.seed!(seed)
    val = [
        random_uniform( interactivity..., Integer(ceil(N/2)));
        random_uniform(-reverse(collect(interactivity))..., Integer(floor(N/2)));
    ]

    rot = _fill!(rot, idx.=>val; kwargs...)

    # Now, select elements of the rotation matrix that are relevant to the given permutation.
    return rot[permutation, permutation]
end

"""
    _fill!(mat, x; kwargs...)
This is a helper function that fills and returns `mat` with the values listed in
`x = idx => val`.
"""
function _fill!(mat, x::Pair; symmetric=false, kwargs...)
    mat[x[1]...] = x[2]

    if symmetric
        mat = Matrix{}(mat)
        mat[reverse(x[1])...] = x[2]
    end

    return mat
end

function _fill!(mat, lst; kwargs...)
    [mat = _fill!(mat, x; kwargs...) for x in lst]
    return mat
end

_fill!(mat::LinearAlgebra.Diagonal, args...; kwargs...) = _fill!(Matrix{}(mat), args...; kwargs...)

"""
    random_permutation(rng::UnitRange, nperm::Int; kwargs...)
This function returns an `nperm`-element list of unique permutations of `rng`.

# Keyword Arguments
- `seed=SEED`
"""
function random_permutation(rng::UnitRange, nperm::Int; seed=SEED)
    # Ensure the input number of permutations is feasible.
    nperm = min(factorial(last(rng)), nperm)

    Random.seed!(seed)
    idx = StatsBase.sample(1:factorial(rng.stop), nperm; replace=false)
    return [Combinatorics.nthperm(rng, ii) for ii in idx]
end

"""
    random_index(x; kwargs...)
    random_index(x, N; kwargs...)
This function returns a list of random indices for which `x` is defined.

# Arguments
- `x::AbstractArray`: If `x` is a triangular matrix, return only values above/below
    (and on, if `x` is not unitary) the diagonal
- `N::Int`: number of indices to return

# Keyword Arguments
- `seed=SEED` for random number generator

# Returns
- `idx::Vector` of indices

```jldoctest
julia> random_index(lower_triangular(8), 4)
4-element Array{Array{Int64,1},1}:
 [7, 3]
 [5, 1]
 [6, 5]
 [6, 4]
```
"""
random_index(x, N::Bool; kwargs...) = N ? random_index(x; kwargs...) : []

function random_index(x, N::Int; kwargs...)
    idx = random_index(x; kwargs...)
    return idx[1:min(length(idx),N)]
end

random_index(x; kwargs...) = _shuffle(list_index(x); kwargs...)

random_index(D::LinearAlgebra.Diagonal; kwargs...) = _random_index(D, !=; kwargs...)
random_index(L::LinearAlgebra.UnitLowerTriangular; kwargs...) = _random_index(L, <; kwargs...)
random_index(U::LinearAlgebra.UnitUpperTriangular; kwargs...) = _random_index(U, >; kwargs...)
# random_index(L::LinearAlgebra.LowerTriangular; kwargs...) = _random_index(L, <=; kwargs...)
# random_index(U::LinearAlgebra.UpperTriangular; kwargs...) = _random_index(U, >=; kwargs...)

"""
This is a helper function to apply `random_index` to a matrix with a constraint
on its index (ex: triangular)
"""
function _random_index(mat, constraint; kwargs...)
    idx = _constrain_index(list_index(mat), constraint)
    return _shuffle(idx; kwargs...)
end

"This is a helper function to apply a constraint to a list of matrix indices.
At the moment, it only works for ``2\\times 2`` matrices."
function _constrain_index(idx::Vector{Vector{Int}}, fun)
    return idx[broadcast(fun, getindex.(idx,2), getindex.(idx,1))]
end

"This is a helper function to set a seed and shuffle the index in one step."
function _shuffle(idx; seed=SEED, kwargs...)
    Random.seed!(seed)
    return Random.shuffle(idx)
end

"""
    list_index(x::T) where T<:AbstractArray
This function returns a sorted list of all vector or matrix indices.

```jldoctest
julia> list_index(lower_triangular(3))
9-element Array{Array{Int64,1},1}:
 [1, 1]
 [1, 2]
 [1, 3]
 [2, 1]
 [2, 2]
 [2, 3]
 [3, 1]
 [3, 2]
 [3, 3]
```
"""
list_index(mat::T) where T<:AbstractMatrix = sort(permute(UnitRange.(1, size(mat))...))
list_index(vec::T) where T<:AbstractVector = collect(1:length(vec))

"This function returns a `len`-digit integer equal to the first four significant digits of `x`"
function make_seed(x::Float64, len=4)
    return convert(Integer, round(parse(Float64,
        Printf.@sprintf("%.10e", x)[1:len+1]) * 10^(len-1)))
end


text:
# ---
# title: 1293. Shortest Path in a Grid with Obstacles Elimination
# id: problem1293
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Breadth-first Search
# link: <https://leetcode.com/problems/shortest-path-in-a-grid-with-obstacles-elimination/description/>
# hidden: true
# ---
#
# Given a `m * n` grid, where each cell is either `0` (empty) or `1` (obstacle).
# In one step, you can move up, down, left or right from and to an empty cell.
#
# Return the minimum number of steps to walk from the upper left corner `(0, 0)`
# to the lower right corner `(m-1, n-1)` given that you can eliminate  **at
# most** `k` obstacles. If it is not possible to find such walk return -1.
#
#
#
# **Example 1:**
#
#
#
#     Input:
#     grid =
#     [[0,0,0],
#      [1,1,0],
#      [0,0,0],
#      [0,1,1],
#      [0,0,0]],
#     k = 1
#     Output: 6
#     Explanation: The shortest path without eliminating any obstacle is 10.
#     The shortest path with one obstacle elimination at position (3,2) is 6. Such path is (0,0) -> (0,1) -> (0,2) -> (1,2) -> (2,2) -> **(3,2)** -> (4,2).
#
#
#
#
# **Example 2:**
#
#
#
#     Input:
#     grid =
#     [[0,1,1],
#      [1,1,1],
#      [1,0,0]],
#     k = 1
#     Output: -1
#     Explanation: We need to eliminate at least two obstacles to find such a walk.
#
#
#
#
# **Constraints:**
#
#   * `grid.length == m`
#   * `grid[0].length == n`
#   * `1 <= m, n <= 40`
#   * `1 <= k <= m*n`
#   * `grid[i][j] == 0 **or** 1`
#   * `grid[0][0] == grid[m-1][n-1] == 0`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
#!/usr/bin/env julia

using Luxor, Colors

function drawbbox(apoly)
    gsave()
    setline(0.3)
    setdash("dotted")
    box(polybbox(apoly), :stroke)
    grestore()
end

function drawpoly(p, x, y, counter)
    gsave()
    translate(x, y)
    sethue("purple")
    psort = polysortbyangle(p)
    drawbbox(psort)
    poly(psort, close=true, :stroke)

    # label points
    for (n, pt) in enumerate(psort)
        circle(pt, 1, :fill)
        text(string(n), pt.x, pt.y)
    end

    sethue("green")
    cp = polycentroid(psort)
    circle(cp, 2, :fill)
    text(string(counter), cp.x, cp.y)

    # highlight cases where centroid isn't inside poly. Might be an error...
    if ! isinside(cp, psort)
        gsave()
        setdash("dotted")
        setopacity(0.2)
        sethue("red")
        circle(cp, 100, :fillstroke)
        grestore()
    end

    grestore()
end

function polycentroidsort(width, height)
    pagetiles = Tiler(width, height, 5, 5, margin=50)
    tilesize = pagetiles.tilewidth/2
    for (pos, n) in pagetiles
      if rand(Bool)
        p = randompointarray(rand(-tilesize:-tilesize), rand(-tilesize:-tilesize), rand(tilesize:tilesize), rand(tilesize:tilesize), rand(5:12))
      else
        p = ngon(0, 0, tilesize, rand(3:12), vertices=true)
      end
      drawpoly(p, pos.x, pos.y, n)
    end
end

width, height = 3000, 3000
fname = "/tmp/polycentroidsort.pdf"
Drawing(width, height, fname)
fontsize(20)
origin()
polycentroidsort(width, height)
finish()
println("finished test: output in $(fname)")


text:
using JuMP: @variable, @constraint, @objective, Model, with_optimizer, optimize!, value
using CPLEX

function find_feasible_point(A::Matrix{T}, b::Vector{T}, r_min::T=zero(T), r_max::T=T(Inf)) where T
	# Attempts to solve the feasibility problem
	# find       x
	# such that  Ax ≤ b
	#            r_min ≤ ‖x‖ ≤ r_max

	# First use CPLEX to find the solution x_min of
	# minimize     x'x
	# subject to   Ax ≤ b
	m, n = size(A)
	model = Model(with_optimizer(CPLEX.Optimizer))
    @variable(model, x[1:n])
    @objective(model, Min, dot(x, x))
	@constraint(model, A*x .<= b)
	optimize!(model)
    x_min = value.(x)
	# Check that x_min is feasible and ‖x_min‖ ≤ r_max
	@assert norm(x_min) <= r_max string("The problem is infeasible. Min radius:", norm(x_min))
	@assert maximum(A*x_min - b) <= 1e-11 string("Min-radius solution returned by Gurobi violates constraints by:", maximum(A*x_min - b))

	# Use GeneralQP to find a solution x
	# maximize     x'x
	# subject to   Ax ≤ b
	# that has radius less than or equal to (r_min + r_max)/2
	x = GeneralQP.solve(Matrix(-one(T)*I, n, n), zeros(T, n), A, b, x_min; r_max=(r_min+r_max)/2, printing_interval=5000)
	# Check that r_min ≤ ‖x‖ ≤ r_max and Ax ≤ b
	@assert norm(x) - r_min >= -1e-9 "Failed to find a feasible point"
	@assert norm(x) - r_max <= 1e-9 "Failed to find a feasible point"
	@assert maximum(A*x - b) <= 1e-9 "Failed to find a feasible point"

	return x
end


text:
using LinearAlgebra
using Plots

datapoints=[2.4 420;
            2.0 350;
            2.1 310;
            1.8 280;
            1.3 75]
A = hcat(ones(size(datapoints,1)),datapoints[:,1])
b = datapoints[:,2]

x = A\b

line = t -> x[1] + x[2]*t

scatter(datapoints[:,1], datapoints[:,2], legend=:none)
plot!(line, legend=:none)
ylabel!("weight, in pounds")
xlabel!("length, in meters")


text:
# call add_procs(n) to add extra processors if desired
# load packages
using KOSHELPER
@everywhere using KOS
using KOS_FIELD_SOLVE
#
# other Julia packages
@everywhere using Interpolations
using PyPlot
#
ion()
# distance function
@everywhere function d_ball_0_1( x::Array{Float64,1} )
  # out=d_ball_0_1(x)
  #
  # Distance of points to the domain boundary
  # Domain is a ball of radius 1 centred at (0,0)
  # x: point
  R=1.
  return R-sqrt(sum(x.^2))
end
#
@everywhere function vec_d_ball_0_1( xs )
  # out=d_ball_0_1(x)
  # same as above
  # points are rows of xs
  R=1.
  return (R-sqrt.(sum(xs.^2,length(size(xs)))))[:]
end
# exterior-data function g (in vector and scalar forms)
@everywhere function fn_ext( x::Array{Float64,1})
   # Exterior values
   tmp=sin((sum(x.^2)))
   return tmp::Float64
end
#
@everywhere function fn_ext( x::Array{Float64,2})
  tmp=sin.((sum(x.^2, 2)))
  return tmp[:]::Array{Float64,1}
end
#
# right-hand side function f (in vector and scalar forms)
@everywhere function fn_rhs( x::Array{Float64,1} )
  # Right-hand side designed to expose slow coupling rates.  #
  # x: point(s) as columns for evaluation
  tmp=2+sum(x.^2)
  return tmp::Float64
end
#
@everywhere function fn_rhs( x::Array{Float64,2} )
  tmp=2+sum(x.^2, 2)
  return tmp[:]::Array{Float64,1}
end
##########################################
##########################################
function mwe(alpha)
# set up example
p=KOS.Params(fn_ext,fn_rhs,vec_d_ball_0_1,alpha)
# set up params
dict=Dict("max_samples"=>Int32(1e5),
  "tol"=>1e-2,
  "no_samps"=>4000,
  "no_levels"=>7,"l0"=>5,
  "report"=>true,
  "solve_method"=>"ml")
# set up grid
R=1.2
G=KOSHELPER.get_tri_mesh([-R,R,-R,R],
                            dict["no_levels"],p.d)
# get the solution
 out,elapsed_time=get_solution(p, dict, G)
# make sure boundary conditions are applied to out
G.get_ext(out,p.fn_ext)
# put solution on a data structure suitable for plotting
Z=G.move_to_TriMesh(out)
#
if true # plot_flag
  do_plot(G,Z,"",string("fig_soln",".pdf") )
end
println("Elapsed time=", elapsed_time)
end
#################
################ end eg_field.jl
####################################


text:
"""
SimQueue type for holding arbitrary objects `O`. It allows objects to be waiting
in a queue with optional maximum queue size. Servers can get objects from the
queue with optional maximum number of requests pending for fulfillment.

Fields:
* `fifo_queue`    if `true` `queue` is FIFO, otherwise LIFO
* `max_queue`     maximum `queue` size
* `queue`         vector of objects in a queue
* `fifo_requests` if `true` `requests` is FIFO, otherwise LIFO
* `max_requests`  maximum `requests` size
* `requests`      vector of request functions

Functions in `requests` must accept two arguments `Scheduler` and `O`. When `O`
arrives to a queue there is a try to immediately feed it to pending requests.
When new request arrives there is a try to immediately provide it with `O`.

Initially an empty `SimQueue` with no requests is constructed.
By default `queue` and `requests` have FIFO policy and are unbounded.
"""
mutable struct SimQueue{O} <: AbstractReservoir
    fifo_queue::Bool
    max_queue::Int
    queue::Vector{O}
    fifo_requests::Bool
    max_requests::Int
    requests::Vector{Function}

    function SimQueue{O}(;fifo_queue::Bool=true, max_queue::Int=typemax(Int),
                         fifo_requests::Bool=true,
                         max_requests::Int=typemax(Int)) where O
        max_queue > 0 || error("max_queue must be positive")
        max_requests > 0 || error("max_requests must be positive")
        new(fifo_queue, max_queue, Vector{O}(),
            fifo_requests, max_requests, Vector{Function}())
    end
end

function dispatch!(s::Scheduler, q::SimQueue)
    # technically could be if as it should never happen that
    # the loop executes more than once, but user might tweak the internals ...
    while !isempty(q.requests) && !isempty(q.queue)
        req = pop!(q.requests)
        obj = pop!(q.queue)
        register!(s, x -> req(x, obj)) # plan to execute request immediately
    end
end

function request!(s::Scheduler, q::SimQueue, request::Function)
    length(q.requests) < q.max_requests || return false
    qend = q.fifo_requests ? pushfirst! : push!
    qend(q.requests, request)
    dispatch!(s, q)
    return true
end

function waive!(q::SimQueue, request::Function)
    idx = findfirst(isequal(request), q.requests)
    isa(idx, Nothing) && return false
    deleteat!(q.requests, idx)
    return true
end

function provide!(s::Scheduler, q::SimQueue{O}, object::O) where O
    if length(q.queue) < q.max_queue
        qend = q.fifo_queue ? pushfirst! : push!
        qend(q.queue, object)
        dispatch!(s, q)
        return true
    end
    return false
end

"""
    withdraw!(q, object)

Allows to remove first occurrence that would be served of `object` from `SimQueue`.

Returns `true` on success and `false` if `object` was not found.
"""
function withdraw!(q::SimQueue{O}, object::O) where O
    idx = findfirst(isequal(object), q.queue)
    isa(idx, Nothing) && return false
    deleteat!(q.queue, idx)
    return true
end

"""
    length(q)

Returns number of objects in the queue.
"""
Base.length(q::SimQueue) = length(q.queue)


text:
#!/usr/bin/julia
#
# Definitions of useful potential functions for use in simulations.
#
using UnitfulAtomic;
export gs_polarisation, radial_gaussian, radial_gaussian_harmonic_freq,
       tweezer_potential_depth, simple_harmonic, translate_trajectory, merge_radial_gaussian_harmonic_freq;

"""
    gs_polarisation(atom::String, λ::Number)

Find the frequency dependent ground state polarisability for Rubidium and
Caesium at the given wavelength.
"""
function gs_polarisation(atom::String, λ::Number)
  # Apply the relevant data.
  # Data from M. S. Safronova et al., Phys. Rev. A 73, 022505 (2006).
  # https://doi-org.ezphost.dur.ac.uk/10.1103/PhysRevA.73.022505
  if atom == "Cs" || atom == "cs"
    dip12 = 4.4890;
    dip32 = 6.3238;
    δE12 = 0.050932;
    δE32 = 0.053456;
    A = 17.35;
  elseif atom == "Rb" || atom == "rb"
    dip12 = 4.231;
    dip32 = 5.977;
    δE12 = 0.057314;
    δE32 = 0.058396;
    A = 10.54;
  else throw(ArgumentError("only data for Rb and Cs has been implemented")); end

  # Calculate frequency in atomic units.
  ω = ustrip(auconvert(2π*Unitful.c0/λ));

  # Return the ground state dynamic polarisability.
  α = 1/3 * ((δE12 * dip12^2)/(δE12^2 - ω^2)
             + (δE32 * dip32^2)/(δE32^2 - ω^2))
      + A;

  # Apply the atomic units.
  return α * unit(1Unitful.u"a0_au^3");
end

"""
    radial_gaussian(ω₀::Number, U₀::Number, r::AbstractVector{T}) where T<:Number

Returns the numerical values for the radial slice of a Gaussian function with
the parameters given defined over the values contained in `r`.
"""
function radial_gaussian(ω₀::Number, U₀::Number, r::AbstractVector{T}) where T<:Number
  return U₀ * (1 .- exp.(-2r.^2 / ω₀^2));
end

"""
    radial_gaussian_harmonic_freq(ω₀::Number, U₀::Number, m::Number)

Calculate the radial harmonic frequency for the Gaussian potential with the
given parameters
"""
function radial_gaussian_harmonic_freq(ω₀::Number, U₀::Number, m::Number)
  (4U₀/(m * ω₀^2))^(1/2)
end

"""
    tweezer_potential_depth(P::Number, α::Number, ω₀::Number)

Calculate the trap depth, U₀, for an optical tweezer generated by a Gaussian
laser beam with power `P`, frequency-dependent atomic ground state
polarisability of `α` (in a₀³) and beam waist `ω₀`.

Uses the SI unit conversion of `α(C m²/V) = 4π ϵ₀ a₀³ α(a₀³)` with the potential
depth given by `U₀ = (P α)/(π c ϵ₀ ω₀²)`.
"""
function tweezer_potential_depth(P::Number, α::Number, ω₀::Number)
  return (4*ustrip(α)*Unitful.uconvert(Unitful.u"m", 1Unitful.u"a0_au")^3*P)/
         (Unitful.c0 * ω₀^2);
end

"""
    simple_harmonic(m::Number, ω::Number, x::AbstractVector{T}) where T<:Number

Returns the numerical values for the simple harmonic function with the
parameters given defined over the values contained in `x`.
"""
function simple_harmonic(m::Number, ω::Number, x::AbstractVector{T}) where T<:Number
  return 0.5 * m * ω^2 * x.^2;
end

"""
    translate_trajectory(r, b, Ua, ω₀a, Ub, ω₀b)

Translate between the intended trajectory for the minimum of the radial
component of two merging Gaussian potentials and the location of the potential.
"""
function translate_trajectory(r, b, Ua, ω₀a, Ub, ω₀b)
  c = Ub/ω₀b^2 * (r-b) * exp(-2(r-b)^2/ω₀b^2);
  arg = -(2c * ω₀a/Ua)^2;
  #println("Argument: ", arg);
  if (abs(arg) > 1.0/exp(1))
    println("WARNING: outside of radius of convergence of series expansion!");
    println("Argument was ", arg);
  end
  #println(arg);

  lambert = 0.0;
  for i ∈ 1:10
    lambert += (-i)^(i-1)/factorial(i) * arg^i;
  end
  #println("Correction: ", 0.5ω₀a*sqrt(-lambert));
  return r + 0.5ω₀a*sqrt(-lambert);
end

"""
    merge_radial_gaussian_harmonic_freq(m, a, Ua, wa, b, Ub, wb)

Get the equivalent harmonic frequency of the trap in the radial direction given
the locations of the merging Gaussian potentials.
"""
function merge_radial_gaussian_harmonic_freq(m, a, Ua, wa, b, Ub, wb)
  # Imaginary values for repulsive potential, real for attractive.
  ω = 2.0*sqrt(complex(1/m*(Ub/wb^2*exp(-2(a-b)^2/wb^2)*(4(a-b)^2/wb^2 - 1) - Ua/wa^2)));
  return abs(ω);
end

function desired_trap_frequency2(γ, t, tfinal, ω0)
  s = t/tfinal;
  #println("s = ", s);
  ρ = 1 + 10*(γ-1)*s^3 - 15*(γ-1)*s^4 + 6*(γ-1)*s^5;
  #println("ρ = ", ρ);
  ρdd = (60*(γ-1)*s - 180*(γ-1)*s^2 + 120*(γ-1)*s^3)/Unitful.uconvert(unit(ω0)^-1, tfinal)^2;
  #println("ρdd = ", ρdd);
  return (ω0^2/ρ^4 - ρdd/ρ);
end

function tweezer_power_from_potential(U, w, α)
  return ((Unitful.c0 * w^2 * -U)/(4*ustrip(α)*Unitful.uconvert(Unitful.u"m", 1Unitful.u"a0_au")^3));
end

function Ua_from_ω2(ω2, m, Ub, wa, wb, a, b)
  return (wa^2*(Ub/wb^2 * exp(-2*(a-b)^2/wb^2)*((4*(a-b)^2)/wb^2 - 1) - m*ω2/4));
end

function Ub_from_ω2(ω2, m, Ua, wa, wb, a, b)
  return ((m*ω2/4 + Ua/wa^2) * wb^2/(exp(-2*(a-b)^2/wb^2)*((4*(a-b)^2)/wb^2 - 1)));
end


text:
using DSP

"""
    highlow_butterworth_filter(data,sampling_rate; low_pass=30, high_pass=1, bw_n_pole=5, offset=true)

Applies a high and low-pass filter of butterworth design (n pole 5). For altering the
threshold values for filters, change add keyword arguments low_pass for low pass filter cut-off
(default=30) and high_pass for high-pass cut-off (default=1). To change the nth order
of the butterworth filter, use keyword bw_n_pole (default=5). Offset due to filtering at 0z,
i.e the mean of the time series, can be turned off by setting `offset=false` which
manually adds back the offset

Returns filtered data
"""
function highlow_butterworth_filter(
    data, sampling_rate;
    low_pass=30,
    high_pass=1,
    bw_n_pole=5,
    offset = true,
)

# Setting up low-pass filter properties
    low_threshold = low_pass
    responsetype_low = Lowpass(low_threshold, fs=sampling_rate)
    designmethod = Butterworth(bw_n_pole)
    # Applying filter
    low_pass_filter = digitalfilter(responsetype_low, designmethod)
    lowp_filtered_data = filtfilt(low_pass_filter, data)

    # Setup high-pass filter properties
    high_threshold = high_pass
    responsetype_high = Highpass(high_threshold, fs=sampling_rate)
    # Applying filter
    high_pass_filter = digitalfilter(responsetype_high, designmethod)
    filtered_data = filtfilt(high_pass_filter, lowp_filtered_data)

    if offset == false
        # Resetting the offset
        filtered_data = filtered_data .+ mean(data, dims=1)
    end
        return filtered_data

end

"""
    highlow_butterworth_filter(data::Dict, sampling_rate; low_pass=30, high_pass=1, bw_n_pole=5)

Data contains all conditions and is a Dict (subject). Applies a high and low-pass filter of butterworth design (n pole 5). For altering the
threshold values for filters, change add keyword arguments low_pass for low pass filter cut-off
(default=30) and high_pass for high-pass cut-off (default=1). To change the nth order
of the butterworth filter, use keyword bw_n_pole (default=5).

Returns filtered data as a Dict
"""
function highlow_butterworth_filter(
    data::Dict, sampling_rate;
    low_pass=30,
    high_pass=1,
    bw_n_pole=5,
    offset=true
)

    filtered_data = Dict()
    for (condition, cond_data) in data
        filtered_data[condition] = highlow_butterworth_filter(
            cond_data,
            sampling_rate,
            low_pass=low_pass,
            high_pass=high_pass,
            bw_n_pole=bw_n_pole,
            offset=offset,
        )
    end

    return filtered_data

end

max_stars_repo_name: ElectronicTeaCup/MegTools

max_issues_repo_name: ElectronicTeaCup/MegTools

max_forks_repo_name: ElectronicTeaCup/MegTools


text: using PlotlyJS


text:
const TAIL_MIN_ASSETS = 0.
const TAIL_NUM_YEARS = 5

function performancebyein(data::NCCSData; tailminassets = TAIL_MIN_ASSETS,
  numyears = TAIL_NUM_YEARS, lastyear = 2015,
  outputpath::String = TABLE_PATH,
  identifierfields::Vector{Symbol} = [:ein, :name, :fisyr, :datefiscal],
  return1y::Symbol = :lreturnnetbmsp5001yr,
  wealthfield::Symbol = :adjnetassets,
  sumfields::Vector{Symbol} = [:monthscovered, :lt30ret, :lsp500ret,
    :bmff3, :bmff5, :bmsp500],
  annualizefields::Bool = true)::DataFrame

  #make some assertions on what must be present in the field lists
  @assert :monthscovered ∈ sumfields
  @assert :fisyr ∈ identifierfields

  local firstyear::Int = lastyear - numyears
  local fields2annualize::Vector{Symbol}

  #partition off the part we want
  validfirmsdf = view(data.df, (data.df[wealthfield] .> tailminassets) .&
    (data.df[:fisyr] .≤ lastyear) .& (data.df[:fisyr] .> firstyear),
    [identifierfields; return1y; wealthfield; sumfields;]
    )

  #form the consolidated df

  taildf::DataFrame = DataFrame(ein=validfirmsdf[:ein])
  N::Int = size(taildf, 1)

  #pre-allocate the columns
  for f ∈ setdiff(names(validfirmsdf), [:ein])
    local newtype::Type
    local origtype::Type = eltype(data.df[f])

    if Float64 <: origtype
      newtype = MFloat64
    elseif Int <: origtype
      newtype = MInt
    elseif Date <: origtype
      newtype = MDate
    elseif String <: origtype
      newtype = MString
    else
      newtype = MSymbol
    end

    println("field $f: oldtype= $(origtype) newtype = $newtype")

      taildf[f] = Vector{Union{Missing, newtype}}(missing, N)
  end
  taildf[:N] = Vector{MInt}(undef, N)

  #create the same set of fields but annualized
  #NOTE: Annualized fields are not logged
  if annualizefields
    for f ∈ [return1y; sumfields]
      f_ann::Symbol = Symbol(f, "_ann")
      taildf[f_ann] = similar(taildf[f])
    end

    fields2annualize = setdiff([return1y; sumfields;], [:monthscovered])
  end

  taildfindex::Dict = Dict(taildf[i,:ein] => taildf[i,:] for i::Int  ∈ 1:N)

  #consolidate the eprformacne data
  for subdf ∈ groupby(validfirmsdf, :ein)
    tailrow = taildfindex[subdf[1,:ein]]

    #identifier fields receive the value from the first row
    #NOTE: No verification that identifier fields are the same for a given org
    for f ∈ identifierfields
      tailrow[f] = subdf[1,f]
    end

    tailrow[return1y] = sum(subdf[return1y])

    #subtotal the sum fields
    for f ∈ sumfields
      tailrow[f] = sum(subdf[f])
    end

    #if we want annualized copies
    if annualizefields
      yearscoveredinv::MFloat64 = 1.0 / (tailrow[:monthscovered] / 12.0 )

      for f ∈ fields2annualize
        tailrow[Symbol(f,"_ann")] = exp(tailrow[f] * yearscoveredinv) - 1.0
      end
    end

    #get average assets
    tailrow[wealthfield] = mean(subdf[wealthfield]) ./ 1_000_000

    tailrow[:N] = size(subdf,1)
  end

  taildf = taildf[completecases(taildf[[identifierfields; return1y; wealthfield]]),:]
  taildf = taildf[taildf[:N] .== numyears,:]

  return taildf
end

#writes out the data to the tails
function summarizetails(data::NCCSData;
  numyears = TAIL_NUM_YEARS, lastyear = 2015,
  return1y::Symbol = :lreturnnetbmsp5001yr,
  numberpertail::Int = 200,
  outputpath::String = TABLE_PATH)::Nothing

  firstyear::Int = lastyear - numyears
  outpath = "$outputpath\\tails_$(firstyear)_$(lastyear)_$return1y"

  #get a frame with cumulative returns by year
  taildf::DataFrame = performancebyein(data, numyears=numyears, lastyear=lastyear,
    return1y=return1y)

  #organize the return dataframe
  sort!(taildf, [return1y])
  N::Int = size(taildf, 1)

  taildfbottom::SubDataFrame = view(taildf, 1:numberpertail,:)
  taildftop::SubDataFrame = view(taildf, (N-numberpertail + 1):N,:)

  #finally write out all the tail data
  taildftop |> CSV.write("$(outpath)_top.csv")
  taildfbottom |> CSV.write("$(outpath)_bottom.csv")

  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module Born

import GeoPhyInv.Acquisition
import GeoPhyInv.Models
import GeoPhyInv.Data
using SpecialFunctions
using FFTW
using DSP

mutable struct Param
end

"""
k = wavenumber = 2pif/v0
"""
function G0_homo_acou(x::T, z::T, k::T, ρ0::T) where {T<:AbstractFloat}
	sqrt(x*x + z*z) == 0.0 ? error("distance cannot be zero") : nothing
	k == 0.0 ? error("wavenumber cannot be zero") : nothing
	G0 = -0.25 * ρ0 * im * complex.(besselj0(k*sqrt(x*x + z*z)), -bessely0(k*sqrt(x*x + z*z)))
	return G0
end

function mod(;
	     vp0::Float64=2000.0,
	     ρ0::Float64=2000.0,
	     model_pert::Models.Seismic=nothing,
             born_flag::Bool=false,
	     tgridmod::StepRangeLen=nothing,
	     tgrid::StepRangeLen = tgridmod,
	     acqgeom::Acquisition.Geom=nothing,
	     acqsrc::Acquisition.Src=nothing,
	     src_flag::Int64=2,
		  )
	(getfield(acqgeom,:nss) != getfield(acqsrc,:nss))  ? error("different supersources") : nothing
	(getfield(acqgeom,:ns) != getfield(acqsrc,:ns))  ? error("different sources") : nothing

	nt = (length(tgridmod) == length(acqsrc.tgrid)) ? length(tgrid) : error("acqsrc tgrid")
	np2 = nextpow(2, 2*nt);

	if(born_flag)
		mesh_x=model_pert.mgrid[2]
		mesh_z=model_pert.mgrid[1]
		nz=length(mesh_z)
		nx=length(mesh_x)
		δx = step(mesh_x)
		δz = step(mesh_z)

		δmodtt = Models.Seismic_get(model_pert, :KI) - (vp0 * vp0 * ρ0)^(-1)
		δmodrr = Models.Seismic_get(model_pert, :ρI) - (ρ0)^(-1)
	end

	fnpow2grid = DSP.fftfreq(np2,inv(step(tgridmod)));

	nss = acqgeom.nss
	nr = acqgeom.nr
	ns = acqgeom.ns
	data = Data.TD(
	      [zeros(nt,nr[iss]) for iss=1:nss, ifield=1:1],
	      [:P],
	      tgridmod,acqgeom)

	for ifield = 1:length(data.fields), iss = 1:nss
		sx = acqgeom.sx[iss][:]
		rx = acqgeom.rx[iss][:]
		sz = acqgeom.sz[iss][:]
		rz = acqgeom.rz[iss][:]

		for ir = 1:nr[iss]

		dtemp=zeros(nt)
		dpow2all=complex.(zeros(np2), zeros(np2));
		wpow2=complex.(zeros(np2), zeros(np2));

		for is=1:acqgeom.ns[iss]
			# zero pad wavelet
			for it in 1:nt
				wpow2[it]=complex(acqsrc.wav[iss,ifield][it,is])
			end
			FFTW.fft!(wpow2) # source wavelet in the frequency domain

			x = sx[is] - rx[ir]
			z = sz[is] - rz[ir]
			# analytical expression for every frequency, except zero
			dpow2 = complex.(zeros(np2), zeros(np2));
			dpow2[1] = complex.(0.0, 0.0)
			for iω in 2:np2
				ω = 2. * pi * abs(fnpow2grid[iω])
				k = ω / vp0

				if(born_flag)
					term = complex(0., 0.)
					for ix=1:nx
						@simd for iz=1:nz
							if(δmodtt[iz,ix] ≠ 0.0)
								term += (G0_homo_acou(sx[is]-mesh_x[ix], sz[is]-mesh_z[iz], k, ρ0)[1]
									 * G0_homo_acou(rx[ir]-mesh_x[ix], rz[ir]-mesh_z[iz], k, ρ0)[1] .* ω .* ω .* δmodtt[iz,ix]) * δx * δz
								# factor due to integration
							end
						end
					end

				else
					if(src_flag == 2)
						term = G0_homo_acou(x, z, k, ρ0);
					elseif(src_flag == 1)
						dpow2[iω] = G0_homo_acou(x, z, k, ρ0) * im * abs(fnpow2grid[iω])
					else
						error("invalid src_flag")
					end
				end
				if(fnpow2grid[iω] > 0)
					dpow2[iω] = term;
				elseif(fnpow2grid[iω] < 0)
					dpow2[iω] = conj(term);
				end

			end
			# convolution and stack over simultaneous sources
			dpow2all += dpow2 .* wpow2;
		end

		# back to time domain
		FFTW.ifft!(dpow2all)

		# truncate
		for it in 1:nt
			data.d[iss,ifield][it,ir] = real(dpow2all[it])
		end
		end
	end
	return data
end

end # module


text:
# generate examples
import Literate

EXAMPLEDIR = joinpath(@__DIR__, "src", "literate")
GENERATEDDIR = joinpath(@__DIR__, "src", "examples")
mkpath(GENERATEDDIR)
for example in readdir(EXAMPLEDIR)
    if endswith(example, ".jl")
        input = abspath(joinpath(EXAMPLEDIR, example))
        script = Literate.script(input, GENERATEDDIR)
        code = strip(read(script, String))
        mdpost(str) = replace(str, "@__CODE__" => code)
        Literate.markdown(input, GENERATEDDIR, postprocess = mdpost)
        Literate.notebook(input, GENERATEDDIR, execute = true)
    else
        cp(joinpath(EXAMPLEDIR, example), joinpath(GENERATEDDIR, example); force=true)
    end
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# pseudo-definitions to show how everything behaves
#
# throw(label, val) = # throw a value to a dynamically enclosing block
#
# function rethrow(val)
#     global current_exception = val
#     throw(current_handler(), current_exception)
# end
#
# rethrow() = rethrow(current_exception)
#
# function throw(val)
#     global catch_backtrace = backtrace()
#     rethrow(val)
# end

## native julia error handling ##

error(s::AbstractString) = throw(ErrorException(s))
error(s...) = throw(ErrorException(Main.Base.string(s...)))

rethrow() = ccall(:jl_rethrow, Bottom, ())
rethrow(e) = ccall(:jl_rethrow_other, Bottom, (Any,), e)
backtrace() = ccall(:jl_backtrace_from_here, Array{Ptr{Void},1}, (Int32,), false)
catch_backtrace() = ccall(:jl_get_backtrace, Array{Ptr{Void},1}, ())

## keyword arg lowering generates calls to this ##
kwerr(kw, args...) = throw(MethodError(typeof(args[1]).name.mt.kwsorter, (kw,args...)))

## system error handling ##

systemerror(p, b::Bool; extrainfo=nothing) = b ? throw(Main.Base.SystemError(string(p), Libc.errno(), extrainfo)) : nothing

## assertion functions and macros ##

assert(x) = x ? nothing : throw(Main.Base.AssertionError())
macro assert(ex, msgs...)
    msg = isempty(msgs) ? ex : msgs[1]
    if !isempty(msgs) && (isa(msg, Expr) || isa(msg, Symbol))
        # message is an expression needing evaluating
        msg = :(Main.Base.string($(esc(msg))))
    elseif isdefined(Main, :Base) && isdefined(Main.Base, :string)
        msg = Main.Base.string(msg)
    else
        # string() might not be defined during bootstrap
        msg = :(Main.Base.string($(Expr(:quote,msg))))
    end
    :($(esc(ex)) ? $(nothing) : throw(Main.Base.AssertionError($msg)))
end

# NOTE: Please keep the constant values specified below in sync with the doc string
const DEFAULT_RETRY_N = 1
const DEFAULT_RETRY_ON = e->true
const DEFAULT_RETRY_MAX_DELAY = 10.0

"""
    retry(f, [retry_on]; n=1, max_delay=10.0) -> Function

Returns a lambda that retries function `f` up to `n` times in the
event of an exception. If `retry_on` is a `Type` then retry only
for exceptions of that type. If `retry_on` is a function
`test_error(::Exception) -> Bool` then retry only if it is true.

The first retry happens after a gap of 50 milliseconds or `max_delay`,
whichever is lower. Subsequently, the delays between retries are
exponentially increased with a random factor up to `max_delay`.

**Examples**
```julia
retry(http_get, e -> e.status == "503")(url)
retry(read, UVError)(io)
```
"""
function retry(f::Function, retry_on::Function=DEFAULT_RETRY_ON; n=DEFAULT_RETRY_N, max_delay=DEFAULT_RETRY_MAX_DELAY)
    (args...) -> begin
        delay = min(0.05, max_delay)
        for i = 1:n+1
            try
                return f(args...)
            catch e
                if i > n || try retry_on(e) end !== true
                    rethrow(e)
                end
            end
            delay = min(max_delay, delay)
            sleep(delay * (0.8 + (rand() * 0.2)))
            delay = delay * 5
        end
    end
end

retry(f::Function, t::Type; kw...) = retry(f, e->isa(e, t); kw...)

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
module Jutils

include("Transforms.jl")
include("Elements.jl")
include("Functions.jl")
include("Topologies.jl")
include("Mesh.jl")
include("Integration.jl")

end # module

max_stars_repo_name: TheBB/Jutils

max_issues_repo_name: TheBB/Jutils

max_forks_repo_name: TheBB/Jutils


text:
module PairwiseAlignment

export
    pairalign,
    # types
    AlignmentResult,
    # scoring systems
    AbstractSubstitutionMatrix,
    SubstitutionMatrix,
    AffineGapScoreModel,
    CostModel,
    AbstractSubstitutionMatrix,
    # alignment types
    GlobalAlignment,
    LocalAlignment,
    SemiGlobalAlignment,
    # distances
    EditDistance,
    HammingDistance,
    LevenshteinDistance,
    # gapped sequence
    GappedSequence,
    n_chars,
    n_gaps,
    leading_chars,
    leading_gaps,
    trailing_chars,
    trailing_gaps,
    counts,
    reversed_counts,
    push_chars!,
    push_gaps!,
    gapmap

# alignment types
# ---------------

# alignment
immutable GlobalAlignment end
immutable LocalAlignment end
immutable SemiGlobalAlignment end

# distances
immutable EditDistance end
immutable HammingDistance end
immutable LevenshteinDistance end

include("gapped_sequence.jl")
include("model.jl")
include("result.jl")
include("pairalign.jl")
include("algorithm/common.jl")
include("algorithm/affinegap_global_align.jl")
include("algorithm/affinegap_local_align.jl")
include("algorithm/affinegap_semiglobal_align.jl")
include("algorithm/affinegap_banded_global_align.jl")
include("algorithm/edit_distance.jl")
include("algorithm/hamming_distance.jl")

end # module


text:
########################################################################
function ppca(X; Q=1)
########################################################################

    D, N = size(X)

    μ = vec(mean(X, dims=2))

    JITTER = 1e-8

    @printf("Running PPCA for %d number of data items of dimension %d, projected to %d components\n", N, D, Q)

    @printf("\t Fixing jitter to %f\n", JITTER)

    #---------------------------------------------
    function marginalLogLikelihood(W, σ)
    #---------------------------------------------

        sum(logpdf(MvNormal(μ, W*W' + σ*σ*I + JITTER*I), X))

    end

    #---------------------------------------------
    function unpack(param)
    #---------------------------------------------

        local W = reshape(param[1:end-1], D, Q)

        local σ = exp(param[end]) + 1e-8

        return W, σ

    end

    #---------------------------------------------
    function objective(param)
    #---------------------------------------------

        @assert(length(param) == D*Q + 1)

        local W, σ = unpack(param)

        try

            return -1.0 * marginalLogLikelihood(W, σ)

        catch err

            if isa(err, PosDefException)

                @warn("covariance is not positive definite, returning Inf")

                return Inf

            else

                throw(err)

            end

        end

    end

    #---------------------------------------------
    # Run optimiser
    #---------------------------------------------

    opt    = Optim.Options(show_trace = true, iterations = 100_000)

    result = optimize(objective, [randn(D*Q); randn()*3], NelderMead())# LBFGS(), opt, autodiff=:forward)

    W, σ   = unpack(result.minimizer)

    # @show W, μ, σ
    #
    # #---------------------------------------------
    # # Define projections
    # #---------------------------------------------
    #
    # M = W'*W  + σ*σ*I # eq. (12.41) in Bishop
    #
    # proj(x::Array{T,1} where T<:Real) = vec(M \ (W'*(x - μ)))[1]
    #
    # fwd(z::Real) = vec(W*z + μ)

    return W, μ, σ#, proj, fwd

end


text:
using DataStructures

const data = parse.(Int, reduce(hcat, collect.(readlines())))

neighbours(pos; expand = 1) =
    filter([
        CartesianIndex(pos.I .- n)
        for n ∈ ((-1, 0), (1, 0), (0, -1), (0, 1))
    ]) do I
        1 <= I.I[1] <= size(data, 1) * expand && 1 <= I.I[2] <= size(data, 2) * expand
    end

div1mod1(x, y) = (div(x, y) - (mod(x, y) == 0 ? 1 : 0), mod1(x, y))

function get_val(data, pos::CartesianIndex)
    qx, rx = div1mod1(pos[1], size(data, 1))
    qy, ry = div1mod1(pos[2], size(data, 2))

    reduced_pos = CartesianIndex(rx, ry)
    mod1(data[reduced_pos] + qx + qy, 9)
end

function dijkstra(start, finish; expand = 1)
    q = PriorityQueue{CartesianIndex,Int}(Base.Order.Forward)
    seen = Set{CartesianIndex}()

    enqueue!(q, start => 0)

    while !isempty(q)
        pos, d = dequeue_pair!(q)
        if pos ∈ seen
            continue
        else
            push!(seen, pos)
        end
        if pos == finish
            return d
        end
        for nei in neighbours(pos; expand)
            if nei ∉ seen
                if nei ∉ keys(q)
                    enqueue!(q, nei => d + get_val(data, nei))
                else
                    q[nei] = min(q[nei], d + get_val(data, nei))
                end
            end
        end
    end
    return nothing
end

# Part One
dijkstra(CartesianIndex(1, 1), CartesianIndex(size(data)...)) |> println

# Part Two
dijkstra(CartesianIndex(1, 1), CartesianIndex(size(data)...) * 5; expand = 5) |> println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
using ITensors

N = 10
s = siteinds("Qubit", N)

# Make all of the operators
X = ops(s, [("X", n) for n in 1:N])
H = ops(s, [("H", n) for n in 1:N])
CX = ops(s, [("CX", n, m) for n in 1:N, m in 1:N])

# Start with the state |0000...⟩
ψ0 = productMPS(s, "0")

# Change to the state |1010...⟩
gates = [X[n] for n in 1:2:N]
ψ = apply(gates, ψ0; cutoff=1e-15)
@assert inner(ψ, productMPS(s, n -> isodd(n) ? "1" : "0")) ≈ 1

# Change to the state |10111011...⟩
append!(gates, [CX[n, n + 3] for n in 1:4:(N - 3)])
ψ = apply(gates, ψ0; cutoff=1e-15)
@assert inner(ψ, productMPS(s, ["1", "0", "1", "1", "1", "0", "1", "1", "1", "0"])) ≈ 1

# Change the state |10111011...⟩ to the (|+⟩, |-⟩) basis
append!(gates, [H[n] for n in 1:N])
ψ = apply(gates, ψ0; cutoff=1e-15)


text:
"""
    Term_Check(x::BnBSolver,y::BnBModel,k_int::Int64)

Checks if algorithm continues for the branch and bound problem. Checks that:
* Nodes are remaining on the stack.
* The iteration limit (`x.max_iter`) is not exceeded if iteration limit
flag (`x.iter_lim`) is true.
* Maximum nodes (`x.max_nodes`) on stack are not exceeded.
* Absolute global tolerance (`x.BnB_atol`) is not reached.
* Relative global tolerance (`x.BnB_rtol`) is not reached.
* Target upper bound (`x.target_upper`) is not reached.
Outputs a description of termination check if it returns false.
"""
function Term_Check(x::BnBSolver,y::BnBModel,k_int::Int64)
  t1 = length(y.LBD)>0
  t2 = (x.iter_lim ? k_int<x.max_iter : true)
  t3 = (length(y.LBD)<x.max_nodes)
  t4 = (y.UBDg-y.LBDg) > x.BnB_atol
  t5 = (y.UBDg-y.LBDg) > abs(y.LBDg)*x.BnB_rtol || ~(y.LBDg > -Inf)
  t6 = x.target_upper <= y.UBDg
  if t1 & t2 & t3 & t4 & t5 & t6
    return true
  else
    if (x.Verbosity=="Normal"||x.Verbosity=="Full")
      if ~(length(y.LBD)>0)
        if (y.first_num>0)
          println("Empty Stack: Exhaustive Search Finished")
        else
          println("Empty Stack: Infeasible")
        end
      elseif ~(length(y.LBD)<x.max_nodes)
        println("Node Limit Exceeded")
      elseif ~(x.iter_lim ? k_int<x.max_iter : true)
        println("Maximum Iteration Exceeded")
      else
        println("Convergence Tolerance Reached")
      end
    end
    return false
  end
end

"""
    Conv_Check(x::BnBSolver,ubd::Float64,lbd::Float64)

Checks if convergence tolerance is reach for the branch and bound problem.
Checks that:
* Absolute global tolerance (`x.BnB_atol`) is reached.
* Relative global tolerance (`x.BnB_rtol`) is reached.
Outputs a description of termination check if it returns false.
"""
function Conv_Check(x::BnBSolver,ubd::Float64,lbd::Float64)
  return ((abs(ubd-lbd) <= x.BnB_atol) || (abs(ubd-lbd)/(min(abs(lbd),abs(ubd))) <= x.BnB_rtol))
end

"""
    Repeat_Node_Default(x::BnBSolver,y::BnBModel,Xin::Vector{Interval{Float64}},
                                 Xout::Vector{Interval{Float64}})

Default check for repeating a node. Always returns false.
"""
function Repeat_Node_Default(x::BnBSolver,y::BnBModel{Interval{T}},
                             Xin::Vector{Interval{T}},
                             Xout::Vector{Interval{T}}) where {T<:AbstractFloat}
  return false
end

function Repeat_Node_Default(x::BnBSolver,y::BnBModel{MCInterval{T}},
                             Xin::Vector{MCInterval{T}},
                             Xout::Vector{MCInterval{T}}) where {T<:AbstractFloat}
  return false
end


text:
using Test
using Algorithm101
const ALGO = Algorithm101

@testset "Monte Carlo Integral" begin
    f = x -> x
    @test isapprox(ALGO.mc_integral(f, 0:0.00001:1), 0.5, atol=1e-2)

    f = x -> x^2
    @test isapprox(ALGO.mc_integral(f, 0:0.00001:1), 0.333, atol=1e-2)

    f = x -> x * (atan(x))^2
    @test isapprox(ALGO.mc_integral(f, 0:0.00001:1), 0.173, atol=1e-2)

    f = x ->  1 / (1 + (tan(x))^4)
    @test isapprox(ALGO.mc_integral(f, 0:0.00001:1), 0.769, atol=1e-2)
end


text:
@testset "Trends" begin
  Random.seed!(123)

  # constant trend
  d = georef((z=rand(100),), CartesianGrid(100))
  z̄ = trend(d, :z)[:z]
  @test all(abs.(diff(z̄)) .< 0.01)

  # linear trend
  μ = range(0,stop=1,length=100)
  ϵ = 0.1rand(100)
  d = georef((z=μ+ϵ,), CartesianGrid(100))
  z̄ = trend(d, :z)[:z]
  @test all([abs(z̄[i] - μ[i]) < 0.1 for i in 1:length(z̄)])

  # quadratic trend
  r = range(-1,stop=1,length=100)
  μ = [x^2 + y^2 for x in r, y in r]
  ϵ = 0.1rand(100,100)
  d = georef((z=μ+ϵ,))
  d̄ = trend(d, :z, degree=2)
  z̄ = reshape(d̄[:z], 100, 100)
  @test all([abs(z̄[i] - μ[i]) < 0.1 for i in 1:length(z̄)])

  if visualtests
    p₁ = heatmap(μ+ϵ, title="z")
    p₂ = heatmap(z̄, title="z trend")
    plt = plot(p₁, p₂, size=(900,300))
    @test_reference "data/trends.png" plt
  end
end


text:
# Autogenerated wrapper script for XPA_jll for x86_64-w64-mingw32
export libxpa, xpaaccess, xpaget, xpainfo, xpamb, xpans, xpaset

JLLWrappers.@generate_wrapper_header("XPA")
JLLWrappers.@declare_library_product(libxpa, "libxpa.dll")
JLLWrappers.@declare_executable_product(xpaaccess)
JLLWrappers.@declare_executable_product(xpaget)
JLLWrappers.@declare_executable_product(xpainfo)
JLLWrappers.@declare_executable_product(xpamb)
JLLWrappers.@declare_executable_product(xpans)
JLLWrappers.@declare_executable_product(xpaset)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libxpa,
        "bin\\libxpa.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        xpaaccess,
        "bin\\xpaaccess.exe",
    )

    JLLWrappers.@init_executable_product(
        xpaget,
        "bin\\xpaget.exe",
    )

    JLLWrappers.@init_executable_product(
        xpainfo,
        "bin\\xpainfo.exe",
    )

    JLLWrappers.@init_executable_product(
        xpamb,
        "bin\\xpamb.exe",
    )

    JLLWrappers.@init_executable_product(
        xpans,
        "bin\\xpans.exe",
    )

    JLLWrappers.@init_executable_product(
        xpaset,
        "bin\\xpaset.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
s = get_simulation("hopper_3D_quaternion", "flat_3D_nc", "flat_nc")

# time
h = 0.1
T = 50

# Rn + quaternion space
rq_space = rn_quaternion_space(num_var(s.model, s.env) - 1, x -> Gz_func(s.model, s.env, x),
			collect([(1:3)..., (8:num_var(s.model, s.env))...]),
			collect([(1:3)..., (7:num_var(s.model, s.env)-1)...]),
			collect((4:7)),
			collect((4:6)))

# initial conditions
p0 = [0.0; 0.0; 0.5]
v0 = [0.0; 0.0; 0.0]
quat0 = [1.0; 0.0; 0.0; 0.0]
ω0 = [0.0; 1.0; 0.0]

q0 = SVector{s.model.dim.q}([p0; quat0; 0.5])
q1 = SVector{s.model.dim.q}([p0 + v0 * h; 0.5 * h * L_multiply(quat0) * [sqrt((2.0 / h)^2.0 - ω0' * ω0); ω0]; 0.5])

@assert norm(q0[4:7]) ≈ 1.0
@assert norm(q1[4:7]) ≈ 1.0

# simulator
sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
	space = rq_space,
	ip_opts = ContactImplicitMPC.InteriorPointOptions(
		r_tol = 1.0e-6, κ_tol = 1.0e-6,
		diff_sol = true,
		solver = :lu_solver),
	sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = false))

# simulate
@time status = ContactImplicitMPC.simulate!(sim)
@test status

include(joinpath(pwd(), "src/dynamics/hopper_3D_quaternion/visuals.jl"))
vis = Visualizer()
render(vis)
# open(vis)
visualize_robot!(vis, s.model, sim.traj.q, h = h)

@assert all([norm(q[4:7]) ≈ 1.0 for q in sim.traj.q])


text:
"""
    initMat!(matrix::SparseMatrixCSC{Tv, Ti}, diag_l::Ti, diag_u::Ti, size::Ti; scale::Real = 1.0, shift::Real = 0.0, sparsity::Real = 0.9) where {Tv<:Complex, Ti<:Integer}

Initialization of matrix to be generated. This function is for non-Hermitian case, in which the entries of matrix are complex scalars.

In this function, the parameters `diag_l` and `diag_u` determine the range (`[diag_l, diag_u]`) of lower triangular part to be filled. The two parameters should also be negatif, which
refer to the offfsets of diagonals in the lower triangular part. The parameter `size` defines the size of matrix to be generated.

The non-zero elements of the initialized matrices are randomly generated in the interval (0,1). Two optionals parameters `shift` and `scale` allow shifting and scaling this interval.

Another optional parameter `sparsity` is more important. It determines the possility of the elements within the band of diagonal determined by `diag_l` and `diag_u` to be non-zeros.
In other words, the initalized matrix would be more sparse with a smaller number of `sparisty` parameter.

## Examples
```julia
julia> Am=spzeros(ComplexF64, 50, 50);
julia> initMat!(Am, -20, -10, 50)
julia> Am
50×50 SparseMatrixCSC{ComplexF64, Int64} with 335 stored entries:

⣄
⣟⡳⣄
⣿⣿⣢⣵⣄
⠙⢟⣾⣿⣽⢲⣄
⠀⠀⠙⢽⣿⣛⣿⣵⣄
⠀⠀⠀⠀⠙⠽⣻⡾⣿⣷⣄
⠀⠀⠀⠀⠀⠀⠙⢟⣻⢿⢻⣗⣄
⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣽⣯⣻⣳⣄
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣛⢿⡾⣷⣄
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢯⡽⣯⣜⣷⢄
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠚⠋⠛⠚⠓

```

## Examples
```julia
julia> Am=spzeros(ComplexF64, 50, 50);
julia> initMat!(Am, -20, -10, 50, sparsity=0.1)
julia> Am
50×50 SparseMatrixCSC{ComplexF64, Int64} with 36 stored entries:

⠀⠠⢀
⠉⠂
⠀⠀⠀⠂
⠀⠀⠀⠀⢂⠀⠈
⠀⠀⠀⠀⠀⠀⠂⠈⠀⠤
⠀⠀⠀⠀⠀⠀⠀⢀⢀⠅⠂⠄
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠠⠈⠰
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠈⠀⠡⢀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠄⢄⡆⢀

```

"""
function initMat!(
    matrix::SparseMatrixCSC{Tv, Ti},
    diag_l::Ti,
    diag_u::Ti,
    size::Ti;
    scale::Real = 1.0,
    shift::Real = 0.0,
    sparsity::Real = 0.9
) where {Tv<:Complex, Ti<:Integer}

    diag_l = abs(diag_l)
    diag_u = abs(diag_u)

    if diag_l < diag_u
        error(
            "for initialisationg of matrix, please ensure abs(diag_l) < abas(diag_u) ",
            )
        end

    if diag_u >= size
        error(
            "for initialisationg of matrix, please ensure abs(diag_u) < size ",
            )
    end

    cnt::Ti = (diag_l - diag_u + 1) * (2 * size - diag_l - diag_u) / 2

    rnd = rand(MersenneTwister(1234), Float64, cnt)

    for i = 1:cnt
        if real(rnd[i]) < 1.0 - sparsity
            rnd[i] = 0.0
        end
    end

    idx = 0
    for i = diag_u + 1:size
        for j = max(1, i - diag_l) : i - diag_u
            idx += 1
            matrix[i, j] = Tv(scale * (rnd[idx]) + shift)
        end
    end
end

"""
    initMat!(matrix::SparseMatrixCSC{Tv, Ti}, diag_l::Ti, diag_u::Ti, size::Ti; scale::Real = 1.0, shift::Real = 0.0, sparsity::Real = 0.9) where {Tv<:Real, Ti<:Integer}

Initialization of matrix to be generated. This function is for non-Symmetric case, in which the entries of matrix are real scalars.

The usage of this function is quite similar as the one for non-Hermitian case, the only difference is that the entries of matrices are real scalars.

## Examples
```julia
julia> Am=spzeros(Float64, 50, 50);
julia> initMat!(Am, -20, -10, 50, sparsity=0.5)
julia> Am
50×50 SparseMatrixCSC{Float64, Int64} with 178 stored entries:

⣄
⠌⡱⣄
⠍⠫⢠⡀⡀
⠈⠏⣌⠺⣑⠰⢀
⠀⠀⠀⠙⣣⢑⡿⢡⢀
⠀⠀⠀⠀⠀⠅⣊⡼⠖⠷⡄
⠀⠀⠀⠀⠀⠀⠉⢀⣫⠇⢛⡆
⠀⠀⠀⠀⠀⠀⠀⠀⠑⠟⡕⣭⠉⣱⡄
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⢐⡚⢘⠔⣡⢄
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠪⡌⡆⣜⣆⢀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠒⠉⠂⠘⠁

```

"""
function initMat!(
    matrix::SparseMatrixCSC{Tv, Ti},
    diag_l::Ti,
    diag_u::Ti,
    size::Ti;
    scale::Real = 1.0,
    shift::Real = 0.0,
    sparsity::Real = 0.9
) where {Tv<:Real, Ti<:Integer}

    diag_l = abs(diag_l)
    diag_u = abs(diag_u)

    if diag_u < 2
        error(
            "for initialisationg of non-symmetric matrix, please ensure abs(diag_u) >= 2 ",
            )
    end

    if diag_l < diag_u
        error(
            "for initialisationg of matrix, please ensure abs(diag_l) < abs(diag_u) ",
            )
    end

    if diag_u >= size
        error(
            "for initialisationg of matrix, please ensure abs(diag_u) < size ",
            )
    end

    cnt::Ti = (diag_l - diag_u + 1) * (2 * size - diag_l - diag_u) / 2

    rnd = rand(MersenneTwister(1234), Float64, cnt)

    for i = 1:cnt
        if real(rnd[i]) < 1.0 - sparsity
            rnd[i] = 0.0
        end
    end

    idx = 0
    for i = diag_u + 1:size
        for j = max(1, i - diag_l) : i - diag_u
            idx += 1
            matrix[i, j] = Tv(scale * (rnd[idx]) + shift)
        end
    end

end


text:
@polly function kernel_gemver(alpha, beta, A, u1, v1, u2, v2, w, x, y, z)
    n = size(u1,1)
    for i = 1:n, j = 1:n
        A[i,j] = A[i,j] + u1[i] * v1[j] + u2[i] * v2[j]
    end

    for i = 1:n, j = 1:n
        x[i] = x[i] + beta * A[j,i] * y[j]
    end

    for i = 1:n
        x[i] = x[i] + z[i]
    end

    for i = 1:n, j = 1:n
        w[i] = w[i] + alpha * A[i,j] * x[j]
    end
end

let
    n = 2000

    A = zeros(Float32, n, n)
    u1 = zeros(Float32, n)
    v1 = zeros(Float32, n)
    u2 = zeros(Float32, n)
    v2 = zeros(Float32, n)
    w = zeros(Float32, n)
    x = zeros(Float32, n)
    y = zeros(Float32, n)
    z = zeros(Float32, n)

    alpha = Float32(1.5)
    beta = Float32(1.2)

    fn = Float32(n)

    for i = 1:n
        u1[i] = i;
        u2[i] = ((i+1)/fn)/2.0
        v1[i] = ((i+1)/fn)/4.0
        v2[i] = ((i+1)/fn)/6.0
        y[i] = ((i+1)/fn)/8.0
        z[i] = ((i+1)/fn)/9.0
        x[i] = 0.0
        w[i] = 0.0
        for j = 1:n
            A[i,j] = (i*j % n) / n
        end
    end

    SUITE["gemver"] = @benchmarkable kernel_gemver($alpha, $beta, A, u1, v1, u2, v2, w, x, y, z) setup = (A = copy($A); u1 = copy($u1); v1 = copy($v1); u2 = copy($u2); v2 = copy($v2); w = copy($w); x = copy($x); y = copy($y); z = copy($z))
end


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 6a65913e-2745-11eb-2b15-9f16d977a8e3
begin
	import Pkg
	#Pkg.activate(mktempdir())
	Pkg.activate(".")
 	Pkg.add(["PlutoUI","Plots"])
	using PlutoUI
	using Plots
	plotly()
end

# ╔═╡ 0ba73c90-26a9-11eb-152d-9bbf1422ced7
html"<button onclick=present()>Present</button>"

# ╔═╡ 1a2960de-26a9-11eb-121b-2d0e2f2ec64c
md"""## யாதும் ஊரே யாவறும் கேளிர்!
### To us all towns are our own, everyone our kin!"""

# ╔═╡ 26a15d50-26a9-11eb-0da8-bb8870bbae00
  md"""### 1. கிடைத்தளத்தில் எறிதல்"""

# ╔═╡ 872b74c0-26aa-11eb-3c2e-8125dd5fdb85
md""" ###### 'u' திசைவேகத்தில் 't' நேரத்தில் எறிபொருள் கடந்த தொலைவு 'x'
      ###### எனில்.    x= ut or t=x/u"""

# ╔═╡ 10302340-26ad-11eb-1055-cf1e5bf798de
 md"""##### y= (1/2)gt^2.
      y=(1/2)g*(x/u)^2 """

# ╔═╡ 07545bf2-2753-11eb-17e3-5b7f80ce3107
begin
	md"""
	 கிடைத்தள திசைவேகம்:
	$(@bind u PlutoUI.Slider(1:9,default=2, show_value=true))
	 உயரம்:
	$(@bind H1 PlutoUI.Slider(100:100:1000,default=2, show_value=true))"""

end

# ╔═╡ 7ec82f3e-2744-11eb-07a1-df646ddeb880
begin
#u=10; # 10 m/s
g= 1.62;
md"""*ஆரம்ப உயரம்--> $H1*  &	 *ஈர்ப்பு விசை--> $g (m/s)*
	"""
end

# ╔═╡ 27d2d830-2748-11eb-3308-d3a8dca46f75
let
	x1=1:100;
	y=@. H1-(1/2)*g*(x1/u)^2
	indexArray = findall( x -> x <= 0, y );
	if(length(findall( x -> x <= 0, y )) ==0)
		indexArray=x1[end]
	end
	y[isless.(y, 0)] .=0
	#x1[isless.(y, 0)] .=0
	plot(x1,y,xlims=(0,x1[indexArray[1]]+10),c=:green,leg=:false)
	scatter!(x1,y)
	xaxis!("நெடுக்கம்")
	yaxis!("உயரம்")

end

# ╔═╡ 30d5c7c0-276c-11eb-0bcb-ad29cd883af7
md""" #### கணக்குகள்
"""

# ╔═╡ a7440c90-274a-11eb-3976-47e25458fac4
begin
md""" ##### வெவ்வேறு கோள்கள் ஈர்ப்பு விசை கொண்டு நெடுக்கத்தை கணக்கீடுக ????"""
end

# ╔═╡ 03205060-276d-11eb-1d73-2de0b39d7a70
HTML("""<img src="https://upload.wikimedia.org/wikipedia/commons/b/b7/Solar_System_size_to_scale_tamil.svg">""")

# ╔═╡ Cell order:
# ╟─0ba73c90-26a9-11eb-152d-9bbf1422ced7
# ╠═1a2960de-26a9-11eb-121b-2d0e2f2ec64c
# ╟─6a65913e-2745-11eb-2b15-9f16d977a8e3
# ╟─26a15d50-26a9-11eb-0da8-bb8870bbae00
# ╟─872b74c0-26aa-11eb-3c2e-8125dd5fdb85
# ╟─10302340-26ad-11eb-1055-cf1e5bf798de
# ╠═7ec82f3e-2744-11eb-07a1-df646ddeb880
# ╟─07545bf2-2753-11eb-17e3-5b7f80ce3107
# ╠═27d2d830-2748-11eb-3308-d3a8dca46f75
# ╟─30d5c7c0-276c-11eb-0bcb-ad29cd883af7
# ╟─a7440c90-274a-11eb-3976-47e25458fac4
# ╠═03205060-276d-11eb-1d73-2de0b39d7a70

max_stars_repo_name: ResearchersHelp/LearningMotherTongue

max_issues_repo_name: ResearchersHelp/LearningMotherTongue

max_forks_repo_name: ResearchersHelp/LearningMotherTongue


text:
const AVec = AbstractVector
const Rydberg2eV = 13.605698066

struct ElecHam
   basic_data
   ws_rvecs
   rvec_idx
   hopping
end

#struct LatticeIFC
#   basic_data
#   ws_rvecs
#   rvec_idx
#   forceconst
#end


text:
using Stheno: get_f, get_y, Observation, merge, GPC
using BlockArrays, LinearAlgebra

function abs_rel_errs(x, y)
    δ = abs.(vec(x) .- vec(y))
    return [δ δ ./ vec(y)]
end

@timedtestset "conditioning" begin
    @timedtestset "Observation" begin
        rng, N, N′, D = MersenneTwister(123456), 5, 6, 2
        X, X′ = ColVecs(randn(rng, D, N)), ColVecs(randn(rng, D, N′))
        y, y′ = randn(rng, N), randn(rng, N′)
        f = GP(1, EQ(), GPC())

        fX, fX′ = f(X), f(X′)
        c1, c2 = fX←y, fX′←y′
        @test Observation(fX, y) == c1
        @test get_f(c1) === fX && get_f(c2) === fX′
        @test get_y(c1) === y && get_y(c2) === y′

        c = merge((c1, c2))
        @test get_y(c) == BlockVector([y, y′])
        @test get_y(Obs((c1, c2))) == BlockVector([y, y′])
    end
    @timedtestset "condition once" begin
        rng, N, N′, D = MersenneTwister(123456), 10, 3, 2
        x = collect(range(-3.0, stop=3.0, length=N))
        f = GP(1, EQ(), GPC())
        y = rand(rng, f(x))

        # Test mechanics for finite conditioned process with single conditioning.
        f′ = f | (f(x, 1e-9)←y)
        @test maximum(abs.(rand(rng, f′(x, 1e-9)) - y)) < 1e-3
        @test maximum(abs.(mean(f′(x)) - y)) < 1e-3
        @test all(abs.(cov(f′(x))) .< 1e-6)
    end
    @timedtestset "condition repeatedly" begin
        rng, N, N′ = MersenneTwister(123456), 5, 7
        xx′ = collect(range(-3.0, stop=3.0, length=N+N′))
        idx = randperm(rng, length(xx′))[1:N]
        idx_1, idx_2 = idx, setdiff(1:length(xx′), idx)
        x, x′ = xx′[idx_1], xx′[idx_2]

        f = GP(1, EQ(), GPC())
        y = rand(rng, f(xx′))
        y1, y2 = y[idx_1], y[idx_2]

        # Construct posterior using one conditioning operation.
        f′ = f | (f(xx′, 0.1) ← y)

        # Construct posterior using two conditioning operations.
        f′1 = f | (f(x, 0.1) ← y1)
        f′2 = f′1 | (f′1(x′, 0.1) ← y2)

        # Check that conditioning twice yields the same answer.
        @test mean(f′(xx′)) ≈ mean(f′2(xx′))
        @test cov(f′(xx′)) ≈ cov(f′2(xx′))
        @test cov(f′(x), f′(x′)) ≈ cov(f′2(x), f′2(x′))
    end
    @timedtestset "Consistency Tests" begin
        rng, N, P, Q = MersenneTwister(123456), 11, 5, 3
        f = GP(sin, EQ(), GPC())
        x = collect(range(-1.0, 1.0; length=N))
        y = randn(rng, N)
        f′ = f | (f(x) ← y)
        x0, x1, x2, x3 = randn(rng, P), randn(rng, Q), randn(rng, Q), randn(rng, P)
        abstractgp_interface_tests(f′, f, x0, x1, x2, x3)

        # Multiple things to condition on
        P′ = 7
        x′ = randn(rng, P′)
        y, y′ = rand(rng, [f(x, 1e-9), f(x′, 1e-9)])
        f′ = f | (f(x, 1e-9)←y, f(x′, 1e-9)←y′)
        abstractgp_interface_tests(f′, f, x0, x1, x2, x3)

        # Condition multiple things
        f′, f′′ = (f, f) | (f(x, 1e-9)←y)
        abstractgp_interface_tests(f′, f′′, x0, x1, x2, x3)

        # Condition multiple things on multiple things
        f′, f′′ = (f, f) | (f(x, 1e-9)←y, f(x′, 1e-9)←y′)
        abstractgp_interface_tests(f′, f′′, x0, x1, x2, x3)
    end
    @timedtestset "Diff Tests" begin
        rng, N, P, Q = MersenneTwister(123456), 11, 5, 3
        x_obs, A_obs = collect(range(-5.0, 5.0; length=N)), randn(rng, N, N)
        y_obs = rand(rng, GP(sin, EQ(), GPC())(x_obs, _to_psd(A_obs)))
        standard_1D_tests(
            MersenneTwister(123456),
            Dict(:x=>x_obs, :y=>y_obs, :A=>A_obs),
            θ->begin
                f = GP(sin, EQ(), GPC())
                f′ = f | (f(θ[:x], _to_psd(θ[:A]))←θ[:y])
                return f′, f′
            end,
            collect(range(-4.5, 4.5; length=P)),
            collect(range(-4.0, 5.0; length=Q)),
        )
    end
    @timedtestset "Conditioning with cross" begin
        rng, N, N′, σ² = MersenneTwister(123456), 3, 7, 1.0
        xx′ = collect(range(-3.0, stop=3.0, length=N+N′))
        xp = collect(range(-4.0, stop=4.0, length=N+N′+10))
        xp′ = collect(range(-4.0, stop=4.0, length=N+N′+11))
        f = GP(sin, EQ(), GPC())
        yy′ = rand(rng, f(xx′, σ²))

        # Chop up into blocks.
        idx = randperm(rng, length(xx′))[1:N]
        idx_1, idx_2 = idx, setdiff(1:length(xx′), idx)
        x, x′ = xx′[idx_1], xx′[idx_2]
        y, y′ = yy′[idx_1], yy′[idx_2]

        f′ = f | (f(xx′, σ²)←yy′)
        fb′ = f | (cross([f, f])(BlockData([x, x′]), σ²)←vcat(y, y′))
        fmc′ = f | (f(x, σ²)←y, f(x′, σ²)←y′)
        f′_obs = f | Obs(f(x, σ²)←y, f(x′, σ²)←y′)

        @test mean(f′(xp)) ≈ mean(fb′(xp))
        @test mean(f′(xp)) ≈ mean(fmc′(xp))
        @test mean(f′(xp)) ≈ mean(f′_obs(xp))

        @test cov(f′(xp)) ≈ cov(fb′(xp))
        @test cov(f′(xp)) ≈ cov(fmc′(xp))
        @test cov(f′(xp)) ≈ cov(f′_obs(xp))

        @test cov(f′(xp), f′(xp)) ≈ cov(fb′(xp), fb′(xp))
        @test cov(f′(xp), f′(xp)) ≈ cov(fmc′(xp), fmc′(xp))
        @test cov(f′(xp), f′(xp)) ≈ cov(f′_obs(xp), f′_obs(xp))

        @test cov(f′(xp), f′(xp′)) ≈ cov(fb′(xp), fb′(xp′))
        @test cov(f′(xp), f′(xp′)) ≈ cov(fmc′(xp), fmc′(xp′))
        @test cov(f′(xp), f′(xp′)) ≈ cov(f′_obs(xp), f′_obs(xp′))
    end
end


text:
mutable struct SigmoidalFeature
    _means::AbstractArray{Float64,1}
    _coeff::Float64
    _n_kernels::Int64 # cache the number of kernels
    function SigmoidalFeature(means::AbstractArray{Float64,1}, coeffs::Float64)
        new(means, coeffs, size(means)[1])
    end
end

function _sigmoid(mean::Float64, coeff::Float64, x::AbstractArray{Float64,1})
    return tanh.((x .- mean) .* (coeff * 0.5)) .* 0.5 .+ 0.5
end

function transform(feature::SigmoidalFeature, x::AbstractArray{Float64,1})
    # x is an array of size n_samples
    # returns Phi = [phi_1(x),,,, phi_N(x)] where phi_i(x) is a sigmoidal centered around means[i]
    n_samples = size(x)[1]
    n_kernels = feature._n_kernels
    Phi = zeros(n_kernels, n_samples)
    for i = 1:n_kernels
        mean = feature._means[i]
        Phi[i, :] = _sigmoid(mean, feature._coeff, x)
    end

    return Phi
end

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
using FinEtools
File = "snap.py"
filetag(v) = Float64(Int(round(v*100000))) / 100000
filename(i) = begin
	file = "snap"
	if i < 10
		file = file * "000$(i)"
	elseif i < 100
		file = file * "00$(i)"
	end
end
lambdas = collect(linearspace(0.0, 1.0, 49))
for i in 1:10
	push!(lambdas, lambdas[end])
end
for (i, lambda) in enumerate(lambdas[2:end])
	VTKfile = "bertoldi_compression-$(filetag(lambda)).vtk"
	@show s = stat(VTKfile)
if s.mode == 0o000000
	break
end
open(File, "w") do f
	print(f, """
#### import the simple module from the paraview
from paraview.simple import *
#### disable automatic camera reset on 'Show'
paraview.simple._DisableFirstRenderCameraReset()

# create a new 'Legacy VTK Reader'
bertoldi_compressionboundaryvtk = LegacyVTKReader(FileNames=['bertoldi_compression-boundary.vtk'])

# get active view
renderView1 = GetActiveViewOrCreate('RenderView')
# uncomment following to set a specific view size
renderView1.ViewSize = [966, 966]

# show data in view
bertoldi_compressionboundaryvtkDisplay = Show(bertoldi_compressionboundaryvtk, renderView1)

# trace defaults for the display properties.
bertoldi_compressionboundaryvtkDisplay.Representation = 'Surface'
bertoldi_compressionboundaryvtkDisplay.ColorArrayName = [None, '']
bertoldi_compressionboundaryvtkDisplay.OSPRayScaleFunction = 'PiecewiseFunction'
bertoldi_compressionboundaryvtkDisplay.SelectOrientationVectors = 'None'
bertoldi_compressionboundaryvtkDisplay.ScaleFactor = 0.0101258232648
bertoldi_compressionboundaryvtkDisplay.SelectScaleArray = 'None'
bertoldi_compressionboundaryvtkDisplay.GlyphType = 'Arrow'
bertoldi_compressionboundaryvtkDisplay.GlyphTableIndexArray = 'None'
bertoldi_compressionboundaryvtkDisplay.GaussianRadius = 0.00050629116324
bertoldi_compressionboundaryvtkDisplay.SetScaleArray = [None, '']
bertoldi_compressionboundaryvtkDisplay.ScaleTransferFunction = 'PiecewiseFunction'
bertoldi_compressionboundaryvtkDisplay.OpacityArray = [None, '']
bertoldi_compressionboundaryvtkDisplay.OpacityTransferFunction = 'PiecewiseFunction'
bertoldi_compressionboundaryvtkDisplay.DataAxesGrid = 'GridAxesRepresentation'
bertoldi_compressionboundaryvtkDisplay.SelectionCellLabelFontFile = ''
bertoldi_compressionboundaryvtkDisplay.SelectionPointLabelFontFile = ''
bertoldi_compressionboundaryvtkDisplay.PolarAxes = 'PolarAxesRepresentation'
bertoldi_compressionboundaryvtkDisplay.ScalarOpacityUnitDistance = 0.011328756643810481

# init the 'GridAxesRepresentation' selected for 'DataAxesGrid'
bertoldi_compressionboundaryvtkDisplay.DataAxesGrid.XTitleFontFile = ''
bertoldi_compressionboundaryvtkDisplay.DataAxesGrid.YTitleFontFile = ''
bertoldi_compressionboundaryvtkDisplay.DataAxesGrid.ZTitleFontFile = ''
bertoldi_compressionboundaryvtkDisplay.DataAxesGrid.XLabelFontFile = ''
bertoldi_compressionboundaryvtkDisplay.DataAxesGrid.YLabelFontFile = ''
bertoldi_compressionboundaryvtkDisplay.DataAxesGrid.ZLabelFontFile = ''

# init the 'PolarAxesRepresentation' selected for 'PolarAxes'
bertoldi_compressionboundaryvtkDisplay.PolarAxes.PolarAxisTitleFontFile = ''
bertoldi_compressionboundaryvtkDisplay.PolarAxes.PolarAxisLabelFontFile = ''
bertoldi_compressionboundaryvtkDisplay.PolarAxes.LastRadialAxisTextFontFile = ''
bertoldi_compressionboundaryvtkDisplay.PolarAxes.SecondaryRadialAxesTextFontFile = ''

# reset view to fit data
renderView1.ResetCamera()

#changing interaction mode based on data extents
renderView1.CameraPosition = [0.0005937481759999974, 0.0504437499, 10000.0]
renderView1.CameraFocalPoint = [0.0005937481759999974, 0.0504437499, 0.0]

# get the material library
materialLibrary1 = GetMaterialLibrary()

# update the view to ensure updated data information
renderView1.Update()

# reset view to fit data
renderView1.ResetCamera()

# create a new 'Legacy VTK Reader'
snapshotvtk = LegacyVTKReader(FileNames=['$(VTKfile)'])

# show data in view
snapshotvtkDisplay = Show(snapshotvtk, renderView1)

# trace defaults for the display properties.
snapshotvtkDisplay.Representation = 'Surface'
snapshotvtkDisplay.ColorArrayName = [None, '']
snapshotvtkDisplay.OSPRayScaleArray = 'u'
snapshotvtkDisplay.OSPRayScaleFunction = 'PiecewiseFunction'
snapshotvtkDisplay.SelectOrientationVectors = 'u'
snapshotvtkDisplay.ScaleFactor = 0.0101258232648
snapshotvtkDisplay.SelectScaleArray = 'None'
snapshotvtkDisplay.GlyphType = 'Arrow'
snapshotvtkDisplay.GlyphTableIndexArray = 'None'
snapshotvtkDisplay.GaussianRadius = 0.00050629116324
snapshotvtkDisplay.SetScaleArray = ['POINTS', 'u']
snapshotvtkDisplay.ScaleTransferFunction = 'PiecewiseFunction'
snapshotvtkDisplay.OpacityArray = ['POINTS', 'u']
snapshotvtkDisplay.OpacityTransferFunction = 'PiecewiseFunction'
snapshotvtkDisplay.DataAxesGrid = 'GridAxesRepresentation'
snapshotvtkDisplay.SelectionCellLabelFontFile = ''
snapshotvtkDisplay.SelectionPointLabelFontFile = ''
snapshotvtkDisplay.PolarAxes = 'PolarAxesRepresentation'
snapshotvtkDisplay.ScalarOpacityUnitDistance = 0.009992410134288721

# init the 'GridAxesRepresentation' selected for 'DataAxesGrid'
snapshotvtkDisplay.DataAxesGrid.XTitleFontFile = ''
snapshotvtkDisplay.DataAxesGrid.YTitleFontFile = ''
snapshotvtkDisplay.DataAxesGrid.ZTitleFontFile = ''
snapshotvtkDisplay.DataAxesGrid.XLabelFontFile = ''
snapshotvtkDisplay.DataAxesGrid.YLabelFontFile = ''
snapshotvtkDisplay.DataAxesGrid.ZLabelFontFile = ''

# init the 'PolarAxesRepresentation' selected for 'PolarAxes'
snapshotvtkDisplay.PolarAxes.PolarAxisTitleFontFile = ''
snapshotvtkDisplay.PolarAxes.PolarAxisLabelFontFile = ''
snapshotvtkDisplay.PolarAxes.LastRadialAxisTextFontFile = ''
snapshotvtkDisplay.PolarAxes.SecondaryRadialAxesTextFontFile = ''

# update the view to ensure updated data information
renderView1.Update()

# create a new 'Warp By Vector'
warpByVector1 = WarpByVector(Input=snapshotvtk)
warpByVector1.Vectors = ['POINTS', 'u']

# show data in view
warpByVector1Display = Show(warpByVector1, renderView1)

# trace defaults for the display properties.
warpByVector1Display.Representation = 'Surface'
warpByVector1Display.ColorArrayName = [None, '']
warpByVector1Display.OSPRayScaleArray = 'u'
warpByVector1Display.OSPRayScaleFunction = 'PiecewiseFunction'
warpByVector1Display.SelectOrientationVectors = 'u'
warpByVector1Display.ScaleFactor = 0.010142636760398721
warpByVector1Display.SelectScaleArray = 'None'
warpByVector1Display.GlyphType = 'Arrow'
warpByVector1Display.GlyphTableIndexArray = 'None'
warpByVector1Display.GaussianRadius = 0.0005071318380199361
warpByVector1Display.SetScaleArray = ['POINTS', 'u']
warpByVector1Display.ScaleTransferFunction = 'PiecewiseFunction'
warpByVector1Display.OpacityArray = ['POINTS', 'u']
warpByVector1Display.OpacityTransferFunction = 'PiecewiseFunction'
warpByVector1Display.DataAxesGrid = 'GridAxesRepresentation'
warpByVector1Display.SelectionCellLabelFontFile = ''
warpByVector1Display.SelectionPointLabelFontFile = ''
warpByVector1Display.PolarAxes = 'PolarAxesRepresentation'
warpByVector1Display.ScalarOpacityUnitDistance = 0.009938686024079718

# init the 'GridAxesRepresentation' selected for 'DataAxesGrid'
warpByVector1Display.DataAxesGrid.XTitleFontFile = ''
warpByVector1Display.DataAxesGrid.YTitleFontFile = ''
warpByVector1Display.DataAxesGrid.ZTitleFontFile = ''
warpByVector1Display.DataAxesGrid.XLabelFontFile = ''
warpByVector1Display.DataAxesGrid.YLabelFontFile = ''
warpByVector1Display.DataAxesGrid.ZLabelFontFile = ''

# init the 'PolarAxesRepresentation' selected for 'PolarAxes'
warpByVector1Display.PolarAxes.PolarAxisTitleFontFile = ''
warpByVector1Display.PolarAxes.PolarAxisLabelFontFile = ''
warpByVector1Display.PolarAxes.LastRadialAxisTextFontFile = ''
warpByVector1Display.PolarAxes.SecondaryRadialAxesTextFontFile = ''

# hide data in view
Hide(snapshotvtk, renderView1)

# update the view to ensure updated data information
renderView1.Update()

# change representation type
warpByVector1Display.SetRepresentationType('Surface With Edges')

# current camera placement for renderView1
renderView1.InteractionMode = '2D'
renderView1.CameraPosition = [0.0012350857876053747, 0.050203248295647984, 0.27629011483940236]
renderView1.CameraFocalPoint = [0.0012350857876053747, 0.050203248295647984, 0.0]
renderView1.CameraParallelScale = 0.05252248397463676

RenderAllViews()
# save screenshot
SaveScreenshot('$(filename(i)).png', renderView1, ImageResolution=[966, 966])

#### saving camera placements for all active views

# current camera placement for renderView1
renderView1.InteractionMode = '2D'
renderView1.CameraPosition = [0.0012350857876053747, 0.050203248295647984, 0.27629011483940236]
renderView1.CameraFocalPoint = [0.0012350857876053747, 0.050203248295647984, 0.0]
renderView1.CameraParallelScale = 0.05252248397463676

#### uncomment the following to render all views

# alternatively, if you want to write images, you can use SaveScreenshot(...).
quit()
""")
end

run(`"pvpython" $File`)

end

run(`"C:\\Program Files\\ImageMagick-7.0.8-Q16\\magick.exe" -delay 20 snap0*.png anim.gif`)


text:
using Test
using LanguageServer.URIs2

@testset "URIs2" begin
    @testset "filepath2uri to string" begin
        # TODO Remove this Windows flag later, it is not in the original, but we need support for relative paths first
        if Sys.iswindows()
            @test filepath2uri("c:/win/path") |> string == "file:///c%3A/win/path"
            @test filepath2uri("C:/win/path") |> string == "file:///c%3A/win/path"
            @test filepath2uri("c:/win/path/") |> string == "file:///c%3A/win/path/"
        end
        @test filepath2uri("/c:/win/path") |> string == "file:///c%3A/win/path"
    end

    @testset "filepath2uri to string - Windows special" begin
        if Sys.iswindows()
            @test filepath2uri("c:\\win\\path") |> string == "file:///c%3A/win/path"
            @test filepath2uri("c:\\win/path") |> string == "file:///c%3A/win/path"
        # else TODO Put this else back in once we support these paths on Unix
            # @test filepath2uri("c:\\win\\path") |> string == "file:///c%3A%5Cwin%5Cpath"
            # @test filepath2uri("c:\\win/path") |> string == "file:///c%3A%5Cwin/path"
        end
    end

    @testset "uri2filepath to string - Windows special" begin
        if Sys.iswindows()
            @test uri2filepath(filepath2uri("c:\\win\\path")) == "c:\\win\\path"
            @test uri2filepath(filepath2uri("c:\\win/path")) == "c:\\win\\path"

            @test uri2filepath(filepath2uri("c:/win/path")) == "c:\\win\\path"
            @test uri2filepath(filepath2uri("c:/win/path/")) == "c:\\win\\path\\"
            @test uri2filepath(filepath2uri("C:/win/path")) == "c:\\win\\path"
            @test uri2filepath(filepath2uri("/c:/win/path")) == "c:\\win\\path"
            @test_broken uri2filepath(filepath2uri("./c/win/path")) == "\\.\\c\\win\\path"
         # else TODO Put this else back in once we support relative paths
            # @test uri2filepath(filepath2uri("c:/win/path")) == "c:/win/path"
            # @test uri2filepath(filepath2uri("c:/win/path/")) == "c:/win/path/"
            # @test uri2filepath(filepath2uri("C:/win/path")) == "c:/win/path"
            # @test uri2filepath(filepath2uri("/c:/win/path")) == "c:/win/path"
            # @test_broken uri2filepath(filepath2uri("./c/win/path")) == "/./c/win/path"
        end
    end

    @testset "uri2filepath - no `uri2filepath` when no `path`" begin
        value = URI("file://%2Fhome%2Fticino%2Fdesktop%2Fcpluscplus%2Ftest.cpp")

        @test value.authority == "/home/ticino/desktop/cpluscplus/test.cpp"
        @test_broken value.path == "/"
        if Sys.iswindows()
            @test_broken uri2filepath(value) == "\\"
        else
            @test_broken uri2filepath(value) == "/"
        end
    end

    @testset "parse" begin
        value = URI("http:/api/files/test.me?t=1234")
        @test value.scheme == "http"
        @test_broken value.authority == ""
        @test value.path == "/api/files/test.me"
        @test value.query == "t=1234"
        @test value.fragment === nothing

        value = URI("http://api/files/test.me?t=1234")
        @test value.scheme == "http"
        @test value.authority == "api"
        @test value.path == "/files/test.me"
        @test value.query == "t=1234"
        @test value.fragment === nothing

        value = URI("file:///c:/test/me")
        @test value.scheme == "file"
        @test value.authority == ""
        @test value.path == "/c:/test/me"
        @test value.fragment === nothing
        @test value.query=== nothing
        @test uri2filepath(value) == (Sys.iswindows() ? "c:\\test\\me" : "c:/test/me")

        value = URI("file://shares/files/c%23/p.cs")
        @test value.scheme == "file"
        @test value.authority == "shares"
        @test value.path == "/files/c#/p.cs"
        @test value.fragment === nothing
        @test value.query === nothing
        @test uri2filepath(value) == (Sys.iswindows() ? "\\\\shares\\files\\c#\\p.cs" : "//shares/files/c#/p.cs")

        value = URI("file:///c:/Source/Z%C3%BCrich%20or%20Zurich%20(%CB%88zj%CA%8A%C9%99r%C9%AAk,/Code/resources/app/plugins/c%23/plugin.json")
        @test value.scheme == "file"
        @test value.authority == ""
        @test value.path == "/c:/Source/Zürich or Zurich (ˈzjʊərɪk,/Code/resources/app/plugins/c#/plugin.json"
        @test value.fragment === nothing
        @test value.query === nothing

        value = URI("file:///c:/test %25/path")
        @test value.scheme == "file"
        @test value.authority == ""
        @test value.path == "/c:/test %/path"
        @test value.fragment === nothing
        @test value.query === nothing

        value = URI("inmemory:")
        @test value.scheme == "inmemory"
        @test value.authority === nothing
        @test value.path === ""
        @test value.query === nothing
        @test value.fragment === nothing

        value = URI("foo:api/files/test")
        @test value.scheme == "foo"
        @test value.authority === nothing
        @test value.path == "api/files/test"
        @test value.query === nothing
        @test value.fragment === nothing

        value = URI("file:?q")
        @test value.scheme == "file"
        @test_broken value.authority == ""
        @test_broken value.path == "/"
        @test value.query == "q"
        @test value.fragment === nothing

        value = URI("file:#d")
        @test value.scheme == "file"
        @test value.authority === nothing
        @test_broken value.path == "/"
        @test value.query === nothing
        @test value.fragment == "d"

        value = URI("f3ile:#d")
        @test value.scheme == "f3ile"
        @test value.authority === nothing
        @test value.path === ""
        @test value.query === nothing
        @test value.fragment == "d"

        value = URI("foo+bar:path")
        @test value.scheme == "foo+bar"
        @test value.authority === nothing
        @test value.path == "path"
        @test value.query === nothing
        @test value.fragment === nothing

        value = URI("foo-bar:path")
        @test value.scheme == "foo-bar"
        @test value.authority === nothing
        @test value.path == "path"
        @test value.query === nothing
        @test value.fragment === nothing

        value = URI("foo.bar:path")
        @test value.scheme == "foo.bar"
        @test value.authority === nothing
        @test value.path == "path"
        @test value.query === nothing
        @test value.fragment === nothing
    end
end


text:
# """
#     MOI.add_constrained_variable(model::Optimizer, set::MOI.Interval{Int})

# Add a variable to the model, enforcing to give an interval with it as you cannot create variables without bounds in SeaPearl.
# """

function MOI.add_variable(model::Optimizer)
    id = length(keys(model.moimodel.variables)) + 1

    vi = MOI.VariableIndex(id)

    moiv = MOIVariable("", nothing, nothing, vi)
    push!(model.moimodel.variables, moiv)

    return vi
end


text:
using Test
using CUDA
using CUTENSORNET
import CUTENSORNET: CuTensorNetwork, rehearse_contraction, perform_contraction!, AutoTune, NoAutoTune

@testset "CUTENSORNET" begin
    n = 8
    m = 16
    k = 32
    #@testset for elty in [Float16, Float32, Float64, ComplexF16, ComplexF32, ComplexF64, Int8, Int32, UInt8, UInt32]
    @testset for elty in [Float16, Float32]
        # make the simplest TN: two matrices
        A = rand(elty, n, k)
        B = rand(elty, k, m)
        A_modes = Int32[1, 2]
        B_modes = Int32[2, 3]
        A_extents = collect(size(A))
        B_extents = collect(size(B))
        A_strides = collect(strides(A))
        B_strides = collect(strides(B))
        A_aligns  = Int32(256)
        B_aligns  = Int32(256)

        C = zeros(elty, n, m)
        C_modes = Int32[1, 3]
        C_extents = collect(size(C))
        C_strides = collect(strides(C))
        C_aligns  = Int32(256)
        @testset for max_ws_size in [2^10, 2^20] # test that slicing works
            @testset for tuning in [AutoTune(), NoAutoTune()]
                ctn = CuTensorNetwork(elty, [A_modes, B_modes], [A_extents, B_extents], [A_strides, B_strides], [A_aligns, B_aligns], C_modes, C_extents, C_strides, C_aligns)
                info, plan = rehearse_contraction(ctn, max_ws_size)
                ctn.input_arrs = CuArray.([A, B])
                ctn.output_arr = CuArray(C)
                ctn = perform_contraction!(ctn, info, plan, max_ws_size, tuning)
                dC = ctn.output_arr
                @test collect(dC) ≈ A*B
            end
        end
    end
end


text:
using BenchmarkTools, StaticArrays
abstract type System end
@inline function eom_towel(x)
  x1, x2, x3 = x[1], x[2], x[3]
  SVector( 3.8*x1*(1-x1) - 0.05*(x2+0.35)*(1-2*x3),
  0.1*( (x2+0.35)*(1-2*x3) - 1 )*(1 - 1.9*x1),
  3.78*x3*(1-x3)+0.2*x2 )
end
@inline function jacob_towel(x)
  @SMatrix [3.8*(1 - 2x[1]) -0.05*(1-2x[3]) 0.1*(x[2] + 0.35);
  -0.19((x[2] + 0.35)*(1-2x[3]) - 1)  0.1*(1-2x[3])*(1-1.9x[1])  -0.2*(x[2] + 0.35)*(1-1.9x[1]);
  0.0  0.2  3.78(1-2x[3]) ]
end
u0 = @SVector [0.085, -0.121, 0.075]
struct DiscreteDS{D, T<:Real, F, J} <: System
  state::SVector{D,T}
  eom::F
  jacob::J
end

mutable struct DiscreteDSmut{D, T<:Real, F, J} <: System
  state::SVector{D,T}
  eom::F
  jacob::J
end

function evolve(state, ds::System, N::Int = 1)
  f = ds.eom
  for i in 1:N
    state = f(state)
  end
  return state
end

function evolve!(ds::DiscreteDS, N::Int = 1)
  st = ds.state
  st = evolve(st, ds, N)
  return DiscreteDS(st, ds.eom, ds.jacob)
end

function evolve!(ds::DiscreteDSmut, N::Int = 1)
  st = ds.state
  st = evolve(st, ds, N)
  ds.state = st
  return ds
end

ds = DiscreteDS(u0, eom_towel, jacob_towel)
dsmut = DiscreteDSmut(u0, eom_towel, jacob_towel)

ds = evolve!(ds, 100)
evolve!(dsmut, 100)

a = @btime evolve!($ds, 10000)
b = @btime evolve!($dsmut, 10000)


text:
"""
Uses `Plots.jl` to provide position and velocity
plots for R2BP, CR3BP, and NBP orbits, as well
as lagrange plots and zero velocity curves.
"""
module OrbitPlots

# Module Exports
export plotpositions, plotpositions!
export plotvelocities, plotvelocities!
export lagrangeplot, lagrangeplot!
export plotbodies, plotbodies!
export zerovelocityplot, zerovelocityplot!

# Module Dependencies
using Plots
using LinearAlgebra
using Unitful, UnitfulAstro, UnitfulAngles

using ..AstrodynamicsCore

# Source Files
include("R2BP/R2BPPlots.jl")
include("CR3BP/CR3BPPlots.jl")
include("CR3BP/LagrangePlot.jl")
include("CR3BP/ZeroVelocityPlot.jl")
end


text:
module DelayNetTools

using Random
# GLMakie, AbstractPlotting

include("types.jl")
export Spike, Synapse, Neuron, SimpleNeuronType, SimpleNeuron, SRMNeuron, DelayLine
export ModelParams, TrialParams, Experiment, Model, ExperimentOutput
export CategoricalSpiker

include("signalgen.jl")
export denseperiodic, sparseperiodic
export densepoisson, sparsepoisson, sparserefractorypoisson
export channeldup, channelscatter, stochasticblock, genrandinput

include("victor.jl")

# include("auditory.jl") # needs auditory filters

include("neuron.jl")
export izhupdate_euler, izhupdate_rk4
export resetsrmneuron, srmneuronupdate
export clearbuffer, pushevent, advance  # revise so have !

include("netgen.jl")
export analyzeblock, delayhistvals, sampletargets, genpatch, gensyn

include("io.jl")

include("util.jl")
export runexperiment, stageexperiment, pgsearch

include("analysis.jl")
export populationrate, populationpst, imagefromspikes
export sparseac, denseac, densehalfac, corranalysis
export pstspikes #, psthistograms

include("delayml.jl")
export trainclassifier, test, trainmulticlassifier, testmulti

end

max_stars_repo_name: analogouscircuit/DelayNetTools

max_issues_repo_name: analogouscircuit/DelayNetTools

max_forks_repo_name: analogouscircuit/DelayNetTools


text:
function iex_gibbs_insert_delete_update!(
    S_curr::InteractionSequence,
    S_prop::InteractionSequence,
    γ_curr::Float64,
    i::Int,
    posterior::SIM,
    mcmc::SimMcmcInsertDeleteGibbs,
    P::CumCondProbMatrix,
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool
    )
    V = posterior.V

    @inbounds I_tmp = S_prop[i]
    n = length(I_tmp)

    δ = rand(1:mcmc.ν_ed)

    d = rand(0:min(n,δ))
    m = n + δ_tmp - 2*d

    ind_del = view(mcmc.ind_del, 1:d)
    ind_add = view(mcmc.ind_add, 1:(δ_tmp-d))
    vals_del = view(mcmc.vals, 1:d)

    # Sample indexing info and new entries (all in-place)
    StatsBase.seqsample_a!(1:n, ind_del)
    StatsBase.seqsample_a!(1:m, ind_add)

    # *** HERE IS DIFFERENT FROM MODEL SAMPLER ***
    # The delete_insert_informed() function does the sampling + editing
    log_ratio = delete_insert_informed!(
        I_tmp,
        ind_del, ind_add, vals_del,
        P
    )

    log_ratio += log(min(n, δ)+1) - log(min(m, δ)+1)

    aux_model = SIM(
        S_prop, γ_curr,
        dist,
        V,
        K_inner,
        K_outer
    )

    if aux_init_at_prev
        @inbounds tmp = deepcopy(aux_data[end])
        draw_sample!(aux_data, aux_mcmc, aux_model, init=tmp)
    else
        draw_sample!(aux_data, aux_mcmc, aux_model)
    end

    aux_log_lik_ratio = -γ_curr * (
        sum(x -> dist(x, S_curr) for x in  aux_data)
        - sum(x -> dist(x, S_prop) for x in aux_data)
    )

    suff_stat_prop = sum(x -> dist(x, S_prop) for x in data)
    log_lik_ratio = -γ_curr * (
        suff_stat_prop - suff_stat_curr
    )

    log_multinom_ratio_term = log_multinomial_ratio(S_curr, S_prop)

    log_prior_ratio = -γ_prior * (
        dist(S_prop, mode_prior) - dist(S_curr, mode_prior)
    )

    # Log acceptance probability
    log_α = (
        log_lik_ratio +
        log_prior_ratio +
        aux_log_lik_ratio +
        log_ratio +
        log_multinom_ratio_term
    )

    # Accept-reject step. Use info in mcmc.ind_update to know which interaction are to be copied over
    if log(rand()) < log_α
        copy!(S_curr[i], I_tmp)
        return 1, suff_stat_prop
    else
        copy!(I_tmp, S_curr[i])
        return 0, suff_stat_curr
    end

end

function iex_gibbs_flip_update!(
    S_curr::InteractionSequence,
    S_prop::InteractionSequence,
    γ_curr::Float64,
    i::Int,
    posterior::SIM,
    mcmc::SimMcmcInsertDeleteGibbs,
    P::CumCondProbMatrix,
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool
    )
    V = posterior.V

    @inbounds I_tmp = S_prop[i]
    n = length(I_tmp)

    δ = rand(1:min(mcmc.ν_ed, n))

    ind = view(mcmc.ind_add, 1:δ)

    # Sample entries to flip
    StatsBase.seqsample_a!(1:n, ind)

    log_ratio = flip_informed_excl!(
        I_tmp,
        ind,
        P
    )

    aux_model = SIM(
        S_prop, γ_curr,
        dist,
        V,
        K_inner,
        K_outer
        )

    if aux_init_at_prev
        @inbounds tmp = deepcopy(aux_data[end])
        draw_sample!(aux_data, aux_mcmc, aux_model, init=tmp)
    else
        draw_sample!(aux_data, aux_mcmc, aux_model)
    end

    aux_log_lik_ratio = -γ_curr * (
        sum(x -> dist(x, S_curr) for x in  aux_data)
        - sum(x -> dist(x, S_prop) for x in aux_data)
    )

    suff_stat_prop = sum(x -> dist(x, S_prop) for x in data)
    log_lik_ratio = -γ_curr * (
        suff_stat_prop - suff_stat_curr
    )

    log_multinom_ratio_term = log_multinomial_ratio(S_curr, S_prop)

    log_prior_ratio = -γ_prior * (
        dist(S_prop, mode_prior) - dist(S_curr, mode_prior)
    )

    # Log acceptance probability
    log_α = (
        log_lik_ratio +
        log_prior_ratio +
        aux_log_lik_ratio +
        log_ratio +
        log_multinom_ratio_term
    )

    # Accept-reject step. Use info in mcmc.ind_update to know which interaction are to be copied over
    if log(rand()) < log_α
        copy!(S_curr[i], I_tmp)
        return 1, suff_stat_prop
    else
        copy!(I_tmp, S_curr[i])
        return 0, suff_stat_curr
    end

end

function iex_gibbs_scan!(
    S_curr::InteractionSequence,
    S_prop::InteractionSequence,
    γ_curr::Float64,
    posterior::SIM,
    mcmc::SimMcmcInsertDeleteGibbs,
    P::CumCondProbMatrix,
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool,
    count::Vector{Int},
    acc_count::Vector{Int}
    )

    for i in 1:length(S_curr)
        if rand() < 0.5
            was_acc, suff_stat_curr = iex_gibbs_insert_delete_update!(
                S_curr, S_prop,
                γ_curr, i,
                posterior, mcmc,
                P, aux_data,
                suff_stat_curr,
                aux_init_at_prev
            )
            count[1] += 1
            acc_count[1] += was_acc
        else
            was_acc, suff_stat_curr = iex_gibbs_insert_delete_update!(
                S_curr, S_prop,
                γ_curr, i,
                posterior, mcmc,
                P, aux_data,
                suff_stat_curr,
                aux_init_at_prev
            )
            count[2] += 1
            acc_count[2] += was_acc
        end
    end
end

function draw_sample_mode!(
    sample_out::Union{InteractionSequenceSample{Int}, SubArray},
    mcmc::Union{SimIexInsertDelete,SimIexInsertDeleteProportional},
    posterior::SimPosterior,
    γ_fixed::Float64;
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::InteractionSequence{Int}=sample_frechet_mean(posterior.data, posterior.dist),
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    )

    if loading_bar
        iter = Progress(
            length(sample_out) * lag + burn_in, # How many iters
            1,  # At which granularity to update loading bar
            "Chain for γ = $(γ_fixed) and n = $(posterior.sample_size) (mode conditional)....")  # Loading bar. Minimum update interval: 1 second
    end

    # Define aliases for pointers to the storage of current vals and proposals
    curr_pointers = mcmc.curr_pointers
    prop_pointers = mcmc.prop_pointers
    β = mcmc.β
    aux_mcmc = mcmc.aux_mcmc

    S_curr = InteractionSequence{Int}()
    S_prop = InteractionSequence{Int}()
    for i in 1:length(S_init)
        migrate!(S_curr, curr_pointers, i, 1)
        migrate!(S_prop, prop_pointers, i, 1)
        copy!(S_curr[i], S_init[i])
        copy!(S_prop[i], S_init[i])
    end

    γ_curr = γ_fixed

    sample_count = 1 # Keeps which sample to be stored we are working to get
    i = 0 # Keeps track all samples (included lags and burn_ins)

    tr_dim_count = 0
    tr_dim_acc_count = 0
    ed_count = 0
    ed_acc_count = 0
    flp_count = 0
    flp_acc_count = 0

    aux_data = [[Int[]] for i in 1:posterior.sample_size]
    # Initialise the aux_data
    aux_model = SIM(
        S_curr, γ_curr,
        posterior.dist,
        posterior.V,
        posterior.K_inner,
        posterior.K_outer)
    draw_sample!(aux_data, aux_mcmc, aux_model, burn_in=10000)

    # Evaluate sufficient statistic
    suff_stat_curr = mapreduce(
        x -> posterior.dist(S_curr, x),
        +,
        posterior.data
    )
    suff_stats = Float64[suff_stat_curr] # Storage for all sufficient stats (for diagnostics)
    P, vmap, vmap_inv = get_informed_proposal_matrix(posterior, mcmc.α)
    p_ins = get_informed_insertion_dist(posterior, mcmc.α)
    while sample_count ≤ length(sample_out)
        i += 1
        # Store value
        if (i > burn_in) & (((i - burn_in - 1) % lag)==0)
            @inbounds sample_out[sample_count] = deepcopy(S_curr)
            push!(suff_stats, suff_stat_curr)
            sample_count += 1
        end
        # W.P. do update move (accept-reject done internally by function call)
        if rand() < β
            if rand() < 0.5
                was_acc, suff_stat_curr = double_iex_multinomial_edit_accept_reject!(
                    S_curr, S_prop,
                    posterior, γ_curr,
                    mcmc, P,
                    aux_data,
                    suff_stat_curr,
                    aux_init_at_prev
                )
                ed_acc_count += was_acc
                ed_count += 1
            else
                was_acc, suff_stat_curr = double_iex_flip_accept_reject!(
                    S_curr, S_prop,
                    posterior, γ_curr,
                    mcmc, P,
                    aux_data,
                    suff_stat_curr,
                    aux_init_at_prev
                )
                flp_acc_count += was_acc
                flp_count += 1
            end
        # Else do trans-dim move. We will do accept-reject move here
        else
            was_acc, suff_stat_curr = double_iex_trans_dim_informed_accept_reject!(
                S_curr, S_prop,
                posterior, γ_curr,
                mcmc,
                p_ins,
                aux_data,
                suff_stat_curr,
                aux_init_at_prev
            )
            tr_dim_acc_count += was_acc
            tr_dim_count += 1
        end
        if loading_bar
            next!(iter)
        end
    end
    for i in 1:length(S_curr)
        migrate!(curr_pointers, S_curr, 1, 1)
        migrate!(prop_pointers, S_prop, 1, 1)
    end
    return (
        ed_count, ed_acc_count,
        flp_count, flp_acc_count,
        tr_dim_count, tr_dim_acc_count,
        suff_stats
    )
end

function accept_reject_mode!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    posterior::SimPosterior,
    γ_curr::Float64,
    mcmc::Union{SimIexInsertDelete,SimIexInsertDeleteProportional},
    P::CumCondProbMatrix,
    p_ins::Categorical,
    aux_data::InteractionSequenceSample{Int},
    acc_count::Vector{Int},
    count::Vector{Int},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool=false
    )

    β = mcmc.β
    if rand() < β
        if rand() < 0.5
            was_accepted, suff_stat_curr = double_iex_multinomial_edit_accept_reject!(
                S_curr, S_prop,
                posterior, γ_curr,
                mcmc, P,
                aux_data,
                suff_stat_curr,
                aux_init_at_prev
            )
            acc_count[1] += was_accepted
            count[1] += 1
        else
            was_accepted, suff_stat_curr = double_iex_flip_accept_reject!(
                S_curr, S_prop,
                posterior, γ_curr,
                mcmc, P,
                aux_data,
                suff_stat_curr,
                aux_init_at_prev
            )
            acc_count[2] += was_accepted
            count[2] += 1
        end
    else
        was_accepted, suff_stat_curr = double_iex_trans_dim_informed_accept_reject!(
            S_curr, S_prop,
            posterior, γ_curr,
            mcmc, p_ins,
            aux_data,
            suff_stat_curr,
            aux_init_at_prev
        )
        acc_count[3] += was_accepted
        count[3] += 1
    end
    return suff_stat_curr
end