text:
# --------------------------------------------------------------------------- #
# Labels -- a list of cell position coordinates

struct Label
  data :: Vector{Tuple{Int, Int}}
end

function Base.convert(::Type{Label}, mat::Matrix{Int})
  @assert (size(mat, 1) == 2)
  return Label([(mat[:,i]...) for i in 1:size(mat, 2)])
end

function Base.convert(::Type{Matrix{Int}}, lab::Label)
  return hcat([[coords...] for coords in lab.data]...)
end

Label() = Label(Tuple{Int, Int}[])
Label(mat::Matrix{Int}) = convert(Label, mat)

# --------------------------------------------------------------------------- #
# See tuples of images and labels as "labeled Images"

const LabeledGreyscaleImage = Tuple{GreyscaleImage, Label}
const LabeledRGBImage = Tuple{RGBImage, Label}
const LabeledImage = Union{LabeledGreyscaleImage, LabeledRGBImage}

# --------------------------------------------------------------------------- #
# Proximity maps
# Proxymaps are "washed out" 2d arrays with peaks at given label positions

function proxymap(width, height, lbls :: Vector{Label};
                  kernelsize :: Int = 7,
                  stddev :: Real = kernelsize/4,
                  peakheight :: Real = 100)

  # Create proximity maps from labels
  # Each point in the label gets converted to a gaussian activation

  n = length(lbls)

  # Pixel map
  prmap = zeros(Float32, width, height, 1, n)
  for i in 1:n
    for (x, y) in lbls[i]
      prmap[y, x, 1, i] = 1.
    end
  end

  # Convolve the pixel-map with a suitable kernel
  c = ceil(Int, kernelsize/2)
  kernel = Float32[ exp(- ((i - c)^2 + (j - c)^2) / (2stddev^2))
                    for i in 1:kernelsize, j in 1:kernelsize ]
  kernel = reshape(kernel, size(kernel)..., 1, 1)

  # Scale the output such that deviations from real labels are
  # punished more than deviations from the background
  pad = floor(Int, kernelsize/2)
  prmap = reshape(conv4(kernel, prmap, padding=pad), width, height, n)

  return peakheight * prmap
end

function proxymap(width, height, lbl :: Label; kwargs...)
  return reshape(proxymap(width, height, [lbl]; kwargs...), width, height)
end

# --------------------------------------------------------------------------- #
# Adjacency
# Calculates a list of smallest distances between the spots in two labels

function adjacency(dlbl, tlbl)
  ddata, tdata = convert.(Matrix{Int}, (dlbl, tlbl))
  if isempty(ddata) && isempty(tdata)
    return 0.
  elseif isempty(ddata) || isempty(tdata)
    return Inf
  end

  n, m = length(tlbl), length(dlbl)
  vals = zeros(n + m)
  for i in 1:n
    dists = sqrt.(sum(abs2, ddata .- tdata[:,i], 1))[1,:]
    vals[i] = minimum(dists)
  end
  for i in 1:m
    dists = sqrt.(sum(abs2, ddata[:,i] .- tdata, 1))[1,:]
    vals[n+i] = minimum(dists)
  end
  return vals
end

function adjacency(model, img :: Image, lbl :: Label; kwargs...)
  return adjacency(label(model, img; kwargs...), lbl)
end

# --------------------------------------------------------------------------- #
# Convenience functions

Base.size(l :: Label, args...; kwargs...) = size(l.data, args..., kwargs...)
Base.length(l :: Label, args...; kwargs...) = length(l.data, args..., kwargs...)
Base.getindex(l :: Label, k :: Integer) = l.data[k]
Base.setindex!(l :: Label, v, k :: Integer) = setindex!(l.data, (v...), k)
Base.endof(l :: Label) = length(l)

Base.start(l::Label) = start(l.data)
Base.next(l::Label, state) = next(l.data, state)
Base.done(l::Label, state) = done(l.data, state)
Base.eltype(::Type{Label}) = Tuple{Int, Int}

Base.push!(l :: Label, coords) = push!(l.data, coords)
Base.pop!(l :: Label) = pop!(l.data)
Base.deleteat!(l :: Label, idx) = deleteat!(l.data, idx)
Base.join(l1 :: Label, l2 :: Label) = Label([l1.data; l2.data])

Base.isempty(l :: Label) = isempty(l.data)

function crop(lbl :: Label, x :: Int, y :: Int, w :: Int, h :: Int)
  o = (x-1, y-1)
  data = [z .- o for z in lbl.data if x <= z[1] < x+w && y <= z[2] < y+h]
  return Label(data)
end

crop(lmg :: LabeledImage, x, y, w, h) = crop.(lmg, x, y, w, h)

function lblscale(lbl :: Label, s :: Real)
  return Label([round.(Int, s .* c) for c in lbl])
end

max_stars_repo_name: tscode/DCellC

max_issues_repo_name: tscode/DCellC


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using JuliaFEM, FEMBase, LinearAlgebra, Materials, DelimitedFiles

include("continuum.jl")

X = Dict(
    1 => [0.0, 0.0, 0.0],
    2 => [1.0, 0.0, 0.0],
    3 => [1.0, 1.0, 0.0],
    4 => [0.0, 1.0, 0.0],
    5 => [0.0, 0.0, 1.0],
    6 => [1.0, 0.0, 1.0],
    7 => [1.0, 1.0, 1.0],
    8 => [0.0, 1.0, 1.0])

body_element = Element(Hex8, (1, 2, 3, 4, 5, 6, 7, 8))
body_elements = [body_element]
update!(body_elements, "geometry", X)
update!(body_elements, "youngs modulus", 200.0e3)
update!(body_elements, "poissons ratio", 0.3)
update!(body_elements, "yield stress", 100.0)
update!(body_elements, "K_n", 100.0)
update!(body_elements, "n_n", 10.0)
update!(body_elements, "C_1", 10000.0)
update!(body_elements, "D_1", 100.0)
update!(body_elements, "C_2", 50000.0)
update!(body_elements, "D_2", 1000.0)
update!(body_elements, "Q", 50.0)
update!(body_elements, "b", 0.1)

bc_element_1 = Element(Poi1, (1,))
bc_element_2 = Element(Poi1, (2,))
bc_element_3 = Element(Poi1, (3,))
bc_element_4 = Element(Poi1, (4,))
bc_element_5 = Element(Poi1, (5,))
bc_element_6 = Element(Poi1, (6,))
bc_element_7 = Element(Poi1, (7,))
bc_element_8 = Element(Poi1, (8,))

bc_elements = [bc_element_1, bc_element_2, bc_element_3, bc_element_4,
               bc_element_5, bc_element_6, bc_element_7, bc_element_8]
update!(bc_elements, "geometry", X)
for element in (bc_element_1, bc_element_2, bc_element_3, bc_element_4)
    update!(element, "displacement 3", 0.0)
end
for element in (bc_element_5, bc_element_6, bc_element_7, bc_element_8)
    update!(element, "displacement 3", 0.0 => 0.0)
    update!(element, "displacement 3", 1.0 => 5.0e-3)
    update!(element, "displacement 3", 3.0 => -5.0e-3)
    update!(element, "displacement 3", 5.0 => 5.0e-3)
    update!(element, "displacement 3", 7.0 => -5.0e-3)
    update!(element, "displacement 3", 9.0 => 5.0e-3)
    update!(element, "displacement 3", 10.0 => 0.0)
end

update!(bc_element_1, "displacement 1", 0.0)
update!(bc_element_1, "displacement 2", 0.0)
update!(bc_element_2, "displacement 2", 0.0)
update!(bc_element_4, "displacement 1", 0.0)

update!(bc_element_5, "displacement 1", 0.0)
update!(bc_element_5, "displacement 2", 0.0)
update!(bc_element_6, "displacement 2", 0.0)
update!(bc_element_8, "displacement 1", 0.0)

#update!(bc_element_5, "displacement 1", 0.0)
#update!(bc_element_5, "displacement 2", 0.0)
#update!(bc_element_5, "displacement 3", 0.0 => 0.0)
#update!(bc_element_5, "displacement 3", 1.0 => 1.0e-3)

# Initialize material model to integration points
for ip in get_integration_points(body_element)
    mat = Material(Chaboche, tuple())
    mat.dtime = 0.05
    Materials.initialize!(mat, body_element, ip, 0.0)
    ip.fields["material"] = field(mat)
end

body = Problem(Continuum3D, "1 element problem", 3)
bc = Problem(Dirichlet, "fix displacement", 3, "displacement")
add_elements!(body, body_elements)
add_elements!(bc, bc_elements)
analysis = Analysis(Nonlinear, "solve problem")
# xdmf = Xdmf("results"; overwrite=true)
# add_results_writer!(analysis, xdmf)
add_problems!(analysis, body, bc)
# time_end = 1.0
time_end = 10.0
dtime = 0.05

for problem in get_problems(analysis)
    FEMBase.initialize!(problem, analysis.properties.time)
end

while analysis.properties.time < time_end
    analysis.properties.time += dtime
    update!(body_element, "displacement", analysis.properties.time => Dict(j => zeros(3) for j in 1:8))
    @info("time = $(analysis.properties.time)")
    for element in body_elements
        for ip in get_integration_points(element)
            material = ip("material", analysis.properties.time)
            preprocess_analysis!(material, element, ip, analysis.properties.time)
        end
    end
    run!(analysis)
    for element in body_elements
        for ip in get_integration_points(element)
            material = ip("material", analysis.properties.time)
            postprocess_analysis!(material, element, ip, analysis.properties.time)
        end
    end
    # update material internal parameters
end

# close(xdmf)

using Plots
if true
    ip1 = first(get_integration_points(body_element))
    t = range(0, stop=time_end, length=Int(time_end/dtime)+1)
    s11(t) = ip1("stress", t)[1]
    s22(t) = ip1("stress", t)[2]
    s33(t) = ip1("stress", t)[3]
    s12(t) = ip1("stress", t)[4]
    s23(t) = ip1("stress", t)[5]
    s31(t) = ip1("stress", t)[6]
    e11(t) = ip1("strain", t)[1]
    e22(t) = ip1("strain", t)[2]
    e33(t) = ip1("strain", t)[3]
    s(t) = ip1("stress", t)
    function vmis(t)
        s11, s22, s33, s12, s23, s31 = ip1("stress", t)
        return sqrt(1/2*((s11-s22)^2 + (s22-s33)^2 + (s33-s11)^2 + 6*(s12^2+s23^2+s31^2)))
    end
    path = joinpath("one_elem_disp_chaboche", "unitelement_results.rpt")
    data = readdlm(path, Float64; skipstart=4)
    t_ = data[:,1]
    s11_ = data[:,2]
    s12_ = data[:,3]
    s13_ = data[:,4]
    s22_ = data[:,5]
    s23_ = data[:,6]
    s33_ = data[:,7]
    e11_ = data[:,8]
    e12_ = data[:,9]
    e13_ = data[:,10]
    e22_ = data[:,11]
    e23_ = data[:,12]
    e33_ = data[:,13]
    plot(e11.(t), s11.(t), label="\$\\sigma_{11}\$", legend=:topleft,
        fg_legend=:transparent, bg_legend=:transparent)
    plot!(e22.(t), s22.(t), label="\$\\sigma_{22}\$")
    plot!(e33.(t), s33.(t), linecolor=:red, label="\$\\sigma_{33}\$")
    plot!(e11_, s11_, ls=:dash, label="\$\\sigma_{11} \\quad \\mathrm{Commercial}\$")
    plot!(e22_, s22_, ls=:dash, label="\$\\sigma_{22} \\quad \\mathrm{Commercial}\$")
    plot!(e33_, s33_, linecolor=:black, lw=1, ls=:dash,
        label="\$\\sigma_{33} \\quad \\mathrm{Commercial}\$")
    title!("Chaboche plasticity model\nOne element model with uniaxial stress")
    # xlabel!("\$\\varepsilon\$")
    # ylabel!("\$\\sigma\$")
    # labels = ["s11" "s22" "s33" "s12" "s23" "s31"]
    # plot(t, s11, title="stress at integration point 1", label="s11")
    # plot!(t, s22, label="s22")
    # plot!(t, s33, label="s33")
    # plot!(t, s12, label="s12")
    # plot!(t, s23, label="s23")
    # plot!(t, s31, label="s31")
end

ext: jl


text:
module LinearCovarianceModels

export LCModel, dim,
    # families
    generic_subspace, generic_diagonal, toeplitz, tree, trees,
    # ML degree
    ml_degree_witness, MLDegreeWitness,
    model, parameters, solutions, is_dual, ml_degree, verify,
    # solve specific instance
    mle, critical_points, covariance_matrix, logl, gradient_logl, hessian_logl, classify_point,
    # MLE helper
    mle_system, dual_mle_system, mle_system_and_start_pair, dual_mle_system_and_start_pair,
    # helpers
    vec_to_sym, sym_to_vec

using LinearAlgebra

import HomotopyContinuation
import DynamicPolynomials
import Distributions: Normal

const HC = HomotopyContinuation
const DP = DynamicPolynomials

include("tree_data.jl")

outer(A) = A*A'

"""
    rand_pos_def(n)

Create a random positive definite `n × n` matrix. The matrix is generated
by first creating a `n × n` matrix `X` where each entry is independently drawn
from the `Normal(μ=0, σ²=1.0)` distribution. Then `X*X'./n` is returned.
"""
rand_pos_def(n) = outer(rand(Normal(0, 1.0), n, n)) ./ n

n_vec_to_sym(k) = div(-1 + round(Int, sqrt(1+8k)), 2)
n_sym_to_vec(n) = binomial(n+1,2)

"""
    vec_to_sym(v)

Converts a vector `v` to a symmetrix matrix by filling the lower triangular
part columnwise.

### Example
```
julia> v = [1,2,3, 4, 5, 6];
julia> vec_to_sym(v)
3×3 Array{Int64,2}:
 1  2  3
 2  4  5
 3  5  6
 ```
"""
function vec_to_sym(s)
    n = n_vec_to_sym(length(s))
    S = Matrix{eltype(s)}(undef, n, n)
    l = 1
    for i in 1:n, j in i:n
        S[i,j] = S[j,i] = s[l]
        l += 1
    end
    S
end

"""
    sym_to_vec(S)

Converts a symmetric matrix `S` to a vector by filling the vector with lower triangular
part iterating columnwise.
"""
sym_to_vec(S) = (n = size(S, 1); [S[i,j] for i in 1:n for j in i:n])

"""
    LCModel(Σ::Matrix{<:DP.AbstractPolynomialLike})

Create a linear covariance model from the parameterization `Σ`.
This uses as input a matrix of polynomials created by the [`DynamicPolynomials`](https://github.com/JuliaAlgebra/DynamicPolynomials.jl) package.

## Example

```
using DynamicPolynomials # load polynomials package

# use DynamicPolynomials to create variables θ₁, θ₂, θ₃.
@polyvar θ[1:3]

# create our model as matrix of DynamicPolynomials
Σ = [θ[1] θ[2] θ[3]; θ[2] θ[1] θ[2]; θ[3] θ[2] θ[1]]

# create model
model = LCModel(Σ)
```
"""
struct LCModel{T1<:DP.AbstractPolynomialLike, T2<:Number}
    Σ::Matrix{T1}
    B::Vector{Matrix{T2}}

    function LCModel(Σ::Matrix{T1}, B::Vector{Matrix{T2}}) where {T1,T2}
        all(DP.maxdegree.(vec(Σ)) .<= 1) || throw(ArgumentError("Input is not a linear covariance model"))
        issymmetric(Σ) || throw(ArgumentError("Input is not a symmetric matrix!"))
        new{T1,T2}(Σ, B)
    end
end
LCModel(Σ::Matrix) = LCModel(Σ, get_basis(Σ))
LCModel(Σ::AbstractMatrix) = LCModel(Matrix(Σ .+ false))

function get_basis(Σ)
    vars = DP.variables(vec(Σ))
    map(1:length(vars)) do i
        [p(vars[i] => 1,
           vars[1:i-1]=>zeros(Int, max(i-1,0)),
           vars[i+1:end]=>zeros(Int, max(length(vars)-i,0))) for p in Σ]
    end
end

Base.size(M::LCModel) = (size(M.Σ, 1), length(M.B))
Base.size(M::LCModel, i::Int) = size(M)[i]
function Base.show(io::IO, M::LCModel)
    println(io, "$(dim(M))-dimensional LCModel:")
    Base.print_matrix(io, M.Σ)
end
Base.broadcastable(M::LCModel) = Ref(M)

"""
    dim(M::LCModel)

Returns the dimension of the model.
"""
dim(M::LCModel) = length(M.B)

"""
    toeplitz(n::Integer)

Returns a symmetric `n×n` toeplitz matrix.
"""
function toeplitz(n::Integer)
    DP.@polyvar θ[1:n]
    sum(0:n-1) do i
        if i == 0
            θ[1] .* diagm(0 => ones(n))
        else
            θ[i+1] .* (diagm(i => ones(n-i)) + diagm(-i => ones(n-i)))
        end
    end |> LCModel
end

"""
    tree(n, id::String)

Get the covariance matrix corresponding to the tree with the given `id` on `n` leaves.
Returns `nothing` if the tree was not found.

## Example
```
julia> tree(4, "{{1, 2}, {3, 4}}")
4×4 Array{PolyVar{true},2}:
 t₁  t₅  t₇  t₇
 t₅  t₂  t₇  t₇
 t₇  t₇  t₃  t₆
 t₇  t₇  t₆  t₄
 ```
"""
function tree(n::Integer, id::String)
    4 ≤ n ≤ 7 || throw(ArgumentError("Only trees with 4 to 7 leaves are supported."))
    for data in TREE_DATA
        if data.n == n && data.id == id
            return make_tree(data.tree)
        end
    end
    nothing
end

function make_tree(tree::Matrix{Symbol})
    var_names = sort(unique(vec(tree)))
    D = Dict(map(v -> (v, DP.PolyVar{true}(String(v))), var_names))
    LCModel(map(v -> D[v], tree))
end

"""
    trees(n)

Return all trees with `n` leaves as a tuple (id, tree).
"""
function trees(n::Int)
    4 ≤ n ≤ 7 || throw(ArgumentError("Only trees with 4 to 7 leaves are supported."))
    map(d -> (id=d.id, tree=make_tree(d.tree)), filter(d -> d.n == n, TREE_DATA))
end

"""
    generic_subspace(n::Integer, m::Integer); pos_def::Bool=true)

Generate a generic family of symmetric ``n×n`` matrices living in an ``m``-dimensional
subspace. If `pos_def` is `true` then positive definite matrices are used as a basis.
"""
function generic_subspace(n::Integer, m::Integer; pos_def::Bool=true)
    m ≤ binomial(n+1,2) || throw(ArgumentError("`m=$m` is larger than the dimension of the space."))
    DP.@polyvar θ[1:m]
    if pos_def
        LCModel(sum(θᵢ .* rand_pos_def(n) for θᵢ in θ))
    else
        LCModel(sum(θᵢ .* Symmetric(randn(n,n)) for θᵢ in θ))
    end
end

"""
    generic_diagonal(n::Integer, m::Integer)

Generate a generic family of ``n×n`` diagonal matrices living in an ``m``-dimensional
subspace.
"""
function generic_diagonal(n::Integer, m::Integer)
    m ≤ n || throw(ArgumentError("`m=$m` is larger than the dimension of the space."))
    DP.@polyvar θ[1:m]
    LCModel(sum(θᵢ .* diagm(0 => randn(n)) for θᵢ in θ))
end

"""
    mle_system(M::LCModel)

Generate the MLE system corresponding to the family of covariances matrices
parameterized by `Σ`.
Returns the named tuple `(system, variables, parameters)`.
"""
function mle_system(M::LCModel)
    Σ = M.Σ
    θ = DP.variables(vec(Σ))
    m = DP.nvariables(θ)
    n = size(Σ, 1)
    N = binomial(n+1,2)

    DP.@polyvar k[1:N] s[1:N]

    K, S = vec_to_sym(k), vec_to_sym(s)
    l = -tr(K * Σ) + tr(S * K * Σ * K)
    ∇l = DP.differentiate(l, θ)
    KΣ_I = vec(K * Σ - Matrix(I, n,n))
    (system=[∇l; KΣ_I], variables=[θ; k], parameters=s)
end

"""
    dual_mle_system(M::LCModel)

Generate the dual MLE system corresponding to the family of covariances matrices
parameterized by `Σ`.
Returns the named tuple `(system, variables, parameters)`.
"""
function dual_mle_system(M::LCModel)
    Σ = M.Σ
    θ = DP.variables(vec(Σ))
    m = DP.nvariables(θ)
    n = size(Σ, 1)
    N = binomial(n+1,2)

    DP.@polyvar k[1:N] s[1:N]

    K, S = vec_to_sym(k), vec_to_sym(s)
    l = -tr(K * Σ) + tr(S * Σ)
    ∇l = DP.differentiate(l, θ)
    KΣ_I = vec(K * Σ - Matrix(I, n,n))
    (system=[∇l; KΣ_I], variables=[θ; k], parameters=s)
end

"""
    mle_system_and_start_pair(M::LCModel)

Generate the mle_system and a corresponding start pair `(x₀,p₀)`.
"""
function mle_system_and_start_pair(M::LCModel)
    system, vars, params = mle_system(M)
    θ = DP.variables(vec(M.Σ))
    θ₀ = randn(ComplexF64, length(θ))
    Σ₀ = [p(θ => θ₀) for p in M.Σ]
    K₀ = inv(Σ₀)
    x₀ = [θ₀; sym_to_vec(K₀)]
    A, b = HC.linear_system(DP.subs.(system[1:length(x₀)], Ref(vars => x₀)), params)
    p₀ = A \ b

    (system=system, x₀=x₀, p₀=p₀, variables=vars, parameters=params)
end

"""
    dual_mle_system_and_start_pair(M::LCModel)

Generate the dual MLE system and a corresponding start pair `(x₀,p₀)`.
"""
function dual_mle_system_and_start_pair(M::LCModel)
    system, vars, params = dual_mle_system(M)
    θ = DP.variables(vec(M.Σ))
    θ₀ = randn(ComplexF64, length(θ))
    Σ₀ = [p(θ => θ₀) for p in M.Σ]
    K₀ = inv(Σ₀)
    x₀ = [θ₀; sym_to_vec(K₀)]
    A, b = HC.linear_system(DP.subs.(system[1:length(x₀)], Ref(vars => x₀)), params)
    p₀ = A \ b

    (system=system, x₀=x₀, p₀=p₀, variables=vars, parameters=params)
end

"""
    MLDegreeWitness

Data structure holding an MLE model. This also holds a set of solutions for a generic instance,
which we call a witness.
"""
struct MLDegreeWitness{T1, T2, V<:AbstractVector}
    model::LCModel{T1,T2}
    solutions::Vector{V}
    p::Vector{ComplexF64}
    dual::Bool
end

function MLDegreeWitness(Σ::AbstractMatrix, solutions, p, dual)
    MLDegreeWitness(LCModel(Σ), solutions, p, dual)
end

function Base.show(io::IO, R::MLDegreeWitness)
    println(io, "MLDegreeWitness:")
    println(io, " • ML degree → $(length(R.solutions))")
    println(io, " • model dimension → $(dim(model(R)))")
    println(io, " • dual → $(R.dual)")
end

"""
    model(W::MLDegreeWitness)

Obtain the model corresponding to the `MLDegreeWitness` `W`.
"""
model(R::MLDegreeWitness) = R.model

"""
    solutions(W::MLDegreeWitness)

Obtain the witness solutions corresponding to the `MLDegreeWitness` `W`
with given parameters.
"""
solutions(W::MLDegreeWitness) = W.solutions

"""
    parameters(W::MLDegreeWitness)

Obtain the parameters of the `MLDegreeWitness` `W`.
"""
parameters(W::MLDegreeWitness) = W.p

"""
    is_dual(W::MLDegreeWitness)

Indicates whether `W` is a witness for the dual MLE.
"""
is_dual(W::MLDegreeWitness) = W.dual

"""
    ml_degree(W::MLDegreeWitness)

Returns the ML degree.
"""
ml_degree(W::MLDegreeWitness) = length(solutions(W))

"""
    ml_degree_witness(Σ::LCModel; ml_degree=nothing, max_tries=5, dual=false)

Compute a [`MLDegreeWitness`](@ref) for a given model Σ. If the ML degree is already
known it can be provided to stop the computations early. The stopping criterion is based
on a heuristic, `max_tries` indicates how many different parameters are tried a most until
an agreement is found.
"""
ml_degree_witness(Σ; kwargs...) = ml_degree_witness(LCModel(Σ); kwargs...)
function ml_degree_witness(M::LCModel; ml_degree=nothing, max_tries = 5, dual=false)
    if dual
        F, x₀, p₀, x, p = dual_mle_system_and_start_pair(M)
    else
        F, x₀, p₀, x, p = mle_system_and_start_pair(M)
    end
    result = HC.monodromy_solve(F, x₀, p₀; target_solutions_count=ml_degree,
                                            parameters=p, max_loops_no_progress=5)
    if HC.nsolutions(result) == ml_degree
        return MLDegreeWitness(M, HC.solutions(result), result.parameters, dual)
    end

    best_result = result
    result_agreed = false
    for i in 1:max_tries
        q₀ = randn(ComplexF64, length(p₀))
        S_q₀ = HC.solutions(HC.solve(F, HC.solutions(result); parameters=p, start_parameters=p₀, target_parameters=q₀))
        new_result = HC.monodromy_solve(F, S_q₀, q₀; parameters=p, max_loops_no_progress=3)
        if HC.nsolutions(new_result) == HC.nsolutions(best_result)
            result_agreed = true
            break
        elseif HC.nsolutions(new_result) > HC.nsolutions(best_result)
            best_result = new_result
        end
    end
    MLDegreeWitness(M, HC.solutions(best_result), best_result.parameters, dual)
end

"""
    verify(W::MLDegreeWitness; trace_tol=1e-5, options...)

Tries to verify that the computed ML degree witness is complete, i.e., that
the ML degree is correct. This uses the [`verify_solution_completeness`](https://www.juliahomotopycontinuation.org/HomotopyContinuation.jl/stable/monodromy/#HomotopyContinuation.verify_solution_completeness)
of HomotopyContinuation.jl. All caveats mentioned there apply.
The `options` are also passed to `verify_solution_completeness`.
"""
function verify(W::MLDegreeWitness; trace_tol=1e-5, kwargs...)
    if W.dual
        F, var, params = dual_mle_system(model(W))
    else
        F, var, params = mle_system(model(W))
    end
    HC.verify_solution_completeness(F, solutions(W), parameters(W); parameters=params, trace_tol=trace_tol, kwargs...)
end

"""
    critical_points(W::MLDegreeWitness, S::AbstractMatrix;
            only_positive_definite=true, only_non_negative=false,
            options...)

Compute all critical points to the MLE problem of `W` for the given sample covariance matrix
`S`. If `only_positive_definite` is `true` only positive definite solutions are considered.
If `only_non_negative` is `true` only non-negative solutions are considered.
The `options` are argument passed to the [`solve`](https://www.juliahomotopycontinuation.org/HomotopyContinuation.jl/stable/solving/#HomotopyContinuation.solve) routine from `HomotopyContinuation.jl`.
"""
function critical_points(W::MLDegreeWitness, S::AbstractMatrix;
               only_positive_definite=true, only_non_negative=false, kwargs...)
    issymmetric(S) || throw("Sample covariance matrix `S` is not symmetric. Consider wrapping it in `Symmetric(S)` to enforce symmetry.")
    if W.dual
        F, var, params = dual_mle_system(model(W))
    else
        F, var, params = mle_system(model(W))
    end
    result = HC.solve(F, solutions(W); parameters=params,
                         start_parameters=W.p,
                         target_parameters=sym_to_vec(S),
                         kwargs...)

    m = size(model(W), 2)
    θs = map(s -> s[1:m], HC.real_solutions(result))
    if only_positive_definite
        filter!(θs) do θ
            isposdef(covariance_matrix(model(W), θ))
        end
    end

    if only_non_negative
        filter!(θs) do θ
            all(covariance_matrix(model(W), θ) .≥ 0)
        end
    end

    res = map(θs) do θ
        (θ, logl(model(W), θ, S), classify_point(model(W), θ, S))
    end

    if !isempty(res)
        best_val = maximum(θs)
        for i in 1:length(res)
            if first(res[i]) == best_val && res[i][3] == :local_maximum
                res[i] = (res[i][1], res[i][2], :global_maximum)
            end
        end
    end
    res
end

"""
    covariance_matrix(M::LCModel, θ)

Compute the covariance matrix corresponding to the value of `θ` and the given model
`M`.
"""
covariance_matrix(W::MLDegreeWitness, θ) = covariance_matrix(model(W), θ)
covariance_matrix(M::LCModel, θ) = sum(θ[i] * M.B[i] for i in 1:size(M,2))

"""
    logl(M::LCModel, θ, S::AbstractMatrix)

Evaluate the log-likelihood ``log(det(Σ⁻¹)) - tr(SΣ⁻¹)`` of the MLE problem.
"""
function logl(M::LCModel, θ, S::AbstractMatrix)
    logl(covariance_matrix(M, θ), S)
end
logl(Σ::AbstractMatrix, S::AbstractMatrix) = -logdet(Σ) - tr(S*inv(Σ))

"""
    gradient_logl(M::LCModel, θ, S::AbstractMatrix)

Evaluate the gradient of the log-likelihood ``log(det(Σ⁻¹)) - tr(SΣ⁻¹)`` of the MLE problem.
"""
gradient_logl(M::LCModel, θ, S::AbstractMatrix) = gradient_logl(M.B, θ, S)
function gradient_logl(B::Vector{<:Matrix}, θ, S::AbstractMatrix)
    Σ = sum(θ[i] * B[i] for i in 1:length(B))
    Σ⁻¹ = inv(Σ)
    map(1:length(B)) do i
        -tr(Σ⁻¹ * B[i]) + tr(S *  Σ⁻¹ * B[i] * Σ⁻¹)
    end
end

"""
    hessian_logl(M::LCModel, θ, S::AbstractMatrix)

Evaluate the hessian of the log-likelihood ``log(det(Σ⁻¹)) - tr(SΣ⁻¹)`` of the MLE problem.
"""
hessian_logl(M::LCModel, θ, S::AbstractMatrix) = hessian_logl(M.B, θ, S)
function hessian_logl(B::Vector{<:Matrix}, θ, S::AbstractMatrix)
    m = length(B)
    Σ = sum(θ[i] * B[i] for i in 1:m)
    Σ⁻¹ = inv(Σ)
    H = zeros(eltype(Σ), m, m)
    for i in 1:m, j in i:m
        kernel = Σ⁻¹ * B[i] * Σ⁻¹ * B[j]
        H[i,j] = H[j,i] = tr(kernel) - 2tr(S * kernel * Σ⁻¹)
    end
    Symmetric(H)
end

"""
    classify_point(M::LCModel, θ, S::AbstractMatrix)

Classify the critical point `θ` of the log-likelihood function.
"""
function classify_point(M::LCModel, θ, S::AbstractMatrix)
    H = hessian_logl(M, θ, S)
    emin, emax = extrema(eigvals(H))
    if emin < 0 && emax < 0
        :local_maximum
    elseif emin > 0 && emax > 0
        :local_minimum
    else
        :saddle_point
    end
end

"""
    mle(W::MLDegreeWitness, S::AbstractMatrix; only_positive_definite=true, only_positive=false)

Compute the MLE for the matrix `S` using the MLDegreeWitness `W`.
Returns the parameters for the MLE covariance matrix or `nothing` if no solution was found
satisfying the constraints (see options below).

## Options

* `only_positive_definite`: controls whether only positive definite
covariance matrices should be considered.
* `only_positive`: controls whether only (entrywise) positive covariance matrices
should be considered.
"""
function mle(W::MLDegreeWitness, S::AbstractMatrix; only_positive_definite=true, only_positive=false, kwargs...)
    is_dual(W) && throw(ArgumentError("`mle` is currently only supported for MLE not dual MLE."))
    results = critical_points(W, S; kwargs...)
    ind = findfirst(r -> r[3] == :global_maximum, results)
    isnothing(ind) ? nothing : results[ind][1]
end

end # module


text:
using ArgParse

s = ArgParseSettings()

@add_arg_table s begin
  "model"
    help = "Path to AdaGram model file"
    arg_type = String
    required = true
end

args = parse_args(ARGS, s)

println("started with $args")

import JSON
import AdaGram

println("loading model...")
vm, dict = AdaGram.load_model(args["model"])
println("done.")

all_words = [
    "альбом",
    "байка",
    "баян",
    "билет",
    "блок",
    "бомба",
    "борщ",
    "бухгалтер",
    "бык",
    "вата",
    "вешалка",
    "вилка",
    "винт",
    "воск",
    "горшок",
    ]

output = Dict()

for word  in all_words
    println(word)
    if !haskey(dict.word2id, word)
        println("missing")
        continue
    end
    output[word] = Any[]
    for (idx, weight) in enumerate(AdaGram.expected_pi(vm, dict.word2id[word]))
        if weight > 0.001
            sense = Dict()
            push!(output[word], sense)
            sense["weight"] = weight
            neighbours = AdaGram.nearest_neighbors(vm, dict, word, idx, 5)
            sense["neighbours"] = neighbours
            println("$idx: $weight")
            println("$neighbours")
        end
    end
end

JSON.print(open("adagram.json", "w"), output)

max_stars_repo_name: lopuhin/WSI-LDA

max_issues_repo_name: lopuhin/WSI-LDA

max_forks_repo_name: lopuhin/WSI-LDA


text:
@testset "lessorqual.jl" begin
    @testset "LessOrEqualConstant()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 6, "x", trailer)

        constraint = SeaPearl.LessOrEqualConstant(x, 3, trailer)

        @test constraint in x.onDomainChange
        @test constraint.active.value
    end
    @testset "propagate!(::LessOrEqualConstant)" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 6, "x", trailer)

        constraint = SeaPearl.LessOrEqualConstant(x, 3, trailer)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test length(x.domain) == 2
        @test 2 in x.domain
        @test !(4 in x.domain)
        @test prunedDomains == SeaPearl.CPModification("x" => [4, 5, 6])

        y = SeaPearl.IntVar(2, 2, "y", trailer)

        cons2 = SeaPearl.LessOrEqualConstant(y, 1, trailer)

        @test !SeaPearl.propagate!(cons2, toPropagate, prunedDomains)

        @test isempty(y.domain)

    end
    @testset "LessOrEqual()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 6, "x", trailer)
        y = SeaPearl.IntVar(2, 6, "y", trailer)

        constraint = SeaPearl.LessOrEqual(x, y, trailer)

        @test constraint in x.onDomainChange
        @test constraint in y.onDomainChange
        @test constraint.active.value
    end
    @testset "propagate!(::LessOrEqual)" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 6, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)

        constraint = SeaPearl.LessOrEqual(x, y, trailer)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test length(x.domain) == 2
        @test !(5 in x.domain)
        @test 3 in x.domain
        @test prunedDomains == SeaPearl.CPModification("x" => [4, 5, 6])

        z = SeaPearl.IntVar(1, 1, "z", trailer)
        constraint2 = SeaPearl.LessOrEqual(x, z, trailer)

        @test !SeaPearl.propagate!(constraint2, toPropagate, prunedDomains)

    end
end


text:
#!/usr/local/bin/julia

using Miners;
using ReinforcementLearning;

@time learn(Gatherers);

max_stars_repo_name: tudor-berariu/phd

max_issues_repo_name: tudor-berariu/phd

max_forks_repo_name: tudor-berariu/phd


text:
using TuringModels

Turing.setadbackend(:reverse_diff);
Turing.turnprogress(false);

# In Rethinking the model actually has priors that are Uniform[-Inf, Inf], or as the Stan manual (2.17.0, pp. 127) tells us:

# "A parameter declared without constraints is thus given a uniform prior
# on (−∞, ∞) ..."

# Setting -Inf and Inf in Turing doesn't work. So the below works of
# course better since we've restrained it to [-1,1]

@model m8_2(y) = begin
    σ ~ Uniform(0, Inf)
    α ~ Uniform(-Inf, Inf)

    for i ∈ 1:length(y)
        y[i] ~ Normal(α, σ)
    end
end

y = [-1,1];

# Sample

chns = sample(m8_2(y), NUTS(0.65), 1000)

# Results rethinking

m82rethinking = "
       mean   sd  5.5% 94.5% n_eff Rhat
a      9.22 0.15  9.00  9.46   205    1
bR    -0.20 0.08 -0.34 -0.07   192    1
bA    -1.95 0.24 -2.36 -1.59   203    1
bAR    0.40 0.14  0.19  0.63   186    1
sigma  0.95 0.05  0.88  1.04   361    1
";

# Describe the posterior samples

describe(chns)

# End of `08/m8.2t.jl`


text:
#using NeuralVerification: ReLU, Id, Network, Layer
using OVERT: symbols, expand
using SymEngine
# include("../nv/utils/activation.jl")
# include("../nv/utils/network.jl")
# include("../nv/utils/problem.jl")
# include("../nv/utils/util.jl")

# helper functions for the conversion to smt2
act_dict = Dict(ReLU()=> :relu)
function network2expr(nn::Network, state::Array{Symbol})
    exprs = symbols.(state) # starts out as [x1, x2] -> [w1*x1+b1, w2*x2+b2]
    for l in nn.layers
        println("activation is: $(l.activation)")
        if l.activation == Id()
            exprs = l.weights*exprs + l.bias
        else
            exprs = elementwise_apply(act_dict[l.activation], l.weights*exprs + l.bias)
        end
    end
    return Basic2Expr(exprs)
end

# --------stolen from certificateverification repo
# Generate a metric network which has ReLU activation on its ouput layer.
# the output will be interpreted as the diagonal of our M matrix
function random_network(layer_sizes, activations)
    layers = Layer[]
    for i = 1:(length(layer_sizes)-1)
        layer_in = layer_sizes[i]
        layer_out = layer_sizes[i+1]
        weights = randn(layer_out, layer_in)
        bias = randn(layer_out)
        # All ReLU including the last layer because we want
        # the output to be >= 0
        push!(layers, Layer(weights, bias, activations[i]))
    end
    return Network(layers)
end

# Apply a symbolic function elementwise to a list of Basics
function elementwise_apply(symbol::Symbol, items::Array)
    # Convert to a list of expressions
    expressions = Basic2Expr(items)
    # Wrap each expression in the function
    wrapped_expressions = [:($symbol($expr)) for expr in expressions]
    # Convert back to a list of basics
    return Basic.(wrapped_expressions)
end

function Basic2Expr(e)
    return Meta.parse.(string.(expand.(e)))
end
# --------------------------------------------------

function test_network2expr()
    network = random_network([2,2,2], [ReLU(), Id()])
    state = [:x1, :x2]
    exprs = network2expr(network, state)
    x_rand = rand(1,2)
    # check equivalency by subbing in values to Expr?
    map = Dict(zip(state, x_rand))
    eval_expr = [eval(substitute(e, map)) for e in exprs]
end

function timestamp_expr(e::Expr, N)
    """
    Convert x to x_$N and x_tp1 to x_$(N+1)
    """
end

function substitute(expression, map)
    # replace each key with each value in the expression
    e_string = string(expression)
    for (k,v) in map
        e_string = replace(e_string, string(k)=>string(v))
    end
    return Meta.parse(e_string)
end

function add_controller(u::Array{Symbol}, u_expr::Array{Expr}, x::Array{Symbol}, formula::SMTLibFormula, N::Int)
    # u is like [:u_1, :u_2]
    # x is like [:x_1, :x_2]
    # u_expr is like [:(W12*relu(W_11*x_1 + b1) + b2)] where W11, b1, etc. are constants like 3.453

    # u_var = add controller() [separate function]
    # timestamp inputs and outputs (maybe handle LHS and RHS separately...)
    # time inputs:
    u_timed = [Meta.parse("$(v)_$N") for v in u]
    x_timed =  [Meta.parse("$(v)_$N") for v in x]
    map = Dict(zip([u...,x...], [u_timed...,x_timed...]))
    # time exprs
    u_expr_timed = [substitute(e, map) for e in u_expr]
    # assertion = convert_any_constraint(Expr)
    for (i,u_i) in enumerate(u_timed)
        e = :($(u_i) == $(u_expr_timed[i]))
        #println("c.e = $e")
        assertion = assert_literal(e, formula.stats)
        push!(formula.formula, assertion) # add assertion to formula
    end
    return u_timed
end

function test_add_controller()
    u = [:u1, :u2]
    u_exprs = [:(x1 - W*x2 + relu(x3)), :(x1^2 - relu(x2))]
    x = [:x1, :x2, :x3]
    formula = SMTLibFormula()
    N = 2
    u_timed = add_controller(u, u_exprs, x, formula, N)
end

function add_dynamics(x::Array{Symbol}, u::Array{Symbol}, N::Int, dynamics_map::Dict, dt::T where T <: Real, formula::SMTLibFormula)
    """
    inputs:
        # u is like [:u1, :u2]
        # x is like [:x1, :x2]
        # dynamics map is like: Dict(:x1 => :x2, :x2 => :(u1^2 - 5*x1), )

    # state_vars = add dynamics [separate function]
    # timestamp inputs and outputs
    # construct RHS integration
    # assertion = convert_any_constraints(Expr)
    #add assertion to formula
    """
    # timestamp inputs
    u_timed = [Meta.parse("$(v)_$N") for v in u]
    x_timed =  [Meta.parse("$(v)_$N") for v in x]
    map = Dict(zip([u...,x...], [u_timed...,x_timed...]))
    for x_i in x
        dx = substitute(dynamics_map[x_i], map)
        x_i_n = Meta.parse("$(x_i)_$N")
        x_i_np1 = Meta.parse("$(x_i)_$(N+1)")
        # x_i_tp1 = x_t + xdot_t * dt
        e = :($(x_i_np1) == $(x_i_n) + $dx*$dt)
        println("d.e = $e")
        assertion = assert_literal(e, formula.stats)
        push!(formula.formula, assertion)
    end
    return [Meta.parse("$(v)_$(N+1)") for v in x]
end

# function define_relu()
#     return ["(define-fun relu ((arg Real)) Real (max arg 0))"]
# end

function gen_full_formula(formula::SMTLibFormula)
    pushfirst!(formula.formula, declare_reals(formula.stats)...)
    pushfirst!(formula.formula, define_relu())
    pushfirst!(formula.formula, header()...)
    push!(formula.formula, footer()...)
    return formula
end

function add_output_constraints!(formula::SMTLibFormula, c::Expr, x::Array{Symbol}, N::T where T <: Real)
    x_N = [Meta.parse("$(v)_$N") for v in x]
    map = Dict(zip([x...], [x_N...]))
    c_timed = substitute(c, map)
    push!(formula.formula, assert_literal(c_timed, formula.stats))
end

function run_dreal(full_fname::String; δ=0.001, jobs=1, dreal_path="/opt/dreal/4.20.12.1/bin/dreal")
    # returns "true" if property is unsat (holds)
    println("running dreal on file $(full_fname) with precision $δ and $jobs jobs.")
    result = read(`$dreal_path $full_fname --precision $δ --jobs $jobs`, String)
    @debug("result: ", result)
    # write result file
    println("Writing result to: $(full_fname)")
    write_result(full_fname, result)
    return occursin("unsat", result)
end

function add_output_constraints_and_check_property(formula::SMTLibFormula, output_constraints::Array{Expr}, x::Array{Symbol}, t::T where T<:Real; δ=0.001, jobs=1, dreal_path="/opt/dreal/4.21.06.1/bin/dreal", experiment_name="dreal_test", dirname="examples/jmlr/comparisons/smtlibfiles/")
    unsat = true
    # add output constraints
    ## create a separate file for each output constraint
    t_dreal_alone = 0
    for (i,c) in enumerate(output_constraints)
        formula_i = deepcopy(formula)
        add_output_constraints!(formula_i, c, x, t)
        # write to file
        formula_i = gen_full_formula(formula_i::SMTLibFormula)
        full_fname = write_to_file(formula_i::SMTLibFormula, experiment_name*"_constraint$(i)_time$t.smt2"; dirname=dirname*"smtlibfiles/")
        # call dreal on file
        t_dreal_start = time()
        unsat &= run_dreal(full_fname; δ=δ, jobs=jobs, dreal_path=dreal_path)
        dt_dreal = time() - t_dreal_start
        t_dreal_alone += dt_dreal
    end
    return unsat, t_dreal_alone
end


text:
# Reproduce Figure 2 of Ruigrok et al., J Seismol, 2017

import Pkg
Pkg.activate(joinpath(@__DIR__, ".."))

using Beamforming, Seis, Plots

"Compute all unique pairs of stations"
station_pairs(stas) = [(stas[i], stas[j]) for i in 1:length(stas) for j in i+1:length(stas)]

# Plotting defaults
default(fontfamily="Helvetica")

# Three ARCES stations.
# Coordinates from https://www.fdsn.org/media/meetings/2015/NORSAR_Network_Report.pdf
ara1 = Station(sta="ARA1", lon=25.5071, lat=69.5363)
ara2 = Station(sta="ARA2", lon=25.5078, lat=69.5338)
arb2 = Station(sta="ARB2", lon=25.5134, lat=69.5357)

stas = [ara1, ara2, arb2]
n = length(stas)
t = [Trace(0, 1, 0) for _ in 1:n]
t.sta = stas

# Station cartesian coordinates and mean longitude and latitude
x, y, z, lon, lat = Beamforming.array_geometry(t)

# Cartesian stations
cstas = [CartStation(x=1000xx, y=1000yy, z=1000zz, sta=sta)
         for (xx, yy, zz, sta) in zip(x, y, z, stas.sta)]

# Beamforming parameters for all
sx = sy = -1:0.01:1
smax = max(maximum.((sx, sy))...)
ds = step(sx)
ps = 0
θs = 0

# Plotting arguments
cpt = cgrad([:white, :blue, :cyan, :yellow, :orange, :red, :darkred])
kwargs = (c=cpt, clim=(-10,0), s_km=true, powscale=:dB)

# Single 5 Hz wave, absolute power
f = 5
df = 1
P_5 = Beamforming.crosscorrelation_array_response(stas, smax, ds, f, f, df, powf=abs)
arf_5 = plot(P_5; kwargs...)

# 3 – 7 Hz, absolute power
f1 = 3
f2 = 7
df = 0.1
P_3_7 = Beamforming.crosscorrelation_array_response(stas, smax, ds, f1, f2, df, powf=abs)
arf_3_7 = plot(P_3_7; kwargs...)

# 3 – 7 Hz, squared power
P_3_7² = Beamforming.crosscorrelation_array_response(stas, smax, ds, f1, f2, df, powf=abs2)
arf_3_7² = plot(P_3_7²; kwargs...)

# Map of array
pairs = station_pairs(cstas)
map = plot([first.(pairs).x last.(pairs).x]', [first.(pairs).y last.(pairs).y]',
    lc=:black, grid=false, framestyle=:box, aspect_ratio=:equal, primary=false,
    xlabel="Easting (m)", ylabel="Northing (m)")
scatter!(map, cstas.x, cstas.y, m=(:dtriangle, :green), label="")
annotate!(map, cstas.x .+ 30, cstas.y, text.(cstas.sta, 8))

# Plot of azimuths and distances
X, Y = [], []
for (θ, h2) in zip(Beamforming._azimuth_distance_ccbf(x, y)...)
    h = h2/2
    push!(X, [0, h*cosd(θ), NaN, 0, -h*cosd(θ)])
    push!(Y, [0, h*sind(θ), NaN, 0, -h*sind(θ)])
end
azmap = plot(X, Y, primary=false, aspect_ratio=:equal, framestyle=:box, grid=false, arrow=:arrow,
    xlabel="East distance (km)", ylabel="North distance (km)")

# Colorbar
cbar = heatmap([0, 2], -10:0.01:0, [collect(-10:0.01:0) collect(-10:0.01:0)], c=cpt, xticks=nothing, cbar=false,
    aspect_ratio=8, xlim=(0,2), ylabel="Beam power (dB)")

# All plots assembled together
allplots = plot(map, azmap, cbar, arf_5, arf_3_7, arf_3_7²,
    layout=(2,3), size=(900,500), label="", cbar=false, titlefontsize=12)

# Save plot if run as a script, otherwise return (and show)
if !isinteractive() && abspath(@__FILE__) == abspath(PROGRAM_FILE)
    savefig(allplots, joinpath(@__DIR__, "Ruigrok_et_al_2017_Fig2.pdf"))
else
    allplots
end

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming


text:
"""
System to be coupled with Sesam3 that provides constant inputs, i.e. drivers.
Need to specify parameters
- `i_L0`: litter input rate
- `β_Ni0`, `β_Pi0`: C:N and C:P ratio of litter input
- `i_IN0`, `i_IP0`: external inputs of inorganic N and P (N deposition, weathering)
optional:
- `u_PlantNmax0`, `u_PlantPmax0`: maximum plant uptake of inorganic N and P
- `k_PlantN0`, `k_PlantP0`: plant uptake rate per inorganic pool of N and P

Defaults of plant uptake are set so that uptake rate is very high (1e3) so that
the maximum is taken up, i.e. wins competition with microbes.
The maximum plan uptake is set to match the input by litter.
"""
function plant_const(;name)
    # @parameters t [unit = uT]
    # D = Differential(t)
    # @parameters i_L0 [unit = uQ/uT] i_IN0 [unit = uQ/uT] β_Ni0
    # @parameters u_PlantNmax0 [unit = uQ/uT] k_PlantN0 [unit = uT^-1]
    # @variables (begin
    #     i_L(t),[unit = uQ/uT], β_Ni(t), i_IN(t),[unit = uQ/uT],
    #     u_PlantNmax(t),[unit = uQ/uT], k_PlantN(t),[unit = uT^-1]
    # end)
    @parameters t
    D = Differential(t)
    @parameters i_L0  i_IN0
    @parameters u_PlantNmax0  k_PlantN0 β_Ni0 i_IP0  β_Pi0 u_PlantPmax0  k_PlantP0 s_EP0
    @variables (begin
        i_L(t),
        β_Ni(t), i_IN(t),
        u_PlantNmax(t), k_PlantN(t),
        β_Pi(t), i_IP(t),
        u_PlantPmax(t), k_PlantP(t),
        s_EP(t)
    end)
    eqs = [
        i_L ~ i_L0,
        β_Ni ~ β_Ni0,
        i_IN ~ i_IN0,
        u_PlantNmax ~ u_PlantNmax0,
        k_PlantN ~ k_PlantN0,
        β_Pi ~ β_Pi0,
        i_IP ~ i_IP0,
        u_PlantPmax ~ u_PlantPmax0,
        k_PlantP ~ k_PlantP0,
        s_EP ~ s_EP0,
    ]
    defaults=Pair{Num, Any}[
        # rate not constrained, to make u_PlantNmax0 effective
        k_PlantN0 => 1e3,
        # take as much N as provided with litter
        u_PlantNmax0 => i_L0/β_Ni0,
        # rate not constrained
        k_PlantP0 => 1e3,
        # take as much P as provided with litter
        u_PlantPmax0 => i_L0/β_Pi0,
        # defaults in order to work with the CN-only version
        β_Pi0 => 50.0, i_IP0 => 0.0, s_EP0 => 0.0,
    ]
    ODESystem(eqs; name, defaults)
end


text:
"""
    MultiScaleSkipDeepEquilibriumNetwork(main_layers::Tuple, mapping_layers::Matrix, shortcut_layers::Tuple,
                                         solver; post_fuse_layers::Union{Tuple,Nothing}=nothing, p=nothing,
                                         sensealg=SteadyStateAdjoint(0.1f0, 0.1f0, 10), kwargs...)
    MultiScaleSkipDeepEquilibriumNetwork(main_layers::Tuple, mapping_layers::Matrix, solver;
                                         post_fuse_layers::Union{Tuple,Nothing}=nothing, p=nothing,
                                         sensealg=SteadyStateAdjoint(0.1f0, 0.1f0, 10), kwargs...)

Multiscale Deep Equilibrium Network as proposed in [baimultiscale2020](@cite)

## Arguments

* `main_layers`: Tuple of Explicit Neural Networks. The first network needs to take 2 inputs, the other ones only take 1 input
* `mapping_layers`: Matrix of Explicit Neural Networks. The ``(i, j)^{th}`` network takes the output of ``i^{th}`` `main_layer`
                    and passes it to the ``j^{th}`` `main_layer`
* `shortcut_layers`: Shortcuts for the network (If not given, then we create SkipDEQV2)
* `solver`: Solver for the optimization problem (See: [`ContinuousDEQSolver`](@ref) & [`DiscreteDEQSolver`](@ref))
* `post_fuse_layers`: Tuple of Explicit Neural Networks. Applied after the `mapping_layers` (Default: `nothing`)
* `p`: Optional parameters for the `model`
* `sensealg`: See [`SteadyStateAdjoint`](@ref)
* `kwargs`: Additional Parameters that are directly passed to `solve`

## Example

```julia
# MSkipDEQ
model = MultiScaleSkipDeepEquilibriumNetwork(
    (
        Parallel(+, Dense(4, 4, tanh_fast), Dense(4, 4, tanh_fast)),
        Dense(3, 3, tanh_fast), Dense(2, 2, tanh_fast),
        Dense(1, 1, tanh_fast)
    ),
    [
        NoOpLayer() Dense(4, 3, tanh_fast) Dense(4, 2, tanh_fast) Dense(4, 1, tanh_fast);
        Dense(3, 4, tanh_fast) NoOpLayer() Dense(3, 2, tanh_fast) Dense(3, 1, tanh_fast);
        Dense(2, 4, tanh_fast) Dense(2, 3, tanh_fast) NoOpLayer() Dense(2, 1, tanh_fast);
        Dense(1, 4, tanh_fast) Dense(1, 3, tanh_fast) Dense(1, 2, tanh_fast) NoOpLayer()
    ],
    (
        Dense(4, 4, tanh_fast),
        Dense(4, 3, tanh_fast),
        Dense(4, 2, tanh_fast),
        Dense(4, 1, tanh_fast)
    ),
    ContinuousDEQSolver(VCABM3(); abstol=0.01f0, reltol=0.01f0),
)

model(rand(Float32, 4, 1))

# MSkipDEQV2
model = MultiScaleSkipDeepEquilibriumNetwork(
    (
        Parallel(+, Dense(4, 4, tanh_fast), Dense(4, 4, tanh_fast)),
        Dense(3, 3, tanh_fast), Dense(2, 2, tanh_fast),
        Dense(1, 1, tanh_fast)
    ),
    [
        NoOpLayer() Dense(4, 3, tanh_fast) Dense(4, 2, tanh_fast) Dense(4, 1, tanh_fast);
        Dense(3, 4, tanh_fast) NoOpLayer() Dense(3, 2, tanh_fast) Dense(3, 1, tanh_fast);
        Dense(2, 4, tanh_fast) Dense(2, 3, tanh_fast) NoOpLayer() Dense(2, 1, tanh_fast);
        Dense(1, 4, tanh_fast) Dense(1, 3, tanh_fast) Dense(1, 2, tanh_fast) NoOpLayer()
    ],
    ContinuousDEQSolver(VCABM3(); abstol=0.01f0, reltol=0.01f0),
)

model(rand(Float32, 4, 1))
```

See also: [`DeepEquilibriumNetwork`](@ref), [`SkipDeepEquilibriumNetwork`](@ref), [`MultiScaleDeepEquilibriumNetwork`](@ref)
"""
struct MultiScaleSkipDeepEquilibriumNetwork{M3<:Union{Nothing,Parallel},N,M1<:Parallel,M2<:Union{Chain,FChain},RE1,
                                            RE2,RE3,P,A,K,S} <: AbstractDeepEquilibriumNetwork
    main_layers::M1
    mapping_layers::M2
    shortcut_layers::M3
    main_layers_re::RE1
    mapping_layers_re::RE2
    shortcut_layers_re::RE3
    p::P
    ordered_split_idxs::NTuple{N,Int}
    solver::A
    kwargs::K
    sensealg::S
    stats::DEQTrainingStats

    function MultiScaleSkipDeepEquilibriumNetwork(main_layers::Parallel,
                                                  mapping_layers::Union{Chain,FChain},
                                                  shortcut_layers::Union{Nothing,Parallel}, re1, re2, re3, p,
                                                  ordered_split_idxs, solver::A, kwargs::K, sensealg::S,
                                                  stats) where {A,K,S}
        @assert length(mapping_layers) == 2
        @assert mapping_layers[1] isa MultiParallelNet

        p_main_layers, re_main_layers = destructure_parameters(main_layers)
        p_mapping_layers, re_mapping_layers = destructure_parameters(mapping_layers)
        p_shortcut_layers, re_shortcut_layers = shortcut_layers === nothing ? ([], nothing) :
                                                destructure_parameters(shortcut_layers)

        ordered_split_idxs = tuple(cumsum([0, length(p_main_layers), length(p_mapping_layers),
                                           length(p_shortcut_layers)])...)

        p = p === nothing ? vcat(p_main_layers, p_mapping_layers, p_shortcut_layers) : convert(typeof(p_main_layers), p)

        return new{typeof(shortcut_layers),length(ordered_split_idxs),
                   typeof.((main_layers, mapping_layers, re_main_layers, re_mapping_layers, re_shortcut_layers, p))...,
                   A,K,S}(main_layers, mapping_layers, shortcut_layers, re_main_layers,
                          re_mapping_layers, re_shortcut_layers, p, ordered_split_idxs, solver, kwargs, sensealg, stats)
    end
end

Flux.@functor MultiScaleSkipDeepEquilibriumNetwork

function MultiScaleSkipDeepEquilibriumNetwork(main_layers::Tuple, mapping_layers::Matrix, shortcut_layers::Tuple,
                                              solver; post_fuse_layers::Union{Tuple,Nothing}=nothing, p=nothing,
                                              sensealg=SteadyStateAdjoint(0.1f0, 0.1f0, 10), kwargs...)
    mapping_layers = if post_fuse_layers === nothing
        @assert size(mapping_layers, 1) == size(mapping_layers, 2) == length(main_layers) == length(shortcut_layers)
        Chain(MultiParallelNet(Parallel.(+, map(x -> tuple(x...), eachcol(mapping_layers)))...), NoOpLayer())
    else
        @assert size(mapping_layers, 1) ==
                size(mapping_layers, 2) ==
                length(main_layers) ==
                length(post_fuse_layers) ==
                length(shortcut_layers)
        Chain(MultiParallelNet(Parallel.(+, map(x -> tuple(x...), eachcol(mapping_layers)))...),
              Parallel(flatten_merge, post_fuse_layers...))
    end

    main_layers = Parallel(flatten_merge, main_layers...)
    shortcut_layers = Parallel(flatten_merge, shortcut_layers...)

    return MultiScaleSkipDeepEquilibriumNetwork(main_layers, mapping_layers, shortcut_layers,
                                                nothing, nothing, nothing, p, nothing, solver, kwargs, sensealg,
                                                DEQTrainingStats(0))
end

function MultiScaleSkipDeepEquilibriumNetwork(main_layers::Tuple, mapping_layers::Matrix, solver;
                                              post_fuse_layers::Union{Tuple,Nothing}=nothing, p=nothing,
                                              sensealg=SteadyStateAdjoint(0.1f0, 0.1f0, 10), kwargs...)
    mapping_layers = if post_fuse_layers === nothing
        @assert size(mapping_layers, 1) == size(mapping_layers, 2) == length(main_layers)
        Chain(MultiParallelNet(Parallel.(+, map(x -> tuple(x...), eachcol(mapping_layers)))...), NoOpLayer())
    else
        @assert size(mapping_layers, 1) == size(mapping_layers, 2) == length(main_layers) == length(post_fuse_layers)
        Chain(MultiParallelNet(Parallel.(+, map(x -> tuple(x...), eachcol(mapping_layers)))...),
              Parallel(flatten_merge, post_fuse_layers...))
    end

    main_layers = Parallel(flatten_merge, main_layers...)

    return MultiScaleSkipDeepEquilibriumNetwork(main_layers, mapping_layers, nothing, nothing,
                                                nothing, nothing, p, nothing, solver, kwargs, sensealg,
                                                DEQTrainingStats(0))
end

function (mdeq::MultiScaleSkipDeepEquilibriumNetwork)(x::AbstractArray{T}) where {T}
    current_nfe = mdeq.stats.nfe

    p1, p2, p3 = split_array_by_indices(mdeq.p, mdeq.ordered_split_idxs)
    initial_conditions = mdeq.shortcut_layers_re(p3)(x)
    u_sizes = size.(initial_conditions)
    u_split_idxs = vcat(0, cumsum(length.(initial_conditions) .÷ size(x, ndims(x)))...)
    u0 = Zygote.@ignore vcat(Flux.flatten.(initial_conditions)...)

    N = length(u_sizes)
    update_is_variational_hidden_dropout_mask_reset_allowed(false)

    function dudt_(u, _p)
        mdeq.stats.nfe += 1

        uₛ = split_array_by_indices(u, u_split_idxs)
        p1, p2, _ = split_array_by_indices(_p, mdeq.ordered_split_idxs)

        u_reshaped = ntuple(i -> reshape(uₛ[i], u_sizes[i]), N)

        main_layers_output = mdeq.main_layers_re(p1)((u_reshaped[1], x), u_reshaped[2:end]...)

        return mdeq.mapping_layers_re(p2)(main_layers_output)
    end

    dudt(u, _p, t) = vcat(Flux.flatten.(dudt_(u, _p))...) .- u

    ssprob = SteadyStateProblem(dudt, u0, mdeq.p)
    res = solve(ssprob, mdeq.solver; u0=u0, sensealg=mdeq.sensealg, mdeq.kwargs...).u
    x_ = dudt_(res, mdeq.p)
    residual = Zygote.@ignore Tuple(map((iu) -> reshape(iu[2], u_sizes[iu[1]]),
                                    enumerate(split_array_by_indices(dudt(res, mdeq.p, nothing), u_split_idxs))))

    update_is_variational_hidden_dropout_mask_reset_allowed(true)

    return x_, DeepEquilibriumSolution(x_, initial_conditions, residual, T(0), mdeq.stats.nfe - current_nfe)
end

function (mdeq::MultiScaleSkipDeepEquilibriumNetwork{Nothing})(x::AbstractArray{T}) where {T}
    current_nfe = mdeq.stats.nfe

    p1, p2 = split_array_by_indices(mdeq.p, mdeq.ordered_split_idxs)

    _initial_conditions = Zygote.@ignore [l(x) for l in map(l -> l.layers[1], mdeq.mapping_layers[1].layers)]
    _initial_conditions = mdeq.mapping_layers_re(p2)((x, zero.(_initial_conditions[2:end])...))
    initial_conditions = mdeq.main_layers_re(p1)((zero(_initial_conditions[1]), _initial_conditions[1]),
                                                 _initial_conditions[2:end]...)
    u_sizes = size.(initial_conditions)
    u_split_idxs = vcat(0, cumsum(length.(initial_conditions) .÷ size(x, ndims(x)))...)
    u0 = vcat(Flux.flatten.(initial_conditions)...)

    N = length(u_sizes)
    update_is_variational_hidden_dropout_mask_reset_allowed(false)

    function dudt_(u, _p)
        mdeq.stats.nfe += 1

        uₛ = split_array_by_indices(u, u_split_idxs)
        p1, p2, _ = split_array_by_indices(_p, mdeq.ordered_split_idxs)

        u_reshaped = ntuple(i -> reshape(uₛ[i], u_sizes[i]), N)

        main_layers_output = mdeq.main_layers_re(p1)((u_reshaped[1], x), u_reshaped[2:end]...)

        return mdeq.mapping_layers_re(p2)(main_layers_output)
    end

    dudt(u, _p, t) = vcat(Flux.flatten.(dudt_(u, _p))...) .- u

    ssprob = SteadyStateProblem(dudt, u0, mdeq.p)
    res = solve(ssprob, mdeq.solver; u0=u0, sensealg=mdeq.sensealg, mdeq.kwargs...).u
    x_ = dudt_(res, mdeq.p)
    residual = Zygote.@ignore Tuple(map((iu) -> reshape(iu[2], u_sizes[iu[1]]),
                                    enumerate(split_array_by_indices(dudt(res, mdeq.p, nothing), u_split_idxs))))

    update_is_variational_hidden_dropout_mask_reset_allowed(true)

    return x_, DeepEquilibriumSolution(x_, initial_conditions, residual, T(0), mdeq.stats.nfe - current_nfe)
end


text:
module MuxPkgServer

import HTTP
using URIs: URI
using Base: UUID, SHA1

const uuid_re = raw"[0-9a-z]{8}-[0-9a-z]{4}-[0-9a-z]{4}-[0-9a-z]{4}-[0-9a-z]{12}(?-i)"
const hash_re = raw"[0-9a-f]{40}"
const registry_re = Regex("^/registry/($(uuid_re))/($(hash_re))\$")
const resource_re = Regex("""
                          ^/registry/$(uuid_re)/($(hash_re))\$ |
                          ^/package/$(uuid_re)/($(hash_re))\$ |
                          ^/artifact/($(hash_re))\$
                          """, "x")

function serve_registry(http::HTTP.Stream)
    # TODO: Maybe catch registries here to cache which upstream server it came from...
end

function serve_registries(http::HTTP.Stream, upstreams::Vector{URI})
    regs = Dict{UUID,SHA1}()
    for upstream in upstreams
        url = URI(upstream; path="/registries")
        req = HTTP.request("GET", url; status_exception=false)
        req.status == 200 || continue
        for l in eachline(IOBuffer(req.body))
            m = match(registry_re, l)
            m === nothing && continue
            get!(regs, UUID(m[1]), SHA1(m[2]))
        end
    end
    str = sprint() do io
        for (uuid, sha1) in regs
            println(io, "/registry/$(uuid)/$(sha1)")
        end
    end
    serve_data(http, str; content_type = "text/plain")
    return nothing
end

function find_upstream(http::HTTP.Stream, upstreams::Vector{URI})
    target = http.message.target
    # TODO: Could be done asynchronous
    for upstream in upstreams
        url = URI(upstream; path=target)
        req = HTTP.request("HEAD", url; status_exception=false)
        @info "HEAD request to" url req.status
        if req.status == 200
            return proxy_pass(http, upstream)
            # This could just 301 to take load from this server, but that can possibly
            # make client configuration more difficult if any of the upstream servers
            # are authenticated for example.
            # HTTP.setstatus(http, 302)
            # HTTP.setheader(http, "Location" => url)
            # HTTP.setheader(http, "Content-Length" => "0")
            # HTTP.startwrite(http)
        end
    end
    return fourxx(http, 404)
end

"""
    MuxPkgServer.start(;
        host="0.0.0.0",
        port=8080,
        upstreams=["https://pkg.julialang.org"])
"""
function start(; host="0.0.0.0", port=8080,
                 upstreams = ["https://pkg.julialang.org"],
              )
    upstreams = URI.(upstreams)

    @info "Starting MuxPkgServer" host port upstreams
    # Start the MuxPkgServer instance
    HTTP.listen(host, port) do http::HTTP.Stream
        # Target resource
        resource = http.message.target

        if occursin(resource_re, resource)
            # /registry/${UUID}/${SHA1}, /package/${UUID}/${SHA1} and /artifact/${SHA1}
            return find_upstream(http, upstreams)
        # elseif occursin(registry_re, resource)
        #     return serve_registry(http, gid)
        elseif resource == "/registries"
            return serve_registries(http, upstreams)
        end

        # Resource URL does not exist
        return fourxx(http, 404)
    end # HTTP.listen
end

##############################
## HTTP.jl helper functions ##
##############################

# Proxy a request to upstream and write back the response
function proxy_pass(http::HTTP.Stream, upstream::URI; target::String=http.message.target)
    # Upstream URI
    proxy_uri = URI(upstream; path=target)

    # Pass along all headers except Authorization
    filter!(x -> x.first != "Authorization", http.message.headers)
    filter!(x -> x.first != "Host", http.message.headers) # TODO: see what nginx does

    # Open connection to upstream (TODO: make sure downstream used GET...)
    HTTP.open("GET", proxy_uri, http.message.headers; status_exception=false) do http2
        # Read status code and headers
        req = HTTP.startread(http2)
        @info "In open remote" req.status HTTP.isredirect(http2)
        if HTTP.isredirect(req)
            # TODO: Should this function even be called for redirects?
            #       Have asked HTTP.jl developers.
            return nothing
        end
        # Write response back to http
        @info "writing response" req.status
        HTTP.setstatus(http, req.status)
        foreach(h -> HTTP.setheader(http, h), req.headers)
        HTTP.startwrite(http)
        if http.message.method == "GET"
            HTTP.write(http, http2)
        end
    end
    return nothing
end

# Throw 4XX status codes
function fourxx(http::HTTP.Stream, s::Int, msg::Union{String,Nothing}=nothing)
    HTTP.setstatus(http, s)
    if msg === nothing
        HTTP.setheader(http, "Content-Length" => "0")
        HTTP.startwrite(http)
        HTTP.close(http.stream)
    else
        HTTP.setheader(http, "Content-Type" => "text/plain")
        HTTP.setheader(http, "Content-Length" => string(sizeof(msg)))
        HTTP.startwrite(http)
        HTTP.write(http, msg)
        HTTP.close(http.stream)
    end
    return nothing
end

function serve_data(http::HTTP.Stream, data::String; content_type::String)
    HTTP.setstatus(http, 200)
    HTTP.setheader(http, "Content-Type" => content_type)
    HTTP.setheader(http, "Content-Length" => string(sizeof(data)))
    HTTP.startwrite(http)
    if http.message.method == "GET"
        HTTP.write(http, data)
    end
    return nothing
end

end # module


text:
# Julia/LLVM IR generation and transformation passes

function module_setup(mod::LLVM.Module)
    triple!(mod, Int === Int64 ? "nvptx64-nvidia-cuda" : "nvptx-nvidia-cuda")

    # add debug info metadata
    push!(metadata(mod), "llvm.module.flags",
         MDNode([ConstantInt(Int32(1), JuliaContext()),    # llvm::Module::Error
                 MDString("Debug Info Version"),
                 ConstantInt(DEBUG_METADATA_VERSION(), JuliaContext())]))
end

# make function names safe for PTX
safe_fn(fn::String) = replace(fn, r"[^aA-zZ0-9_]"=>"_")
safe_fn(f::Core.Function) = safe_fn(String(typeof(f).name.mt.name))
safe_fn(f::LLVM.Function) = safe_fn(LLVM.name(f))

# NOTE: we remove `throw_...` functions in the ThrowRemoval pass, but that relies on
#       function names, which is fragile. Remove them here as well, to be safe.
const exception_arguments = Vector{LLVM.Value}()
function raise_exception(insblock::BasicBlock, ex::Value)
    fun = LLVM.parent(insblock)
    mod = LLVM.parent(fun)
    ctx = context(mod)

    trap = if haskey(functions(mod), "llvm.trap")
        functions(mod)["llvm.trap"]
    else
        LLVM.Function(mod, "llvm.trap", LLVM.FunctionType(LLVM.VoidType(ctx)))
    end

    let builder = Builder(ctx)
        position!(builder, insblock)

        cuprintf!(builder, "ERROR: an unknown exception occurred$cuprintf_endline during kernel execution")
        call!(builder, trap)

        dispose(builder)
    end

    # mark arguments that passed to `throw` for removal, as they are often GC-allocated.
    push!(exception_arguments, ex)
end

# generate a pseudo-backtrace from a stack of methods being emitted
function backtrace(ctx::CompilerContext, method_stack::Vector{Core.MethodInstance})
    bt = StackTraces.StackFrame[]
    for method_instance in method_stack
        method = method_instance.def
        frame = StackTraces.StackFrame(method.name, method.file, method.line)
        pushfirst!(bt, frame)
    end
    bt
end

# NOTE: we use an exception to be able to display a stack trace using the logging framework
struct MethodSubstitutionWarning <: Exception
    original::Method
    substitute::Method
end
Base.showerror(io::IO, err::MethodSubstitutionWarning) =
    print(io, "You called $(err.original), maybe you intended to call $(err.substitute) instead?")

function irgen(ctx::CompilerContext)
    # get the method instance
    isa(ctx.f, Core.Builtin) && throw(KernelError(ctx, "function is not a generic function"))
    world = typemax(UInt)
    meth = which(ctx.f, ctx.tt)
    sig = Base.signature_type(ctx.f, ctx.tt)::Type
    (ti, env) = ccall(:jl_type_intersection_with_env, Any,
                      (Any, Any), sig, meth.sig)::Core.SimpleVector
    meth = Base.func_for_method_checked(meth, ti)
    linfo = ccall(:jl_specializations_get_linfo, Ref{Core.MethodInstance},
                  (Any, Any, Any, UInt), meth, ti, env, world)

    # set-up the compiler interface
    function hook_module_setup(ref::Ptr{Cvoid})
        ref = convert(LLVM.API.LLVMModuleRef, ref)
        module_setup(LLVM.Module(ref))
    end
    function hook_raise_exception(insblock::Ptr{Cvoid}, ex::Ptr{Cvoid})
        insblock = convert(LLVM.API.LLVMValueRef, insblock)
        ex = convert(LLVM.API.LLVMValueRef, ex)
        raise_exception(BasicBlock(insblock), Value(ex))
    end
    dependencies = Vector{LLVM.Module}()
    function hook_module_activation(ref::Ptr{Cvoid})
        ref = convert(LLVM.API.LLVMModuleRef, ref)
        push!(dependencies, LLVM.Module(ref))
    end
    method_stack = Vector{Core.MethodInstance}()
    function hook_emit_function(method_instance, code, world)
        push!(method_stack, method_instance)

        # check for recursion
        if method_instance in method_stack[1:end-1]
            throw(KernelError(ctx, "recursion is currently not supported";
                              bt=backtrace(ctx, method_stack)))
        end

        # check for Base methods that exist in CUDAnative too
        # FIXME: this might be too coarse
        method = method_instance.def
        if Base.moduleroot(method.module) == Base &&
           isdefined(CUDAnative, method_instance.def.name)
            substitute_function = getfield(CUDAnative, method.name)
            tt = Tuple{method_instance.specTypes.parameters[2:end]...}
            if hasmethod(substitute_function, tt)
                method′ = which(substitute_function, tt)
                if Base.moduleroot(method′.module) == CUDAnative
                    @warn "calls to Base intrinsics might be GPU incompatible" exception=(MethodSubstitutionWarning(method, method′), backtrace(ctx, method_stack))
                end
            end
        end
    end
    function hook_emitted_function(method, code, world)
        @compiler_assert last(method_stack) == method ctx
        pop!(method_stack)
    end
    params = Base.CodegenParams(cached             = false,
                                track_allocations  = false,
                                code_coverage      = false,
                                static_alloc       = false,
                                prefer_specsig     = true,
                                module_setup       = hook_module_setup,
                                module_activation  = hook_module_activation,
                                raise_exception    = hook_raise_exception,
                                emit_function      = hook_emit_function,
                                emitted_function   = hook_emitted_function)

    # get the code
    mod = let
        ref = ccall(:jl_get_llvmf_defn, LLVM.API.LLVMValueRef,
                    (Any, UInt, Bool, Bool, Base.CodegenParams),
                    linfo, world, #=wrapper=#false, #=optimize=#false, params)
        if ref == C_NULL
            throw(InternalCompilerError(ctx, "the Julia compiler could not generate LLVM IR"))
        end

        llvmf = LLVM.Function(ref)
        LLVM.parent(llvmf)
    end

    # the main module should contain a single jfptr_ function definition,
    # e.g. jfptr_kernel_vadd_62977
    definitions = LLVM.Function[]
    for llvmf in functions(mod)
        if !isdeclaration(llvmf)
            push!(definitions, llvmf)
        end
    end
    wrapper = nothing
    for llvmf in definitions
        if startswith(LLVM.name(llvmf), "jfptr_")
            @compiler_assert wrapper == nothing ctx first=wrapper second=llvmf
            wrapper = llvmf
        end
    end
    @compiler_assert wrapper != nothing ctx

    # the jfptr wrapper function should point us to the actual entry-point,
    # e.g. julia_kernel_vadd_62984
    # FIXME: Julia's globalUnique starting with `-` is probably a bug.
    entry_tag = let
        m = match(r"^jfptr_(.+)_[-\d]+$", LLVM.name(wrapper))
        @compiler_assert m != nothing ctx name=LLVM.name(wrapper)
        m.captures[1]
    end
    unsafe_delete!(mod, wrapper)
    entry = let
        re = Regex("^julia_$(entry_tag)_[-\\d]+\$")
        entrypoints = LLVM.Function[]
        for llvmf in definitions
            if llvmf != wrapper
                llvmfn = LLVM.name(llvmf)
                if occursin(re, llvmfn)
                    push!(entrypoints, llvmf)
                end
            end
        end
        @compiler_assert length(entrypoints) == 1 ctx functions=Tuple(LLVM.name.(definitions)) tag=entry_tag entrypoints=Tuple(LLVM.name.(entrypoints))
        entrypoints[1]
    end

    # link in dependent modules
    for dep in dependencies
        link!(mod, dep)
    end

    # clean up incompatibilities
    for llvmf in functions(mod)
        llvmfn = LLVM.name(llvmf)

        # only occurs in debug builds
        delete!(function_attributes(llvmf), EnumAttribute("sspstrong", 0, JuliaContext()))

        # dependent modules might have brought in other jfptr wrappers, delete them
        if startswith(LLVM.name(llvmf), "jfptr_") && isempty(uses(llvmf))
            unsafe_delete!(mod, llvmf)
            continue
        end

        # llvmcall functions aren't to be called, so mark them internal (cleans up the IR)
        if startswith(llvmfn, "jl_llvmcall")
            linkage!(llvmf, LLVM.API.LLVMInternalLinkage)
            continue
        end

        # rename functions
        if !isdeclaration(llvmf)
            # Julia disambiguates local functions by prefixing with `#\d#`.
            # since we don't use a global function namespace, get rid of those tags.
            if occursin(r"^julia_#\d+#", llvmfn)
                llvmfn′ = replace(llvmfn, r"#\d+#"=>"")
                if !haskey(functions(mod), llvmfn′)
                    LLVM.name!(llvmf, llvmfn′)
                    llvmfn = llvmfn′
                end
            end

            # anonymous functions are just named `#\d`, make that somewhat more readable
            m = match(r"_#(\d+)_", llvmfn)
            if m !== nothing
                llvmfn′ = replace(llvmfn, m.match=>"_anonymous$(m.captures[1])_")
                LLVM.name!(llvmf, llvmfn′)
                llvmfn = llvmfn′
            end

            # finally, make function names safe for ptxas
            # (LLVM should to do this, but fails, see eg. D17738 and D19126)
            llvmfn′ = safe_fn(llvmfn)
            if llvmfn != llvmfn′
                LLVM.name!(llvmf, llvmfn′)
                llvmfn = llvmfn′
            end
        end
    end

    # HACK: remove unused arguments to exception functions
    for val in exception_arguments
        if isa(val, LLVM.Instruction) && isempty(uses(val))
            bb = LLVM.parent(val)
            unsafe_delete!(bb, val)
        end
    end
    empty!(exception_arguments)

    # rename the entry point
    llvmfn = replace(LLVM.name(entry), r"_\d+$"=>"")
    ## append a global unique counter
    global globalUnique
    globalUnique += 1
    llvmfn *= "_$globalUnique"
    LLVM.name!(entry, llvmfn)

    # minimal optimization to get rid of useless generated code (llvmcall, kernel wrapper)
    ModulePassManager() do pm
        add!(pm, ModulePass("ThrowRemoval", remove_throw!))
        add!(pm, FunctionPass("ControlFlowFixup", fixup_controlflow!))
        always_inliner!(pm)
        verifier!(pm)
        run!(pm, mod)
    end

    return mod, entry
end

# HACK: this pass replaces `julia_throw_*` void functions with a simple print & trap
#
# this is necessary even though we already have a `raise_exception` hook that just prints,
# as many `throw_...` functions are now `@noinline` which means their arguments survive and
# may cause GC allocations.
#
# TODO: replace with a Cassette-style overdub of the `throw` builtin
function remove_throw!(mod::LLVM.Module)
    ctx = LLVM.context(mod)

    trap = if haskey(functions(mod), "llvm.trap")
        functions(mod)["llvm.trap"]
    else
        LLVM.Function(mod, "llvm.trap", LLVM.FunctionType(LLVM.VoidType(ctx)))
    end

    changed = false

    # NOTE: module pass, since we delete functions
    for f in collect(functions(mod))
        fn = LLVM.name(f)
        ft = eltype(llvmtype(f))

        for re in [# common exceptions as defined in the runtime
                   r"jl_(bounds_error)_.+",
                   # user-code throw functions
                   # FIXME: this is coarse
                   r"julia_throw_(.+)_\d+"]
            m = match(re, fn)
            if m != nothing && return_type(ft) == LLVM.VoidType(ctx)
                ex = m.captures[1]

                # create a function that prints the exception name, and exits
                fn′ = "ptx_throw_$ex"
                f′ = if haskey(functions(mod), fn′)
                    functions(mod)[fn′]
                else
                    ft′ = LLVM.FunctionType(LLVM.VoidType(ctx))
                    f′ = LLVM.Function(mod, fn′, ft′)
                    let builder = Builder(ctx)
                        entry = BasicBlock(f′, "entry", ctx)
                        position!(builder, entry)

                        desc = replace(ex, '_'=>' ')
                        cuprintf!(builder, "ERROR: a $ex exception occurred$cuprintf_endline during kernel execution")
                        call!(builder, trap)
                        ret!(builder)

                        dispose(builder)
                    end
                    f′
                end

                # replace uses of the original function
                for use in uses(f)
                    call = user(use)
                    @assert isa(call, LLVM.CallInst)
                    let builder = Builder(ctx)
                        position!(builder, call)
                        call!(builder, f′)
                        dispose(builder)
                    end
                    unsafe_delete!(LLVM.parent(call), call)
                end

                # remove the original function or declaration
                @assert isempty(uses(f))
                unsafe_delete!(mod, f)

                changed = true
            end
        end
    end

    return changed
end

# HACK: this pass removes control flow that confuses `ptxas`
#
# basically, we only want a single exit from a function. exiting (ret, trap, exit) from
# divergent branches causes ptxas to emit invalid code.
#
# TODO: do this with structured CFG with LLVM?
function fixup_controlflow!(f::LLVM.Function)
    ctx = LLVM.context(f)

    exit_ft = LLVM.FunctionType(LLVM.VoidType(ctx))
    exit = InlineAsm(exit_ft, "exit;", "", true)

    changed = false

    # remove `noreturn` attributes
    #
    # when calling a `noreturn` function, LLVM places an `unreachable` after the call.
    # this leads to an early `ret` from the function.
    attrs = function_attributes(f)
    delete!(attrs, EnumAttribute("noreturn", 0, ctx))

    for bb in blocks(f)
        # remove calls to `trap`, for obvious reasons.
        for inst in instructions(bb)
            if isa(inst, LLVM.CallInst) && LLVM.name(called_value(inst)) == "llvm.trap"
                # replace with call to `exit`
                # FIXME: this still confuses `ptxas`, `brkpt` seems to work but that's fragile
                #let builder = Builder(ctx)
                #    position!(builder, inst)
                #    call!(builder, exit)
                #    dispose(builder)
                #end

                unsafe_delete!(bb, inst)
                changed = true
            end
        end

        # remove `unreachable `terminators.
        #
        # at this point, the optimizer hasn't run, so these hints are placed there by Julia.
        # this happens e.g. with exceptions, and cannot be controlled by a compiler hook.
        unreachable = terminator(bb)
        if isa(unreachable, LLVM.UnreachableInst)
            unsafe_delete!(bb, unreachable)
            changed = true

            try
                terminator(bb)
                # the basic-block is still terminated properly, nothing to do
                # (this can happen with `ret; unreachable`)
                # TODO: `unreachable; unreachable`
            catch ex
                isa(ex, UndefRefError) || rethrow(ex)

                let builder = Builder(ctx)
                    position!(builder, bb)

                    # find the predecessors to this block
                    predecessors = LLVM.BasicBlock[]
                    for bb′ in blocks(f), inst in instructions(bb′)
                        if isa(inst, LLVM.BrInst)
                            if bb in successors(inst)
                                push!(predecessors, bb′)
                            end
                        end
                    end

                    # find the fall through successors
                    if length(predecessors) == 1
                        predecessor = first(predecessors)
                        br = terminator(predecessor)

                        # find the other successors
                        targets = successors(br)
                        if length(targets) == 2
                            for target in targets
                                if target != bb
                                    br!(builder, target)
                                    break
                                end
                            end
                        else
                            @warn "unreachable control flow with $(length(targets))-way branching predecessor"
                            unreachable!(builder)
                        end
                    elseif length(predecessors) > 1
                        @warn "unreachable control flow with multiple predecessors"
                        unreachable!(builder)
                    else
                        # this block has no predecessors, so we can't fall through
                        #
                        # a block without predecessors is either never called,
                        # or the only block in a function.
                        ft = eltype(llvmtype(f))
                        if return_type(ft) == LLVM.VoidType(ctx)
                            # even though returning can lead to invalid control flow,
                            # it mostly happens with functions that just throw,
                            # and leaving the unreachable there would make the optimizer
                            # place another after the call.
                            ret!(builder)
                        else
                            unreachable!(builder)
                        end
                    end

                    dispose(builder)
                end
            end
        end
    end

    return changed
end


text:
# structs

struct Disp
    F ::Float64
    pairwise_F ::Array{Float64}
    medians ::Union{Tuple,Vector}
    residuals ::Vector{Matrix{Float64}}
    means ::Vector{Float64}
    group  ::Vector
    levels  ::Vector
    group_residuals ::Vector
end
#Helper functions

# calculate distances from medians
function get_residuals(x,c)
    d = x .-c
    sum(d .^2,dims =2)
end

# doube center matrix
function double_center(D)
    A = D .- mean(D,dims = 1)
    return A .- mean(A,dims = 2)
end

# calculate medians
# Get spatial median as a row vector
geo_median(X) = DirectionalStatistics.geometric_median([row for row in eachrow(X)])'
#apply `geo_median` function to get positive_inds and negative_inds median for each group
function geo_median(pco_space, group, positive_inds)

    spMedPos = [geo_median(pco_space[group .== g,positive_inds]) for g in unique(group)]
    spMedNeg = [geo_median(pco_space[group .== g, .!positive_inds]) for g in unique(group)]

    return (spMedPos,spMedNeg )
end
function geo_median(pco_space, group)

    return [geo_median(pco_space[group .== g,:]) for g in unique(group)]
end

function dispersion(D,group;metric = false)
    """
    Quantify the dispersion around the group spatial median for each group in D, according to the given distance metric.
    See Distances.jl for list of available metrics and the process of implementic your own metric. Once X is converted
    to a distance matrix, it is then transformed by principal coordinate analysis before medians and distances are calculated.
    The algorithm is based on Anderson (2006). But we aim to be consistent with the betadisper implementation in the
    R package Vegan, which has implemented some changes since the original paper was published.

    The function accepts:
        D: A dissimiarity or distance matrix
        group: A vector containing group labels
        metric: (key word) boolean value specifying whether or not the distance used used metric

    This function returns a named tuple containing:
        F = Global F-Statistic
        pairwise_F = pairwise F-statistics
        medians = spatial (aka geometric) median of each group in the transformed coordinates
        residuals = vector containing each observations distance from its group median
        means = vector containing the means distance to median for each group
        group = the original grouping vector
        levels = unique items from 'group'

    Anderson, M.J. (2006) Distance-based tests for homogeneity of multivariate dispersions. Biometrics 62(1), 245--253.
    https://github.com/vegandevs/vegan/blob/master/R/betadisper.R
    """
    levels = sort(unique(group))
    # vegdist objects contain lower triangular matrices wheras D is symmetric
    # thus we skip a D + D' step here.
    if metric == false
        # We double center the matrix as in vegan.betadisper. Though this is not identical Anderson (2006).
        A = double_center(D .^2 )
        vals, vecs = eigen(Hermitian(-A/2))
        # Anderson (2006) multiplies eigen vectors corresponding to negative eigen values by sqrt(-1).
        # We ommit this step in keeping with vegan.betadisper
        pco_space = vecs * diagm(sqrt.(abs.(vals)))
        # record indices corresponding to postitive eigenvalues
        positive_inds = real.(vals) .> 0.0
        med = geo_median(pco_space,ones(size(pco_space)[1]),positive_inds)
        medians = geo_median(pco_space, group, positive_inds)
        # calculate distances (to median) for "positive_inds" and "negative_inds" vectors separately. in orer to
        # subtract "negative_inds" from "positive_inds". See Anderson (2006) for details.
        dis_pos = [get_residuals(pco_space[group .== levels[i],positive_inds], medians[1][i]) for i in 1:length(levels)]
        dis_neg = [get_residuals(pco_space[group .== levels[i], .!positive_inds], medians[2][i]) for i in 1:length(levels)]
        # Where dis_neg > dis_pos, set residual = 0 by taking only the real part
        # of the square root of a negative. This was implemented in vegan.betadisper after discussion in
        # issue #306. But this is not done in Anderson (2006)
        residuals = [(real.(sqrt.(Complex.(dis_pos[i] .- dis_neg[i]))) )   for i in 1:length(dis_pos)]

        group_dis_pos = get_residuals(vcat(medians[1]...), med[1][1])
        group_dis_neg = get_residuals(vcat(medians[2]...), med[2][1])
        group_residuals = [(real.(sqrt.(Complex.(group_dis_pos[i] .- group_dis_neg[i]))) )   for i in 1:length(group_dis_pos)]

    else
        medians = geo_median(D, group)
        med = geo_median(D,ones(size(D)[1]))
        residuals = [get_residuals(D[group .== levels[i],:], medians[i]) for i in 1:length(levels)]
        group_residuals = vec(get_residuals(vcat(medians...), med[1]))
    end

    F = f(residuals)
    F_pairs = f_pairs(residuals)
    means = NamedArray(mean.(residuals),string.(levels),"group")
    return Disp(F, F_pairs, medians, residuals, means,group, levels,group_residuals)
end

function dispersion(X,group, distance)
    """
        Quantify the dispersion around the group spatial median for each group in X, according to the given distance metric.
        See Distances.jl for list of available metrics and the process of implementic your own metric. Once X is converted
        to a distance matrix, it is then transformed by principal coordinate analysis before medians and distances are calculated.
        The algorithm is based on Anderson (2006). But we aim to be consistent with the betadisper implementation in the
        R package Vegan, which has implemented some changes since the original paper was published.

        The function accepts:
            X: A numerical array
            group: A vector containing group labels
            distance: A distance/dissimilarity measure implemented in Distances.jl

            This function returns a named tuple containing:
            F = Global F-Statistic
            pairwise_F = pairwise F-statistics
            medians = spatial (aka geometric) median of each group in the transformed coordinates
            residuals = vector containing each observations distance from its group median
            means = vector containing the means distance to median for each group
            group = the original grouping vector
            levels = unique items from 'group'

        Anderson, M.J. (2006) Distance-based tests for homogeneity of multivariate dispersions. Biometrics 62(1), 245--253.
        https://github.com/vegandevs/vegan/blob/master/R/betadisper.R
    """

    D = Distances.pairwise(distance(),X,dims = 1)
    return dispersion(D,group)

end


text:
# TODO: https://openreview.net/forum?id=r1eyceSYPr
function ucd!(
    rbm::RBM, data::AbstractArray;
    batchsize::Int = 1
)
    error("not implemented yet")
end


text:
# Spinnaker.jl: wrapper for FLIR/Point Grey Spinnaker SDK
# Copyright (C) 2019 Samuel Powell

# System.jl: interface to Spinnaker system

"""
  Spinnaker SDK system object.

  System() returns new System object from which interfaces and devices can be discovered.
"""
mutable struct System
  handle::spinSystem

  function System()
    hsystem_ref = Ref(spinSystem(C_NULL))
    spinSystemGetInstance(hsystem_ref)
    @assert hsystem_ref[] != C_NULL
    sys = new(hsystem_ref[])
    finalizer(_release!, sys)
    return sys
  end
end

unsafe_convert(::Type{spinSystem}, sys::System) = sys.handle
unsafe_convert(::Type{Ptr{spinSystem}}, sys::System) = pointer_from_objref(sys)

# Release handle to system
function _release!(sys::System)
  spinSystemReleaseInstance(sys)
  sys.handle = C_NULL
  return nothing
end

"""
  version(::System)-> version, build

  Query the version number of the Spinnaker library in use. Returns a VersionNumber object
  and seperate build number.
"""
function version(sys::System)
  hlibver_ref = Ref(spinLibraryVersion(0,0,0,0))
  spinSystemGetLibraryVersion(sys, hlibver_ref)
  libver = VersionNumber(hlibver_ref[].major, hlibver_ref[].minor, hlibver_ref[].type)
  return libver, hlibver_ref[].build
end


text:
export
    example01,
    example02,
    example03,
    example04,
    #example05,
    example06,
    example07

function example01()
    x = linspace(0, 3pi, 100)
    c = cos(x)
    s = sin(x)

    p = FramedPlot(
            title="title!",
            xlabel="\\Sigma x^2_i",
            ylabel="\\Theta_i")

    add(p, FillBetween(x, c, x, s))
    add(p, Curve(x, c, color="red"))
    add(p, Curve(x, s, color="blue"))

    p
end

function example02()
    srand(42)
    p = FramedPlot(
            aspect_ratio=1,
            xrange=(0,100),
            yrange=(0,100))

    n = 21
    x = linspace(0, 100, n)
    yA = 40 .+ 10*randn(n)
    yB = x + 5*randn(n)

    a = Points(x, yA, kind="circle")
    setattr(a, label="a points")

    b = Points(x, yB)
    setattr(b, label="b points")
    style(b, kind="filled circle")

    s = Slope(1, (0,0), kind="dotted")
    setattr(s, label="slope")

    l = Legend(.1, .9, {a,b,s})

    add(p, s, a, b, l)
    p
end

function example03()
    srand(42)
    p = FramedPlot(
            title="Title",
            xlabel="X axis",
            ylabel="Y axis")

    add(p, Histogram(hist(randn(1000))...))
    add(p, PlotLabel(.5, .5, "Histogram", color=0xcc0000))

    t1 = Table(1, 2)
    t1[1,1] = p
    t1[1,2] = p

    t2 = Table(2, 1)
    t2[1,1] = t1
    t2[2,1] = p
    t2
end

function example04()
    x = linspace(0., 2pi, 40)
    s = sin(x)
    c = cos(x)

    inset = FramedPlot(title="inset")
    setattr(inset.frame, draw_ticks=false)

    add(inset, Curve(x, s, kind="dashed"))

    p = FramedPlot(aspect_ratio=1)
    setattr(p.frame, tickdir=+1, draw_spine=false)

    add(p, SymmetricErrorBarsY(x, s, 0.2*ones(length(x))))
    add(p, Points(x, s, color="red"))
    add(p, PlotInset((.6,.6), (.95,.95), inset))
    p
end

function example05()
    x = linspace(0., 2pi, 30)
    y = sin(x)

    p = FramedArray(2, 2,
            title="title",
            aspect_ratio=0.75,
            xlabel="x label",
            ylabel="y label",
            uniform_limits=true,
            cellspacing=1.)

    add(p, LineY(0, kind="dot"))

    add(p[1,1], Curve(x, .25*y))
    add(p[1,2], Curve(x, .50*y))
    add(p[2,1], Curve(x, .75*y))
    add(p[2,2], Curve(x, y))
    p
end

function example06()
    x = linspace(pi, 3pi, 60)
    c = cos(x)
    s = sin(x)

    p = FramedPlot(aspect_ratio=1)
    setattr(p.frame1, draw_grid=true, tickdir=1)

    setattr(p.x1, label="bottom", subticks=1)
    setattr(p.y1, label="left", draw_spine=false)
    setattr(p.x2, label="top", range=(10,1000), log=true)

    setattr(p.y2, label="right", draw_ticks=false,
        ticklabels=["-1", "-1/2", "0", "1/2", "1"])

    add(p, Curve(x, c, kind="dash"))
    add(p, Curve(x, s))
    p
end

function example07()
    n = 300
    x = linspace(-10., 10., n)
    t = linspace(-1., 1., n)
    z = (3. .+ 4*cosh(2x' .- 8t) .+ cosh(4x' .- 64t)) ./
        (3*cosh(x' .- 28t) + cosh(3x' .- 36t)) .^ 2

    t = Table(2,2)
    t[1,1] = imagesc(z, (minimum(z),0.6maximum(z)))

    p = imagesc(z)
    xlim(0, 1000)
    ylim(400, 0)
    t[2,1] = p

    p = imagesc(z)
    xlim(0, 1000)
    ylim(0, 400)
    t[1,2] = p

    t
end


text:
# stable periodic trajectories for ρ ∈ {350, 160, 100.5, 99.65, 28}
const σ=10
const β=8/3

lorenz(t,v) = [
	σ*(v[2]-v[1])
	v[1]*(v[4]-v[3])-v[2]
	v[1]*v[2]-β*v[3]
	]

lorenz′(t,v) = [
	-σ σ 0 0
	(v[4]-v[3]) -1 -v[1] v[1]
	v[2] v[1] -β 0
	]

# plot the fixpoints and intersections
@noinline function cbPlotSolution(v)
	ρ = v[end]
	d = sqrt(β*(ρ-1))
	scatter3D([d;-d], [d;-d], ρ-1, edgecolors="none", c="r")

	ts,fs,ns = lorenzProjectionInternal(v)
	tmp = transpose(reduce(hcat, fs))
	scatter3D(tmp[:,1],tmp[:,2],tmp[:,3])
end

# takes 2π-periodic functions, n > Nyquist rate
# returns set of scalars to distinguish solutions
#TODO oversampling?
@noinline function lorenzProjectionInternal(V)
	W,ω,ρ = unwrap(V)
	m = size(W,1)-1
	ftmp,N = interps(W), 2m

	d = sqrt(β*(ρ-1))
	A,B = [d;d;ρ-1],[-d;-d;ρ-1]
	const p = 0.08872374069251765
	T = linspace(-pi+p,pi+p,2N)

	pA(t) = ftmp[3](t) - (ρ+7) #ftmp[1](t) + ftmp[2](t) + 2d #ftmp[1](t) - ftmp[2](t) #ftmp[3](t) - ρ + 1 ... interesting
	# pB(t) = ftmp[1](t) + ftmp[2](t) - 2d
	t₀ = T[1]
	a₀ = pA(T[1])
	# b₀ = pB(T[1])

	ts,fs,ns = [],[],[]
	for t in T[2:end]
		tmp = pA(t)
		if a₀ ≤ 0 < tmp
			push!(ts, bisection(pA, t₀, t))
			push!(fs, map(g->g(ts[end]), ftmp))
			push!(ns, norm(fs[end]-A))
		end
		a₀ = tmp

		# tmp = pB(t)
		# if b₀ ≤ 0 < tmp
		# 	push!(ts, bisection(pB, t₀, t))
		# 	push!(fs, map(g->g(ts[end]), ftmp))
		# 	push!(ns, -norm(fs[end]-B))
		# end
		# b₀ = tmp

		t₀ = t
	end

	return ts, fs, sort(ns)
end

@noinline lorenzProjection(V) = lorenzProjectionInternal(V)[3]
#@noinline lorenzProjection(V) = [V[end-1]]

#=
lorenzProjection(v...) = lorenzProjection2(v...) function lorenzProjection2(ftmp,ρ,N)
	T = linspace(-.1,2pi-.1,3N+2)
	plane(t) = ftmp[1](t) - ftmp[2](t)
	t₀,p₀ = T[1], plane(T[1])

	rtn = []
	for t in T[2:end-1]
		tmp = plane(t)
		p₀ ≤ 0 < tmp && push!(rtn, bisection(plane, t₀, t, ϵ=1e-4))
		t₀,p₀ = t,tmp
	end

	return map(x->norm(map(g->g(x), ftmp)), rtn)
end
=#

# N = 2m
# ρ = ℵ
# tmp = sqrt(β*(ρ-1))
# # fixpoints
# A,B = [tmp; tmp; ρ-1], [-tmp; -tmp; ρ-1]
#
# # vector from one fixpoint to the other, used as normal
# n = [ 2tmp; 2tmp; 0 ]
#
# T = linspace(-.1,2pi-.1,N+2)
# pA(t) = dot(n, f(T[1])-A)
# pB(t) = dot(n, f(T[1])-B)
# t₀ = T[1]
# a₀, b₀ = pA(t₀), pB(t₀)
#
# rtn = []
#
# I = T[2:end-1]
# state = start(I)
# t,state = next(I,state)
# tmp = pA(t)
# a₀ ≤ 0 < tmp && push!(rtn, bisection(pA, t₀, t, ϵ=1e-4))
# a₀ = tmp
#
# tmp = pB(t)
# b₀ ≤ 0 < tmp && push!(rtn, bisection(pB, t₀, t, ϵ=1e-4))
# b₀ = tmp
#
# t₀ = t
#
# rtn
#
# return map(f, rtn)

#=
function orthogonal_proj(zfront, zback)
    a = (zfront+zback)/(zfront-zback)
    b = -2*(zfront*zback)/(zfront-zback)
    return [
		1 0 0 0
		0 1 0 0
		0 0 a b
		0 0 0 zback]
end
proj3d.persp_transformation = orthogonal_proj
=#


text:
module OpenDriveManagerLite

using Lint
using Compat

macro lintpragma(s)
end

##############################################################################
##
## Ensure lib can be found
##
##############################################################################
export
    AbstractTrackCoord,

    times!,
	plus!,
    minus!,
    init!

import Base: copy, copy!, deepcopy, show, print, ==, *, +, -

const LIB_ODRMGR = let
    curpath = dirname(@__FILE__() )
    pathtoadd = curpath[1:end-3]*"deps/downloads/"
    find_library(["libodrmanjl"], [pathtoadd])
end

if isempty(LIB_ODRMGR)
    warn("Could not communicate with libodrmanjl. Make sure it exists and that it is on the path.")
end

include("includes/trackcoord.jl")
include("includes/lanecoord.jl")
include("includes/coord.jl")
include("includes/odrmanagerlite.jl")

end # module

max_stars_repo_name: sisl/OpenDriveManagerLite

max_issues_repo_name: sisl/OpenDriveManagerLite

max_forks_repo_name: sisl/OpenDriveManagerLite


text:
using DiffEqSensitivity, Flux, Zygote, DelayDiffEq, Test

## Setup DDE to optimize
function delay_lotka_volterra(du,u,h,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = (α - β*y)*h(p,t-0.1)[1]
  du[2] = dy = (δ*x - γ)*y
end
h(p,t) = ones(eltype(p),2)
prob = DDEProblem(delay_lotka_volterra,[1.0,1.0],h,(0.0,10.0),constant_lags=[0.1])
p = [2.2, 1.0, 2.0, 0.4]
function predict_fd_dde(p)
  solve(prob,MethodOfSteps(Tsit5()),p=p,saveat=0.0:0.1:10.0,reltol=1e-4,sensealg=ForwardDiffSensitivity())[1,:]
end
loss_fd_dde(p) = sum(abs2,x-1 for x in predict_fd_dde(p))
loss_fd_dde(p)
@test !iszero(Zygote.gradient(loss_fd_dde,p)[1])

function predict_rd_dde(p)
  solve(prob,MethodOfSteps(Tsit5()),p=p,saveat=0.1,reltol=1e-4,sensealg=TrackerAdjoint())[1,:]
end
loss_rd_dde(p) = sum(abs2,x-1 for x in predict_rd_dde(p))
loss_rd_dde(p)
@test !iszero(Zygote.gradient(loss_rd_dde,p)[1])

@test Zygote.gradient(loss_fd_dde,p)[1] ≈ Zygote.gradient(loss_rd_dde,p)[1] rtol=1e-2


text:
# This file contans the common interface for LightGraphs.

_NI(m) = error("Not implemented: $m")

"""
    AbstractEdge

An absract type representing a single edge between two vertices of a graph.
"""
abstract type AbstractEdge end

"""
    AbstractEdgeIter

An abstract type representing an edge iterator.
"""
abstract type AbstractEdgeIter end

"""
    AbstractGraph

An abstract type representing a graph.
"""
abstract type AbstractGraph end

@traitdef IsDirected{G<:AbstractGraph}
@traitimpl IsDirected{G} <- is_directed(G)

#
# Interface for AbstractEdges
#

"""
    src(e)

Return the source vertex of edge `e`.
"""
src(e::AbstractEdge) = _NI("src")

"""
    dst(e)

Return the destination vertex of edge `e`.
"""
dst(e::AbstractEdge) = _NI("dst")

Pair(e::AbstractEdge) = _NI("Pair")
Tuple(e::AbstractEdge) = _NI("Tuple")

"""
    reverse(e)

Create a new edge from `e` with source and destination vertices reversed.
"""
reverse(e::AbstractEdge) = _NI("reverse")

==(e1::AbstractEdge, e2::AbstractEdge) = _NI("==")

#
# Interface for AbstractGraphs
#
"""
    edgetype(g)

Return the type of graph `g`'s edge
"""
edgetype(g::AbstractGraph) = _NI("edgetype")

"""
    eltype(g)

Return the type of the graph's vertices (must be <: Integer)
"""
eltype(g::AbstractGraph) = _NI("eltype")

"""
    nv(g)

Return the number of vertices in `g`.
"""
nv(g::AbstractGraph) = _NI("nv")

"""
    ne(g)

Return the number of edges in `g`.
"""
ne(g::AbstractGraph) = _NI("ne")

"""
    vertices(g)

Return (an iterator to or collection of) the vertices of a graph.

### Implementation Notes
A returned iterator is valid for one pass over the edges, and
is invalidated by changes to `g`.

"""
vertices(g::AbstractGraph) = _NI("vertices")

"""
    edges(g)
Return (an iterator to or collection of) the edges of a graph.

### Implementation Notes
A returned iterator is valid for one pass over the edges, and
is invalidated by changes to `g`.
"""
edges(x...) = _NI("edges")

is_directed(x...) = _NI("is_directed")
is_directed{T}(::Type{T}) = _NI("is_directed")
"""
    add_vertex!(g)

Add a new vertex to the graph `g`.
Return true if the vertex was added successfully, false otherwise.
"""
add_vertex!(x...) = _NI("add_vertex!")

"""
    add_edge!(g, e)

Add a new edge `e` to `g`. Return false if add fails
(e.g., if vertices are not in the graph, or edge already exists), true otherwise.
"""
add_edge!(x...) = _NI("add_edge!")

"""
    rem_vertex!(g)

Remove the vertex `v` from graph `g`. Return false if removal fails
(e.g., if vertex is not in the graph), true otherwise.
"""
rem_vertex!(x...) = _NI("rem_vertex!")

"""
    rem_edge!(g, e)

Remove the edge `e` from `g`. Return false if edge removal fails
(e.g., if edge does not exist), true otherwise.
"""
rem_edge!(x...) = _NI("rem_edge!")

"""
    has_vertex(g, v)

Return true if `v` is a vertex of `g`.
"""
has_vertex(x...) = _NI("has_vertex")

"""
    has_edge(g, e)

Return true if the graph `g` has an edge `e`.
"""
has_edge(x...) = _NI("has_edge")

"""
    in_neighbors(g, v)

Return a list of all neighbors connected to vertex `v` by an incoming edge.

### Implementation Notes
Returns a reference, not a copy. Do not modify result.
"""
in_neighbors(x...) = _NI("in_neighbors")

"""
    out_neighbors(g, v)

Return a list of all neighbors connected to vertex `v` by an outgoing edge.

# Implementation Notes
Returns a reference, not a copy. Do not modify result.
"""
out_neighbors(x...) = _NI("out_neighbors")

"""
    zero(g)

Return a zero-vertex, zero-edge version of the same type of graph as `g`.
"""
zero(g::AbstractGraph) = _NI("zero")

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles