text:
const ColorMetric = Colors.DifferenceMetric

"""
    CIEDE2000 <: SemiMetric
    ciede2000(x, y; metric=DE_2000())

The pixel-wise `CIEDE2000` color-difference between image `a` and `b`.
Optionally, a metric may be supplied, chosen among `DE_2000` (the default), `DE_94`,
`DE_JPC79`, `DE_CMC`, `DE_BFD`, `DE_AB`, `DE_DIN99`, `DE_DIN99d`, `DE_DIN99o`.

## References

Sharma, G., Wu, W., and Dalal, E. N., 2005. *The CIEDE2000 color‐difference formula*.
"""
struct CIEDE2000{M<:ColorMetric} <: SemiMetric
    metric::M
end

CIEDE2000() = CIEDE2000(DE_2000())

# colordiff converts image to Lab space, so we don't promote the storage type here
function (d::CIEDE2000)(imgA::GenericImage, imgB::GenericImage)
    sum(abs, colordiff.(imgA, imgB; metric=d.metric))
end

# helper function
@doc (@doc CIEDE2000)
ciede2000(imgA::GenericImage, imgB::GenericImage; metric=DE_2000()) =
    CIEDE2000(metric)(imgA, imgB)


text:
# ____________________________________________________________________
# SpecLineGauss
#
mutable struct SpecLineGauss <: AbstractComponent
    norm::Parameter
    center::Parameter
    fwhm::Parameter
    voff::Parameter
    resolution::Float64

    function SpecLineGauss(center::Number)
        out = new(Parameter(1),
                  Parameter(center),
                  Parameter(3000),
                  Parameter(0),
                  0.)
        @assert center > 0
        out.norm.low = 0
        out.center.low = 0
        out.fwhm.low = 0
        out.voff.low = 0
        out.center.fixed = true
        return out
    end
end

function prepare!(comp::SpecLineGauss, domain::Domain{1})
    return fill(NaN, length(domain))
end

function evaluate!(buffer, comp::SpecLineGauss, x::Domain{1},
                   norm, center, fwhm, voff)
    x0 = center - (voff / 3.e5) * center
    σ_res = comp.resolution / 2.355 / 3.e5 * center
    σ     = fwhm            / 2.355 / 3.e5 * center
    σ = sqrt(σ^2 + σ_res^2)
    X = (x .- x0) ./ σ

    function profile(x)
        (abs(x) > 4)  &&  (return 0.)
        return norm * exp(-x^2 / 2) / sqrt(2pi) / σ
    end
    map!(profile, buffer, X)
end

ext: jl


text:
module CoupledFields

using LinearAlgebra
using Statistics
using StatsBase: sample, zscore
export InputSpace, ModelObj
export KernelParameters, GaussianKP, PolynomialKP
export gradvecfield
export bf, cca, gKCCA
export CVfn, Rsq_adj
export whiten
include("MLKernels.jl")

"""
    KernelParameters: An abstract type.
All KernelParameters types contain certain parameters which are later passed to internal functions `Kf` and `∇Kf`. \n
A KernelParameters type is set using e.g. `PolynomialKP(X::Matrix{Float64})` or `GaussianKP(X::Matrix{Float64})`.
"""
abstract type KernelParameters end

"""
    GaussianKP: For the gaussian kernel
"""
struct GaussianKP <: KernelParameters
    xx::Matrix{Float64}
    varx::Float64
end

function GaussianKP(X::Matrix{Float64})
    xx = kernelmatrix(SquaredDistanceKernel(1.0), X, X)
    varx = Statistics.median(xx[xx.>10.0^-9])
    return GaussianKP(xx, varx)
end

    function Kf(par::Array{T}, X::Matrix{T}, kpars::GaussianKP) where T<:Float64
        sx2 = 2*par[1]*par[1]*kpars.varx
        return exp.(-kpars.xx/sx2)
    end

    function ∇Kf(par::Array{T}, X::Matrix{T}, kpars::GaussianKP) where T<:Float64
        n,p = size(X)
        sx2 = 2*par[1]*par[1]*kpars.varx
        Gx = Kf(par, X, kpars)
        ∇K = Float64[2.0*Gx[i,k]*(X[i,j]-X[k,j])/sx2  for i in 1:n, k in 1:n, j in 1:p ]
        return ∇K
    end

"""
    PolynomialKP: For the polynomial kernel
"""
struct PolynomialKP <: KernelParameters
    xx::Matrix{Float64}
    function PolynomialKP(X::Matrix{Float64})
        xx = kernelmatrix(LinearKernel(1.0, 1.0), X, X)
        return new(xx)
    end
end

    function Kf(par::Array{T}, X::Matrix{T}, kpars::PolynomialKP) where T<:Float64
        return (kpars.xx).^par[1]
    end

    function ∇Kf(par::Array{T}, X::Matrix{T}, kpars::PolynomialKP) where T<:Float64
        n,p = size(X)
        m = par[1]
        ∇K = Float64[m*X[k,j]*kpars.xx[i,k]^(m-1.0)  for i in 1:n, k in 1:n, j in 1:p ]
        return ∇K
    end

"""
    ModelObj: A type to hold statistical model results
Such as the matrices `W, R, A, T`, where `R=XW` and `T=YA`.
"""
struct ModelObj
    W::Matrix{Float64}
    R::Matrix{Float64}
    A::Matrix{Float64}
    T::Matrix{Float64}
    evals::Vector{Float64}
    pars::Array{Float64}
    method::String
end

"""
    InputSpace: A type to hold the `X` and `Y` fields of the Input space
InputSpace(X, Y, d, lat): The fields are whitened if `d=[d1, d2]` is supplied.
Area weighting is applied if `lat` is supplied.
"""
struct InputSpace
    X::Matrix{Float64}
    Y::Matrix{Float64}
    function InputSpace(a::T, b::T) where T<:Matrix{Float64}
        new(zscore(a,1), zscore(b,1))
    end
    function InputSpace(a::T, b::T, d::Vector{Float64}) where T<:Matrix{Float64}
        new(whiten(a, d[1]), whiten(b, d[2]))
    end
    function InputSpace(a::T, b::T, d::V, lat::V) where {V<:Matrix{Float64}, T<:Vector{Float64}}
        new(whiten(a, d[1], lat=lat), whiten(b, d[2], lat=lat))
    end
end

"""
    gradvecfield(par::Array{Float64}, X::T, Y::T, kpars::KernelParameters) where T<:Matrix{Float64}:
Compute the gradient vector or gradient matrix at each instance of the `X` and `Y` fields, by making use of a kernel feature space.
"""
function gradvecfield(par::Array{Float64}, X::T, Y::T, kpars::KernelParameters ) where T<:Matrix{Float64}
    n,p = size(X)
    Gx = Kf(par, X, kpars) + (10.0^par[2])*n*I
    ∇K = ∇Kf(par, X, kpars)
    return [∇K[i,:,:]' * (Gx \ Y) for i in 1:n]
end

"""
    bf(x::Vector{Float64}, df::Int):
Compute a piecewise linear basis matrix for the vector x.
"""
function bf(x::Vector{Float64}, df::Int)
    if df<2 return x[:,1:1] end
    bp = Statistics.quantile(x, range(0, stop=1, length=df+1))
    n = length(x)
    a1 = repeat(x, inner=(1,df))
    a2 = repeat( permutedims(bp[1:df]), inner=(n,1))

    for j in 1:df
        ix1, ix2 = x.<bp[j], x.>bp[j+1]
        a1[ix1,j] .= bp[j]
        a1[ix2,j] .= bp[j+1]
    end

    return a1-a2
end

"""
    cca(v::Array{Float64}, X::T,Y::T) where T<:Matrix{Float64}:
Regularized Canonical Correlation Analysis using SVD.
"""
function cca(v::Array{Float64}, X::T,Y::T) where T<:Matrix{Float64}
#    n,p = size(X)
    q = size(Y, 2)

    Cxx_fac_inv = (Statistics.cov(X)+(10.0^v[1])*I)^-0.5
    Cyy_fac_inv = (q>1) ? (Statistics.cov(Y)+(10.0^v[2])*I)^-0.5 : Statistics.cov(Y)^(-0.5)
    M = Statistics.cov(X*Cxx_fac_inv, Y*Cyy_fac_inv)
    U, D, V = svd(M)

    Wx = Cxx_fac_inv * U
    Wy = Cyy_fac_inv * V
    Tx = X * Wx
    Ty = Y * Wy

    L = D.^2
    return ModelObj(Wx, Tx, Wy, Ty, L, v, "CCA")
end

cca(v::Array{Float64}, X::T,Y::T, kpars::KernelParameters) where {T<:Matrix{Float64}} = cca(v, X, Y)

"""
    gKCCA(par::Array{Float64}, X::Matrix{Float64}, Y::Matrix{Float64}, kpars::KernelParameters):
Compute the projection matrices and components for gKCCA.
"""
function gKCCA(par::Array{Float64}, X::Matrix{Float64}, Y::Matrix{Float64}, kpars::KernelParameters)
    q = size(Y,2)
    n, p = size(X)
    ∇g = gradvecfield(par, X, Y, kpars)
    M = mapreduce(x -> x*x', +, ∇g)/n
    values, vectors = eigen(Symmetric(M))
    Wx = reverse(vectors, dims=2)
    R = X * Wx

    dmin = minimum([p q 3])
    Ay = zeros(q, dmin)
    T = zeros(n, dmin)
    for j in 1:dmin
        ϕRj = bf(R[:,j], Int(par[3]))
        cc = cca([-9.0 -9.0], ϕRj, Y)
        Ay[:,j] = cc.A[:,1]
        T[:,j] = cc.T[:,1]
    end

    return ModelObj(Wx, R, Ay, T, reverse(values), par, "gKCCA")
end

"""
    whiten(x::Matrix{Float64}, d::Float64; lat=nothing): Whiten matrix
`d` (0-1) Percentage variance of components to retain. \n
`lat` Latitudinal area-weighting.
"""
function whiten(X::Matrix{Float64}, d::Float64; lat=nothing)
    m1 = zscore(X, 1)
#    lat != nothing && scale!(m1, 1.0./sqrt(cos(pi*lat/180)) )
    U,S,V = svd(m1)
    l = cumsum(S.^2)/sum(abs2, S)
    U = U[:, l.≤d]
    U ./= Statistics.std(U, dims=1)
    return U
end

"""
    CVfn(parm::T, X::T, Y::T, modelfn::Function, kerneltype::DataType; verbose::Bool=true, dcv::Int64=2) where T<:Matrix{Float64}
Cross-validation function
"""
function CVfn(parm::T, X::T, Y::T, modelfn::Function, kerneltype::DataType; verbose::Bool=true, dcv::Int64=2) where T<:Matrix{Float64}

    # free parameters
    # df = par[3] # Number of segments in PWLM
    # dcv = 2 # Number of components used in CV

    # other params
    nrg = size(parm, 1)
    pargrid = Float64[parm zeros(nrg)]
    n, p = size(X)

    # Setup
    yr = div(0:n-1,12)+1
    yrs = unique(yr)
    ta = findall(in(sample(yrs, div(length(yrs),2)+1, replace=false, ordered=true)), yr)
    tb = setdiff(1:n,ta)

    xm1, xm2 = mean(X[ta,:], dims=1), mean(X[tb,:], dims=1)
    ym1, ym2 = mean(Y[ta,:], dims=1), mean(Y[tb,:], dims=1)

    X1train, Y1train = X[ta,:].-xm1, Y[ta,:].-ym1
    X2train, Y2train = X[tb,:].-xm2, Y[tb,:].-ym2
    X1test, Y1test = X[ta,:].-xm2, Y[ta,:].-ym2
    X2test, Y2test = X[tb,:].-xm1, Y[tb,:].-ym1

    k1pars = kerneltype(X1train)
    k2pars = kerneltype(X2train)

    # Main loop

    for i in 1:nrg
        percent = round(Int, 100*i/nrg)
        if verbose
           print(string(percent,"% "))
        end
        par = parm[i,:]
        model1 = modelfn(par, X1train, Y1train, k1pars)
        model2 = modelfn(par, X2train, Y2train, k2pars)
        W1, W2 = model1.W[:,1:dcv], model2.W[:,1:dcv]
        A1, A2 = model1.A[:,1:dcv], model2.A[:,1:dcv]
        W1 .*= sign.(W1[1:1,:])
        W2 .*= sign.(W2[1:1,:])
        A1 .*= sign.(A1[1:1,:])
        A2 .*= sign.(A2[1:1,:])
        Rcv = [X2test * W1; X1test * W2]
        Tcv = [Y2test * A1; Y1test * A2]
        pargrid[i,end] = Rsq_adj(Rcv, Tcv, Int(par[3]))
    end

    imax = argmax(pargrid[:,end])
    return pargrid[imax,:]
end

"""
    Rsq_adj(Tx::T, Ty::T, df::Int) where T<:Array{Float64}:
Cross-validation metric
"""
function Rsq_adj(Tx::T, Ty::T, df::Int) where {T<:Array{Float64}}
    y = vec(Ty)
    n,p = size(Tx)
    o = ones(n*p)
    mL = [ bf(Tx[:,j], df) for j in 1:p  ]
#    mL = [ [ bf(Tx[:,j], df) o] for j in 1:p  ]
    Xs =  [ cat([1,2], mL...) o]   # Block diagonal
    px = p*df
    nx = p*n
    Σyhat = (y'*Xs / (Xs' * Xs) ) * (Xs' * y)
    Rsq = ((Σyhat) / (y'y))[1]
    return Rsq - (1-Rsq)*px/(nx-px-1)
end

end


text:
export destructure

function _restructure(m, xs)
  i = 0
  fmap(m) do x
    x isa AbstractArray || return x
    x = reshape(xs[i.+(1:length(x))], size(x))
    i += length(x)
    return x
  end
end

"""
    destructure(m)
Flatten a model's parameters into a single weight vector.
    julia> m = Chain(Dense(10, 5, σ), Dense(5, 2), softmax)
    Chain(Dense(10, 5, σ), Dense(5, 2), softmax)
    julia> θ, re = destructure(m);
    julia> θ
    67-element Array{Float32,1}:
    -0.1407104
    ...
The second return value `re` allows you to reconstruct the original network after making
modifications to the weight vector (for example, with a hypernetwork).
    julia> re(θ .* 2)
    Chain(Dense(10, 5, σ), Dense(5, 2), softmax)
"""
function destructure(m)
  xs = Zygote.Buffer([])
  fmap(m) do x
    x isa AbstractArray && push!(xs, x)
    return x
  end
  return vcat(vec.(copy(xs))...), p -> _restructure(m, p)
end


text:
module TestPkg

using DifferentialEquations

export oscillator, solve_ode, relaxation

function oscillator(u,p,t)
  x, y, z = u
  A, B, G, K, N = p
  du1 = B[1]*(A[1] - x) + G[1,3] * (1.0 - x) * abs(z)^N[1,3] / (abs(z)^N[1,3] + K[1,3]^N[1,3])
  du2 = B[2]*(A[2] - y) + G[2,1] * (1.0 - y) * abs(x)^N[2,1] / (abs(x)^N[2,1] + K[2,1]^N[2,1])
  du3 = B[3]*(A[3] - z) - G[3,2] * z * abs(y)^N[3,2] / (abs(y)^N[3,2] + K[3,2]^N[3,2])
  [du1,du2,du3]
end

function relaxation(u,p,t)
  x, y, z = u
  A, B, G, K, N = p
  du1 = B[1]*(A[1] - x) + G[1,3] * (1.0 - x) * abs(z)^N[1,3] / (abs(z)^N[1,3] + K[1,3]^N[1,3])
  du2 = B[2]*(A[2] - y) + G[2,1] * (1.0 - y) * abs(x)^N[2,1] / (abs(x)^N[2,1] + K[2,1]^N[2,1])
  du3 = B[3]*(A[3] - z) - G[3,2] * z * abs(y)^N[3,2] / (abs(y)^N[3,2] + K[3,2]^N[3,2]) + G[3,3] * (1.0 - z) * abs(z)^N[3,3] / (abs(z)^N[3,3] + K[3,3]^N[3,3])
  [du1,du2,du3]
end

function solve_ode(prob)
  sol = solve(prob, isoutofdomain=(y,p,t)->any(x->x<0,y))
  return sol
end

end

max_stars_repo_name: ftavella/TestPkg

max_issues_repo_name: ftavella/TestPkg

max_forks_repo_name: ftavella/TestPkg


text:
################################################################################
## Modified Cressie Read
################################################################################

#=---------------
evaluate
---------------=#
function evaluate{T<:AbstractFloat}(div::MCR, a::AbstractVector{T}, b::AbstractVector{T})

    if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end

    r = zero(T)

    α  =  div.α
    ϑ   = div.ϑ
    f0, f1, f2, uϑ = div.m

    for i = eachindex(a, b)
        @inbounds ai = a[i]
        @inbounds bi = b[i]
        @inbounds ui = ai/bi
        if ui>=uϑ
      		  r += (f0 + f1*(ui-uϑ) + .5*f2*(ui-uϑ)^2)*bi
        elseif ui > 0 && ui < u₀
            r += ( (ui^(1+α)-1)/(α*(1.0 + α)) - (ui-1)/α )*bi
        elseif ui == 0
            r += bi/(1.0 + α)
        else
            r = +Inf
            break
        end
    end
    return r
end

function evaluate{T<:AbstractFloat}(div::MCR, a::AbstractVector{T})
    r = zero(T)
    α  =  div.α
    ϑ   = div.ϑ
    f0, f1, f2, uϑ = div.m

    @inbounds for i = eachindex(a)
        ui = a[i]
        if ui>=uϑ
      		r += f0 + f1*(ui-uϑ) + .5*f2*(ui-uϑ)^2
        elseif ui > 0 && ui < uϑ
            r += ( (ui^(1+α)-1)/(α*(1.0 + α)) - (ui-1)/α )
        elseif ui==0
            r += 1/(1.0 + α)
        else
            r = +Inf
            break
        end
    end
    return r
end

#=---------------
gradient
---------------=#
function gradient{T<:AbstractFloat}(div::MCR, a::T)
    α  =  div.α
    ϑ  = div.ϑ
    f0, f1, f2, uϑ = div.m
    if a>=uϑ
        u = f1 + f2*(a-uϑ)
    else
        u = gradient(div.d, a)
    end
    return u
end

function gradient{T<:AbstractFloat}(div::MCR, a::T, b::T)
    α  =  div.α
    ϑ  = div.ϑ
    f0, f1, f2, uϑ = div.m
    ui   = a/b
    if ui > uϑ
        u = f1 + f2*(ui-uϑ)
    else
        u = gradient(div.d, a, b)
    end
    return u
end

#=---------------
hessian
---------------=#
function hessian{T<:AbstractFloat}(div::MCR, a::T)
    f0, f1, f2, uϑ = div.m
    if a >= uϑ
        u = f2
    else
        u = hessian(div.d, a)
    end
    return u
end

function hessian{T<:AbstractFloat}(div::MCR, a::T, b::T)
    f0, f1, f2, uϑ = div.m
    ui = a/b
    if ui > uϑ
        u = f2/b
    else
        u = hessian(div.d, a, b)
    end
    return u
end

################################################################################
## Fully Modified Cressie Read
################################################################################

#=---------------
evaluate
---------------=#
function evaluate{T<:AbstractFloat}(div::FMCR, a::AbstractVector{T}, b::AbstractVector{T})

    if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end

    r = zero(T)
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m

    for i = eachindex(a, b)
        @inbounds ai = a[i]
        @inbounds bi = b[i]
        @inbounds ui = ai/bi
        if ui >= uϑ
      		  r += (f0 + f1*(ui-uϑ) + .5*f2*(ui-uϑ)^2)*bi
        elseif ui <= uφ
            r += (g0 + g1*(ui-uφ) + .5*g2*(ui-uφ)^2)*bi
        else
            r += ( (ui^(1+α)-1)/(α*(1.0 + α)) - (ui-1)/α )*bi
        end
    end
    return r
end

function evaluate{T<:AbstractFloat}(div::FMCR, a::AbstractVector{T})
    r = zero(T)
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m

    for i = eachindex(a)
        @inbounds ui = a[i]
        if ui >= uϑ
      		  r += (f0 + f1*(ui-uϑ) + .5*f2*(ui-uϑ)^2)
        elseif ui <= uφ
            r += (g0 + g1*(ui-uφ) + .5*g2*(ui-uφ)^2)
        else
            r += ( (ui^(1+α)-1)/(α*(1.0 + α)) - (ui-1)/α )
        end
    end
    return r
end

#=---------------
gradient
---------------=#
function gradient{T<:AbstractFloat}(div::FMCR, a::T)
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m
    if a >= uϑ
        u = f1 + f2*(a-uϑ)
    elseif a <= uφ
        u = g1 + g2*(a-uφ)
    else
        u = gradient(div.d, a)
    end
    return u
end

function gradient{T<:AbstractFloat}(div::FMCR, a::T, b::T)
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m
    ui = a/b
    if ui >= uϑ
        u = f1 + f2*(ui-uϑ)
    elseif ui <= uφ
        u = g1 + g2*(ui-uφ)
    else
        u = gradient(div.d, a, b)
    end
    return u
end

#=---------------
hessian
---------------=#
function hessian{T<:AbstractFloat}(div::FMCR, a::T)
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m
    if a >= uϑ
        u = f2
    elseif a <= uφ
        u = g2
    else
        u = hessian(div.d, a)
    end
    return u
end

function hessian{T<:AbstractFloat}(div::FMCR, a::T, b::T)
    f0, f1, f2, uϑ, g0, g1, g2, uφ = div.m
    ui = a/b
    if ui >= uϑ
        u = f2/b
    elseif ui <= uφ
        u = g2/b
    else
        u = hessian(div.d, a, b)
    end
    return u
end


text:
# FASTA Index
# ===========
#
# Index for random access to FASTA files.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

# http://www.htslib.org/doc/faidx.html
struct Index
    names::Vector{String}
    lengths::Vector{Int}
    offsets::Vector{Int}
    linebases::Vector{Int}
    linewidths::Vector{Int}
end

function Index(filepath::AbstractString)
    return open(read_faidx, filepath)
end

function read_faidx(input::IO)
    names = String[]
    lengths = Int[]
    offsets = Int[]
    linebases = Int[]
    linewidths = Int[]
    for line in eachline(input)
        values = split(chomp(line), '\t')
        name = values[1]
        length = parse(Int, values[2])
        offset = parse(Int, values[3])
        linebase = parse(Int, values[4])
        linewidth = parse(Int, values[5])
        push!(names, name)
        push!(lengths, length)
        push!(offsets, offset)
        push!(linebases, linebase)
        push!(linewidths, linewidth)
    end
    return Index(names, lengths, offsets, linebases, linewidths)
end

# Set the reading position of `input` to the starting position of the record `name`.
function seekrecord(input::IO, index::Index, name::AbstractString)
    i = findfirst(x -> x == name, index.names)
    if i == 0
        throw(ArgumentError("sequence \"$(name)\" is not in the index"))
    end
    offset = index.offsets[i]
    n_back = 100
    @label seekback
    seek(input, max(offset - n_back, 0))
    data = UInt8[]
    while position(input) < offset
        push!(data, read(input, UInt8))
    end
    for j in lastindex(data):-1:1
        if data[j] == UInt8('>') && (offset ≤ n_back || (j ≥ 2 && data[j-1] == UInt8('\n')))
            seek(input, offset - (lastindex(data) - j + 1))
            return
        end
    end
    if n_back ≥ offset
        # reached the starting position of the input
        error("failed to find the starting position of the record")
    end
    n_back *= 2
    @goto seekback
end


text:
# Get position of walker
position(w::Walker1D) = w.x
position(w::Walker2D) = (w.x, w.y)


text:
# return perm[newid]=oldid sorted by idtable
# perm[newid]  = origid
# invp[origid] = newid
function sort_by_dict(idtable::Dict{String,Int},n::Int)
   perm = zeros(Int,n)
   invp = zeros(Int,n)
   newid = 0
   for key in sort(collect(keys(idtable)))
      origid = idtable[key]
      if key!="0" && origid<=n
         newid = newid + 1
         perm[newid] = origid
         invp[origid] = newid
      end
   end
   return perm,invp
end

function get_qb(upgtable)
   nunk = upgtable.count
   nupg = length(unique(values(upgtable)))
   Qb = zeros(nunk,nupg)
   i = 0
   for key in sort(collect(keys(upgtable)))
      i = i + 1
      group = upgtable[key]
      Qb[i,group] = 1.0
   end
   return Qb
end


text:
using AlphaShapes
include("../src/utils.jl")
using Plots
using Main.AlphaShapes
using LinearAlgebra
using Random
Random.seed!(31415926535897)

N = 1000
r1 = 0.75
r2 = 2.0
X = randn(N,2)
ind = [((norm(X[i,:])>r1) + (norm(X[i,:])<r2)) == 2 for i in 1:size(X,1)]
X = X[ind.==true,:]
# normally distributed points with a whole in the middle
p1 = scatter(X[:,1],X[:,2],label="",aspect_ratio=:equal)

T = AlphaShapes.GetDelaunayTriangulation(X);
p1 = scatter(X[:,1],X[:,2],label="",aspect_ratio=:equal)
[DrawTri!(T[i,:,:]) for i in 1:size(T,1)]
plot!(title="Delaunay Triangulation")

A = AlphaShape(X);
p2 = scatter(X[:,1],X[:,2],label="",aspect_ratio=:equal)
[DrawTri!(A[i,:,:]) for i in 1:size(A,1)]
plot!(title="Optimum Alpha Shape Triangulation")

p = plot(p1,p2,size=(900,300))
savefig(p,"Example.png")


text:
using OrdinaryDiffEq, ParameterizedFunctions, ODE, ODEInterfaceDiffEq, LSODA,
      Sundials, DiffEqDevTools

f = @ode_def LotkaVolterra begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d

p = [1.5,1.0,3.0,1.0]
prob = ODEProblem(f,[1.0;1.0],(0.0,10.0),p)

abstols = 1.0 ./ 10.0 .^ (6:13)
reltols = 1.0 ./ 10.0 .^ (3:10);
sol = solve(prob,Vern7(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)
using Plots; gr()

plot(sol)

setups = [Dict(:alg=>DP5())
          #Dict(:alg=>ode45()) # fail
          Dict(:alg=>dopri5())
          Dict(:alg=>Tsit5())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,maxiters=10000,numruns=100)
plot(wp)

setups = [Dict(:alg=>DP5())
          #Dict(:alg=>ode45())
          Dict(:alg=>Tsit5())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,maxiters=10000,error_estimate=:L2,dense_errors=true,numruns=100)
plot(wp)

setups = [Dict(:alg=>DP8())
          #Dict(:alg=>ode78()) # fails
          Dict(:alg=>Vern7())
          Dict(:alg=>Vern8())
          Dict(:alg=>dop853())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,maxiters=1000,numruns=100)
plot(wp)

setups = [Dict(:alg=>odex())
          Dict(:alg=>ddeabm())
          Dict(:alg=>Vern7())
          Dict(:alg=>Vern8())
          Dict(:alg=>CVODE_Adams())
          #Dict(:alg=>lsoda())
          Dict(:alg=>Vern6())
          Dict(:alg=>ARKODE(Sundials.Explicit(),order=6))
          ]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,maxiters=1000,numruns=100)
plot(wp)

setups = [Dict(:alg=>DP8())
          #Dict(:alg=>ode78())
          Dict(:alg=>Vern7())
          Dict(:alg=>Vern8())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,dense=true,maxiters=1000,error_estimate=:L2,numruns=100)
plot(wp)

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
# A semi-separable matrix is defined by
#
# S = triu(uvᵀ, u) + tril(pqᵀ, -bl),
#
# where u, bl >= 0. u and v are n × ru matrices and p and q are n × rl matrices.
# See [1] for the definition.

# [1] Chandrasekaran and Gu, Fast and stable algorithms for banded plus
# Semiseparable systems of linear equations, SIMAX, 25 (2003), pp. 373-384.

## Constructors

struct SemiseparableMatrix{T,LL<:AbstractMatrix{T},UU<:AbstractMatrix{T}} <: LayoutMatrix{T}
    L::LL
    U::UU
    l::Int
    u::Int

    function SemiseparableMatrix{T,LL,UU}(L::LL, U::UU, l::Integer, u::Integer) where {T,LL<:AbstractMatrix{T},UU<:AbstractMatrix{T}}
        Lm, Ln = size(L)
        Um, Un = size(U)
        Um == Un == Lm == Ln || throw(DimensionMismatch())
        (Un >= u && Lm >= l) || throw(ArgumentError("l and u should be compatible with dimensions"))
        -l < u || throw(ArgumentError("bands must not overlap"))
        new{T,LL,UU}(L, U, l, u)
    end
end

SemiseparableMatrix(A::LL, B::UU, l::Integer, u::Integer) where {T,LL<:AbstractMatrix{T},UU<:AbstractMatrix{T}} =
    SemiseparableMatrix{T,LL,UU}(A, B, l, u)

function convert(::Type{Matrix}, S::SemiseparableMatrix)
    return triu(Matrix(S.U), S.u) + tril(Matrix(S.L), -S.l)
end

size(S::SemiseparableMatrix) = size(S.L)

function getindex(S::SemiseparableMatrix{T}, k::Int, j::Int)  where T
    k-j ≥ S.l && return S.L[k,j]
    j-k ≥ S.u && return S.U[k,j]
    return zero(T)
end

# triangular

colsupport(L::SemiseparableMatrix{<:Any,<:Any,<:Zeros}, j) = minimum(j)+L.l:size(L,1)
colsupport(U::SemiseparableMatrix{<:Any,<:Zeros}, j) = 1:maximum(j)-U.u
rowsupport(L::SemiseparableMatrix{<:Any,<:Any,<:Zeros}, k) = 1:maximum(k)-L.l
rowsupport(U::SemiseparableMatrix{<:Any,<:Zeros}, k) = minimum(k)+U.u:size(U,2)


text:
function findClosestValue(function_handler, inputRange, value, conf)

	lb = 1;
	ub = size(inputRange, 1);

	iter = 1;
	closest = 0;

	#%if feval(function_handler, conf, inputRange(ub,:)) < value || feval(function_handler, conf, inputRange(lb,:)) > value
	if function_handler == "cfFlushRateApprox"
		if cfFlushRateApprox(conf, inputRange[ub,:])[1] < value
			closest = Inf;
			return closest
		end
		if cfFlushRateApprox(conf, inputRange[lb,:])[1] > value
			closest = -Inf;
			return closest
		end

		while ub - lb > 100
			closest = round((ub+lb)/2);
			#%fprintf(1,'Checking %d\n', closest);
			iter = iter + 1;
			cv = cfFlushRateApprox(conf, inputRange[closest,:]);
			if cv[1] < value
				lb = closest;
			else
				if cv[1] > value
					ub = closest;
				else # equals!
					lb = closest;
					ub = closest;
				end
			end
		end
	elseif function_handler == "getConcurrencyLevel"
		if getConcurrencyLevel(conf, inputRange[ub,:])[1] < value
			closest = Inf;
			return closest
		end
		if getConcurrencyLevel(conf, inputRange[lb,:])[1] > value
			closest = -Inf;
			return closest
		end

		while ub - lb > 100
			closest = round((ub+lb)/2);
			#%fprintf(1,'Checking %d\n', closest);
			iter = iter + 1;
			cv = getConcurrencyLevel(conf, inputRange[closest,:]);
			if cv[1] < value
				lb = closest;
			else
				if cv[1] > value
					ub = closest;
				else # equals!
					lb = closest;
					ub = closest;
				end
			end
		end
	end
	closest = round((ub+lb)/2);

	println(iter)
	return closest
end #end of function

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
###
### Compute Wigner 3-j symbols recursively, as specified in paper:
###     "Revisiting the Fourier expansion of Mie scattering matrices
###      in generalized spherical functions", Suniti Sanghavi, 2014
###
### Note: Minor typos from the paper are corrected in this implementation:
###     1. Eqn 25:
###       -   sqrt(2*m+3)
###       +   (2*m+3)
###
###     2. Eqn 31:
###       -   sqrt((l-1) * l * (l+1) * (l+2))
###       +   ((l-1) * l * (l+1) * (l+2)) ^ (-0.5)
###

"""
    $(FUNCTIONNAME)(m::Integer, n::Integer, l::Integer, wigner_A::Array{Float64, 3}, wigner_B::Array{Float64, 3})

Compute the wigner 3j symbol for provided m, n, l, (aka j₁, j₂, j₃) where m₁ = -1, m₂ = 1, and m₃ = 0
"""
function wigner3j_m110!(m::Integer, n::Integer, l::Integer,
                      wigner_A::Array{Float64, 3}, wigner_B::Array{Float64, 3})

    # Only compute the value if this condition is met
    if ((m>=abs(n-l))&&(l>=abs(n-m))&&(n>=abs(m-l)))

        if (!(m<=size(wigner_A,1) && n<=size(wigner_A,2) && l<=size(wigner_A,3)))
            return 0.0
        end

        # Only perform this computation if it has not been calculated yet.
        # In these cases, always store the result before returning
        # (Memoization)
        if (wigner_A[m,n,l+1] === NaN)

            # Base Case 1: m is l+n+1, in which case the result is 0
            if m == (l+n+1)

                wigner_A[m,n,l+1] = 0.0
                return 0.0

            # Base Case 2: m is l+n
            elseif (m == n+l)

                # If n is 1, return fixed expression (Eqn. 28)
                if (n == 1)
                    result = (-1)^l * sqrt(((l+1) * (l+2))/((2l+1) * (2l+2) * (2l+3)))
                    wigner_A[m,n,l+1] = result
                    return result

                # If n is not 1, return recursive expression (recursing on n, Eqn. 27)
                else
                    result = -wigner3j_m110!(l+n-1, n-1, l, wigner_A, wigner_B) *
                              sqrt((n * (2n-1) * ((n+l)^2-1))/ ((n+l) * (2(n+l)+1) * (n^2-1)))
                              wigner_A[m,n,l+1] = result
                    return result
                end

            # Main Recursive Case
            else

                # Define variables (Eqn. 26)
                k = m + 1
                D_k_nl = (1/k) * sqrt((k^2 - 1) * (k^2 - (l-n)^2) * ((n+l+1)^2 - k^2))
                M_k = 1 - (n * (n+1) - l * (l+1))/(k * (k+1))
                kP1 = k+1
                D_kP1_nl = (1/kP1) * sqrt((kP1^2 - 1) * (kP1^2 - (l-n)^2) * ((n+l+1)^2 - kP1^2))

                # Perform recursion (Eqn. 25)
                # NOTE: Typo found in paper's equation. Corrected here.
                # Incorrect: sqrt(2*m+3)
                # Correct: 2*m+3
                result = (1/D_k_nl) *
                         (M_k * (2*m+3) * wigner3j_m110!(m+1, n, l, wigner_A, wigner_B)
                             - D_kP1_nl * wigner3j_m110!(m+2, n, l, wigner_A, wigner_B))

                wigner_A[m,n,l+1] = result
                return result
            end

        # This value has already been calculated, so just return it
        else
            return wigner_A[m,n,l+1]
        end

    # Condition is not met: return 0
    else
        return 0.0
    end
end

"""
    $(FUNCTIONNAME)(m::Integer, n::Integer, l::Integer, wigner_A::Array{Float64, 3}, wigner_B::Array{Float64, 3})

Compute the wigner 3j symbol for provided m, n, l, (aka j₁, j₂, j₃) where m₁ = 0, m₂ = 0, and m₃ = 0
"""
function wigner3j_000!(m::Integer, n::Integer, l::Integer,
                     wigner_A::Array{Float64, 3}, wigner_B::Array{Float64, 3})

    # If m = n+l+1, return 0 (does not satisfy condition)
    if (m == n + l + 1)
        return 0.0
    # If m = n+l, return the result specified (Eqn. 30)
    elseif (m == n + l)
        return wigner3j_m110!(n+l, n, l, wigner_A, wigner_B) * 2 *
               sqrt((n+l) * (n+l+1) * n * (n+1))/(l*(l+1) - (l+n)*(l+n+1) - n*(n+1))
    # If m < n+l,
    else

        # If m+n+l is even, return the result specified (Eqn. 29)
        if (iseven(m + n + l == 0))
            return -wigner3j_000!(m+2, n, l, wigner_A, wigner_B) *
                    sqrt(((m+2)^2 - (n-l)^2 ) / ((m+1)^2 - (n-l)^2 )) *
                    sqrt((1-(1/(n+l-m))) * (1+(1/(m+n+l+2))))

        # If m+n+l is odd, return 0
        else
            return 0.0
        end

    end

end

"""
    $(FUNCTIONNAME)(m::Integer, n::Integer, l::Integer, wigner_A::Array{Float64, 3}, wigner_B::Array{Float64, 3})

Compute the wigner 3j symbol for provided m, n, l, (aka j₁, j₂, j₃) where m₁ = -1, m₂ = -1, and m₃ = 2
"""
function wigner_m1m12!(m::Integer, n::Integer, l::Integer,
                       wigner_A::Array{Float64, 3}, wigner_B::Array{Float64, 3})

    # Simply save and return the result specified in Eqn. 31

    result = (-1)^(m+n+l) * (((l-1) * l * (l+1) * (l+2))^(-0.5)) *
                ( (m*(m+1) + (-1)^(m+n+l)*n*(n+1)) *
                wigner3j_m110!(m, n, l, wigner_A, wigner_B) +
                2 * sqrt(m*(m+1)*n*(n+1)) *
                wigner3j_000!(m, n, l, wigner_A, wigner_B) )

    wigner_B[m, n, l+1] = result

    return result
end

"""
    $(FUNCTIONNAME)(j₁::Integer, j₂::Integer, j₃::Integer, m₁::Integer, m₂::Integer, m₃::Integer, wigner_A::Array{Float64,3}, wigner_B::Array{Float64,3})

Compute the Wigner 3j values for all n to nmax and all m from mstart to (n+l+1), using
wigner_A and wigner_B to store computed values for (-1, 1, 0) and (-1, -1, 2) respectively.

"""
function wigner3j!(j₁::Integer, j₂::Integer, j₃::Integer,
                 m₁::Integer, m₂::Integer, m₃::Integer,
                 wigner_A::Array{Float64,3}, wigner_B::Array{Float64,3})

    # Perform a check that the input j's satisfy condition
    if ((j₁<abs(j₂-j₃))||(j₃<abs(j₂-j₁))||(j₂<abs(j₁-j₃))||m₁>j₁||m₁<-j₁||m₂>j₂||m₂<-j₂||m₃>j₃||m₃<-j₃)
        return 0.0
    end

    # Switch between the three different configurations and call appropriate function
    # (m₁, m₂, m₃): (-1, 1, 0), (-1, -1, 2), (0, 0, 0)
    if (m₁ == -1 && m₂ == 1 && m₃ == 0)
        return wigner3j_m110!(j₁, j₂, j₃, wigner_A, wigner_B)

    elseif (m₁ == -1 && m₂ == -1 && m₃ == 2)
        return wigner_m1m12!(j₁, j₂, j₃, wigner_A, wigner_B)

    elseif (m₁ == 0 && m₂ == 0 && m₃ == 0)
        return wigner3j_000!(j₁, j₂, j₃, wigner_A, wigner_B)

    else
        throw(DomainError((m₁, m₂, m₃,),
            "Unsupported configuration. Supported: (-1,1,0), (-1,-1,2), (0,0,0)"))
    end

end

"""
    $(FUNCTIONNAME)(m_max::Integer, n_max::Integer, l_max::Integer)

Compute the Wigner 3j values for all (m, n, l) combinations up to m/n/l max, for
(m1, m2, m3) = (-1, 1, 0) (wigner_A) and (-1, -1, 2) (wigner_B)

"""
function compute_wigner_values(m_max::Integer, n_max::Integer, l_max::Integer)

    # Create the m, n, l ranges
    ms = 1:m_max
    ns = 1:n_max
    ls = 1:l_max

    # Create empty matrices to
    wigner_A = zeros(m_max, n_max, l_max);
    wigner_B = zeros(m_max, n_max, l_max);
    fill!(wigner_A, NaN)
    fill!(wigner_B, NaN)

    # Using nested loop syntax is so nice!!
    @showprogress 1 "Computing Wigner Symbols... " for l in ls, n in ns, m in ms
                wigner_A[m, n, l] = wigner3j!(m, n, l-1, -1, 1, 0, wigner_A, wigner_B)
                wigner_B[m, n, l] = wigner3j!(m, n, l-1, -1, -1, 2, wigner_A, wigner_B)
    end

    return wigner_A, wigner_B
end

"""
    $(FUNCTIONNAME)(N_max)

Shorthand for compute_wigner_values(2N_max + 1, N_max + 1, 2N_max + 1)
"""
compute_wigner_values(N_max) = compute_wigner_values(2N_max + 1, N_max + 1, 2N_max + 1)

"""
    $(FUNCTIONNAME)(filepath, wigner_A, wigner_B)

Save the Wigner A and Wigner B matrices at the given filepath
"""
function save_wigner_values(filepath, wigner_A, wigner_B)
    @save filepath wigner_A wigner_B
end

"""
    $(FUNCTIONNAME)(filename)

Load the Wigner A and Wigner B matrices from the given filepath
"""
function load_wigner_values(filename)
    @load filename wigner_A wigner_B
    return wigner_A, wigner_B
end


text:
include("../src/FiberNlse.jl")
using Plots;gr()
using FFTW;
using Peaks;
using SpecialFunctions;
U = FiberNlse; # Alias

begin

# CONSTANTS
    T = 50000*U.ps;
    Fs = 200*U.GHz;
# Simulation dimension
    Nₜ, Nₗ = (Int(round(Fs*T)),500);

# Fiber properties

    D = -7*U.ps/U.nm/U.km;
    α = 0*0.046/U.km;
    γ =  1.1/U.W/U.km;
    L = 4131;
    fib = U.Fiber(L, α, D, γ)

# Signal properties
λ1 = 1550.0*U.nm; # Wavelength
λ2 = λ1+1*U.pm; # Wavelength
λ0 = 0.5*(λ1+λ2);
Δλ = 0.5*abs(λ1-λ2)
sim, t = U.configure(Nₜ,Nₗ,D, γ, α, L, T, λ0);
Δw = 2pi*U.c*abs(λ1-λ2)/(λ1*λ2)
end

SP(ψ) = fftshift(abs.(fft(ψ)).^2)/sim.Nₜ^2
ν = fftshift(FFTW.fftfreq(sim.Nₜ, 1.0 / (sim.dt)));

function run(P)
    # Input construction
    Ψₒ = @. sqrt(0.5*P)*(exp.(1im*Δw.*t).+exp.(-1im*Δw.*t))
    U.inputSignal(sim,Ψₒ);
    U.simulate(sim, false);
    return (sim.Ψ[1,:],sim.Ψ[end,:])
end

ψi,ψf = run(0.5);
sp = SP(ψf)
idx = Peaks.argmaxima(sp)
idx = idx[findall(ν[idx] .< 4*Δw/(2pi))]
idx = idx[findall(ν[idx] .> 0)]
plot(SP(ψf))
scatter!(idx, sp[idx])

begin
Pl = (30:10:300)*U.mW
r = []
pf = []
for p in Pl
ψi,ψf = run(p);
sp = SP(ψf)
idx = Peaks.argmaxima(sp)
idx = idx[findall(ν[idx] .< 4*Δw/(2pi))]
idx = idx[findall(ν[idx] .> 0)]

if length(idx)>1
    append!(r, 10*log10(sp[idx[1]]/sp[idx[2]]))
    append!(pf, p)
end
end
end
pth = range(30, 300, 1000)*U.mW
a=0.82pi/maximum(pth)

rr(x) = @. (besselj0(0.5x)^2+besselj1(0.5x)^2)/(besselj1(0.5x)^2+besseljx(2,0.5x)^2)
scatter(pf,r)
plot!(pth, 10*log10.(rr.((a*pth))))

γexp = 0.5*a/L
δγ = abs(γexp-γ)/γ*100


text:
using Base.Test
using BlackScholes

tol = 1.0e-6
@test_approx_eq_eps(BSCall(10.0,10.0,0.15,0.0323,7.0,0.02),1.699938,tol)
@test_approx_eq_eps(BSPut(10.0,10.0,0.15,0.0323,7.0,0.02),0.982739,tol)

@test BSImpliedVol(10.0,10.0,0.0323,7.0,0.02,BSCall(10.0,10.0,0.55,0.0323,7.0,0.02),true) == 0.55
@test BSImpliedVol(10.0,10.0,0.0323,7.0,0.02,BSPut(10.0,10.0,0.55,0.0323,7.0,0.02),false) == 0.55

@test BSImpliedVol(10.0,10.0,0.0323,7.0,0.02,BSCall(10.0,10.0,0.55,0.0323,7.0,0.02),true) == 0.55


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestConstraintFlipSign

using Test

using MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_GreaterToLess()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.GreaterToLess{Float64}(mock)
    MOI.Test.runtests(
        bridged_mock,
        config,
        include = [
            "test_basic_VariableIndex_GreaterThan",
            "test_basic_ScalarAffineFunction_GreaterThan",
            "test_basic_ScalarQuadraticFunction_GreaterThan",
        ],
    )
    MOI.empty!(bridged_mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [0.0, 0.0]),
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [100.0, 0.0]),
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [100.0, -100.0]),
    )
    MOI.Test.test_linear_LessThan_and_GreaterThan(bridged_mock, config)
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.ScalarAffineFunction{Float64},
                MOI.GreaterThan{Float64},
            }(),
        ),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        MOI.set(bridged_mock, attr, ci, 2.0)
        @test MOI.get(bridged_mock, attr, ci) ≈ 2.0
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        2,
        ((MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 1),),
    )
    return
end

function test_LessToGreater()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.LessToGreater{Float64}(mock)
    MOI.Test.runtests(
        bridged_mock,
        config,
        include = [
            "test_basic_VariableIndex_LessThan",
            "test_basic_ScalarAffineFunction_LessThan",
            "test_basic_ScalarQuadraticFunction_LessThan",
        ],
    )
    MOI.empty!(bridged_mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [1.0],
        ),
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [2.0]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [1.0],
        ),
    )
    MOI.Test.test_modification_set_scalaraffine_lessthan(bridged_mock, config)
    MOI.empty!(bridged_mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [1.0],
        ),
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [0.5]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [0.5],
        ),
    )
    MOI.Test.test_modification_coef_scalaraffine_lessthan(bridged_mock, config)
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.ScalarAffineFunction{Float64},
                MOI.LessThan{Float64},
            }(),
        ),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        MOI.set(bridged_mock, attr, ci, 2.0)
        @test MOI.get(bridged_mock, attr, ci) ≈ 2.0
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        1,
        ((MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),),
    )
    return
end

function test_NonnegToNonpos()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.NonnegToNonpos{Float64}(mock)
    MOI.Test.runtests(
        bridged_mock,
        config,
        include = [
            "test_basic_VectorOfVariables_Nonnegatives",
            "test_basic_VectorAffineFunction_Nonnegatives",
            "test_basic_VectorQuadraticFunction_Nonnegatives",
        ],
    )
    MOI.empty!(bridged_mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [0.0, 0.0]),
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [100.0, 0.0]),
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [100.0, -100.0]),
    )
    MOI.Test.test_linear_VectorAffineFunction(bridged_mock, config)
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.Nonnegatives,
            }(),
        ),
    )
    _test_delete_bridge(
        bridged_mock,
        ci,
        2,
        ((MOI.VectorAffineFunction{Float64}, MOI.Nonpositives, 1),),
    )
    MOI.empty!(bridged_mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1.0, 0.0, 2.0],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonpositives) =>
                [[0.0, -2.0, 0.0]],
            (MOI.VectorAffineFunction{Float64}, MOI.Zeros) =>
                [[-3.0, -1.0]],
        )
    MOI.Test.test_conic_linear_VectorOfVariables(bridged_mock, config)
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorOfVariables,
                MOI.Nonnegatives,
            }(),
        ),
    )
    func = MOI.get(bridged_mock, MOI.ConstraintFunction(), ci)
    MOI.delete(bridged_mock, func.variables[2])
    new_func = MOI.VectorOfVariables(func.variables[[1, 3]])
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) == new_func
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.Nonnegatives(2)
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        MOI.set(bridged_mock, attr, ci, [1.0, 2.0])
        @test MOI.get(bridged_mock, attr, ci) ≈ [1.0, 2.0]
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        2,
        ((MOI.VectorAffineFunction{Float64}, MOI.Nonpositives, 0),),
    )
    return
end

function test_NonposToNonneg()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.NonposToNonneg{Float64}(mock)
    MOI.Test.runtests(
        bridged_mock,
        config,
        include = [
            "test_basic_VectorOfVariables_Nonpositives",
            "test_basic_VectorAffineFunction_Nonpositives",
            "test_basic_VectorQuadraticFunction_Nonpositives",
        ],
    )
    MOI.empty!(bridged_mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [0.0, 0.0]),
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [100.0, 0.0]),
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [100.0, -100.0]),
    )
    MOI.Test.test_linear_VectorAffineFunction(bridged_mock, config)
    MOI.empty!(bridged_mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [0.0, 0.0]),
        ),
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0, 0.75]),
        ),
    )
    MOI.Test.test_modification_const_vectoraffine_nonpos(bridged_mock, config)
    MOI.empty!(bridged_mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [0.5]),
        ),
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [0.25]),
        ),
    )
    MOI.Test.test_modification_multirow_vectoraffine_nonpos(
        bridged_mock,
        config,
    )
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.Nonpositives,
            }(),
        ),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        MOI.set(bridged_mock, attr, ci, [1.0, 2.0])
        @test MOI.get(bridged_mock, attr, ci) ≈ [1.0, 2.0]
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        1,
        ((MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0),),
    )
    return
end

end  # module

TestConstraintFlipSign.runtests()

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
@testset "CoreTracker" begin
    @testset "General" begin
        F = equations(katsura(5))
        # test construction
        t1, start_sols = coretracker_startsolutions(F)

        test_show_juno(t1)
        @test !isempty(string(t1))
        test_show_juno(t1.options)
        @test !isempty(string(t1.options))
        test_show_juno(t1.state)
        @test !isempty(string(t1.state))

        @test_nowarn current_Δt(t1)

        raccuracy = refinement_accuracy(t1)
        @test_nowarn set_refinement_accuracy!(t1, 5e-5)
        @test refinement_accuracy(t1) == 5e-5
        set_refinement_accuracy!(t1, raccuracy)

        caccuracy = accuracy(t1)
        @test_nowarn set_accuracy!(t1, 5e-5)
        @test accuracy(t1) == 5e-5
        set_accuracy!(t1, caccuracy)

        cmaxiter = max_corrector_iters(t1)
        @test_nowarn set_max_corrector_iters!(t1, 11)
        @test max_corrector_iters(t1) == 11
        set_max_corrector_iters!(t1, cmaxiter)

        rmaxiter = max_refinement_iters(t1)
        @test_nowarn set_max_refinement_iters!(t1, 11)
        @test max_refinement_iters(t1) == 11
        set_max_refinement_iters!(t1, rmaxiter)

        @test t1 isa CoreTracker
        @test isa(current_x(t1), Vector)
        @test length(current_x(t1)) == 6

        t3, start_sols = coretracker_startsolutions(F, predictor=Euler())
        @test t3.predictor isa Euler

        setup!(t1, first(start_sols), 1.0, 0.4)
        @test status(t1) == CoreTrackerStatus.tracking
        @test current_t(t1) == 1.0

        setup!(t1, first(start_sols), 0.5, 0.4)
        @test status(t1) == CoreTrackerStatus.terminated_invalid_startvalue
        @test current_t(t1) == 0.5

        R = track(t1, first(start_sols), 1.0, 0.0)
        @test R isa CoreTrackerResult
        @test R.returncode == CoreTrackerStatus.success
        @test R.accuracy < 1e-7
        @test_nowarn show(devnull, R)

        out = copy(first(start_sols))
        retcode = track!(out, t1, first(start_sols), 1.0, 0.0)
        @test retcode == CoreTrackerStatus.success
        @test out == R.x
    end

    @testset "Affine tracking" begin
        @polyvar x a y b
        F = [x^2-a, x*y-a+b]
        p = [a, b]

        tracker, starts = coretracker_startsolutions(F, [1.0, 1.0 + 0.0*im], parameters=p, p₁=[1, 0], p₀=[2, 4],
                    affine_tracking=true)
        @test affine_tracking(tracker) == true
        res = track(tracker, starts[1], 1.0, 0.0)
        @test res.returncode == CoreTrackerStatus.success
        @test isa(res.x, Vector{ComplexF64})
        @test length(res.x) == 2

        tracker, starts = coretracker_startsolutions(F, [1.0, 1.0 + 0.0*im], parameters=p, p₁=[1, 0], p₀=[2, 4],
                    affine_tracking=true, auto_scaling=false)
        @test affine_tracking(tracker) == true
        res = track(tracker, starts[1], 1.0, 0.0)
        @test res.returncode == CoreTrackerStatus.success
        @test isa(res.x, Vector{ComplexF64})
        @test length(res.x) == 2

        x = @SVector [1.0, 1.0 + 0.0*im]
        tracker, starts = coretracker_startsolutions(F, x, parameters=p, p₁=[1, 0], p₀=[2, 4], affine_tracking=true)
        @test length(starts) == 1
        res = track(tracker, starts[1], 1.0, 0.0)
        @test isa(res.x, Vector{ComplexF64})
        @test length(res.x) == 2
    end

    @testset "Allocations" begin
        F = equations(cyclic(5))
        tracker, start_sols = coretracker_startsolutions(F; seed=12356)
        s = first(start_sols)
        @allocated track!(tracker, s, 1.0, 0.01)
        # The path tracker not using a rank-deficient QR
        # should not allocate after the first tracked path
        @test (@allocated track!(tracker, s, 1.0, 0.01)) == 0
    end

    @testset "LinearSystem" begin
        A = rand(3, 3)
        b = rand(3)

        @polyvar x[1:3]
        F = A * x - b

        tracker, start_sols = coretracker_startsolutions(F; max_steps=10)
        s = first(start_sols)
        result = track(tracker, s, 1.0, 0.0)
        @test result.returncode == CoreTrackerStatus.success
        @test result.t == 0.0
        x = result.x
        @test norm(x - A \ b) < 1e-6

        x_inter = copy(s)
        retcode = track!(x_inter, tracker, s, 1.0, 0.1)
        @test retcode == CoreTrackerStatus.success
        x_final = zero(x_inter)
        retcode = track!(x_final, tracker, x_inter, 0.1, 0.0)
        @test retcode == CoreTrackerStatus.success
        @test iters(tracker) < 3
        x = current_x(tracker)
        @test norm(x - A \ b) < 1e-6
    end

    @testset "FixedPatch continuation" begin
        F = homogenize(equations(katsura(5)))
        # test construction
        tracker, start_sols = coretracker_startsolutions(F, patch=OrthogonalPatch())
        patch = tracker.state.patch
        fixedtracker = coretracker(F, patch=FixedPatch())
        fixedpatch = fixedtracker.state.patch
        x1 = first(start_sols)

        r1 = track(fixedtracker, x1, 1.0, 0.1)
        v1 = copy(fixedpatch.v̄)
        track!(fixedtracker, r1.x, 0.1, 0.05, setup_patch=false)
        @test v1 == fixedpatch.v̄
        track!(fixedtracker, current_x(tracker), 0.05, 0.01, setup_patch=false)
        @test v1 == fixedpatch.v̄
        track!(fixedtracker, r1.x, 0.1, 0.05, setup_patch=true)
        @test v1 != fixedpatch.v̄

        r1 = track(tracker, x1, 1.0, 0.1)
        v1 = copy(patch.v̄)
        track!(tracker, r1.x, 0.1, 0.0)
        @test v1 != patch.v̄
    end

    @testset "coretracker_startsolutions" begin
        @polyvar x y
        tracker, S = coretracker_startsolutions([x^2+y^2-4, x + y - 3])
        result = track(tracker, first(S), 1.0, 0.0)
        @test result.returncode == CoreTrackerStatus.success

        # test type promotion of startsolutions
        tracker, S = coretracker_startsolutions([x^2+y^2-4, x + y - 3])
        result = track(tracker, [1, 1], 1.0, 0.0)
        @test result.returncode == CoreTrackerStatus.success
    end

    @testset "iterator" begin
        @polyvar x a y b
        F = [x^2-a, x*y-a+b]
        p = [a, b]
        s = [1.0, 1.0 + 0.0*im]
        tracker = coretracker(F, parameters=p, p₁=[1, 0], p₀=[2, 4])

        typeof(first(iterator(tracker, s, 1.0, 0.0))) ==
            Tuple{Vector{ComplexF64}, Float64}

        @test max_step_size(tracker) == Inf
        set_max_step_size!(tracker, 0.01)
        @test max_step_size(tracker) == 0.01

        length(collect(iterator(tracker, s, 1.0, 0.0))) == 101

        # Test iteration protocol
        tracker = coretracker(F, parameters=p, p₁=[1, 0], p₀=[2, 4])
        setup!(tracker, s, 1.0, 0.0)
        for tr in tracker
            # nothing
        end
        @test tracker.state.status == CoreTrackerStatus.success
    end

    @testset "precision option" begin
        @polyvar  x y
        for P in [PRECISION_ADAPTIVE, PRECISION_FIXED_64, PRECISION_FIXED_128]
            tracker, S = coretracker_startsolutions([x^2+y^2-4, x + y - 3]; precision=P)
            @test tracker.options.precision == P
            result = track(tracker, first(S), 1.0, 0.0)
            @test result.returncode == CoreTrackerStatus.success
        end
        # test default
        tracker, S = coretracker_startsolutions([x^2-y^2+4, x + y - 3])
        @test tracker.options.precision == PRECISION_FIXED_64
    end

    @testset "cond updates" begin
        @polyvar x y
        tracker, S = coretracker_startsolutions([x^2+y^2-4, x + y - 3]; affine_tracking=false, seed=130793)
        @test cond(tracker) == 1.0
        @test digits_lost(tracker) == 0.0
        # check that checkstartvalue also updates cond and digits_lost
        setup!(tracker, first(S))
        cond(tracker)
        digits_lost(tracker)
        cond_start, digits_lost_start = cond(tracker), digits_lost(tracker)
        @test cond_start != 1.0
        @test digits_lost_start != 0.0
        track(tracker, first(S), 1.0, 0.0)
        @test cond(tracker) != cond_start
        @test digits_lost(tracker) != digits_lost_start
    end

    @testset "path iterator" begin
        @polyvar x a
        f = [x - a]
        ct = coretracker(f, [1.0], parameters=[a], p₁=[1.],
                p₀=[2.], max_step_size=0.1)
        Xs = Vector{ComplexF64}[]
        for (x, t) in iterator(ct, [1.0], 1.0, 0.0)
            push!(Xs, x)
        end
        @test round.(Int, real.(first.(Xs)) .* 10) == collect(10:20)
    end
end


text:
using StaticArrays
using Base.Test

@testset "StaticArrays" begin
    include("SVector.jl")
    include("MVector.jl")
    include("SMatrix.jl")
    include("MMatrix.jl")
    include("SArray.jl")
    include("MArray.jl")
    include("FieldVector.jl")

    include("core.jl")
    include("abstractarray.jl")
    include("indexing.jl")
    include("mapreduce.jl")
    include("arraymath.jl")
    include("linalg.jl")
    include("matrix_multiply.jl")
    include("deque.jl")
end


text:
#######################################################################
function noisyintersectionML(; posterior = posterior, g = g, σ2 = 1e-6*ones(length(g)))
#######################################################################

    # work with normalised vector

    gunit = g / norm(g)

    linemeansigma = marginaline(posterior, σ2)

    function loss(p)

        local μ, Σ = linemeansigma(p[2])

        -1.0 * logpdf(MvNormal(μ, Σ), p[1]*gunit)

    end

    # get parameters of mean line

    D = length(g)

    a = mean(posterior)[1:D]

    b = mean(posterior)[D+1:2*D]

    # find point on g that minimises distance to mean line

    aux1(s) = norm(s*gunit - b - proj(a, s*gunit - b))^2

    s0 = optimize(aux1, -1000.0, 1000.0).minimizer

    # find point on mean line that minimises distance to point on g we just found

    aux2(t) = norm(s0*gunit - (a*t + b))^2

    t0 = optimize(aux2, -1000.0, 1000.0).minimizer

    # find s, t starting from s0, t0 that minimise above log-likelihood

    opt    = Optim.Options(show_trace = false, iterations = 1_000_000)

    result = Optim.optimize(loss, [s0; t0], NelderMead(), opt).minimizer

end

"""
    randx = noisyintersectionvi(; posterior = posterior, g = g, [σ2 = 1e-6*ones(length(g))])

`posterior` is the posterior distribution returned by bpca, of type ```Distributions.MvNormal```.

`g` ∈ ℜⁿ₊ is the candidate direction vector that is being tested, of type ```Array{Float64, 1}```.

The function returns a function `randx()` that draws a sample intersection point.

"""
function noisyintersectionvi(; posterior = posterior, g = g, σ2 = 1e-6*ones(length(g)))

    @info("Scale parameter re-parametrised to enforce positivity")

    # check dimensions of arguments

    @assert(length(g) == round(Int, length(posterior)/2))

    # check that components of vector are positive

    @assert(all(g .> 0.0))

    # work with normalised vector

    gunit = g / norm(g)

    # instantiate posterior of lines

    linemeansigma = marginaline(posterior, σ2)

    # auxiliary functions for making positive parameters

    makepos(x)    = x^2

    invmakepos(x) = sqrt(x)

    # Initialise with ML estimate

    bestsolutionML = noisyintersectionML(; posterior = posterior, g = gunit, σ2 = σ2)

    # log-likleihood function for VI

    function logp(p)

        local μ, Σ = linemeansigma(p[2])

        logpdf(MvNormal(μ, Σ), makepos(p[1])*gunit)

    end

    # Get ML estimate as initial solution

    @printf("ML solution for scale is %f\n", bestsolutionML[1])

    bestsolutionML[1] = invmakepos(max(1e-3, bestsolutionML[1]))

    # Call variational inference

    postp, logevidence = VI(x -> logp(x), [bestsolutionML.+(i/10)*(randn(length(bestsolutionML))) for i=1:30], optimiser = NelderMead(), S = 500, inititerations = 50, iterations = 10_000, show_every = 25)

    # Retrieve results from VI and instantiate sampling function

    pdfscale = Normal(mean(postp)[1], sqrt(cov(postp)[1,1]))

    randscale = () -> makepos(rand(pdfscale))

    # xavg = mean([randscale() for i=1:10_000]) * gunit

    randx = () -> randscale() * gunit

    @printf("Log-evidence is %f (higher is better)\n", logevidence)

    return randx

end


text:
# 3D Object Recognition based on Correspondence Grouping
# http://www.pointclouds.org/documentation/tutorials/correspondence_grouping.php#correspondence-grouping

using PCLCommon
using PCLIO
using PCLFilters
using PCLFeatures
using PCLKDTree
using PCLRecognition
using PCLVisualization
using Cxx

import PCLCommon: CorrespondenceVal

model_pcd = Pkg.dir("PCLIO", "test", "data", "milk.pcd")
scene_pcd = Pkg.dir("PCLIO", "test", "data", "milk_cartoon_all_small_clorox.pcd")

T = PointXYZRGBA
NT = Normal
DT = SHOT352
RT = ReferenceFrame

model_ss = 0.01
scene_ss = 0.03
descr_rad = 0.02
rf_rad = 0.015
cg_size = 0.01
cg_thresh = 5.0
use_hough = true

model = PointCloud{T}(model_pcd)
scene = PointCloud{T}(scene_pcd)
model_keypoints = PointCloud{T}()
scene_keypoints = PointCloud{T}()
model_normals = PointCloud{NT}()
scene_normals = PointCloud{NT}()
model_descriptors = PointCloud{DT}()
scene_descriptors = PointCloud{DT}()

info("Normal estimation")
norm_est = NormalEstimationOMP{T,NT}()
setKSearch(norm_est, 10)
setInputCloud(norm_est, model)
compute(norm_est, model_normals)

setInputCloud(norm_est, scene)
compute(norm_est, scene_normals)

# Downsample
info("Uniform sampling")
uniform_sampling = UniformSampling{T}()
setInputCloud(uniform_sampling, model)
PCLFilters.setRadiusSearch(uniform_sampling, model_ss)
filter(uniform_sampling, model_keypoints)
println("Model total points: $(length(model)); Selectd Key points: $(length(model_keypoints))")

setInputCloud(uniform_sampling, scene)
PCLFilters.setRadiusSearch(uniform_sampling, scene_ss)
filter(uniform_sampling, scene_keypoints)
println("Scene total points: $(length(scene)); Selectd Key points: $(length(scene_keypoints))")

info("Associating a 3D descriptor to each model and scene keypoint")
descr_est = SHOTEstimationOMP{T,NT,DT}()
PCLFeatures.setRadiusSearch(descr_est, descr_rad)
setInputCloud(descr_est, model_keypoints)
setInputNormals(descr_est, model_normals)
setSearchSurface(descr_est, model)
compute(descr_est, model_descriptors)

setInputCloud(descr_est, scene_keypoints)
setInputNormals(descr_est, scene_normals)
setSearchSurface(descr_est, scene)
compute(descr_est, scene_descriptors)

info("Find Model-Scene Correspondences with KdTree")
match_search = KdTreeFLANN{DT}()
setInputCloud(match_search, model_descriptors)

model_scene_corrs = Correspondences()
for i in 0:length(scene_descriptors)-1
    neigh_indices = icxx"std::vector<int>(1);"
    neigh_sqr_dists = icxx"std::vector<float>(1);"
    p = icxx"$(scene_descriptors[i]).descriptor[0];"
    !isfinite(Cfloat(p)) && continue

    found_neighs = nearestKSearch(match_search, scene_descriptors[i],
        1, neigh_indices, neigh_sqr_dists)
    if found_neighs == 1 && (Float32(neigh_sqr_dists[0]) < Float32(0.25))
        corr = CorrespondenceVal(neigh_indices[0], i, neigh_sqr_dists[0])
        push!(model_scene_corrs, corr)
    end
end
@show length(model_scene_corrs)

info("Clustering")
rototranslations = icxx"""
    std::vector<Eigen::Matrix4f, Eigen::aligned_allocator<Eigen::Matrix4f>>();
"""
clustered_corrs = icxx"std::vector<pcl::Correspondences>();"

clusterer = use_hough ? Hough3DGrouping{T,T,RT,RT}() :
    GeometricConsistencyGrouping{T,T}()

if use_hough
    model_rf = PointCloud{RT}()
    scene_rf = PointCloud{RT}()
    rf_est = BOARDLocalReferenceFrameEstimation{T,NT,RT}()
    setFindHoles(rf_est, true)
    PCLFeatures.setRadiusSearch(rf_est, rf_rad)

    setInputCloud(rf_est, model_keypoints)
    setInputNormals(rf_est, model_normals)
    setSearchSurface(rf_est, model)
    compute(rf_est, model_rf)

    setInputCloud(rf_est, scene_keypoints)
    setInputNormals(rf_est, scene_normals)
    setSearchSurface(rf_est, scene)
    compute(rf_est, scene_rf)

    setHoughBinSize(clusterer, cg_size)
    setHoughThreshold(clusterer, cg_thresh)
    setUseInterpolation(clusterer, true)
    setUseDistanceWeight(clusterer, false)
else
    setGCSize(clusterer, cg_size)
    setGCThreshold(clusterer, cg_thresh)
end

setInputCloud(clusterer, model_keypoints)
use_hough && setInputRf(clusterer, model_rf)
setSceneCloud(clusterer, scene_keypoints)
use_hough && setSceneRf(clusterer, scene_rf)
setModelSceneCorrespondences(clusterer, model_scene_corrs)
recognize(clusterer, rototranslations, clustered_corrs)

@show length(rototranslations)

if isdefined(:vis) && vis
    info("Prepare PCL visualizer...")
    viewer = PCLVisualizer("pcl visualizer")
    color_handler = PointCloudColorHandlerRGBField(scene)
    addPointCloud(viewer, scene, color_handler, id="scene")

    for i in 0:length(rototranslations)-1
        rotated_model = PointCloud{T}()
        transformPointCloud(model, rotated_model, rototranslations[i])

        rotated_model_color_handler = PointCloudColorHandlerCustom(scene,
            255, 0, 0)
        addPointCloud(viewer, rotated_model, rotated_model_color_handler,
            id="instance $i")
    end
    spin(viewer)
    close(viewer)
end


text:
# Gurobi environment and other supporting facilities

type Env
    ptr_env::Ptr{Void}

    function Env()
        a = Array(Ptr{Void}, 1)
        ret = @grb_ccall(loadenv, Cint, (Ptr{Ptr{Void}}, Ptr{UInt8}),
            a, C_NULL)
        if ret != 0
            if ret == 10009
                error("Invalid Gurobi license")
            else
                error("Failed to create environment (error $ret).")
            end
        end
        env = new(a[1])
        # finalizer(env, free_env)  ## temporary disable: which tends to sometimes caused warnings
        env
    end
end

Base.unsafe_convert(ty::Type{Ptr{Void}}, env::Env) = env.ptr_env::Ptr{Void}

function is_valid(env::Env)
    env.ptr_env != C_NULL
end

function free_env(env::Env)
    if env.ptr_env != C_NULL
        @grb_ccall(freeenv, Void, (Ptr{Void},), env.ptr_env)
        env.ptr_env = C_NULL
    end
end

function get_error_msg(env::Env)
    @assert env.ptr_env != C_NULL
    sz = @grb_ccall(geterrormsg, Ptr{UInt8}, (Ptr{Void},), env.ptr_env)
    unsafe_string(sz)
end

# error

type GurobiError
    code::Int
    msg::String

    function GurobiError(env::Env, code::Integer)
        new(convert(Int, code), get_error_msg(env))
    end
end


text:
function make_my_day(n=day(today()))
    dir = joinpath("scripts", "Day"*string(n))
    code_file = joinpath(dir, "code.jl")
    template_file = joinpath("template", "code.jl")

    try
        mkdir(dir)
    catch e
        throw(e)
    end

    (tmp_path, tmp_f) = mktemp()

    !isfile(code_file) && open(template_file, read=true, write=true) do f
        for (line_num, line) in enumerate(eachline(f, keep=true))
            if line_num==1
                write(tmp_f, "module Day"*string(n)*"\n")
            else
                write(tmp_f, line)
            end
        end
        close(tmp_f)
        mv(tmp_path, code_file)
    end

    n ≤ day(today()) && download_input(n)

    edit(code_file)

    return nothing
end

"""
    run_day(n; test=true, time=true)
    run_day(mod::Module; test=true, time=true)

Run day `n`'s puzzle and print outputs. Solutions will only be run if test cases pass. To
turn off the test check, set the `test` flag to `false`. If the `time` flag is set to `true`,
it will also print the time it took to run the solution. The positional argument can also be
a Day module directly.
"""
run_day(n; kwargs...) = run_day(eval(Meta.parse("Advent2020.Day$(n)")); kwargs...)

function run_day(mod::Module; test=true, time=true)
    @unpack get_inputs, get_solution1, get_solution2 = mod
    @unpack test_input1, test_output1, test_input2, test_output2, data = get_inputs()

    # Do not pass if tests don't check out
    test && @testset "Input tests" begin
        @testset "Part 1" begin @test get_solution1(test_input1) == test_output1 end
        @testset "Part 2" begin @test get_solution2(test_input2) == test_output2 end
    end

    # Print and time solution 1
    println("Solution 1: $(get_solution1(data))")
    time && @time get_solution1(data)

    # Print and time solution 2
    println("Solution 2: $(get_solution2(data))")
    time && @time get_solution2(data)

    return nothing
end


text:
using Z3

# Optimize model
# Min x + y
# s.t x > 3
#     y > 5

# Model
m = Z3Model()
@assert isa(m,Z3Model)

# Variables
x = variable(m)
y = variable(m,:Bin)

# Constraints
constraint(m, x == 5)
implies(m, x > 0, y)

# Objective
objective(m, :Min, x*10 + y*100)

#status = solve(m)


text:
module RadioPropagation
	# using Plots
	using Printf

	c = 299792458; # m/s
	i = im;

	"""
	Calculate the two-way two-ray propagation between two points above a flat plane.
	Returns the propagation factor F. The one-way power propagation factor is |F|².

	- Radar Systems Engineering Lecture 5 Propagation through the Atmosphere, IEEE New Hampshire Section IEEE AES Society, 2010

	```julia
	# Example code plotting signal strength as function of range, compared to free space power spread
	using Plots

	c = 3e8;

	transmit_height_m	= 2;
	receive_height_m	= 10;
	target_range_m		= LinRange(20, 100, 10000);
	transmit_frequency_hz =3e9;
	λ = c/transmit_frequency_hz;
	Γ1 = -1;
	Γ2 = -0.6;

	F²( Γ ) = RadioPropagation.two_ray_propagation.( target_range_m, transmit_height_m, receive_height_m, transmit_frequency_hz, Γ );

	F²1 = F²( Γ1 );
	F²2 = F²( Γ2 );

	f⁴_db( F² ) = 10*log10.( abs2.( F² ));

	freespace_loss_db = 10 .*log10.(λ^2 ./( (4*π)^3 .*target_range_m .^4 ));
	bias = freespace_loss_db[1]

	plot1 = plot( 	target_range_m, f⁴_db.(F²1)+freespace_loss_db .-bias,
			 	 	xlabel 	= "Range (one-way) [m]",
					ylabel 	= "Two-way loss [dB]",
					title  	= "Free-space versus two-ray propagation.",
					ylims 	= (-30, 6),
					label 	= "|F|⁴+L, Γ=-1",
					legend	= true,
					xaxis	=:log,
					dpi=300)
	plot!( 	target_range_m, f⁴_db.(F²2)+freespace_loss_db .-bias,
	 		label = "|F|⁴+L, Γ=-0.6",
	 		xaxis	=:log )
	plot!( 	target_range_m, freespace_loss_db .-bias,
	 		label = "L",
	 		xaxis	=:log )
	savefig("figures/example_figure")

	```

	# Arguments
	- 'distance'            The distance parallel to the plane earth.
	- 'transmit_height'     The height of the transmitter above the plane earth.
	- 'receive_height'      The height of the receiver/target above the plane earth.
	- 'Γ'                           The reflection coefficient of the medium of the plane earth.
	"""
	function two_ray_propagation( distance, transmit_height, receive_height, frequency_hz, Γ )
		λ = c/frequency_hz;
		Δϕ=  4*π*transmit_height*receive_height/(λ*distance);
		F = 1+abs(Γ)*exp(i*Δϕ);
		return F;
	end

	"""
	Empirical model for RF rain attenuation for frequencies between 1 and 400 GHz, linear polarization.
	Model uses the closest frequency in the underlying data.

	- M. A. Richards and J. A. Scheer and W. A. Holm, Principles of Modern Radar, SciTech Publishing, 2010.

	```julia
	rain_attenuation_db_per_km_circular_pol( 'v' , 30, 20 )
	3.3400000000000003
	```

	# Arguments
	- 'polarization'        The polarization, 'v', 'h'.
	- 'frequency_ghz'     	The frequency of the propagating waves.
	- 'fall_rate_mm_hour'	The rain intensity [mm/h].
	"""
	function rain_attenuation_db_per_km( polarization::Char, frequency_ghz, fall_rate_mm_hour )
		r = fall_rate_mm_hour;

		# Empirical data:
		frequency_table_ghz = [1, 2, 3 ,6, 7, 8, 10, 12, 15, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80, 90, 100, 120, 150, 200, 300, 400 ];
		ah = [ 3.87e-5, 1.54e-4, 6.5e-4, 1.75e-3, 3.01e-3, 4.54e-3, 1.01e-2, 1.88e-2, 3.67e-2, 7.51e-2, .124, .187, .263, .350, .442, .536, .707, .851, .975, 1.06, 1.12, 1.18, 1.31, 1.45, 1.36, 1.32 ];
		av = [ 3.52e-5, 1.38e-4, 5.91e-4, 1.55e-3, 2.65e-3, 3.95e-3, 8.87e-3, 1.68e-2, 3.47e-2, 6.91e-2, .113, .167, .233, .310, .393, .479, .642, .784, .906, .999, 1.06, 1.13, 1.27, 1.42, 1.35, 1.31 ];
		bh = [ .912, .963, 1.121, 1.308, 1.332, 1.327, 1.276, 1.217, 1.154, 1.099, 1.061, 1.021, .979, .939, .903, .873, .826, .793, .769, .753, .743, .731, .710, .689, .688, .683 ];
		bv = [ .88, .923, 1.075, 1.265, 1.312, 1.310, 1.264, 1.2, 1.128, 1.065, 1.030, 1.0, .963, .929, .897, .868, .824, .793, .769, .754, .744, .732, .711, .69, .689, .683 ];

		# TODO implement interpolation vor arbitrary frequency selection.
		index = argmin( abs.(  frequency_ghz.-frequency_table_ghz ) );
		if polarization == 'v'
			a = av[index];
			b = bv[index];
		else
			a = ah[index];
			b = bh[index];
		end

		α = a*r^b;
		return α;
   end

   """
	Empirical model for rain attenuation for frequencies between 1 and 400 GHz, circular polarization.
	Model uses closest frequency in the underlying data.
	One-way attenuation in dB/km.

	- M. A. Richards and J. A. Scheer and W. A. Holm, Principles of Modern Radar, SciTech Publishing, 2010.

	```julia
	rain_attenuation_db_per_km_circular_pol( 30, 20 )
	3.6755027981960815
	```

	# Arguments
	- 'frequency_ghz'     	The frequency of the propagating waves.
	- 'fall_rate_mm_hour'	The rain intensity [mm/h].
	"""
	function rain_attenuation_db_per_km_circular_pol( frequency_ghz, fall_rate_mm_hour )
		α_lin( pol ) = rain_attenuation_db_per_km( pol, frequency_ghz, fall_rate_mm_hour );
		αv = α_lin('v');
		αh = α_lin('h');

		α = 1/sqrt(2) * sqrt(αv^2+αh^2);

		return α;
	end

	"""
	Empirical model for rain attenuation for frequencies above 5 GHz.
	One-way attenuation in dB/km.

	- M. A. Richards and J. A. Scheer and W. A. Holm, Principles of Modern Radar, SciTech Publishing, 2010.

	```julia
	fog_attenuation_db_per_km_circular_pol( 10, 0.8, 23 );
	4.68976
	```

	# Arguments
	- 'frequency_ghz'	The frequency of the propagating waves.
	- 'M'				The watewater concentration in g/m³.
	- 'T_deg'			The air temperature in degree Celsius.
	"""
	function fog_attenuation_db_per_km( frequency_ghz, M, T_deg )
		T = T_deg;
		f = frequency_ghz

		α= M*( -1.347 + 0.66*f+(11.152/f) - 0.022*T )
		return α;
	end

	"""
	Empirical model for amospheric gaseous attenuation for frequencies in the range 1 - 1000 GHz.
	One-way attenuation in dB/km.

	- Rec. ITU-R P.676-12, Attenuation by atmospheric gases and related effects, ITU-R 2020.

	```julia
	using Plots

	frequency_ghz = collect(1:1000);

	attenuation_db_km = RadioPropagation.atmospheric_attenuation_db_per_km.( frequency_ghz );

	plot1 = plot(	frequency_ghz, attenuation_db_km,
					xlabel 	= "Frequency [GHz]",
					ylabel 	= "One-way loss [dB/km]",
					title  	= "Atmostpheric Attenuation.",
					#ylims 	= (),
					label 	= "Standard",
					yaxis=:log);
	```

	# Arguments
	- 'frequency_ghz'				The frequency of the propagating waves.
	- 'T_kelvin'					The absolute temperature in kelvin.
	- 'water_vapour_density_g_m³' 	The water vapur density in g/m³.
	- 'dry_air_pressure_h_pa' 		The dry air pressure in hpa.
	"""
	function atmospheric_attenuation_db_per_km( frequency_ghz, T_kelvin=288.15, water_vapour_density_g_m³=7.5, dry_air_pressure_h_pa=1013.25 )
		f = frequency_ghz;
		T = T_kelvin;
		ρ = water_vapour_density_g_m³;
		p = dry_air_pressure_h_pa;	# Dry air pressure [hPa].
		e = ρ*T/216.7;				# Water vapour partial pressure [hPa].
		θ = 300/T;

		# Spectroscopic data for oxygen attenuation.
		# 		  f0 			a1			a2		a3		a4 		a5 		a6
		tab1 = 	[ 50.474214		0.975 		9.651	6.690 	0.0		2.566 	6.850
				  50.987745		2.529 		8.653	7.170 	0.0		2.246 	6.800
				  51.503360		6.193 		7.709	7.640 	0.0		1.947 	6.729
				  52.021429		14.320 		6.819	8.110 	0.0		1.667 	6.640
				  52.542418		31.240 		5.983	8.580 	0.0		1.388 	6.526
				  53.066934		64.290 		5.201	9.060 	0.0		1.349 	6.206
				  53.595775		124.600 	4.474	9.550 	0.0		2.227 	5.085
				  54.130025		227.300 	3.800	9.960 	0.0		3.170 	3.750
				  54.671180		389.700 	3.182	10.370 	0.0		3.558 	2.654
				  55.221384		627.100 	2.618	10.890	0.0		2.560 	2.952
				  55.783815		945.300 	2.109	11.340	0.0		-1.172 	6.135
				  56.264774		543.400 	0.014	17.030	0.0		3.525 	-0.978
				  56.363399		1331.800 	1.654	11.890	0.0		-2.378	6.547
				  56.968211		1746.600 	1.255	12.230	0.0		-3.545	6.451
				  57.612486		2120.100 	0.910	12.620	0.0		-5.416	6.056
				  58.323877		2363.700 	0.621	12.950	0.0		-1.932	0.436
				  58.446588		1442.100	0.083	14.910	0.0		6.768 	-1.273
				  59.164204		2379.900 	0.387	13.530	0.0		-6.561 	2.309
				  59.590983		2090.700 	0.207	14.080	0.0		6.957 	-0.776
				  60.306056		2103.400 	0.207	14.150	0.0		-6.395 	0.699
				  60.434778		2438.000 	0.386	13.390	0.0		6.342 	-2.825
				  61.150562		2479.500 	0.621	12.920	0.0		1.014 	-0.584
				  61.800158		2275.900 	0.910	12.630	0.0		5.014 	-6.619
				  62.411220		1915.400 	1.255	12.170	0.0		3.029 	-6.759
				  62.486253		1503.000 	0.083	15.130	0.0		-4.499 	0.844
				  62.997984		1490.200 	1.654	11.740	0.0		1.856 	-6.675
				  63.568526		1078.000 	2.108	11.340	0.0		0.658 	-6.139
				  64.127775		728.700 	2.617	10.880	0.0		-3.036 	-2.895
				  64.678910		461.300 	3.181	10.380	0.0		-3.968 	-2.590
				  65.224078		274.000 	3.800	9.960   0.0		-3.528 	-3.680
				  65.764779		153.000 	4.473	9.550   0.0		-2.548 	-5.002
				  66.302096		80.400 		5.200	9.060   0.0		-1.660 	-6.091
				  66.836834		39.800 		5.982	8.580   0.0		-1.680	-6.393
				  67.369601		18.560 		6.818	8.110   0.0		-1.956 	-6.475
				  67.900868		8.172 		7.708	7.640   0.0		-2.216 	-6.545
				  68.431006		3.397 		8.652	7.170   0.0		-2.492	-6.600
				  68.960312		1.334 		9.650	6.690   0.0		-2.773 	-6.650
				  118.750334 	940.300		0.010	16.640	0.0		-0.439	0.079
				  368.498246	67.400		0.048	16.400	0.0		0.0		0.0
				  424.763020	637.700		0.044	16.400	0.0		0.0		0.0
				  487.249273	237.400		0.049	16.000	0.0		0.0		0.0
				  715.392902	98.100		0.145	16.000	0.0		0.0		0.0
				  773.839490	572.300		0.141	16.200	0.0		0.0		0.0
				  834.145546	183.100		0.145	14.700	0.0		0.0		0.0]

		f0a = @view tab1[:,1];
		a1  = @view tab1[:,2];
		a2  = @view tab1[:,3];
		a3  = @view tab1[:,4];
		a4  = @view tab1[:,5];
		a5  = @view tab1[:,6];
		a6  = @view tab1[:,7];

		# Spectroscopic data for water vapour attenuation.
		# 		 f0 		b1		b2		b3		b4 		b5 		b6
		tab2 = [ 22.235080 	.1079	2.144	26.38	.76		5.087	1.00
				 67.803960 	.0011	8.732	28.58	.69		4.930	.82
				 119.995940 .0007	8.353	29.48	.70		4.780	.79
				 183.310087 2.273	.668 	29.06	.77		5.022	.85
				 321.225630 .0470	6.179	24.04	.67		4.398	.54
				 325.152888 1.514	1.541	28.23	.64		4.893	.74
				 336.227764 .0010	9.825	26.93	.69		4.740	.61
				 380.197353 11.67	1.048	28.11	.54		5.063	.89
				 390.134508 .0045	7.347	21.52	.63		4.810	.55
				 437.346667 .0632	5.048	18.45	.60		4.230	.48
				 439.150807 .9098	3.595	20.07	.63		4.483	.52
				 443.018343 .1920	5.048	15.55	.60		5.083	.50
				 448.001085 10.41	1.405	25.64	.66		5.028	.67
				 470.888999 .3254	3.597	21.34	.66		4.506	.65
				 474.689092 1.260	2.379	23.20	.65		4.804	.64
				 488.490108 .2529	2.852	25.86	.69		5.201	.72
				 503.568532 .0372	6.731	16.12	.61		3.980	.43
				 504.482692 .0124	6.731	16.12	.61		4.010	.45
				 547.676440 .9785	.158	26.00	.70		4.500	1.00
				 552.020960 .1840	.158	26.00	.70		4.500	1.00
				 556.935985 497.0	.159	30.86	.69		4.552	1.00
				 620.700807 5.015	2.391	24.38	.71		4.856	.68
				 645.766085 .0067	8.633	18.00	.60		4.000	.50
				 658.005280 .2732	7.816	32.10	.69		4.140	1.00
				 752.033113 243.4	.396 	30.86	.68		4.352	.84
				 841.051732 .0134	8.177	15.90	.33		5.760	.45
				 859.965698 .1325	8.055	30.60	.68		4.090	.84
				 899.303175 .0547	7.914	29.85	.68		4.530	.90
				 902.611085 .0386	8.429	28.65	.70		5.100	.95
				 906.205957 .1836	5.110	24.08	.70		4.700	.53
				 916.171582 8.400	1.441	26.73	.70		5.150	.78
				 923.112692 .0079	10.293 	29.00	.70		5.000	.80
				 970.315022 9.009	1.919 	25.50	.64		4.940	.67
				 987.926764 134.6	.257 	29.85	.68		4.550	.90
				 1780.000000 17506. .952 	196.3	2.00 	24.15	5.00]

		f0b = @view tab2[:,1];
		b1  = @view tab2[:,2];
		b2  = @view tab2[:,3];
		b3  = @view tab2[:,4];
		b4  = @view tab2[:,5];
		b5  = @view tab2[:,6];
		b6  = @view tab2[:,7];

		S_o		= @. a1*10^(-7)*p*θ^3*exp( a2*(1-θ));
		S_w		= @. b1*10^(-1)*e*θ^3.5*exp(b2*(1-θ));

		Δf_o_	= @. a3*10^(-4) *(p*θ^(0.8-a4) + 1.1*e*θ);
		Δf_o	= @. sqrt( Δf_o_^2 + 2.25*10^(-6) );

		Δf_w_	= @. b3*10^(-4) *(p*θ^b4 + b5*e*θ^b6);
		Δf_w	= @. 0.535*Δf_w_ + sqrt( 0.217*Δf_w_^2 + (2.1316*10^-12*f0b^2)/ θ );

		δ_o		= @. (a5+a6*θ)*10^(-4) *(p+e)*θ^(0.8);
		δ_w		= 0;

		F_o		= @. f/f0a *( (Δf_o-δ_o*(f0a-f))/((f0a-f)^2 + Δf_o^2 ) + (Δf_o-δ_o*(f0a+f))/((f0a+f)^2 + Δf_o^2) );
		F_w		= @. f/f0b *( (Δf_w-δ_w*(f0b-f))/((f0b-f)^2 + Δf_w^2 ) + (Δf_w-δ_w*(f0b+f))/((f0b+f)^2 + Δf_w^2) );

		d		= 5.6*10^(-4) * (p+e)*θ^0.8;

		N_D_prime = f*p*θ^2 * ( (6.14*10^(-5) / ( d*(1+(f/d)^2) )) +  (1.4*10^(-12) *p*θ^1.5) / ( 1+1.9*10^(-5) *f^1.5 ) );

		N_d_prime_oxygen		= sum( S_o.*F_o) + N_D_prime;
		N_d_prime_water_vapour	= sum( S_w.*F_w );

		γ = 0.1820*f*(N_d_prime_oxygen+N_d_prime_water_vapour);
		return γ;
	end
end

max_stars_repo_name: ErikBuer/RadioPropagation

max_issues_repo_name: ErikBuer/RadioPropagation

max_forks_repo_name: ErikBuer/RadioPropagation