text:
@doc raw"
    random_point(M::Manifold)

generate a random point on a manifold. By default it uses `random_point(M,:Gaussian)`.
"
random_point(M::Manifold) = random_point(M, Val(:Gaussian))
@doc raw"
    random_point(M::Manifold, s::Symbol, options...)

generate a random point using a noise model given by `s` with its additional `options`
just passed on.
"
random_point(M::Manifold, s::Symbol, options...) = random_point(M, Val(s), options...)

@doc raw"""
    random_point(M::AbstractPowerManifold, options...)

generate a random point on the `AbstractPowerManfold` `M` given `options` that are
passed on.
"""
function random_point(
    M::AbstractPowerManifold{𝔽,Mt,NestedPowerRepresentation}, options...
) where {𝔽,Mt}
    return [random_point(M.manifold, options...) for i in get_iterator(M)]
end
function random_point(
    M::AbstractPowerManifold{𝔽,Mt,ArrayPowerRepresentation}, options...
) where {𝔽,Mt}
    return cat(
        [random_point(M.manifold, options...) for i in get_iterator(M)]...;
        dims=length(representation_size(M.manifold)) + 1,
    )
end

@doc raw"""
    random_point(M::AbstractGroupManifold, options...)

On an abstract group manifold, the random point is taken from the internally stored `M.manifold`.
"""
random_point(M::AbstractGroupManifold, kwargs...) = random_point(M.manifold, kwargs...)
function random_point(M::AbstractGroupManifold, s::Symbol, options...)
    return random_point(M, Val(s), options...)
end

@doc raw"""
    random_point(M::Circle, :Uniform)

return a random point on the [`Circle`](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold) $\mathbb S^1$ by
picking a random element from $[-\pi,\pi)$ uniformly.
"""
random_point(::Circle, ::Val{:Uniform}) = sym_rem(rand() * 2 * π)
random_point(M::Circle) = random_point(M, Val(:Uniform)) # introduce different default

@doc raw"""
    random_point(M::Euclidean[,:Gaussian, σ::Float64=1.0])

generate a random point on the `Euclidean` manifold `M`, where the
optional parameter determines the type of the entries of the
resulting point on the Euclidean space d.
"""
random_point(M::Euclidean) = randn(representation_size(M))
random_point(M::Euclidean, ::Val{:Gaussian}, σ=1.0) = σ * randn(manifold_dimension(M))

@doc raw"""
    random_point(M::FixedRankMatrices, options...)

return a random point on the FixedRankMatrices manifold.
The orthogonal matrices are sampled from the Stiefel manifold
and the singular values are sampled uniformly at random.
"""
function random_point(::FixedRankMatrices{m,n,k}, options...) where {m,n,k}
    U = random_point(Stiefel(m, k), options...)
    S = sort(rand(k); rev=true)
    V = random_point(Stiefel(n, k), options...)
    return SVDMPoint(U, S, V')
end

@doc raw"""
    random_point(M::Grassmannian, :Gaussian [, σ=1.0])

return a random point `x` on `Grassmannian` manifold `M` by
generating a random (Gaussian) matrix with standard deviation `σ` in matching
size, which is orthonormal.
"""
function random_point(::Grassmann{n,k,𝔽}, ::Val{:Gaussian}, σ::Float64=1.0) where {n,k,𝔽}
    V = σ * randn(𝔽 === ℝ ? Float64 : ComplexF64, (n, k))
    A = qr(V).Q[:, 1:k]
    return A
end

@doc raw"""
    random_point(M::ProductManifold, options...)

return a random point `x` on `Grassmannian` manifold `M` by
generating a random (Gaussian) matrix with standard deviation `σ` in matching
size, which is orthonormal.
"""
function random_point(M::ProductManifold, o...)
    return ProductRepr([random_point(N, o...) for N in M.manifolds]...)
end

@doc raw"""
    random_point(M::Rotations, :Gaussian [, σ=1.0])

return a random point `p` on the manifold `Rotations`
by generating a (Gaussian) random orthogonal matrix with determinant $+1$. Let $QR = A$
be the QR decomposition of a random matrix $A$, then the formula reads $p = QD$
where $D$ is a diagonal matrix with the signs of the diagonal entries of $R$, i.e.
````math
D_{ij}=\begin{cases}
\operatorname{sgn}(R_{ij}) & \text{if} \; i=j \\
0 & \, \text{otherwise.}
\end{cases}
````
It can happen that the matrix gets -1 as a determinant. In this case, the first
and second columns are swapped.
"""
function random_point(M::Rotations, ::Val{:Gaussian}, σ::Real=1.0)
    d = Int(0.5 + sqrt(0.25 + 2 * manifold_dimension(M)))
    if d == 1
        return ones(1, 1)
    else
        A = randn(Float64, d, d)
        s = diag(sign.(qr(A).R))
        D = Diagonal(s)
        C = qr(A).Q * D
        if det(C) < 0
            C[:, [1, 2]] = C[:, [2, 1]]
        end
        return C
    end
end

@doc raw"""
    random_point(M::SymmetricPositiveDefinite, :Gaussian[, σ=1.0])

gerenate a random symmetric positive definite matrix on the
`SymmetricPositiveDefinite` manifold `M`.
"""
function random_point(
    ::SymmetricPositiveDefinite{N}, ::Val{:Gaussian}, σ::Float64=1.0
) where {N}
    D = Diagonal(1 .+ rand(N)) # random diagonal matrix
    s = qr(σ * randn(N, N)) # random q
    return Matrix(Symmetric(s.Q * D * transpose(s.Q)))
end

@doc raw"""
    random_point(M::Stiefel, :Gaussian[, σ=1.0])

return a random (Gaussian) point `x` on the `Stiefel` manifold `M` by generating a (Gaussian)
matrix with standard deviation `σ` and return the orthogonalized version, i.e. return the Q
component of the QR decomposition of the random matrix of size $n×k$.
"""
function random_point(::Stiefel{n,k,𝔽}, ::Val{:Gaussian}, σ::Float64=1.0) where {n,k,𝔽}
    A = σ * randn(𝔽 === ℝ ? Float64 : ComplexF64, n, k)
    return Matrix(qr(A).Q)
end

@doc raw"""
    random_point(M::Sphere, :Gaussian[, σ=1.0])
return a random point on the Sphere by projecting a normal distirbuted vector
from within the embedding to the sphere.
"""
function random_point(M::Sphere, ::Val{:Gaussian}, σ::Float64=1.0)
    return project(M, σ * randn(manifold_dimension(M) + 1))
end

@doc raw"""
    random_point(M::TangentBundle, options...)

generate a random point on the tangent bundle by calling a [`random_point`](@ref) and a
[`random_tangent`](@ref) with the given `options...`
"""
function random_point(M::TangentBundle, options...)
    p = random_point(M.manifold, options...)
    X = random_tangent(M.manifold, p, options...)
    return ProductRepr(p, X)
end

@doc raw"""
    random_tangent(M::AbstractGroupManifold, p, options...)

On an abstract group manifold, the random tangent is taken from the internally stored `M.manifold`s tangent space at `p`.
"""
function random_tangent(M::AbstractGroupManifold, p, kwargs...)
    return random_tangent(M.manifold, p, kwargs...)
end
function random_tangent(M::AbstractGroupManifold, p, s::Symbol, options...)
    return random_tangent(M, p, Val(s), options...)
end

@doc raw"""
    random_tangent(M, p, options...)

generate a random tangent vector in the tangent space of `p` on `M`. By default
this is a `:Gaussian` distribution.
"""
random_tangent(M::Manifold, p, options...) = random_tangent(M, p, :Gaussian, options...)
function random_tangent(M::Manifold, p, s::Symbol, options...)
    return random_tangent(M, p, Val(s), options...)
end

@doc raw"""
    random_tangent(M::Circle, p [, :Gaussian, σ=1.0])

return a random tangent vector from the tangent space of the point `p` on the
[Circle](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold) $\mathbb S^1$ by using a normal distribution with
mean 0 and standard deviation 1.
"""
random_tangent(::Circle, p, ::Val{:Gaussian}, σ::Real=1.0) = σ * randn()

function random_tangent(M::Euclidean, p, ::Val{:Gaussian}, σ::Float64=1.0)
    return σ * randn(manifold_dimension(M))
end

@doc raw"""
    random_tangent(M::Grassmann, p[,type=:Gaussian, σ=1.0])

return a (Gaussian) random vector from the tangent space $T_p\mathrm{Gr}(n,k)$ with mean
zero and standard deviation `σ` by projecting a random Matrix onto the  `p`.
"""
function random_tangent(M::Grassmann, p, ::Val{:Gaussian}, σ::Float64=1.0)
    Z = σ * randn(eltype(p), size(p))
    X = project(M, p, Z)
    X = X ./ norm(X)
    return X
end

@doc raw"""
    random_tangent(M::FixedRankMatrices, p, options...)

generate a random tangent vector in the tangent space of the point `p` on the
`FixedRankMatrices` manifold `M`.
"""
function random_tangent(::FixedRankMatrices{m,n,k}, p, options...) where {m,n,k}
    Up = randn(m, k, options...)
    Vp = randn(n, k, options...)
    A = randn(k, k, options...)
    return UMVTVector(Up - p.U * p.U' * Up, A, Vp' - Vp' * p.Vt' * p.Vt)
end

@doc raw"""
    random_tangent(M::Hyperbolic, p, :Gaussian [, σ=1.0])

generate a random point on the Hyperbolic manifold by projecting a point from the embedding
with respect to the Minkowsky metric.
"""
function random_tangent(M::Hyperbolic, p, ::Val{:Gaussian}, σ=1.0)
    Y = σ * randn(eltype(p), size(p))
    X = project(M, p, Y)
    return X
end

function random_tangent(M::PowerManifold, p, options...)
    rep_size = representation_size(M.manifold)
    X = zero_tangent_vector(M, p)
    for i in get_iterator(M)
        X[M, i] = random_tangent(M.manifold, p[M, i], options...)
    end
    return X
end

@doc raw"""
    random_tangent(M::ProductManifold, p)

generate a random tangent vector in the tangent space of the point `p` on the
`ProductManifold` `M`.
"""
function random_tangent(M::ProductManifold, p, options...)
    X = map(
        (m, p) -> random_tangent(m, p, options...),
        M.manifolds,
        submanifold_components(M, p),
    )
    return ProductRepr(X...)
end

@doc raw"""
    random_tangent(M::Rotations, p[, type=:Gaussian, σ=1.0])

return a random tangent vector in the tangent space
$T_x\mathrm{SO}(n)$ of the point `x` on the `Rotations` manifold `M` by generating
a random skew-symmetric matrix. The function takes the real upper triangular matrix of a
(Gaussian) random matrix $A$ with dimension $n\times n$ and subtracts its transposed matrix.
Finally, the matrix is normalized.
"""
function random_tangent(M::Rotations, p, ::Val{:Gaussian}, σ::Real=1.0)
    d = manifold_dimension(M)
    if d == 0
        return zeros(1, 1)
    else
        A = σ * randn(Float64, d, d)
        A = triu(A, 1) - transpose(triu(A, 1))
        A = (1 / norm(A)) * A
        return A
    end
end

@doc raw"""
    random_tangent(M::Sphere, p[, :Gaussian, σ=1.0])

return a random tangent vector in the tangent space of `p` on the `Sphere` `M`.
"""
function random_tangent(M::Sphere, p, ::Val{:Gaussian}, σ::Float64=1.0)
    n = σ * randn(size(p)) # Gaussian in embedding
    return project(M, p, n) #project to TpM (keeps Gaussianness)
end

@doc raw"""
    random_tangent(M::Stiefel, p[,type=:Gaussian, σ=1.0])

return a (Gaussian) random vector from the tangent space $T_p\mathrm{St}(n,k)$ with mean
zero and standard deviation `σ` by projecting a random Matrix onto the  `p`.
"""
function random_tangent(M::Stiefel, p, ::Val{:Gaussian}, σ::Float64=1.0)
    Z = σ * randn(eltype(p), size(p))
    X = project(M, p, Z)
    X = X ./ norm(X)
    return X
end

@doc raw"""
    random_tangent(M::SymmetricPositiveDefinite, p[, :Gaussian, σ = 1.0])

generate a random tangent vector in the tangent space of the point `p` on the
`SymmetricPositiveDefinite` manifold `M` by using a Gaussian distribution
with standard deviation `σ` on an ONB of the tangent space.
"""
function random_tangent(M::SymmetricPositiveDefinite, p, ::Val{:Gaussian}, σ::Float64=0.01)
    # generate ONB in TxM
    I = one(p)
    B = get_basis(M, p, DiagonalizingOrthonormalBasis(I))
    Ξ = get_vectors(M, p, B)
    Ξx = vector_transport_to.(Ref(M), Ref(I), Ξ, Ref(p), Ref(ParallelTransport()))
    return sum(σ * randn(length(Ξx)) .* Ξx)
end

@doc raw"""
    random_tangent(M::SymmetricPositiveDefinite, p, :Rician [,σ = 0.01])

generate a random tangent vector in the tangent space of `p` on
the `SymmetricPositiveDefinite` manifold `M` by using a Rician distribution
with standard deviation `σ`.
"""
function random_tangent(::SymmetricPositiveDefinite, p, ::Val{:Rician}, σ::Real=0.01)
    # Rician
    C = cholesky(Hermitian(p))
    R = C.L + sqrt(σ) * triu(randn(size(p, 1), size(p, 2)), 0)
    return R * R'
end

@doc raw"""
    random_tangent(M::TangentBundle, p, options...)

generate a random tangent vector at `p` on the tangent bundle by calling
[`random_tangent`](@ref) with the given `options...` twice.
"""
function random_tangent(M::TangentBundle, p, options...)
    X = random_tangent(M.manifold, p[M, :point], options...)
    Y = random_tangent(M.manifold, p[M, :point], options...)
    return ProductRepr(X, Y)
end


text:
# look at SimpleExample.jl for TCP based usage example

nprocs() > 1 ? thxl = 2 : nothing

USEREADY=1
include("SlamConvenienceFunctions.jl")

function tcpStringSLAMServer(;slamdata=nothing,port::Int=60001)
  println("Empty slam object created")
  if slamdata == nothing
    slamdata = SLAMWrapper(emptyFactorGraph(), nothing, 0)
  end

  println("Listenting on $(port)")
  server = listen(port)
  loop = true
  while loop
   sock = accept(server)
     while isopen(sock)
         loop, retstr = parseTCP!(slamdata, readline(sock)[1:(end-1)])
         loop ? (isopen(sock) ? println(sock, retstr) : nothing) : close(sock)
     end
     println("connection lost")
  end
  !loop ? close(server) : nothing
  return slamdata
end

ext: jl


text:
function gaussian(A::AbstractMatrix, b::AbstractVector)
  #= Solve linear systems using Gaussian elimination

    Linear systems should be in the form of Ax=b where A is a matrix and b & x
    are column vectors.

    Keyword arguments:
      A -- A matrix containing the coefficients of the variables to solve for
      b -- A column vector containing the RHS of the linear system
    Returns:
      A column vector containing the solutions for x1..xn
  =#
  N, M = size(A);
  x = zeros(M);

  # Copies of A & b to avoid mutation
  cpA = copy(A);
  cpb = copy(b);

  for j = 1:M
    # Maximise the pivot point by looking for larger pivot points in the rows
    # beneath the current point's row.
    maxPivotRowIndex = indmax(abs(cpA[j:end, j])) + (j - 1);
    cpA[maxPivotRowIndex, :], cpA[j, :] = cpA[j, :], cpA[maxPivotRowIndex, :];
    cpb[maxPivotRowIndex, :], cpb[j, :] = cpb[j, :], cpb[maxPivotRowIndex, :];

    # Iterate through the columns then rows of the linear system
    for i = j+1:N
      # If the current element is 0 then it is already reduced
      if cpA[i, j] == 0
        continue
      end

      # Perform row reduction on the current element's row using the pivot
      # point.
      multFactor = cpA[i,j] / cpA[j,j];
      cpA[i, :] = cpA[i, :] - multFactor * cpA[j, :];
      cpb[i] = cpb[i] - multFactor*cpb[j];
    end
  end

  # Perform a back substitution for x1, x2, ..., xn.
  x[end] = cpb[end] / cpA[end, end]; # Calculate xn
  for i = N:-1:1 # Calculate xn-1...x1
    x[i] = (cpb[i] - sum(cpA[i, i+1:N] * x[i+1:N])) / cpA[i,i];
  end
  return x;
end

max_stars_repo_path: Julia/Linear Algebra/gaussian.jl

max_stars_repo_name: alexjohnj/numerical-methods

max_issues_repo_path: Julia/Linear Algebra/gaussian.jl

max_issues_repo_name: alexjohnj/numerical-methods

max_forks_repo_path: Julia/Linear Algebra/gaussian.jl

max_forks_repo_name: alexjohnj/numerical-methods


text:
include("..\\libs\\FunctionsRobot.jl")

function chessDesk!(r)
    back_path = BackPath(r, (Nord, West))
    R = ChessRobot(r, iseven(numSteps(back_path)))
    marks!(R, (Sud, Ost))
    BackPath(R, (Nord, West))
    back!(R, back_path)

end


text:
using SpgpGpu
using Documenter

DocMeta.setdocmeta!(SpgpGpu, :DocTestSetup, :(using SpgpGpu); recursive=true)

makedocs(;
    modules=[SpgpGpu],
    authors="Deepak Selvan <dselvan@gmail.com> and contributors",
    repo="https://github.com/dselvan/SpgpGpu.jl/blob/{commit}{path}#{line}",
    sitename="SpgpGpu.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://dselvan.github.io/SpgpGpu.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/dselvan/SpgpGpu.jl",
    devbranch="main"
)


text:
using RCall

@test RCall.render("x = 1")[3] == 1

@test RCall.render("x = 'α'")[3] == 1

@test RCall.render("x = \$y")[3] == 1

@test RCall.render("x = \$(rand(1))")[3] == 1

@test RCall.render("x = \$α")[2]["α"] == :α

@test RCall.render("x = ]")[4] == "unexpected ']'"

@test RCall.render("x = \$(begin")[3] == 2

@test RCall.render("x = \$(begin)")[3] == 3

@test RCall.render("x = ")[4] == "unexpected end of input"

@test rcopy(R"sum($[7,1,3])") == sum([7,1,3])

iris = rcopy(:iris)
model =  R"lm(Sepal.Length ~ Sepal.Width,data=$iris)"
@test rcopy(RCall.getclass(model)) == "lm"
@test isapprox(rcopy(R"sum($iris$Sepal.Length)"), sum(iris[Symbol("Sepal.Length")]), rtol=4*eps())


text:
push!(LOAD_PATH, pwd())
const assets = "assets"
using Images, ImageMagick, Colors, Photopolarimetry
import Tk.ChooseDirectory
#include("Photopolar.jl")
include("PhotopolarGUIfunctions.jl")

path = Input(ChooseDirectory())
flip = Input(false)
flop = Input(false)
torun = Input{Any}(leftbutton)
rotate = Input(0)
cropright = Input(100)
cropbottom = Input(100)
cropleft = Input(0)
croptop = Input(0)
scale = Input(20)
topI = Input(100)
topdolp = Input(100)
topdolcp = Input(100)
topdorcp = Input(100)
botprofile = Input(0)
topprofile = Input(100)
Sp = Input(5)
ellipse_window = Input(0.5)
minpol = Input(0.2)
maxpol = Input(1.0)
eccentricity = Input(0.01)

torun2 = Input{Any}(leftbutton)

fnames = lift(getfnames,path)
SZ = lift(x -> getsz(x[1]),fnames)
lift(x -> write_originalRGB(x[1]),fnames)

row1 = ["Flip", checkbox(false) >>> flip, "Flop", checkbox(false) >>> flop, button("Done!") >>> torun]
row2 = ["Rotate", slider(-90:90,value=0) >>> rotate, "Scale", slider(1:100,value=20) >>> scale]
row3 = ["Crop top", slider(0:99,value=0) >>> croptop, "Crop bottom", slider(1:100,value=100) >>> cropbottom]
row4 = ["Crop left", slider(0:99,value=0) >>> cropleft, "Crop right", slider(1:100,value=100) >>> cropright]

mymin(u,l) = u > l ? l : u - 1
botprofile2 = lift(mymin, topprofile, botprofile)

s = lift(_ -> path2stokes(value(path),value(flip),value(flop),value(rotate),value(cropright),value(cropbottom),value(cropleft),value(croptop),value(scale), angleoffset = deg2rad(60)), torun, typ=Any, init=empty)
polarimgs = lift(calculatepolar, s, typ=Any, init=empty)

I = lift((x,t) -> recolor_I(x[1],t),polarimgs,topI,typ=Any, init=empty)
dolp = lift((x,t) -> recolor_dolp(x[2],t),polarimgs,topdolp,typ=Any, init=empty)
docp = lift((x,b,t) -> recolor_docp(x[4],b,t),polarimgs,topdolcp,topdorcp,typ=Any, init=empty)

profileimg = lift(stokes2profile, s, Sp, typ=Any, init=empty)
ellipse_img_name = lift(latexellipse, s, ellipse_window, maxpol, minpol, eccentricity, typ=Any, init=empty)

orientationtab = vbox(map(pad(0.5em),[
	hbox(map(pad(0.5em), row1)),
	hbox(map(pad(0.5em), row2)),
	hbox(map(pad(0.5em), row3)),
	hbox(map(pad(0.5em), row4)),
	consume(flip,flop,rotate,cropright,cropbottom,cropleft,croptop,scale) do flp, flP, rtt, crprght, crpbttm, crplft, crptp, scl
		@async cleanassets()
		name = updateRGB(flp,flP,rtt,crprght,crpbttm,crplft,crptp,scl)
		image("$assets/$name.jpg")
	end,
	]))

intensitytab = vbox(map(pad(0.5em),[
	hbox("Intensity", slider(0:100, value = 100) >>> topI),
	hbox("DoLP", slider(0:100, value = 100) >>> topdolp),
	hbox("DoLCP", slider(0:100, value = 100) >>> topdolcp),
	hbox("DoRCP", slider(0:100, value = 100) >>> topdorcp),
	consume(I, dolp, docp, polarimgs, typ=Any, init=empty) do i, dlp, dcp, p
		buff = Image(ones(RGB{FixedPointNumbers.UFixed{UInt8,8}},10,size(i,2)))
		cat(1,i, buff, dlp, buff, p[3], buff, dcp)
	end
	]))

profiletab = vbox(map(pad(0.5em),[
	hbox(map(pad(0.5em),["upper", slider(2:100, value = 100) >>> topprofile])),
	consume(topprofile) do x
		hbox(map(pad(0.5em),["Lower", slider(0:(x - 1), value = 0) >>> botprofile]))
	end,
	hbox(map(pad(0.5em),[tex("S_p"), slider(2:10, value = 5) >>> Sp])),
iiimg = lift(recolor_profile, profileimg, botprofile2, topprofile, typ=Any, init=empty)
]))

lift(_ -> println(join([splitdir(value(path))[2],value(flip),value(flop),value(rotate),value(cropright),value(cropbottom),value(cropleft),value(croptop),value(scale),value(ellipse_window), value(minpol), value(maxpol), value(eccentricity)],",")), torun2, typ=Any, init=empty)

ellipsetab = vbox(map(pad(0.5em),[
	hbox(map(pad(0.5em),["Cell size (cm)", slider(0.1:0.05:2, value = 0.5) >>> ellipse_window])),
	hbox(map(pad(0.5em),["Min pol.", slider(0:.01:1, value = 0.2) >>> minpol])),
	hbox(map(pad(0.5em),["Max pol.", slider(0:.01:1, value = 1.0) >>> maxpol])),
	hbox(map(pad(0.5em),["Eccentricity", slider(0:.01:1, value = 0.01) >>> eccentricity])),
	button("Done!") >>> torun2,
	lift(image,ellipse_img_name, typ=Any, init=empty)
]))

#=ok = Input{Any}(leftbutton)
donetab = vbox(map(pad(0.5em),[
	button("Print", raised = true) >>> ok,
	consume(ok, typ=Any, init=empty) do o
		colorimg = load("$assets/RGB.jpg")
		aopimg = convert(Image, value(X)[3])
		Iimg = convert(Image, value(I))
		dolpimg = convert(Image, value(dolp))
		docpimg = convert(Image, value(docp))
		buff = round(Int,0.05*mean(widthheight(colorimg)))
		buffimg = Iimg[1:buff,:]
		buffimg[:] = RGB{U8}(1,1,1)
		img = cat(1,colorimg, buffimg, Iimg, buffimg, dolpimg, buffimg, aopimg, buffimg, docpimg)
		name = string(abs(rand(Int)))
		save( "$assets/$name.jpg",img)
		image("assets/$name.jpg")
	end
	]))=#

function main(window)
	push!(window.assets, "layout2")
	push!(window.assets, "widgets")
	push!(window.assets, "tex")

	tabbar = tabs([ hbox("Orientation"), hbox("Photopolarimetry"), hbox("Interneuron"), hbox("Ellipse")])
	tabcontent = pages([ orientationtab, intensitytab, profiletab, ellipsetab])

	t, p = wire(tabbar, tabcontent, :tab_channel, :selected)
	vbox(t, p)

end


text:
"""
$(TYPEDEF)
`Annealing` type defines the evolution of a time-dependent Hamiltonian in both closed-system and open-system settings. It is called `Annealing` because HOQST started as a toolbox for quantum annealing.

# Fields
$(FIELDS)
"""
mutable struct Annealing{hType,uType} <: AbstractAnnealing{hType,uType}
    "Hamiltonian for the annealing."
    H::hType
    "Initial state for the annealing."
    u0::uType
    "Function of annealing parameter s wrt to t"
    annealing_parameter::Any
    "A system bath interaction set."
    interactions::Union{InteractionSet,Nothing}
end

Evolution = Annealing

function Annealing(
    H::AbstractHamiltonian{T},
    u0;
    coupling = nothing,
    bath = nothing,
    interactions = nothing,
    annealing_parameter = (tf, t) -> t / tf,
) where {T}
    if !isnothing(coupling) && !isnothing(bath)
        if !isnothing(interactions)
            throw(ArgumentError("Both interactions and coupling/bath are specified. Please merge coupling/bath into interactions."))
        end
        interactions = InteractionSet(Interaction(coupling, bath))
    end
    if !(T<:Complex)
    end
    Annealing(H, u0, annealing_parameter, interactions)
end

set_u0!(A::Annealing, u0) = A.u0 = u0
set_annealing_parameter!(A::Annealing, param) = A.annealing_parameter = param

"""
$(TYPEDEF)
Defines a complete set of ODE parameters, which includes Hamiltonian, total annealing time, open system and control objects.
# Fields
$(FIELDS)
"""
struct ODEParams
    # H and opensys may be move to dedicated OpenSysOp in the future
    "Hamiltonian"
    L::Any
    "Total annealing time"
    tf::Real
    "Function to convert physical time to annealing parameter"
    annealing_parameter::Function
    "Annealing control object"
    control::Any
end

ODEParams(L, tf::Real, annealing_param; control = nothing) =
    ODEParams(L, tf, annealing_param, control)
(P::ODEParams)(t::Real) = P.annealing_parameter(P.tf, t)


text:
# Plotting functions

using PyPlot
using Distributions

function create_correlation_figure(data, labels, prior_distributions)

    function plt_corr_vals(ax, m,n,correlation)
        ax.set_axis_bgcolor("r")
        ax.text(0.5, 0.5, "$(round(correlation[m,n],digits=3))", horizontalalignment="center", verticalalignment="center")
    end

    function plt_hist(ax, column, data)
        ax.axis("on")
        col_data = data[:,column]

        ax.hist(col_data, bins=7, density= false)
        ax.locator_params(nbins=1, axis="x")
        ax.locator_params(nbins=1, axis="y")
        ax.get_xaxis().set_visible(true)
    end

    function plt_corr_data(ax, m, n, data)
        ax.axis("on")
        ax.locator_params(nbins=2, axis="x")
        ax.locator_params(nbins=2, axis="y")
        ax.scatter(data[:,n], data[:,m])
        for s in ["top", "right"]
            ax.spines[s].set_visible(false)
        end
        # bound the x and y limits to [0,1]
        for f in [(:get_xlim, :set_xlim), (:get_ylim, :set_ylim)]
            (minval, maxval) =  ax[f[1]]()
            minval = maximum([0., minval])
            maxval = minimum([1., maxval])
            ax[f[2]]([minval, maxval])
        end
    end

    function plt_prior(ax, prior_distrib)
        # plot prior
        xl =  ax.get_xlim()
        xVals = range(xl[1], stop=xl[2], length=30)
        y = [pdf(prior_distrib,x) for x in xVals]
        ax.plot(xVals,y,linewidth=2.0, "g-")
    end

    correlation_matrix = cor(data)
    nr_cols = length(labels)+1
    fig, axes = subplots(nrows=nr_cols, ncols=nr_cols, figsize=(8,nr_cols*1.))
    for m = 1:(nr_cols)
       for n = 1:(nr_cols)
            ax = axes[m,n]
            ax.axis("off")

            if m==1 && n==1
                # nothing
            elseif m==1 && n>1
               ax.text(0.5, 0.5, "$(labels[n-1])", horizontalalignment="center", verticalalignment="baseline")
            elseif n==1 && m>1
               ax.text(0.5, 0.5, "$(labels[m-1])", horizontalalignment="right", verticalalignment="center")
            elseif m==n
                plt_hist(ax, m-1, data)
                plt_prior(ax, prior_distributions[m-1])
            elseif m>n
                plt_corr_data(ax, m-1, n-1, data)
            else
                plt_corr_vals(ax, m-1, n-1,correlation_matrix)
            end
        end
    end
    fig.tight_layout()
end

function plot_result_compact_hist(data::Array{Float64,2}, params, sup_title;
                                nr_bins = 100, add_MAP=true, plot_prior=true,
                                plot_prior_mode = true, print_MAP_in_title = true,
                                title_font_size = 6)
  plot_result(data, params, sup_title;
            nr_bins = nr_bins, add_MAP=add_MAP, plot_prior=plot_prior,
            plot_prior_mode = plot_prior_mode, nr_cols = length(params),
            print_MAP_in_title = print_MAP_in_title)
end
# ---- Plot: re-call get_MAP().
# Attention: set top_n_fraction according to what you want!
# Cause otherwise the values in the .csv files and the values in the figures
# are slightly different.
function plot_result(data::Array{Float64,2}, params, sup_title;
                    nr_bins = 50, add_MAP=true, plot_prior=true, plot_prior_mode = true,
                    nr_cols = 2, print_MAP_in_title = true, title_font_size = 10,
                    top_n_fraction = .01, conditioning_interval=0.08,
                    hist_density = true)

    nr_params = length(params)
    nr_rows = Int(ceil(nr_params/nr_cols))
    fig, axes = subplots(nrows=nr_rows, ncols=nr_cols, figsize=(9.,nr_rows*1.8))
    nr_map_samples = Int(ceil(size(data)[1] * top_n_fraction))

    top_MAP, top_n_avgMAP, top_n_min, top_n_max, top_n_stdMAP, top_n_samples = get_MAP(data, nr_map_samples)
        for i=1:nr_params

        axes[i].set_title(params[i].name,fontsize=title_font_size)
        axes[i].locator_params(nbins=4)

        (n, bins, patches) = axes[i].hist(data[:,i], nr_bins, density = hist_density, edgecolor = "none");

        plot_padding = 0.05 * (maximum(params[i].prior_Distribution) - minimum(params[i].prior_Distribution))
        axes[i].set_xlim([minimum(params[i].prior_Distribution)-plot_padding,maximum(params[i].prior_Distribution)+plot_padding])
        if plot_prior
            xl =  axes[i].get_xlim()
            min_x_prior = max(minimum(params[i].prior_Distribution), (xl[1]+plot_padding))
            max_x_prior = min(maximum(params[i].prior_Distribution), (xl[2]-plot_padding))

            xVals = range(min_x_prior, stop=max_x_prior, length=100)
            y = [pdf(params[i].prior_Distribution,x) for x in xVals]
            axes[i].plot(xVals,y,linewidth=3.0, "g-")

            if plot_prior_mode
                priorD = params[i].prior_Distribution
                priorMode = mode(priorD)
                pdfMode = pdf(priorD, priorMode)
                axes[i].plot([priorMode, priorMode],[0, pdfMode], "g--", linewidth=4.0)
            end

        end
        if add_MAP || print_MAP_in_title

            if add_MAP axes[i].plot(top_MAP[i], 0, "ro", markersize=10) end

            if print_MAP_in_title axes[i].set_title("$(params[i].name) ($(round(top_MAP[i],digits=2)))", fontsize=title_font_size) end
        end
    end
    fig.suptitle(sup_title, y=1.)
    fig.tight_layout()
    return fig
end
# ---- Plot: provide top_MAP, do not re-call get_MAP()
function plot_result(data::Array{Float64,2}, params, sup_title,
                    top_MAP;
                    nr_bins = 50, add_MAP=true, plot_prior=true, plot_prior_mode = true,
                    nr_cols = 2, print_MAP_in_title = true, title_font_size = 10,
                    top_n_fraction = .01, conditioning_interval=0.08,
                    hist_density = true)
    nr_params = length(params)
    nr_rows = Int(ceil(nr_params/nr_cols))
    fig, axes = subplots(nrows=nr_rows, ncols=nr_cols, figsize=(9.,nr_rows*1.8))
    for i=1:nr_params
        axes[i].set_title(params[i].name,fontsize=title_font_size)
        axes[i].locator_params(nbins=4)
        (n, bins, patches) = axes[i].hist(data[:,i], nr_bins, density = hist_density, edgecolor = "none");

        plot_padding = 0.05 * (maximum(params[i].prior_Distribution) - minimum(params[i].prior_Distribution))
        axes[i].set_xlim([minimum(params[i].prior_Distribution)-plot_padding,maximum(params[i].prior_Distribution)+plot_padding])
        if plot_prior
            xl =  axes[i].get_xlim()
            min_x_prior = max(minimum(params[i].prior_Distribution), (xl[1]+plot_padding))
            max_x_prior = min(maximum(params[i].prior_Distribution), (xl[2]-plot_padding))

            xVals = range(min_x_prior, stop=max_x_prior, length=100)
            y = [pdf(params[i].prior_Distribution,x) for x in xVals]
            axes[i].plot(xVals,y,linewidth=3.0, "g-")

            if plot_prior_mode
                priorD = params[i].prior_Distribution
                priorMode = mode(priorD)
                pdfMode = pdf(priorD, priorMode)
                axes[i].plot([priorMode, priorMode],[0, pdfMode], "g--", linewidth=4.0)
            end

        end
        if add_MAP || print_MAP_in_title
            if add_MAP axes[i].plot(top_MAP[i], 0, "ro", markersize=10) end
            if print_MAP_in_title axes[i].set_title("$(params[i].name) ($(round(top_MAP[i],digits=2)))", fontsize=title_font_size) end
        end
    end
    fig.suptitle(sup_title, y=1.)
    fig.tight_layout()
    return fig
end
# ---- Scatterplot LL of top samples: re-call get_MAP().
# Attention: set top_n_fraction according to what you want!
# Cause otherwise the values in the .csv files and the values in the figures
# are slightly different.
function plot_result_logLL(data::Array{Float64,2}, params, sup_title;
                        nr_cols = 2, print_MAP_in_title = true,
                        title_font_size = 10, top_n_fraction = .01
                        )
    nr_params = length(params)
    nr_rows = Int(ceil(nr_params/nr_cols))

    fig, axes = subplots(nrows=nr_rows, ncols=nr_cols, figsize=(9.,nr_rows*1.8))
    nr_map_samples = Int(ceil(size(data)[1] * top_n_fraction))

    top_MAP, top_n_avgMAP, top_n_min, top_n_max, top_n_stdMAP, top_n_samples = get_MAP(data, nr_map_samples)
        for i=1:nr_params
        axes[i].set_title(params[i].name,fontsize=title_font_size)
        axes[i].locator_params(nbins=4)
        plot_padding = 0.05 * (maximum(params[i].prior_Distribution) - minimum(params[i].prior_Distribution))
        axes[i].set_xlim([minimum(params[i].prior_Distribution)-plot_padding,maximum(params[i].prior_Distribution)+plot_padding])
        if print_MAP_in_title axes[i].set_title("$(params[i].name) ($(round(top_MAP[i],digits=2)))", fontsize=title_font_size) end

        axes[i].scatter(top_n_samples[:,i], top_n_samples[:,end], s=3)
    end
    fig.suptitle(sup_title, y=1.)
    fig.tight_layout()
    return fig
end
# ---- Scatterplot LL of top samples:
# provide top_MAP and top_n_samples, do not re-call get_MAP()
function plot_result_logLL(data::Array{Float64,2}, params, sup_title,
                        top_MAP, top_n_samples;
                        nr_cols = 2, print_MAP_in_title = true,
                        title_font_size = 10, top_n_fraction = .01
                        )

    nr_params = length(params)
    nr_rows = Int(ceil(nr_params/nr_cols))
    fig, axes = subplots(nrows=nr_rows, ncols=nr_cols, figsize=(9.,nr_rows*1.8))
        for i=1:nr_params
        axes[i].set_title(params[i].name,fontsize=title_font_size)
        axes[i].locator_params(nbins=4)
        plot_padding = 0.05 * (maximum(params[i].prior_Distribution) - minimum(params[i].prior_Distribution))
        axes[i].set_xlim([minimum(params[i].prior_Distribution)-plot_padding,maximum(params[i].prior_Distribution)+plot_padding])
        if print_MAP_in_title axes[i].set_title("$(params[i].name) ($(round(top_MAP[i],digits=2)))", fontsize=title_font_size) end

        axes[i].scatter(top_n_samples[:,i], top_n_samples[:,end], s=3)
    end
    fig.suptitle(sup_title, y=1.)
    fig.tight_layout()
    return fig
end
function plot_subj_stats(subj_stats::Step_by_step_signals, title_text_part)
        plt.figure(figsize=(9.,1.1))
       LL_subj = round(subj_stats.LL,3)
    probs = subj_stats.action_selection_prob

    legendTxt = []
    if length(subj_stats.RPE) > 0
        plt.plot(subj_stats.RPE, "xr")
        legendTxt = [legendTxt; "RPE"]
    end

    if length(subj_stats.SPE) > 0
        plt.plot(subj_stats.SPE, "ob")
        legendTxt = [legendTxt; "SPE"]
    end

    markerline, stemlines, baseline = plt.stem(probs, "g")
    legendTxt = [legendTxt; "p(a)"]

    plt.legend(legendTxt)
    meanProb = round(mean(probs),digits=3)
    plt.title("subj $(subj_stats.subject_id), $title_text_part, LL=$LL_subj, mean selection prob: $meanProb", fontsize=10)
    plt.locator_params(axis='y',nbins=2)

end

function compare_SarsaL_RPE(data:: Array{Int64,2}, param_population, param_subject)
    println("TODO: move the function compare_SarsaL_RPE")
    LL, signalsPopulation = getLL_SarsaLambda(data, param_population[1], param_population[2], param_population[3], param_population[4], param_population[5])
    LL, signalsSubject = getLL_SarsaLambda(data, param_subject[1], param_subject[2], param_subject[3], param_subject[4], param_subject[5])

    plt.figure(figsize=(9.,1.8))
    markerline, stemlines, baseline = plt.stem(signalsPopulation.RPE, "g")
    plt.setp(markerline, "markerfacecolor", "g", "markeredgecolor", "g", "markersize", 4)
    markerline, stemlines, baseline = plt.stem(signalsSubject.RPE, "r--")
    plt.setp(markerline, "markerfacecolor", "r", "markeredgecolor", "r", "markersize", 4)
    plt.setp(baseline, "color", "k", "linewidth", 1)

    plt.legend(["RPE prior", "RPE subj"], fontsize =7)
    dot_prod= dot(signalsPopulation.RPE, signalsSubject.RPE)
    prod_norms = norm(signalsPopulation.RPE) * norm(signalsSubject.RPE)
    cos_Theta = dot_prod/prod_norms
    plt.title("cos(pop, subj)=$(round(cos_Theta,digits=3))", fontsize=10)

end


text:
# Methods for the Conic interface

abstract type AbstractConicModel <: AbstractMathProgModel end
export AbstractConicModel

@define_interface begin
    ConicModel
    getdual
    getvardual
    supportedcones
    setbvec!
end


text:
export CMAEvolutionStrategyOpt
struct CMAEvolutionStrategyOpt end

function __map_optimizer_args(prob::OptimizationProblem, opt::CMAEvolutionStrategyOpt;
    cb=nothing,
    maxiters::Union{Number, Nothing}=nothing,
    maxtime::Union{Number, Nothing}=nothing,
    abstol::Union{Number, Nothing}=nothing,
    reltol::Union{Number, Nothing}=nothing,
    kwargs...)

    if !isnothing(reltol)
        @warn "common reltol is currently not used by $(opt)"
    end

    mapped_args = (;lower = prob.lb,
    upper = prob.ub,
    kwargs...)

    if !isnothing(maxiters)
        mapped_args = (; mapped_args..., maxiter=maxiters)
    end

    if !isnothing(maxtime)
        mapped_args = (; mapped_args..., maxtime=maxtime)
    end

    if !isnothing(abstol)
        mapped_args = (; mapped_args..., ftol=abstol)
    end

    return mapped_args
end

function __solve(prob::OptimizationProblem, opt::CMAEvolutionStrategyOpt, data = DEFAULT_DATA;
                 cb = (args...) -> (false),
                 maxiters::Union{Number, Nothing} = nothing,
                 maxtime::Union{Number, Nothing} = nothing,
                 abstol::Union{Number, Nothing}=nothing,
                 reltol::Union{Number, Nothing}=nothing,
                 progress = false, kwargs...)
    local x, cur, state

    if data != DEFAULT_DATA
        maxiters = length(data)
    end

    cur, state = iterate(data)

    function _cb(trace)
        cb_call = cb(decompose_trace(trace).metadata["x"],trace.value...)
        if !(typeof(cb_call) <: Bool)
            error("The callback should return a boolean `halt` for whether to stop the optimization process.")
        end
        cur, state = iterate(data, state)
        cb_call
    end

    maxiters = _check_and_convert_maxiters(maxiters)
    maxtime = _check_and_convert_maxtime(maxtime)

    _loss = function(θ)
        x = prob.f(θ, prob.p, cur...)
        return first(x)
    end

    opt_args = _map_optimizer_args(prob, opt, cb=_cb, maxiters=maxiters, maxtime=maxtime,abstol=abstol, reltol=reltol; kwargs...)

    t0 = time()
    opt_res = CMAEvolutionStrategy.minimize(_loss, prob.u0, 0.1; opt_args...)
    t1 = time()

    opt_ret = opt_res.stop.reason

    SciMLBase.build_solution(prob, opt, opt_res.logger.xbest[end], opt_res.logger.fbest[end]; original=opt_res, retcode=opt_ret)
end


text:
using DataFrames
using Effects
using GLM
using LinearAlgebra
using StableRNGs
using StandardizedPredictors
using Statistics
using StatsModels
using Test

b0, b1, b2, b1_2 = beta = [0.0, 1.0, 1.0, -1.0]

@testset "simple" begin
    x = collect(-10:10)
    dat = (; :x => x, :y => x .* b1 .+ b0 + randn(StableRNG(42), length(x)))
    model = lm(@formula(y ~ x), dat)

    design = Dict(:x => 1:20)
    eff = effects(design, model)
    # test effect
    @test eff.y ≈ eff.x .* last(coef(model)) .+ first(coef(model))
    # test error
    pred = [1 first(eff.x)]
    # if we drop support for Julia < 1.4, this first can become only
    @test first(eff.err) ≈ first(sqrt(pred * vcov(model) * pred'))
    # test CI
    @test eff.lower ≈ eff.y - eff.err
    @test eff.upper ≈ eff.y + eff.err

    @testset "contrasts" begin
        # need to check this works when the response has also
        # had a name change
        contrasts = Dict(:x => Center(10), :y => ZScore())
        model_centered = lm(@formula(y ~ x), dat; contrasts=contrasts)

        eff_centered = effects(design, model_centered)

        μ = mean(dat.y)
        σ = std(dat.y)
        # different contrasts shouldn't affect the predictions/effects
        @test eff_centered[!, "y(centered: 0.21 scaled: 6.27)"] ≈ zscore(eff.y, µ, σ)
        @test eff_centered.lower ≈ zscore(eff.lower, µ, σ)
        @test eff_centered.upper ≈ zscore(eff.upper, µ, σ)
    end
end

@testset "multiple" begin
    x = collect(-10:10)
    z = @. sqrt(x + 100)
    dat = DataFrame(; x=x, z=z, y=randn(StableRNG(42), length(x)))

    @testset "additive" begin
        form = @formula(y ~ 1 + x + z)
        dat[!, :y] += modelmatrix(form, dat) * beta[1:(end - 1)]
        model = lm(form, dat)

        design = Dict(:x => 1:20)
        maximin(v) = mean(extrema(v))
        for typical in (maximin, mean)
            eff = effects(design, model; typical=typical)
            # test effect
            intercept = ones(length(design[:x]))
            z_typical = typical(z) * intercept
            bs = coef(model)
            @test eff.y ≈ @. bs[1] + eff.x * bs[2] + bs[3] * z_typical
            # test error
            pred = [intercept eff.x z_typical]
            @test eff.err ≈ sqrt.(diag(pred * vcov(model) * pred'))
            # test CI
            @test eff.lower ≈ eff.y - eff.err
            @test eff.upper ≈ eff.y + eff.err
        end
    end

    @testset "multiplicative" begin
        form = @formula(y ~ 1 + x * z)
        dat[!, :y] += modelmatrix(form, dat) * beta
        model = lm(form, dat)

        design = Dict(:x => 1:20)
        maximin(v) = mean(extrema(v))
        for typical in (maximin, mean)
            eff = effects(design, model; typical=typical)
            # test effect
            intercept = ones(length(design[:x]))
            z_typical = typical(z) * intercept
            zx_typical = z_typical .* design[:x]
            bs = coef(model)
            @test eff.y ≈ @. bs[1] + eff.x * bs[2] + bs[3] * z_typical + bs[4] * zx_typical
            # test error
            pred = [intercept eff.x z_typical zx_typical]
            @test eff.err ≈ sqrt.(diag(pred * vcov(model) * pred'))
            # test CI
            @test eff.lower ≈ eff.y - eff.err
            @test eff.upper ≈ eff.y + eff.err
        end
    end
end


text:
# new exports from this file, WIP

export calcDyadScaleAdaptive
export covEllipseParameterized, plotCovEllipseLayer
export plotSLAM2D_KeyAndRef

## Src

global DISABLESTBPRTLINES = false

function togglePrtStbLines()
  global DISABLESTBPRTLINES
  DISABLESTBPRTLINES = !DISABLESTBPRTLINES
end

function plotLsrScanFeats(br::Array{Float64,2})
  Cart = zeros(size(br))
  Cart[:,1] = br[:,2].*cos(br[:,1])
  Cart[:,2] = br[:,2].*sin(br[:,1])
  plot(x=Cart[:,1],y=Cart[:,2],Geom.point,
  Guide.xticks(ticks=collect(-60:10:60)),
  Guide.yticks(ticks=collect(0:10:80)))
end

function plotFeatTrackers(trkrs::Dict{Int64,Feature}, bfts::Array{Float64,2})
  musX = Float64[]
  varX = Float64[]
  musY = Float64[]
  varY = Float64[]
  allPtsX = Float64[]
  allPtsY = Float64[]

  for ftr in trkrs
    pts = getPoints(ftr[2].bel)
    allPtsX = [allPtsX; vec(pts[1,:])]
    allPtsY = [allPtsY; vec(pts[2,:])]

    push!(musX, Statistics.mean(vec(pts[1,:])))
    push!(varX, Statistics.std(vec(pts[1,:])))
    push!(musY, Statistics.mean(vec(pts[2,:])))
    push!(varY, Statistics.std(vec(pts[2,:])))
  end

  X = Float64[]
  Y = Float64[]

  if size(bfts,2) > 0
    if bfts[1,1] != 0.0 && bfts[2,1] != 0.0 && bfts[3,1] != 0.0
      for i in 1:size(bfts,2)
          u, R = p2c(vec(bfts[:,i]))
          push!(X, u[1])
          push!(Y, u[2])
      end
    end
  end

  # Guide.yticks(ticks=collect(-60:10:60)),
  # Guide.xticks(ticks=collect(0:10:80))
  p = plot(layer(x=musX, y=musY, Geom.point, Theme(default_color=colorant"red")),
  layer(x=allPtsX, y=allPtsY, Geom.histogram2d),
  Guide.yticks(ticks=collect(-70:10:70)),
  Guide.xticks(ticks=collect(-40:10:80)))
  for i in 1:length(X)
    push!(p.layers, Gadfly.layer(x=[0.0;X[i]], y=[0.0;Y[i]], Geom.line, Gadfly.Theme(default_color=colorant"magenta"))[1])
  end
  p
end

function saveImgSeq(d::Dict{Int64,Array{Float64,2}}; from::Int=1,to::Int=10,step::Int=1)
  for i in from:step:to
    p = plotLsrScanFeats(lsrBR(d[i]));
    Gadfly.draw(PNG(string("imgs/img",i,".png"),25cm,25cm),p)
  end
  nothing
end

"""
    $SIGNATURES

Calculate the ellipse from covariance matrix and return as lambda function.

Notes
- https://cookierobotics.com/007/

Related

plotCovEllipseLayer
"""
function covEllipseParameterized(distr::MvNormal; meanOffset::Bool=true)
  println(round.(cov(distr), digits=3) )
  a = cov(distr)[1,1]
  b = cov(distr)[1,2] + cov(distr)[2,1]
  b *= 0.5
  c = cov(distr)[2,2]

  λ1 = 0.5*(a+c) + sqrt(0.25*(a-c)^2 + b^2)
  λ2 = 0.5*(a+c) - sqrt(0.25*(a-c)^2 + b^2)
  θ = if b ≈ 0 && a >= c
    0
  elseif b ≈ 0 && a < c
    pi/2
  else
    atan(λ1-a, b)
  end
  sλ1 = sqrt(λ1)
  sλ2 = sqrt(λ2)
  sθ =  sin(θ)
  cθ =  cos(θ)
  ox = meanOffset ? distr.μ[1] : 0
  oy = meanOffset ? distr.μ[2] : 0
  (t) -> [sλ1*cθ*cos(t)-sλ2*sθ*sin(t)+ox; sλ1*sθ*cos(t)+sλ2*cθ*sin(t)+oy]
end

covEllipseParameterized(pts::Array{Float64,2}; meanOffset::Bool=true) = covEllipseParameterized( fit(MvNormal, pts), meanOffset=meanOffset )
covEllipseParameterized(X::Union{<:BallTreeDensity,<:ManifoldKernelDensity}; meanOffset::Bool=true) = covEllipseParameterized( getPoints(X), meanOffset=meanOffset )
covEllipseParameterized(dfg::AbstractDFG, sym::Symbol; meanOffset::Bool=true, solveKey::Symbol=:default) = covEllipseParameterized( getKDE(dfg, sym, solveKey), meanOffset=meanOffset, solveKey=solveKey )

"""
    $SIGNATURES

Plotting tool to draw Gadfly layers of ellipses of 2D covariance fitted to the belief of factor graph variable nonparametric points.

Related

covEllipseParameterized, plotSLAM2DPosesLandms
"""
function plotCovEllipseLayer( dfg::AbstractDFG,
                              vsym::Symbol;
                              solveKey::Symbol=:default,
                              drawPoints::Bool=true,
                              ellipseColor::AbstractString="gray30",
                              pointsColor::AbstractString="gray30",
                              drawEllipse::Bool=true  )
  #
  PL = []

  if !(drawEllipse || drawPoints)
    return PL
  end

  # points to work from
  bel = getBelief(dfg, vsym, solveKey)
  pp__ = getPoints(bel)
  pp_ = if bel.manifold isa SpecialEuclidean
    # assume Pose2
    (x->x.parts[1]).(pp__)
  else
    # assume Point2
    pp__
  end

  @cast pp[i,j] := pp_[j][i]

  if drawEllipse
    # get ellipse function
    eX = covEllipseParameterized(pp[1:2,:], meanOffset=false)
    vEl = eX.(0:0.02:2pi)
    el = [(x->x[1]).(vEl) (x->x[2]).(vEl)]
    # add suggested PPE mean offset
    el[:,1] .+= getVariablePPE(dfg, vsym, solveKey).suggested[1]
    el[:,2] .+= getVariablePPE(dfg, vsym, solveKey).suggested[2]

    # add the ellipse layers
    plelX2 = Gadfly.layer(x=el[:,1], y=el[:,2], Geom.path, Theme(default_color=parse(Colorant, ellipseColor)))
    push!(PL, plelX2[1])
  end

  # add the points layer if needed
  if drawPoints
    plelX1 = Gadfly.layer(x=pp[1,:],
                          y=pp[2,:],
                          Geom.point,
                          Theme(default_color=parse(Colorant, pointsColor),
                          point_size=1pt))
    #
    push!(PL, plelX1[1])
  end

  return PL
end

"""
    $SIGNATURES

Plot trajectory of Array{,2} with rows as consecutive entries and columns as x,y,theta.
"""
function plotTrajectoryArrayPose2(arr::AbstractMatrix__{<:Real};
                                  spscale::Real=0.5,
                                  triadStride::Int=50)
  #
  @assert size(arr,2)==3
  trajPlt = Gadfly.plot(x=arr[:,1], y=arr[:,2], Geom.path, Coord.cartesian(fixed=true, aspect_ratio=1))

  if triadStride != -1
    Xpp = arr[1:triadStride:end,1]
    Ypp = arr[1:triadStride:end,2]
    Thpp = arr[1:triadStride:end,3]
    addXYLineLayers!(trajPlt, Xpp, Ypp, Thpp, l=spscale)
  end
  return trajPlt
end

## TODO -- you were here with port starboard lines
function stbPrtLineLayers!(pl, Xpp, Ypp, Thpp; l::Real=5.0)
    if DISABLESTBPRTLINES
      return nothing
    end

    lnstpr = [0.0;l;0.0]
    lnstpg = [0.0;-l;0.0]

    Rd  =SE2(lnstpr)
    Gr = SE2(lnstpg)

    for i in 1:length(Xpp)
      lnstt = [Xpp[i];Ypp[i];Thpp[i]]
      Ps = SE2(lnstt)
      lnr = se2vee(Ps*Rd)
      lng = se2vee(Ps*Gr)
      xsr = [Xpp[i];lnr[1]]
      ysr = [Ypp[i];lnr[2]]
      xsg = [Xpp[i];lng[1]]
      ysg = [Ypp[i];lng[2]]

      push!(pl.layers, layer(x=xsr, y=ysr, Geom.path(), Gadfly.Theme(default_color=colorant"red", line_width=1.5pt))[1] )
      push!(pl.layers, layer(x=xsg, y=ysg, Geom.path(), Gadfly.Theme(default_color=colorant"green", line_width=1.5pt))[1] )
    end
    nothing
end

# draw the reference frame as a red-green dyad
function addXYLineLayers!(pl, Xpp, Ypp, Thpp; l::Real=1.0, manualColor::Union{Nothing, AbstractString}=nothing  )
    lnstpr = [l;0.0;0.0]
    lnstpg = [0.0;l;0.0]

    Rd  =SE2(lnstpr)
    Gr = SE2(lnstpg)

    for i in 1:length(Xpp)
      lnstt = [Xpp[i];Ypp[i];Thpp[i]]
      Ps = SE2(lnstt)
      lnr = se2vee(Ps*Rd)
      lng = se2vee(Ps*Gr)
      xsr = [Xpp[i];lnr[1]]
      ysr = [Ypp[i];lnr[2]]
      xsg = [Xpp[i];lng[1]]
      ysg = [Ypp[i];lng[2]]

      push!(pl.layers, layer(x=xsr, y=ysr, Geom.path(), Gadfly.Theme(default_color=parse(Colorant, manualColor === nothing ? "red" : manualColor), line_width=1.5pt))[1] )
      push!(pl.layers, layer(x=xsg, y=ysg, Geom.path(), Gadfly.Theme(default_color=parse(Colorant, manualColor === nothing ? "green" : manualColor), line_width=1.5pt))[1] )
    end
    nothing
end

"""
    $SIGNATURES
Adaptively size the `dyadScale` value for plotSLAM2DPose.
"""
function calcDyadScaleAdaptive( Xpp::AbstractVector{<:Real},
                                Ypp::AbstractVector{<:Real};
                                scaling::Real=0.25)

  dists = diff(Xpp).^2 + diff(Ypp).^2 |> vec .|> sqrt
  scaling*Statistics.mean(dists)
end

"""
    $(SIGNATURES)

2D plot of all poses, assuming poses are labeled from ``::Symbol` type `:x0, :x1, ..., :xn`.  Use `to` and `from` to limit the range of numbers `n` to be drawn.  The underlying histogram can be enabled or disabled, and the size of maximum-point belief estimate cursors can be controlled with `spscale`.

Future:
- Relax to user defined pose labeling scheme, for example `:p1, :p2, ...`
"""
function plotSLAM2DPoses( fg::AbstractDFG;
                          solveKey::Symbol=:default,
                          regexPoses=r"x\d",
                          from::Int=0,
                          to::Int=(2^(Sys.WORD_SIZE-1)-1),
                          variableList::AbstractVector{Symbol}=getVariablesLabelsWithinRange(fg, regexPoses, from=from, to=to),
                          meanmax=:null,
                          ppe=:suggested,
                          recalcPPEs::Bool=false,
                          lbls=true,
                          scale::Real=1,
                          x_off::Real=0,
                          y_off::Real=0,
                          drawhist=false,
                          spscale::Union{Nothing, <:Real}=nothing,
                          dyadScale::Union{Nothing, <:Real}=nothing,
                          drawTriads::Bool=true,
                          drawContour::Bool=true, levels::Int=1,
                          contour::Union{Nothing, Bool}=nothing,
                          line_width=1pt,
                          drawPoints::Bool=true,
                          pointsColor::AbstractString="gray30",
                          drawEllipse::Bool=false,
                          ellipseColor::AbstractString="gray30",
                          manualColor=nothing  )
    #
    # deprecations
    if meanmax != :null
      @warn "plotSLAM2DPoses meanmax keyword is deprecated, use ppe=:suggested instead."
      ppe = meanmax
    end
    !(spscale isa Nothing) ? (dyadScale=spscale; @warn("keyword spscale is deprecated, use dyadScale instead")) : nothing
    !(contour isa Nothing) ? (drawContour=contour; @warn("keyword contour is being deprecated, use drawContour instead")) : nothing

    ## Use PPE.suggested

    Ppes = if recalcPPEs
      map(x->calcVariablePPE(fg, x, solveKey=solveKey), variableList)
    else
      getPPE.(fg, variableList, solveKey)
    end
    mask = Ppes .!== nothing
    variableList = variableList[mask]
    suggPpes = (x->getfield(x,ppe)).(Ppes[mask])

    Xpp::Vector{Float64}  = Float64.( (x->x[1]).(suggPpes) )
    Ypp::Vector{Float64}  = Float64.( (x->x[2]).(suggPpes) )
    Thpp::Vector{Float64} = Float64.( (x->x[3]).(suggPpes) )
    LBLS = string.(variableList)

    Xpp .+= x_off
    Ypp .+= y_off

    Xpp .*= scale
    Ypp .*= scale

    # adaptively scale dyad size
    dyadScale = dyadScale isa Nothing ? calcDyadScaleAdaptive(Xpp, Ypp) : dyadScale

    # lbls = lblsFromTo(1,length(Xpp))
    psplt = Union{}
    thm = manualColor === nothing ? Theme(line_width=1pt) : Theme(line_width=line_width, default_color=parse(Colorant, manualColor))
    if lbls
      psplt = Gadfly.plot(
        Gadfly.layer(x=Xpp,y=Ypp,label=LBLS,Geom.path(), thm, Geom.label),
        Coord.cartesian(fixed=true)
      )
    else
      psplt = Gadfly.plot(
        Gadfly.layer(x=Xpp,y=Ypp,Geom.path(), thm),Coord.cartesian(fixed=true),
        Coord.cartesian(fixed=true)
      )
    end
    # return psplt
    drawTriads && addXYLineLayers!(psplt, Xpp, Ypp, Thpp, l=dyadScale, manualColor=manualColor)
    if drawhist
      Xp,Yp = get2DPoseSamples(fg, from=from, to=to)
      push!(psplt.layers,  Gadfly.layer(x=Xp, y=Yp, Geom.histogram2d)[1] )#(xbincount=100, ybincount=100))
    end
    # add contours to pose estimates
      # varsyms = Symbol.(LBLS)
    if drawContour
      for vsym in variableList
        pln = plotKDE(fg, vsym, solveKey=solveKey, dims=[1;2], levels=levels, c=[(manualColor === nothing ? "gray90" : manualColor);])
        union!(psplt.layers, pln.layers)
      end
    end
    # drawEllipse
    for vsym in variableList
      pln = plotCovEllipseLayer(fg, vsym, solveKey=solveKey, drawEllipse=drawEllipse,drawPoints=drawPoints,ellipseColor=ellipseColor,pointsColor=pointsColor)
      union!(psplt.layers, pln)
    end
    return psplt
end

"""
    $(SIGNATURES)

2D plot of landmarks, assuming `:l1, :l2, ... :ln`.  Use `from` and `to` to control the range of landmarks `n` to include.
"""
function plotSLAM2DLandmarks( fg::AbstractDFG;
                              solveKey::Symbol=:default,
                              regexLandmark::Regex=r"l",
                              from::Int=0, to::Int=(2^(Sys.WORD_SIZE-1)-1),
                              minnei::Int=0,
                              variableList::AbstractVector{Symbol}=getVariablesLabelsWithinRange(fg, regexLandmark, from=from, to=to),
                              meanmax=:null,
                              ppe::Symbol=:suggested,
                              recalcPPEs::Bool=false,
                              lbls=true,showmm=false,
                              scale::Real=1,
                              x_off::Real=0,
                              y_off::Real=0,
                              drawhist=false,
                              drawContour::Bool=true, levels::Int=1,
                              contour::Union{Nothing, Bool}=nothing,
                              manualColor=nothing,
                              c= manualColor===nothing ? "red" : manualColor,
                              MM::Dict{Int,T}=Dict{Int,Int}(),
                              point_size=1pt,
                              drawPoints::Bool=true,
                              pointsColor::AbstractString="gray30",
                              drawEllipse::Bool=false,
                              ellipseColor::AbstractString="gray30",
                              resampleGaussianFit::Int=0  ) where T
    #
    if meanmax != :null
      @warn "plotSLAM2DPoses meanmax keyword is deprecated, use ppe instead."
      ppe = meanmax
    end

    !(contour isa Nothing) ? (drawContour=contour; @warn("keyword contour is being deprecated, use drawContour instead")) : nothing

    ## Use PPE.suggested

    Ppes = if recalcPPEs
      map(x->calcVariablePPE(fg, x, solveKey=solveKey), variableList)
    else
      getPPE.(fg, variableList, solveKey)
    end
    mask = Ppes .!== nothing
    variableList = variableList[mask]
    suggPpes = (x->getfield(x,ppe)).(Ppes[mask])

    Xpp  = (x->x[1]).(suggPpes)
    Ypp  = (x->x[2]).(suggPpes)
    lbltags = string.(variableList)

    Xpp .+= x_off
    Ypp .+= y_off

    Xpp .*= scale
    Ypp .*= scale

    # Xp,Yp = get2DLandmSamples(fg, from=from, to=to)
    # Xpp = Float64[]; Ypp=Float64[]; Thpp=Float64[]; lblstags=String[];
    # # TODO transition to new PPE.suggested
    # if meanmax==:mean
    #   Xpp,Ypp, t, lbltags = get2DLandmMeans(fg, from=from, to=to, regexLandmark=regexLandmark)
    # elseif meanmax==:max
    #   Xpp,Ypp, t, lbltags = get2DLandmMax(fg, from=from, to=to,showmm=showmm,MM=MM, regexLandmark=regexLandmark)
    # end

    psplt = if lbls
      Gadfly.plot(
        Gadfly.layer(x=Xpp,y=Ypp, label=lbltags, Geom.point, Theme(line_width=1pt, default_color=parse(Colorant,c), point_size=point_size), Geom.label),
        Coord.cartesian(fixed=true)
        # ,Gadfly.layer(x=Xp, y=Yp, Geom.histogram2d)#(xbincount=100, ybincount=100)
      )
    else
      Gadfly.plot(
        Gadfly.layer(x=Xpp,y=Ypp, Geom.point, Theme(line_width=1pt, default_color=parse(Colorant,c), point_size=point_size)),
        Coord.cartesian(fixed=true)
      )
    end

    if drawhist
      push!(psplt.layers, Gadfly.layer(x=Xpp, y=Ypp, Geom.histogram2d)[1]) #(xbincount=100, ybincount=100)
    end

    if drawContour
      # make pretty near Gaussian contours?
      if resampleGaussianFit != 0
        #
      end
      varsyms = Symbol.(lbltags)
      for vsym in varsyms
        pln = plotKDE(fg, vsym, solveKey=solveKey, dims=[1;2], levels=levels, c=[(manualColor===nothing ? "gray90" : manualColor);])
        union!(psplt.layers, pln.layers)
      end
    end

    # if drawEllipse
    for vsym in variableList
      pln = plotCovEllipseLayer(fg, vsym, solveKey=solveKey, drawEllipse=drawEllipse,drawPoints=drawPoints,ellipseColor=ellipseColor,pointsColor=pointsColor)
      union!(psplt.layers, pln)
    end

    psplt
end

"""
    $(SIGNATURES)

2D plot of both poses and landmarks contained in factor graph.  Assuming poses and landmarks
are labeled `:x1, :x2, ...` and `:l0, :l1, ...`, respectively.  The range of numbers to
include can be controlled with `from` and `to` along with other keyword functionality for
manipulating the plot.

Notes
- Assumes `:l1`, `:l2`, ... for landmarks --
- Can increase default Gadfly plot size (for JSSVG in browser): `Gadfly.set_default_plot_size(35cm,20cm)`.
- Enable or disable features such as the covariance ellipse with keyword `drawEllipse=true`.

DevNotes
- TODO update to use e.g. `tags=[:LANDMARK]`,
- TODO fix `drawHist`,
- TODO deprecate, `showmm`, `spscale`.

Examples:
```julia
fg = generateCanonicalFG_Hexagonal()
plotSLAM2D(fg)
plotSLAM2D(fg, drawPoints=false)
plotSLAM2D(fg, contour=false, drawEllipse=true)
plotSLAM2D(fg, contour=false, title="SLAM result 1")

# or load a factor graph
fg_ = loadDFG("somewhere.tar.gz")
plotSLAM2D(fg_)
```

Related

[`plotSLAM2DPoses`](@ref), [`plotSLAM2DLandmarks`](@ref), [`plotPose`](@ref), [`plotKDE`](@ref)
"""
function plotSLAM2D(fgl::AbstractDFG;
                    solveKey::Symbol=:default,
                    from::Int=0, to::Int=(2^(Sys.WORD_SIZE-1)-1),
                    minnei::Int=0,
                    meanmax=:null,
                    posesPPE=:suggested,
                    landmsPPE=:suggested,
                    recalcPPEs::Bool=false,
                    lbls=true,
                    scale::Real=1,
                    x_off::Real=0,
                    y_off::Real=0,
                    drawTriads::Bool=true,
                    spscale::Union{Nothing, <:Real}=nothing,
                    dyadScale::Union{Nothing, <:Real}=nothing,
                    levels::Int=1,
                    drawhist=false, MM::Dict{Int,T}=Dict{Int,Int}(),
                    xmin=nothing, xmax=nothing, ymin=nothing, ymax=nothing,
                    showmm=true,
                    window::Union{Nothing, Tuple{Symbol, Real}}=nothing,
                    point_size=4pt,
                    line_width=1pt,
                    regexLandmark=r"l\d+",
                    regexPoses=r"x\d+",
                    manualColor=nothing,
                    drawPoints::Bool=true,
                    pointsColor::AbstractString="gray30",
                    drawContour::Bool=true,
                    drawEllipse::Bool=false,
                    ellipseColor::AbstractString="gray30",
                    contour::Union{Nothing,Bool}=nothing,
                    title::AbstractString=""  ) where T
  #
  # deprecations
  if meanmax != :null
    @warn "plotSLAM2DPosesLandms meanmax keyword is deprecated, use posesPPE or landmsPPE instead."
    posesPPE = meanmax
    landmsPPE = meanmax
  end
  !(spscale isa Nothing) ? (dyadScale=spscale; @warn("keyword spscale is being deprecated, use dyadScale instead")) : nothing
  !(contour isa Nothing) ? (drawContour=contour; @warn("keyword contour is being deprecated, use drawContour instead")) : nothing

  #
  xmin !== nothing && xmax !== nothing && xmin == xmax ? error("xmin must be less than xmax") : nothing
  ymin !== nothing && ymax !== nothing && ymin == ymax ? error("ymin must be less than ymax") : nothing
  ll = listVariables(fgl, regexLandmark)
  p = plotSLAM2DPoses(fgl,
                      solveKey=solveKey,
                      from=from,
                      to=to,
                      regexPoses=regexPoses,
                      ppe=posesPPE,
                      lbls=lbls,
                      scale=scale,
                      x_off=x_off,
                      y_off=y_off,
                      drawhist=drawhist,
                      dyadScale=dyadScale,
                      drawContour=drawContour,
                      drawTriads=drawTriads,
                      manualColor=manualColor,
                      line_width=line_width,
                      drawPoints=drawPoints,
                      ellipseColor=ellipseColor,
                      pointsColor=pointsColor,
                      drawEllipse=drawEllipse,
                      recalcPPEs=recalcPPEs  )
  #
  if length(ll) > 0
    pl = plotSLAM2DLandmarks( fgl,
                              solveKey=solveKey,
                              from=from,
                              to=to,
                              ppe=landmsPPE,
                              minnei=minnei,
                              lbls=lbls,
                              scale=scale,
                              x_off=x_off,
                              y_off=y_off,
                              drawhist=drawhist,
                              MM=MM,
                              showmm=showmm,
                              point_size=point_size,
                              drawContour=drawContour,
                              manualColor=manualColor,
                              drawPoints=drawPoints,
                              ellipseColor=ellipseColor,
                              pointsColor=pointsColor,
                              drawEllipse=drawEllipse,
                              recalcPPEs=recalcPPEs  )
    #
    for l in pl.layers
      push!(p.layers, l)
    end
  end
  if window !== nothing
    focusX = getPPE(fgl, window[1], solveKey).max # getKDEMax( getBelief(getVariable(fgl,window[1]),solveKey) )
    pwind = window[2]
    p.coord = Coord.cartesian(xmin=focusX[1]-pwind,xmax=focusX[1]+pwind,ymin=focusX[2]-pwind,ymax=focusX[2]+pwind)
  end
  co = Coord.Cartesian(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax)
  p.coord = co
  if title != ""
    push!(p.guides, Guide.title(title))
  end
  return p
end

function plotSLAM2DSubmaps( fgl::AbstractDFG, fromto::Array{Int,2};
                            m1hist=false, m2hist=false, m3hist=false,
                            showmm=false, MM::Dict{Int,T} = Dict{Int,Any}(),
                            xmin=nothing, xmax=nothing, ymin=nothing, ymax=nothing ) where T
  #
  p = plotSLAM2DLandmarks(fgl, from=fromto[1,1], to=fromto[1,2], drawhist=m1hist, showmm=showmm, MM=MM)
  if size(fromto,1) >1
    p2 = plotSLAM2DLandmarks(fgl, from=fromto[2,1], to=fromto[2,2], drawhist=m2hist,c="blue", showmm=showmm, MM=MM)
    for l in p2.layers
      push!(p.layers, l)
    end
  end
  if size(fromto,1) >2
    p3 = plotSLAM2DLandmarks(fgl, from=fromto[3,1], to=fromto[3,2], drawhist=m3hist,c="magenta", showmm=showmm, MM=MM)
    for l in p3.layers
      push!(p.layers, l)
    end
  end
  co = Coord.Cartesian(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax)
  p.coord = co
  return p
end

function plotSLAM2DSubmaps( fgl::G, fromto::Array{Int,1}; spread::Int=25,
                            m1hist=false, m2hist=false, m3hist=false,
                            showmm=false, MM::Dict{Int,T}=Dict{Int,Any}(),
                            xmin=nothing, xmax=nothing, ymin=nothing, ymax=nothing ) where {G <: AbstractDFG, T}
  #
  ft = zeros(Int,length(fromto),2)
  for i in 1:length(fromto)
    ft[i,1] = fromto[i]-spread; ft[i,2] = fromto[i]+spread;
  end
  plotSLAM2DSubmaps(fgl, ft, m1hist=m1hist, m2hist=m2hist, m3hist=m3hist, showmm=showmm, MM=MM, xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax)
end

function  plotSLAM2D_KeyAndRef(fg, solveKey=:default; refkey=:simulated, recalcPPEs=true)
  #
  pl_ = plotSLAM2D(fg, solveKey=refkey, drawContour=false, drawPoints=false, drawEllipse=false, manualColor="black", drawTriads=false);
  pl1 = plotSLAM2D(fg, solveKey=solveKey, recalcPPEs=recalcPPEs, drawContour=false, drawPoints=false, drawEllipse=false);

  union!(pl1.layers, pl_.layers);

  pl1
end

# import RoMEPlotting: drawMarginalContour

function plotMarginalContour( fgl::AbstractDFG, lbl::String;
                              solveKey::Symbol=:default,
                              xmin=-150,xmax=150,ymin=-150,ymax=150,
                              n::Int=200 )
  #
  p = getKDE(getVariable(fgl,Symbol(lbl)), solveKey)  # p = getKDE(getVert(fgl,lbl))
  Gadfly.plot(z=(x,y)->evaluateDualTree(p,vectoarr2([x,y]))[1],
    x=collect(range(xmin,stop=xmax,length=n)),
    y=collect(range(ymin,stop=ymax,length=n)),
    Geom.contour,
    Coord.Cartesian(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax),
    Guide.title(lbl)
  )
end

function accumulateMarginalContours(fgl, order;
                                    solveKey::Symbol=:default,
                                    xmin=-150,xmax=150,ymin=-150,ymax=150,n=200 )
  #
  pl = plotMarginalContour(fgl, order[1],solveKey=solveKey, xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,n=n)
  pl2 = nothing
  PL = []
  for or in order[1:end]
    pl2 = plotMarginalContour(fgl, or, solveKey=solveKey, xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,n=n)
    push!(PL, pl2)
    push!(pl.layers, pl2.layers[1])
  end
  return pl, PL
end

#


text:
include("utils.jl")

@testset "Euclidean" begin
    E = Manifolds.Euclidean(3)
    EM = Manifolds.MetricManifold(E,Manifolds.EuclideanMetric())
    @test is_default_metric(EM) == Val{true}()
    @test is_default_metric(E,Manifolds.EuclideanMetric()) == Val{true}()
    x = zeros(3)
    @test det_local_metric(EM,x) == one(eltype(x))
    @test log_local_metric_density(EM,x) == zero(eltype(x))
    @test project_point!(E,x) == x

    manifolds = [ E, EM ]
    types = [Vector{Float64},
             SizedVector{3, Float64},
             MVector{3, Float64},
             Vector{Float32},
             SizedVector{3, Float32},
             MVector{3, Float32},
             Vector{Double64},
             MVector{3, Double64},
             SizedVector{3, Double64}]
    for M in manifolds
        for T in types
            @testset "$M Type $T" begin
                pts = [convert(T, [1.0, 0.0, 0.0]),
                       convert(T, [0.0, 1.0, 0.0]),
                       convert(T, [0.0, 0.0, 1.0])]
                test_manifold(M,
                              pts,
                              test_reverse_diff = isa(T, Vector),
                              test_project_tangent = true,
                              test_musical_isomorphisms = true,
                              test_vector_transport = true,
                              test_mutating_rand = isa(T, Vector),
                              point_distributions = [Manifolds.projected_distribution(M, Distributions.MvNormal(zero(pts[1]), 1.0))],
                              tvector_distributions = [Manifolds.normal_tvector_distribution(M, pts[1], 1.0)])
            end
        end
    end
end


text:
using LPVSpectral
using Test, LinearAlgebra, Statistics, Random, StatsBase, ProximalOperators
using Plots, DSP
Random.seed!(0)
# write your own tests here
function generate_signal(f,w,N, modphase=false)
    x = sort(10rand(N)) # Sample points
    v = range(0, stop=1, length=N) # Scheduling variable

    # generate output signal
    # phase_matrix
    dependence_matrix = Float64[f[(i-1)%length(f)+1].(v) for v in v, i in eachindex(w)] # N x nw
    frequency_matrix = [cos(w*x -0.5modphase*(dependence_matrix[i,j])) for (i,x) in enumerate(x), (j,w) in enumerate(w)] # N x nw
    y = sum(dependence_matrix.*frequency_matrix,dims=2)[:] # Sum over all frequencies
    y += 0.1randn(size(y))
    y,v,x,frequency_matrix, dependence_matrix
end

@testset "LPVSpectral" begin
    @info "Testing LPVSpectral"

    @testset "Windows" begin
        @info "Testing Windows"

        y = 1:100
        t = 1:100
        W = Windows2(y,t,10,0)
        @test length(W) == 10
        @test first(W) == (1:10,1:10)
        res = mapwindows(W) do (y,t)
            -y
        end
        @test res == -1:-1:-100

        W = Windows2(y,t,10,1)
        @test length(W) == 11
        cW = collect(W)
        @test cW[1] == (1:10,1:10)
        @test cW[2] == (10:19,10:19)

        res = mapwindows(W) do (y,t)
            -y
        end
        @test res == -1:-1:-100

        y = 1:100
        t = 1:100
        W = Windows3(y,t,t,10,0)
        @test length(W) == 10
        @test first(W) == (1:10,1:10,1:10)

        W = Windows3(y,t,t,10,1)
        @test length(W) == 11
        cW = collect(W)
        @test cW[1] == (1:10,1:10,1:10)
        @test cW[2] == (10:19,10:19,10:19)

    end

    @testset "Mel" begin
        @info "Testing Mel"

        M = mel(1,256)
        @test size(M,1) == 128
        @test size(M,2) == 256÷2+1
        M = mel(1000,256, fmin=100)
        sum(M[:,1:26]) == 0

        y = randn(1000)
        M = melspectrogram(y)
        @test length(freq(M)) == 128
        @test size(M.power) == (128,14)
        @test length(time(M)) == 14
        plot(M)

        M = mfcc(y)
        @test length(freq(M)) == 20
        @test size(M.mfcc) == (20,14)
        @test length(time(M)) == 14

    end

    @testset "LPV methods" begin
        @info "testing LPV methods"

        N      = 500 # Number of training data points
        f      = [v->2v^2, v->2/(5v+1), v->3exp(-10*(v-0.5)^2),] # Functional dependences on the scheduling variable
        w      = 2pi*[2,10,20] # Frequency vector
        w_test = 2π*collect(2:2:25)
        Y,V,X,frequency_matrix, dependence_matrix = generate_signal(f,w,N,true)

        λ      = 0.02 # Regularization parmater
        normal = true # Use normalized basis functions
        Nv     = 50 # Number of basis functions

        se = ls_spectral_lpv(Y,X,V,w_test,Nv; λ = λ, normalize = normal) # Perform LPV spectral estimation
        @show plot(se, phase=true)
        windowpsd = ls_windowpsd_lpv(Y,X,V,w_test,Nv; λ = λ, normalize = normal)

        spectrum_lpv   = psd(se)

        si = sortperm(spectrum_lpv[:],rev=true)
        @test Set(si[1:3]) == Set([1,5,10])

        si = sortperm(windowpsd[:],rev=true)
        @test Set(si[1:3]) == Set([1,5,10])
    end

    @testset "detrend" begin
        @info "testing detrend"
        tre = [1,2,3]
        @test detrend(tre) == [-1,0,1]
        @test tre == [1,2,3] # Shouldn't have changed
        detrend!(tre)
        @test tre == [-1,0,1] # Should have changed
    end

    ## Test ComplexNormal
    @testset "Complex Normal" begin
        @info "testing Complex Normal"
        n  = 10
        n2 = 5
        a  = randn(n)
        A  = randn(n,n)
        A  = A'A + I
        b  = randn(n2)
        B  = randn(n2,n2)
        B  = B'B + I
        X  = randn(n,n2)
        Y  = randn(n,n2)
        cn = ComplexNormal(X,Y)
        @test size(cn.m) == (n2,)
        @test size(cn.Γ) == (n2,n2)
        @test size(cn.C) == (n2,n2)
        # affine_transform(cn, B, b)
        pdf(cn,b)
        pdf(cn,im*b)

        @test isa(cn_Vxx(A,A), Cholesky)
        @test isa(cn_fVxx(A,A), Matrix)
        @test issymmetric(cn_fVxx(A,A))
        cn_V(A,0.1A)

        cn  = ComplexNormal(a,A)
        cn  = ComplexNormal(im*b,A)

        A   = randn(4,4);
        A   = A'A
        x   = randn(1000,3)
        y   = randn(1000,3)
        cn  = ComplexNormal(x,y)
        z   = rand(cn,1000000);
        cn2 = ComplexNormal(z)
        @test norm(Matrix(cn.Γ)-Matrix(cn2.Γ)) < 0.01
        @test norm(Matrix(cn.C)-Matrix(cn2.C)) < 0.01

    end

    @testset "ls methods" begin
        @info "testing ls methods"
        T = 100
        t = 0:0.1:T-0.1
        f = LPVSpectral.default_freqs(t)
        @test f[1] == 0
        @test f[end] == 5
        @test length(f) == 10T÷2+1

        # f2 = LPVSpectral.default_freqs(t,10)
        # @test f2[1] == 0
        # @test f2[end] == 0.5-1/length(t)/2
        # @test length(f2) == T÷2÷10

        @test LPVSpectral.check_freq(f) == 1
        @test_throws ArgumentError LPVSpectral.check_freq([1,0,2])
        A, z = LPVSpectral.get_fourier_regressor(t,f)
        @test size(A) == (10T,2length(f)-1)

        Base.isapprox(t1::Tuple{Float64,Int64}, t2::Tuple{Float64,Int64}; atol) = all(t -> isapprox(t[1],t[2],atol=atol), zip(t1,t2))
        y = sin.(2pi .* t)
        x,freqs = ls_spectral(y,t)
        @test findmax(abs2.(x)) ≈ (2.0length(freqs), 101) atol=1e-4

        W = ones(length(y))
        x,_ = ls_spectral(y,t,f,W)
        @test findmax(abs2.(x)) ≈ (2.0length(freqs), 101) atol=1e-4

        x,freqs = tls_spectral(y,t)
        @test findmax(abs2.(x)) ≈ (2.0length(freqs), 101) atol=1e-4

        x,freqs = ls_windowpsd(y,t,noverlap=0)
        @test findmax(x)[2] == 13

        x,freqs = ls_windowpsd(y,t,nw=16,noverlap=0)
        @test findmax(abs.(x))[2] ==  7

        x,freqs = ls_windowcsd(y,y,t,noverlap=0)
        @test findmax(abs.(x)) ≈ (2.0length(freqs), 11) atol=1e-4

        x,_ = ls_cohere(y,y,t)
        @test all(x .== 1)

        x,freqs = ls_cohere(y,y .+ 0.5 .*randn.(),t,nw=8,noverlap=-1)
        @show mean(x)
        @test findmax(x) ≈ (1.0, 14) atol=0.15
        @test mean(x) < 0.25

    end

    @testset "plots" begin
        @info "testing plots"

        y = randn(1000)
        @show plot(periodogram(y))
        @show plot(spectrogram(y))
        @show plot(melspectrogram(y))
        @show plot(periodogram(filtfilt(ones(4), [4], y)))
        @show plot(welch_pgram(y))

    end

    @testset "Lasso" begin
        @info "Testing Lasso"
        include("test_lasso.jl")
    end

    @testset "autocov" begin
        @info "Testing autocov"

        y = repeat([1.,0.,-1.],100)
        τ,acf = autocov(1:length(y), y, Inf)
        acf0 = autocov(y, demean=false)
        acfh = [mean(acf[τ.==i]) for i = 0:length(acf0)-1]
        @test acfh ≈ acf0 rtol=0.01
        @test norm(acfh - acf0) < 0.01

        τ,acf = autocor(1:length(y), y, Inf)
        acf0 = autocor(y, demean=false)
        acfh = [mean(acf[τ.==i]) for i = 0:length(acf0)-1]
        @test acfh ≈ acf0 rtol=0.01
        @test norm(acfh - acf0) < 0.1

        y = randn(100)
        τ,acf = autocor(1:length(y), y, Inf)
        acf0 = autocor(y, demean=false)
        acfh = [mean(acf[τ.==i]) for i = 0:length(acf0)-1]
        @test acfh ≈ acf0 rtol=0.01
        @test norm(acfh - acf0) < 0.1

        y = randn(10)
        τ,acf = autocov(1:length(y), y, Inf)
        acf0 = autocov(y, demean=false)
        acfh = [mean(acf[τ.==i]) for i = 0:length(acf0)-1]
        @test acfh ≈ acf0 rtol=0.01
        @test norm(acfh - acf0) < 0.2

        y = randn(10)
        τ,acf = autocor(1:length(y), y, Inf)
        acf0 = autocor(y, demean=false)
        acfh = [mean(acf[τ.==i]) for i = 0:length(acf0)-1]
        @test acfh ≈ acf0 rtol=0.03
        @test norm(acfh - acf0) < 0.2

        y = [randn(10) for _ in 1:10]
        t = reshape(1.0:100,10,10)
        t = [t[:,i] for i in 1:size(t,2)]
        τ,acf = autocov(t, y, Inf)
        acf0 = mean(autocov.(y, demean=false))
        acfh = [mean(acf[τ.==i]) for i = 0:length(acf0)-1]
        @test acfh ≈ acf0 rtol=0.01
        @test norm(acfh - acf0) < 0.2

        τ,acf = autocor(t, y, Inf)
        acf0 = mean(autocor.(y, demean=false))
        acfh = [mean(acf[τ.==i]) for i = 0:length(acf0)-1]
        @test acfh ≈ acf0 rtol=0.01
        @test norm(acfh - acf0) < 0.2

        # See what happens if one series is contant
        y = zeros(10)
        τ,acf = autocor(1:10, y, Inf)
        @test all(acf .== 1)

        using LPVSpectral: _autocov, _autocor, isequidistant

        @test isequidistant(1:5)
        @test isequidistant(1:2:10)
        @test !isequidistant(reverse(1:2:10))
        @test isequidistant(collect(1:5))
        @test isequidistant(collect(1:2:10))
        @test isequidistant(collect(1:0.33:10))

        res = map(1:10) do _
            t = 100rand(100)
            @test !isequidistant(t)
            t0 = 0:99
            y = sin.(0.05 .* t)
            y0 = sin.(0.05 .* t0)
            τ0, acf0 = autocor(t0, y0, Inf, normalize=true)
            τ,acf = autocor(t, y, Inf)
            acff = filtfilt(ones(200),[200], acf)
            @test count(τ .== 0) == length(y)
            # plot(τ[1:10:end],acf[1:10:end])
            # plot!(τ0[1:10:end],acf0[1:10:end])
            # plot!(τ0[1:10:end],acff[1:10:end])
            mean(abs2,acf0-acff) < 0.05
        end
        @test mean(res) > 0.7

        res = map(1:10) do _
            t = 100rand(100)
            @test !isequidistant(t)
            t0 = 0:99
            y = sin.(0.05 .* t)
            y0 = sin.(0.05 .* t0)
            τ0, acf0 = autocov(t0, y0, Inf, normalize=true)
            τ,acf = autocov(t, y, Inf)
            acff = filtfilt(ones(200),[200], acf)
            @test count(τ .== 0) == length(y)
            # plot(τ[1:10:end],acf[1:10:end])
            # plot!(τ0[1:10:end],acf0[1:10:end])
            # plot!(τ0[1:10:end],acff[1:10:end])
            mean(abs2,acf0-acff) < 0.025
        end
        @test mean(res) > 0.7

    end

end


text:
using ThreeBodyDecay
using Test

@testset "Xic2pKpi" begin
    mp = 0.938; mK = 0.49367; mpi  = 0.13957; mXic = 2.46867
    tbs = ThreeBodySystem(mp,mK,mpi; m0=mXic)
    sample = flatDalitzPlotSample(tbs.ms)
    @test true
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
# Data hand processed from https://github.com/dengwirda/inpoly by @miguelraz on 28/03/2020 under terms of fair use.

include("coast.jl")
include("coast-edges.jl")
include("lakes.jl")
include("lakes-edges.jl")
include("river.jl")
include("river-edges.jl")


text:
# This file is a part of SIS3316Digitizers.jl, licensed under the MIT License (MIT).

@enum VMEAddressSpace begin
    VME_A16 = 16
    VME_A24 = 24
    VME_A32 = 32
    VME_A64 = 64
end
export VMEAddressSpace

@enum VMEDataWidth begin
    VME_D8  =  8
    VME_D16 = 16
    VME_D32 = 32
    VME_D64 = 64
end
export VMEDataWidth

@enum VMECycle begin
    VME_SCT         = 0
    VME_BLT         = 1
    VME_MBLT        = 2
    VME_eeVME       = 3
    VME_eeSST160    = 4
    VME_eeSST267    = 5
    VME_eeSST320    = 6
end
export VMECycle

struct VMEMode
    space::VMEAddressSpace
    width::VMEDataWidth
    cycle::VMECycle
end
export VMEMode

const VME_A32_D32_SCT = VMEMode(VME_A32, VME_D32, VME_SCT)
const VME_A32_D32_BLT = VMEMode(VME_A32, VME_D32, VME_SCT)
const VME_A32_D64_MBLT = VMEMode(VME_A32, VME_D64, VME_MBLT)
const VME_A32_D64_2eVME = VMEMode(VME_A32, VME_D64, VME_eeVME)


text:
using GR

mwidth, mheight, width, height = GR.inqdspsize()

w = 800
h = 600

setcolorrep(1, 1.0, 1.0, 1.0)
setcolorrep(2, 0.0, 0.0, 0.0)
setlinecolorind(2)
setlinewidth(1)
setlinetype(GR.LINETYPE_SOLID)
setfillcolorind(2)
setfillintstyle(GR.INTSTYLE_SOLID)

while true
    clearws()

    #w = 400 + round(Int, rand() * 200)
    #h = 400 + round(Int, rand() * 200)
    if w >= h
        ratio = float(h) / w
        msize = mwidth * w / width
        GR.setwsviewport(0, msize, 0, msize * ratio)
        GR.setwswindow(0, 1, 0, ratio)
        setviewport(0, 1, 0, ratio)
    else
        ratio = float(w) / h
        msize = mheight * h / height
        GR.setwsviewport(0, msize * ratio, 0, msize)
        GR.setwswindow(0, ratio, 0, 1)
        setviewport(0, ratio, 0, 1)
    end

    setwindow(1, w, 1, h)

    #drawrect(1, w, 1, h)
    fillrect(1, w, 1, h)
    #polyline([1,w],[1,h])
    #polyline([1,w],[h,1])
    text(0.05, 0.05, "$w x $h")

    updatews()
    sleep(0.5)
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
# Beyn contour integral approach

#using QuadGK # Disabled to lessen requirements
using Distributed
using LinearAlgebra
using Random

export contour_beyn

"""
    λv,V=contour_beyn([eltype,] nep [,mintegrator];[tol,][logger,][σ,][radius,][linsolvercreator,][N,][neigs,][k])

The function computes eigenvalues using Beyn's contour integral approach,
using an ellipse centered at `σ` with radii given in `radius`, or if only one `radius` is given,
the contour is a circle. The numerical quadrature method is specified in `mintegrator`,
which is a type inheriting from `MatrixIntegrator`, by default
`MatrixTrapezoidal`. For a parallell implementation of the
integrator use `MatrixTrapezoidalParallel`.
 The integer `k`
specifies size of the probe subspace. `N` corresponds to the
number of quadrature points. Ellipses are the only supported contours. The
`linsolvercreator` must create a linsolver that can handle (rectangular) matrices
as right-hand sides, not only vectors. We integrate in complex arithmetic so
`eltype` must be complex type.

The kwargs `neigs` specifies the number of wanted eigvals, and `k` is the number
of columns in the matrix to be integrated (default `k=neigs+1`). If you give the `k` parameter and set `neigs=typemax(Int)` all found eigenvalues will be returned. The kwarg `sanity_check` decides if sorting and checking (and removal) of eigpairs should be done.
If disabled, the method
returns `k` (potentially inaccurate) eigpairs. The parameters `errmeasure` and
`tol` and `rank_drop_tol` are used for the sanity check,
to extract accurate eigenvalues.

# Example
```julia-repl
julia> using LinearAlgebra
julia> nep=nep_gallery("dep0");
julia> # Look for two eigvals in unit disk
julia> λv,V=contour_beyn(nep,radius=1.1,neigs=3);
julia> norm(compute_Mlincomb(nep,λv[1],V[:,1])) # Eigenpair 1
1.7462847531404259e-15
julia> norm(compute_Mlincomb(nep,λv[2],V[:,2])) # Eigenpair 2
7.69695692032292e-15
```
# References
* Wolf-Jürgen Beyn, An integral method for solving nonlinear eigenvalue problems, Linear Algebra and its Applications 436 (2012) 3839–3863
"""
contour_beyn(nep::NEP;params...)=contour_beyn(ComplexF64,nep;params...)
contour_beyn(nep::NEP,MIntegrator;params...)=contour_beyn(ComplexF64,nep,MIntegrator;params...)
function contour_beyn(::Type{T},
                      nep::NEP,
                      ::Type{MIntegrator}=MatrixTrapezoidal;
                      tol::Real=sqrt(eps(real(T))), # Note tol is quite high for this method
                      σ::Number=zero(complex(T)),
                      logger=0,
                      linsolvercreator=BackslashLinSolverCreator(),
                      neigs::Integer=2, # Number of wanted eigvals
                      k::Integer=neigs+1, # Columns in matrix to integrate
                      radius::Union{Real,Tuple,Array}=1, # integration radius
                      N::Integer=1000,  # Nof quadrature nodes
                      errmeasure::ErrmeasureType = DefaultErrmeasure(nep),
                      sanity_check=true,
                      rank_drop_tol=tol # Used in sanity checking
                      )where{T<:Number, MIntegrator<:MatrixIntegrator}

    @parse_logger_param!(logger)

    # Geometry
    length(radius)==1 ? radius=(radius,radius) : nothing
    g(t) = complex(radius[1]*cos(t),radius[2]*sin(t)) # ellipse
    gp(t) = complex(-radius[1]*sin(t),radius[2]*cos(t)) # derivative

    n=size(nep,1);

    if (k>n)
        error("Cannot compute more eigenvalues than the size of the NEP with contour_beyn() k=",k," n=",n);

    end
    if (k<=0)
        error("k must be positive, k=",k,
              neigs==typemax(Int) ? ". The kwarg k must be set if you use neigs=typemax" : ".")
    end

    Random.seed!(10); # Reproducability
    Vh=Array{T,2}(randn(real(T),n,k)) # randn only works for real

    function local_linsolve(λ::TT,V::Matrix{TT}) where {TT<:Number}
        local M0inv::LinSolver = create_linsolver(linsolvercreator, nep, λ+σ);
        # This requires that lin_solve can handle rectangular
        # matrices as the RHS
        return lin_solve(M0inv,V);
    end

    # Constructing integrands
    Tv(λ) = local_linsolve(T(λ),Vh)
    f(t) = Tv(g(t))*gp(t)
    push_info!(logger,"Computing integrals")

    local A0,A1
    AA=integrate_interval(MIntegrator, ComplexF64, f,
                        [ s-> one(s), g],0,2*pi,N,logger )
    A0=AA[:,:,1];
    A1=AA[:,:,2];

    # Don't forget scaling
    A0[:,:] = A0 ./(2im*pi);
    A1[:,:] = A1 ./(2im*pi);

    push_info!(logger,"Computing SVD prepare for eigenvalue extraction ",continues=true)
    V,S,W = svd(A0)
    p = count( S/S[1] .> rank_drop_tol);
    push_info!(logger," p=$p");

    V0 = V[:,1:p]
    W0 = W[:,1:p]
    B = (copy(V0')*A1*W0) * Diagonal(1 ./ S[1:p])

    # Extract eigenval and eigvec approximations according to
    # step 6 on page 3849 in the reference
    push_info!(logger,"Computing eigenvalues ")
    λ,VB=eigen(B)
    λ[:] = λ .+ σ
    push_info!(logger,"Computing eigenvectors ")
    V = V0 * VB;
    for i = 1:p
        normalize!(V[:,i]);
    end

    if (!sanity_check)
        sorted_index = sortperm(map(x->abs(σ-x), λ));
        inside_bool = (real(λ[sorted_index].-σ)/radius[1]).^2 + (imag(λ[sorted_index].-σ)/radius[2]).^2 .≤ 1
        inside_perm = sortperm(.!inside_bool)
        return (λ[sorted_index[inside_perm]],V[:,sorted_index[inside_perm]])
    end

    # Compute all the errors
    errmeasures=zeros(real(T),p);
    for i = 1:p
        errmeasures[i]=estimate_error(errmeasure,λ[i],V[:,i]);
    end

    good_index=findall(errmeasures .< tol);

    # Index vector for sorted to distance from σ
    sorted_good_index=
       good_index[sortperm(map(x->abs(σ-x), λ[good_index]))];

    # check that eigenvalues are inside contour, move them to the end if they are outside
   inside_bool = (real(λ[sorted_good_index].-σ)/radius[1]).^2 + (imag(λ[sorted_good_index].-σ)/radius[2]).^2 .≤ 1
   if any(.!inside_bool)
       if neigs ≤ sum(inside_bool)
           @warn "found $(sum(.!inside_bool)) evals outside contour, $p inside. all $neigs returned evals inside, but possibly inaccurate. try increasing N, decreasing tol, or changing radius"
       else
           @warn "found $(sum(.!inside_bool)) evals outside contour, $p inside. last $(neigs-sum(inside_bool)) returned evals outside contour. possible inaccuracy. try increasing N, decreasing tol, or changing radius"
       end
   end
   sorted_good_inside_perm = sortperm(.!inside_bool)

    # Remove all eigpairs not sufficiently accurate
    # and potentially remove eigenvalues if more than neigs.
    local Vgood,λgood
    if( size(sorted_good_index,1) > neigs)
        found_evals=size(sorted_good_index,1);
        push_info!(logger,"Removing unwanted eigvals: neigs=$neigs<$found_evals=found_eigvals")
        Vgood=V[:,sorted_good_index[sorted_good_inside_perm][1:neigs]];
        λgood=λ[sorted_good_index[sorted_good_inside_perm][1:neigs]];
    else
        Vgood=V[:,sorted_good_index[sorted_good_inside_perm]];
        λgood=λ[sorted_good_index[sorted_good_inside_perm]];
    end

    if (p==k)
       @warn "Rank-drop not detected, your eigvals may be correct, but the algorithm cannot verify. Try to increase k. This warning can be disabled with `sanity_check=false`." S
    end

    if (size(λgood,1)<neigs  && neigs < typemax(Int))
       @warn "We found fewer eigvals than requested. Try increasing domain, or decreasing `tol`. This warning can be disabled with `sanity_check=false`." S
    end

    return (λgood,Vgood)
end


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

@doc raw"""SubjectRulesReviewStatus contains the result of a rules check. This check can be incomplete depending on the set of authorizers the server is configured with and any errors experienced during evaluation. Because authorization rules are additive, if a rule appears in a list it's safe to assume the subject has that permission, even if that list is incomplete.

    IoK8sApiAuthorizationV1SubjectRulesReviewStatus(;
        evaluationError=nothing,
        incomplete=nothing,
        nonResourceRules=nothing,
        resourceRules=nothing,
    )

    - evaluationError::String : EvaluationError can appear in combination with Rules. It indicates an error occurred during rule evaluation, such as an authorizer that doesn't support rule evaluation, and that ResourceRules and/or NonResourceRules may be incomplete.
    - incomplete::Bool : Incomplete is true when the rules returned by this call are incomplete. This is most commonly encountered when an authorizer, such as an external authorizer, doesn't support rules evaluation.
    - nonResourceRules::Vector{IoK8sApiAuthorizationV1NonResourceRule} : NonResourceRules is the list of actions the subject is allowed to perform on non-resources. The list ordering isn't significant, may contain duplicates, and possibly be incomplete.
    - resourceRules::Vector{IoK8sApiAuthorizationV1ResourceRule} : ResourceRules is the list of actions the subject is allowed to perform on resources. The list ordering isn't significant, may contain duplicates, and possibly be incomplete.
"""
mutable struct IoK8sApiAuthorizationV1SubjectRulesReviewStatus <: SwaggerModel
    evaluationError::Any # spec type: Union{ Nothing, String } # spec name: evaluationError
    incomplete::Any # spec type: Union{ Nothing, Bool } # spec name: incomplete
    nonResourceRules::Any # spec type: Union{ Nothing, Vector{IoK8sApiAuthorizationV1NonResourceRule} } # spec name: nonResourceRules
    resourceRules::Any # spec type: Union{ Nothing, Vector{IoK8sApiAuthorizationV1ResourceRule} } # spec name: resourceRules

    function IoK8sApiAuthorizationV1SubjectRulesReviewStatus(;evaluationError=nothing, incomplete=nothing, nonResourceRules=nothing, resourceRules=nothing)
        o = new()
        validate_property(IoK8sApiAuthorizationV1SubjectRulesReviewStatus, Symbol("evaluationError"), evaluationError)
        setfield!(o, Symbol("evaluationError"), evaluationError)
        validate_property(IoK8sApiAuthorizationV1SubjectRulesReviewStatus, Symbol("incomplete"), incomplete)
        setfield!(o, Symbol("incomplete"), incomplete)
        validate_property(IoK8sApiAuthorizationV1SubjectRulesReviewStatus, Symbol("nonResourceRules"), nonResourceRules)
        setfield!(o, Symbol("nonResourceRules"), nonResourceRules)
        validate_property(IoK8sApiAuthorizationV1SubjectRulesReviewStatus, Symbol("resourceRules"), resourceRules)
        setfield!(o, Symbol("resourceRules"), resourceRules)
        o
    end
end # type IoK8sApiAuthorizationV1SubjectRulesReviewStatus

const _property_map_IoK8sApiAuthorizationV1SubjectRulesReviewStatus = Dict{Symbol,Symbol}(Symbol("evaluationError")=>Symbol("evaluationError"), Symbol("incomplete")=>Symbol("incomplete"), Symbol("nonResourceRules")=>Symbol("nonResourceRules"), Symbol("resourceRules")=>Symbol("resourceRules"))
const _property_types_IoK8sApiAuthorizationV1SubjectRulesReviewStatus = Dict{Symbol,String}(Symbol("evaluationError")=>"String", Symbol("incomplete")=>"Bool", Symbol("nonResourceRules")=>"Vector{IoK8sApiAuthorizationV1NonResourceRule}", Symbol("resourceRules")=>"Vector{IoK8sApiAuthorizationV1ResourceRule}")
Base.propertynames(::Type{ IoK8sApiAuthorizationV1SubjectRulesReviewStatus }) = collect(keys(_property_map_IoK8sApiAuthorizationV1SubjectRulesReviewStatus))
Swagger.property_type(::Type{ IoK8sApiAuthorizationV1SubjectRulesReviewStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAuthorizationV1SubjectRulesReviewStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthorizationV1SubjectRulesReviewStatus }, property_name::Symbol) =  _property_map_IoK8sApiAuthorizationV1SubjectRulesReviewStatus[property_name]

function check_required(o::IoK8sApiAuthorizationV1SubjectRulesReviewStatus)
    (getproperty(o, Symbol("incomplete")) === nothing) && (return false)
    (getproperty(o, Symbol("nonResourceRules")) === nothing) && (return false)
    (getproperty(o, Symbol("resourceRules")) === nothing) && (return false)
    true
end

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


text:
module QMPS
using Yao, Yao.Blocks
using Yao.Intrinsics: packbits
import Yao.Registers: probs, nqubits

using StatsBase
using StatsBase: mean

using LinearAlgebra
using KrylovKit

include("Adam.jl")
include("Core.jl")
include("AbstractModel.jl")
include("gradient.jl")
include("correlation.jl")
include("ansatz/ansatz.jl")
end

max_stars_repo_name: Roger-luo/QuantumMPS

max_issues_repo_name: Roger-luo/QuantumMPS

max_forks_repo_name: Roger-luo/QuantumMPS