text:
module PFC

function generate_solver(f₀, f; fₘᵢₙ = 0, fₘₐₓ = 1)
    function ϵ⁺(f, g)
        if f < g
            return min(g - f, 2.0*(f - fₘᵢₙ))
        else
            return max(g - f, -2.0*(fₘₐₓ - f))
        end
    end

    function ϵ⁻(f, g)
        if f < g
            return max(f - g, -2.0*(f - fₘᵢₙ))
        else
            return min(f - g, 2.0*(fₘₐₓ - f))
        end
    end

    Φ⁺(f, i, i⁻, i⁺, c) = c*(f[i] + (1.0 - c)/3.0*(
                    ϵ⁺(f[i], f[i⁺])/2.0*(2.0 - c) +
                    ϵ⁻(f[i], f[i⁻])/2.0*(1.0 + c)))

    Φ⁻(f, i, i⁻, i⁺, c) = c*(f[i] - (1.0 + c)/3.0*(
                    ϵ⁺(f[i], f[i⁺])/2.0*(1.0 - c) +
                    ϵ⁻(f[i], f[i⁻])/2.0*(2.0 + c)))

    function solve!(c)
        if c > 0
            Φ = Φ⁺(f₀, 1, length(f), 2, c)
            f[1] = f₀[1] - Φ
            f[2] = f₀[2] + Φ
            for i in 2:length(f)-1
                Φ = Φ⁺(f₀, i, i-1, i+1, c)
                f[i] = f[i] - Φ
                f[i+1] = f₀[i+1] + Φ
            end
            Φ = Φ⁺(f₀, length(f), length(f)-1, 1, c)
            f[end] = f[end] - Φ
            f[1] = f[1] + Φ
        else
            Φ = Φ⁻(f₀, 1, length(f), 2, c)
            f[1] = f₀[1] + Φ
            f[end] = f₀[end] - Φ
            for i in 2:length(f)-1
                Φ = Φ⁻(f₀, i, i-1, i+1, c)
                f[i] = f₀[i] + Φ
                f[i-1] = f[i-1] - Φ
            end
            Φ = Φ⁻(f₀, length(f), length(f)-1, 1, c)
            f[end] = f[end] + Φ
            f[end-1] = f[end-1] - Φ
        end
    end

    function solve!(h, h₀, c)
        if c > 0
            Φ = Φ⁺(h₀, 1, length(f), 2, c)
            h[1] = h₀[1] - Φ
            h[2] = h₀[2] + Φ
            for i in 2:length(f)-1
                Φ = Φ⁺(h₀, i, i-1, i+1, c)
                h[i] = h[i] - Φ
                h[i+1] = h₀[i+1] + Φ
            end
            Φ = Φ⁺(h₀, length(f), length(f)-1, 1, c)
            h[end] = h[end] - Φ
            h[1] = h[1] + Φ
        else
            Φ = Φ⁻(h₀, 1, length(f), 2, c)
            h[1] = h₀[1] + Φ
            h[end] = h₀[end] - Φ
            for i in 2:length(f)-1
                Φ = Φ⁻(h₀, i, i-1, i+1, c)
                h[i] = h₀[i] + Φ
                h[i-1] = h[i-1] - Φ
            end
            Φ = Φ⁻(h₀, length(f), length(f)-1, 1, c)
            h[end] = h[end] + Φ
            h[end-1] = h[end-1] - Φ
        end
    end

    return solve!
end

function generate_solver(f₀, f, c)
    fₘₐₓ = maximum(f₀)

    function ϵ⁺(f, g)
        if f < g
            return min(g - f, 2.0*f)
        else
            return max(g - f, -2.0*(fₘₐₓ - f))
        end
    end

    function ϵ⁻(f, g)
        if f < g
            return max(f - g, -2.0*f)
        else
            return min(f - g, 2.0*(fₘₐₓ - f))
        end
    end

    Φ⁺(f, i, i⁻, i⁺, c) = c*(f[i] + (1.0 - c)/3.0*(
                    ϵ⁺(f[i], f[i⁺])/2.0*(2.0 - c) +
                    ϵ⁻(f[i], f[i⁻])/2.0*(1.0 + c)))

    Φ⁻(f, i, i⁻, i⁺, c) = c*(f[i] - (1.0 + c)/3.0*(
                    ϵ⁺(f[i], f[i⁺])/2.0*(1.0 - c) +
                    ϵ⁻(f[i], f[i⁻])/2.0*(2.0 + c)))

    function solve⁺!()
        Φ = Φ⁺(f₀, 1, length(f), 2, c)
        f[1] = f₀[1] - Φ
        f[2] = f₀[2] + Φ
        for i in 2:length(f)-1
            Φ = Φ⁺(f₀, i, i-1, i+1, c)
            f[i] = f[i] - Φ
            f[i+1] = f₀[i+1] + Φ
        end
        Φ = Φ⁺(f₀, length(f), length(f)-1, 1, c)
        f[end] = f[end] - Φ
        f[1] = f[1] + Φ
    end

    function solve⁻!()
        Φ = Φ⁻(f₀, 1, length(f), 2, c)
        f[1] = f₀[1] + Φ
        f[end] = f₀[end] - Φ
        for i in 2:length(f)-1
            Φ = Φ⁻(f₀, i, i-1, i+1, c)
            f[i] = f₀[i] + Φ
            f[i-1] = f[i-1] - Φ
        end
        Φ = Φ⁻(f₀, length(f), length(f)-1, 1, c)
        f[end] = f[end] + Φ
        f[end-1] = f[end-1] - Φ
    end

    if c > 0
        return solve⁺!
    else
        return solve⁻!
    end
end

end  # module PFC


text:
module MultivariateCalculusTest

using VarianceComponentModels, LinearAlgebra, Random, Test

Random.seed!(123)

# test vech
@testset "vech" begin
  a = reshape(1:9, 3, 3)
  vecha = vech(a)
  @test vecha[1] == a[1, 1]
  @test vecha[2] == a[2, 1]
  @test vecha[3] == a[3, 1]
  @test vecha[4] == a[2, 2]
  @test vecha[5] == a[3, 2]
  @test vecha[6] == a[3, 3]
  a = 1.0
  @test vech(a) == a
end

# test trilind
@testset "trilind" begin
  n = 3
  A = randn(n, n)
  @test norm(vech(A) - A[trilind(A)]) ≈ 0.0
  @test norm([2, 3, 6] - trilind(A, -1)) ≈ 0.0
end

# test triuind
@testset "triuind" begin
  n = 2
  idx = triuind(n, n)
  @test length(idx) == binomial(n + 1, 2)
  @test idx[1] == 1
  @test idx[2] == 3
  @test idx[3] == 4
  idx = triuind(n)
  @test length(idx) == binomial(n + 1, 2)
  @test idx[1] == 1
  @test idx[2] == 3
  @test idx[3] == 4
  idx = triuind(randn(n, n))
  @test length(idx) == binomial(n + 1, 2)
  @test idx[1] == 1
  @test idx[2] == 3
  @test idx[3] == 4
  idx = triuind(randn(n, n), 1)
  @test length(idx) == binomial(n, 2)
  @test idx[1] == 3
end

# test commutation
@testset "commutation" begin
  m, n = 3, 2
  A = randn(m, n)
  @test norm(commutation(m, n) * vec(A) - vec(A')) ≈ 0.0
  @test norm(commutation(A) * vec(A) - vec(A')) ≈ 0.0
  @test norm(spcommutation(m, n) * vec(A) - vec(A')) ≈ 0.0
  @test norm(spcommutation(A) * vec(A) - vec(A')) ≈ 0.0
  B = randn(m, m)
  @test norm(commutation(m) * vec(B) - vec(B')) ≈ 0.0
  @test norm(commutation(eltype(B), m) * vec(B) - vec(B')) ≈ 0.0
  @test norm(spcommutation(eltype(B), m) * vec(B) - vec(B')) ≈ 0.0
end

# test duplication
@testset "duplication" begin
  n = 3
  A = randn(n, n)
  A = 0.5(A + A') # symmetrize
  @test norm(duplication(n) * vech(A) - vec(A)) ≈ 0.0
  @test norm(duplication(A) * vech(A) - vec(A)) ≈ 0.0
  @test norm(spduplication(n) * vech(A) - vec(A)) ≈ 0.0
  @test norm(spduplication(A) * vech(A) - vec(A)) ≈ 0.0
end

# test chol_gradient
@testset "chol_gradient" begin
  n = 3
  A = randn(n, n)
  A = 0.5(A + A') # symmetrize
  L = tril(randn(n, n))
  # calculate gradient wrt L using chol_gradient
  dL1 = chol_gradient(vec(A), L)
  # alternative way to calculate gradient wrt L
  dL2 = 2.0vech((A * L) + L' * A' - Matrix(Diagonal(A * L)))
  @test norm(dL1 - dL2) ≈ 0.0
end

# test kron_gradient
@testset "kron_gradient" begin
  n, q, p, r = 2, 3, 4, 5
  X = randn(n, q)
  Y = randn(p, r)
  M = kron(X, Y)
  # calculate gradient wrt X using kron_gradient
  dX1 = kron_gradient(vec(M), Y, n, q)
  # alternative way to calculate gradient wrt X
  dX2 = zero(X)
  for j = 1:q
    for i = 1:n
      dX2[i, j] = dot(Y, M[(i-1)*p+1:i*p, (j-1)*r+1:j*r])
    end
  end
  @test norm(dX1 - vec(dX2)) < 1e-8
end

# test duplication
@testset "kronaxpy" begin
  m, n, p, q = 3, 4, 5, 6
  A = randn(m, n)
  X = randn(p, q)
  Y = zeros(m * p, n * q)
  kronaxpy!(A, X, Y)
  @test norm(Y - kron(A, X)) ≈ 0.0
end

# test bump_diagonal
@testset "bump_diag" begin
  A = randn(3, 3)
  B = copy(A)
  bump_diagonal!(B, 1e-3)
  @test B[1, 1] - A[1, 1] ≈ 1e-3
  @test B[2, 2] - A[2, 2] ≈ 1e-3
  @test B[3, 3] - A[3, 3] ≈ 1e-3
end

# test clamp_diagonal
@testset "clamp_diag" begin
  n = 3
  A = randn(n, n)
  B = copy(A)
  clamp_diagonal!(B, 0.0, 1.0)
  @test all(0.0 .≤ diag(B) .≤ 1.0)
end

end # module MultivariateCalculusTest

ext: jl

max_stars_repo_name: Hua-Zhou/VarianceComponentModels

max_issues_repo_name: Hua-Zhou/VarianceComponentModels

max_forks_repo_name: Hua-Zhou/VarianceComponentModels


text:
# ---
# title: 552. Student Attendance Record II
# id: problem552
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/student-attendance-record-ii/description/>
# hidden: true
# ---
#
# Given a positive integer **n** , return the number of all possible attendance
# records with length n, which will be regarded as rewardable. The answer may be
# very large, return it after mod 109 \+ 7.
#
# A student attendance record is a string that only contains the following three
# characters:
#
#   1. **'A'** : Absent.
#   2. **'L'** : Late.
#   3. **'P'** : Present.
#
# A record is regarded as rewardable if it doesn't contain **more than one 'A'
# (absent)** or **more than two continuous 'L' (late)**.
#
# **Example 1:**
#
#
#
#     Input: n = 2
#     Output: 8
#     Explanation:
#     There are 8 records with length 2 will be regarded as rewardable:
#     "PP" , "AP", "PA", "LP", "PL", "AL", "LA", "LL"
#     Only "AA" won't be regarded as rewardable owing to more than one absent times.
#
#
# **Note:** The value of **n** won't exceed 100,000.
#
#
## @lc code=start
using LeetCode

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


text:
# Infeasible methods
function InfeasibleProblem(prob::Problem, Z0::SampledTrajectory, R_inf::Real)
    @assert !isnan(sum(sum.(states(Z0))))

    nx,nu = dims(prob)  # original sizes
    N = TO.horizonlength(prob)

    # Create model with augmented controls
    model_inf = InfeasibleModel.(prob.model)

    # Get a trajectory that is dynamically feasible for the augmented problem
    #   and matches the states and controls of the original guess
    Z = infeasible_trajectory(model_inf, Z0)

    # Convert constraints so that they accept new dimensions
    conSet = TO.change_dimension(get_constraints(prob), nx[1], nu[1]+nx[1], 1:nx[1], 1:nu[1])

    # Constrain additional controls to be zero
    inf = InfeasibleConstraint(model_inf[1])  # TODO: handle multiple sizes
    TO.add_constraint!(conSet, inf, 1:N-1)

    # Infeasible Objective
    obj = infeasible_objective(prob.obj, R_inf)

    # Create new problem
    Problem(model_inf, obj, conSet, prob.x0, prob.xf, Z, N, TO.get_initial_time(prob),
            TO.get_final_time(prob))
end

function infeasible_objective(obj::Objective, regularizer)
    n,m = TO.state_dim(obj.cost[1]), TO.control_dim(obj.cost[1])
    Rd = [@SVector zeros(m); @SVector fill(regularizer,n)]
    R = Diagonal(Rd)
    cost_inf = TO.DiagonalCost(Diagonal(@SVector zeros(n)), R, checks=false)
    costs = map(obj.cost) do cost
        cost_idx = TO.change_dimension(cost, n, n+m, 1:n, 1:m)
        cost_idx + cost_inf
    end
    TO.Objective(costs)
end

############################################################################################
#                               INFEASIBLE MODELS                                          #
############################################################################################
RD.@autodiff struct InfeasibleModel{Nx,Nu,D<:DiscreteDynamics} <: DiscreteDynamics
    model::D
    _u::SVector{Nu,Int}  # inds to original controls
    _ui::SVector{Nx,Int} # inds to infeasible controls
    function InfeasibleModel{Nx,Nu}(model::D) where {Nx,Nu,D<:DiscreteDynamics}
        _u = SVector{Nu}(1:Nu)
        _ui = SVector{Nx}((1:Nx) .+ Nu)
        new{Nx,Nu,D}(model, _u, _ui)
    end
end

"""
    InfeasibleModel

An infeasible model is an augmented dynamics model that makes the system artifically fully
actuated by augmenting the control vector with `n` additional controls. The dynamics are
handled explicitly in discrete time:

``x_{k+1} = f(x_k,u_k,dt) + w_k``

where ``w_k`` are the additional `n`-dimensional controls. In practice, these are constrained
to be zero by the end of the solve.

# Constructors
```julia
InfeasibleModel(model::AbstractModel)
```
"""
InfeasibleModel

function InfeasibleModel(model::AbstractModel)
    n,m = RD.dims(model)
    InfeasibleModel{n,m}(model)
end

RD.statevectortype(::Type{<:InfeasibleModel{<:Any,<:Any,D}}) where D = RD.statevectortype(D)
RD.LieState(model::InfeasibleModel) = RD.LieState(model.model)
@inline RD.rotation_type(model::InfeasibleModel) where D = rotation_type(model.model)
RD.default_diffmethod(model::InfeasibleModel) = RD.default_diffmethod(model.model)
RD.default_signature(model::InfeasibleModel) = RD.default_signature(model.model)

# Generic Infeasible Methods
@inline RD.state_dim(model::InfeasibleModel{n}) where n = n
@inline RD.control_dim(model::InfeasibleModel{n,m}) where {n,m} = n+m
@inline RD.errstate_dim(model::InfeasibleModel) = RD.errstate_dim(model.model)

function RD.discrete_dynamics(model::InfeasibleModel,
        x, u, t, dt)
    u0 = u[model._u]
    ui = u[model._ui]
    RobotDynamics.discrete_dynamics(model.model, x, u0, t, dt) + ui
end

function RD.discrete_dynamics!(model::InfeasibleModel{Nx,Nu}, xn,
        x, u, t, dt) where {Nx,Nu}
    u0 = view(u, 1:Nu)
    ui = view(u, Nu+1:Nx+Nu)
    RobotDynamics.discrete_dynamics!(model.model, xn, x, u0, t, dt)
    xn .+= ui
    return
end

function RD.jacobian!(::RD.InPlace, diffmethod::DiffMethod, model::InfeasibleModel{Nx,Nu},
                      J, y, z) where {Nx,Nu}
    n,m = RD.dims(model.model)

    # Evaluate original Jacobian
    J0 = view(J, :, 1:n+m)
    z0 = RD.StaticKnotPoint{Nx,Nu}(n,m, view(z.z, 1:n+m), z.t, z.dt)
    RD.jacobian!(RD.InPlace(), diffmethod, model.model, J0, y, z0)

    # Fill in the rest manually
    for i = 1:n
        J[i, n+m+i] = 1
    end
    return
end

@inline RD.state_diff(model::InfeasibleModel, x::SVector, x0::SVector) =
    RD.state_diff(model.model, x, x0)

@inline RobotDynamics.state_diff_jacobian!(G, model::InfeasibleModel, Z::SampledTrajectory) =
	RobotDynamics.state_diff_jacobian!(G, model.model, Z)

@inline RobotDynamics.∇²differential!(∇G, model::InfeasibleModel, x::SVector, dx::SVector) =
    ∇²differential!(∇G, model.model, x, dx)

Base.position(model::InfeasibleModel, x::SVector) = position(model.model, x)

RobotDynamics.orientation(model::InfeasibleModel, x::SVector) = orientation(model.model, x)

"Calculate a dynamically feasible initial trajectory for an infeasible problem, given a
desired trajectory"
function infeasible_trajectory(models::Vector{<:InfeasibleModel}, Z0::SampledTrajectory{Nx,Nu}) where {Nx,Nu}
    if Nx != Any && Nu != Any
        Nū = Nx + Nu
    else
        Nū = Any
    end
    Z = map(Z0) do z
        n = state_dim(z)
        m = control_dim(z)
        ui = @SVector zeros(n)
        KnotPoint{Nx,Nū}(n, m + n, [RD.getdata(z); ui], z.t, z.dt)
    end
    N = length(Z0)
    for k = 1:N-1
        n = state_dim(models[k])
        RD.propagate_dynamics!(RD.default_signature(models[k]), models[k], Z[k+1], Z[k])
        x′ = state(Z[k+1])
        u_slack = state(Z0[k+1]) - x′
        u = [control(Z0[k]); u_slack]
        RD.setcontrol!(Z[k], u)
        RD.setstate!(Z[k+1], x′ + u_slack)
    end
    return SampledTrajectory(Z)
end

############################################################################################
#  								INFEASIBLE CONSTRAINT 									   #
############################################################################################
"""
    InfeasibleConstraint

Constraints additional ``infeasible'' controls to be zero.
Constructors: ```julia
InfeasibleConstraint(model::InfeasibleModel)
InfeasibleConstraint(n,m)
```
"""
struct InfeasibleConstraint{Nx,Nu} <: TO.ControlConstraint
	ui::SVector{Nx,Int}
	function InfeasibleConstraint{Nx,Nu}() where {Nx,Nu}
		ui = SVector{Nx}((1:Nx) .+ Nu)
		new{Nx,Nu}(ui)
	end
end

InfeasibleConstraint(model::InfeasibleModel{n,m}) where {n,m} = InfeasibleConstraint{n,m}()
RobotDynamics.state_dim(con::InfeasibleConstraint{n,m}) where {n,m} = n
RobotDynamics.control_dim(con::InfeasibleConstraint{n,m}) where {n,m} = n+m
@inline TO.sense(::InfeasibleConstraint) = TO.Equality()
@inline RD.output_dim(::InfeasibleConstraint{n}) where n = n

RD.evaluate(con::InfeasibleConstraint, u::RD.DataVector) = u[con.ui] # infeasible controls
function RD.evaluate!(con::InfeasibleConstraint, c, u::RD.DataVector)
    for (i,j) in enumerate(con.ui)
        c[i] = u[j]
    end
    return
end

function RD.jacobian!(con::InfeasibleConstraint{Nx}, ∇c, c, u::RD.DataVector) where {Nx}
	for (i,j) in enumerate(con.ui)
		∇c[i,j] = 1
    end
    return
end


text:
##==============================================================================
## Accessors
##==============================================================================
##==============================================================================
## PointParametricEst
##==============================================================================
"$(SIGNATURES)"
getMaxPPE(est::AbstractPointParametricEst) = est.max
"$(SIGNATURES)"
getMeanPPE(est::AbstractPointParametricEst) = est.mean
"$(SIGNATURES)"
getSuggestedPPE(est::AbstractPointParametricEst) = est.suggested
"$(SIGNATURES)"
getLastUpdatedTimestamp(est::AbstractPointParametricEst) = est.lastUpdatedTimestamp

##==============================================================================
## Variable Node Data
##==============================================================================

## COMMON
# getSolveInProgress
# isSolveInProgress

##------------------------------------------------------------------------------
## softtype
##------------------------------------------------------------------------------
"""
   $(SIGNATURES)

Variable nodes softtype information holding a variety of meta data associated with the type of variable stored in that node of the factor graph.

Related

getVariableType
"""
function getSofttype(vnd::VariableNodeData)
  return vnd.softtype
end

getSofttype(v::DFGVariable, solvekey::Symbol=:default) = getSofttype(getSolverData(v, solvekey))

getSofttype(dfg::AbstractDFG, lbl::Symbol, solvekey::Symbol=:default) = getSofttype(getVariable(dfg,lbl), solvekey)

"""
    getVariableType

Alias for [`getSofttype`](@ref).
"""
getVariableType(args...) = getSofttype(args...)

##------------------------------------------------------------------------------
## solvedCount
##------------------------------------------------------------------------------

"""
    $SIGNATURES

Get the number of times a variable has been inferred -- i.e. `solvedCount`.

Related

isSolved, setSolvedCount!
"""
getSolvedCount(v::VariableNodeData) = v.solvedCount
getSolvedCount(v::VariableDataLevel2, solveKey::Symbol=:default) = getSolverData(v, solveKey) |> getSolvedCount
getSolvedCount(dfg::AbstractDFG, sym::Symbol, solveKey::Symbol=:default) = getSolvedCount(getVariable(dfg, sym), solveKey)

"""
    $SIGNATURES

Update/set the `solveCount` value.

Related

getSolved, isSolved
"""
setSolvedCount!(v::VariableNodeData, val::Int) = v.solvedCount = val
setSolvedCount!(v::VariableDataLevel2, val::Int, solveKey::Symbol=:default) = setSolvedCount!(getSolverData(v, solveKey), val)
setSolvedCount!(dfg::AbstractDFG, sym::Symbol, val::Int, solveKey::Symbol=:default) = setSolvedCount!(getVariable(dfg, sym), val, solveKey)

"""
    $SIGNATURES

Boolean on whether the variable has been solved.

Related

getSolved, setSolved!
"""
isSolved(v::VariableNodeData) = 0 < v.solvedCount
isSolved(v::VariableDataLevel2, solveKey::Symbol=:default) = getSolverData(v, solveKey) |> isSolved
isSolved(dfg::AbstractDFG, sym::Symbol, solveKey::Symbol=:default) = isSolved(getVariable(dfg, sym), solveKey)

##------------------------------------------------------------------------------
## initialized
##------------------------------------------------------------------------------
"""
    $SIGNATURES

Returns state of vertex data `.initialized` flag.

Notes:
- used by both factor graph variable and Bayes tree clique logic.
"""
function isInitialized(var::DFGVariable, key::Symbol=:default)::Bool
      data = getSolverData(var, key)
      if data == nothing
          #TODO we still have a mixture of 2 error behaviours
        return false
      else
          return data.initialized
    end
end

function isInitialized(dfg::AbstractDFG, label::Symbol, key::Symbol=:default)::Bool
  return isInitialized(getVariable(dfg, label), key)
end

##==============================================================================
## Variables
##==============================================================================
#
# |                     | label | tags | timestamp | ppe | softtypename | solvable | solverData | smallData | bigData |
# |---------------------|:-----:|:----:|:---------:|:---:|:------------:|:--------:|:----------:|:---------:|:-------:|
# | SkeletonDFGVariable |   X   |   X  |           |     |              |          |            |           |         |
# | DFGVariableSummary  |   X   |   X  |     X     |  X  |       X      |          |            |           |    X    |
# | DFGVariable         |   X   |   X  |     x     |  X  |              |     X    |      X     |     X     |    X    |
#
##------------------------------------------------------------------------------

##------------------------------------------------------------------------------
## label
##------------------------------------------------------------------------------

## COMMON
# getLabel

##------------------------------------------------------------------------------
## tags
##------------------------------------------------------------------------------

## COMMON
# getTags
# setTags!

##------------------------------------------------------------------------------
## timestamp
##------------------------------------------------------------------------------

## COMMON
# getTimestamp

"""
    $SIGNATURES

Set the timestamp of a DFGVariable object returning a new DFGVariable.
Note:
Since the `timestamp` field is not mutable `setTimestamp` returns a new variable with the updated timestamp (note the absence of `!`).
Use [`updateVariable!`](@ref) on the returened variable to update it in the factor graph if needed. Alternatively use [`setTimestamp!`](@ref).
See issue #315.
"""
function setTimestamp(v::DFGVariable, ts::DateTime; verbose::Bool=true)
    if verbose
        @warn "verbose=true: setTimestamp(::DFGVariable,...) creates a returns a new immutable DFGVariable object (and didn't change a distributed factor graph object), make sure you are using the right pointers: getVariable(...).  See setTimestamp!(...) and note suggested use is at addVariable!(..., [timestamp=...]).  See DFG #315 for explanation."
    end
    return DFGVariable(v.label, ts, v.nstime, v.tags, v.ppeDict, v.solverDataDict, v.smallData, v.bigData, Ref(v.solvable))
end

function setTimestamp(v::DFGVariableSummary, ts::DateTime; verbose::Bool=true)
    if verbose
        @warn "verbose=true: setTimestamp(::DFGVariableSummary,...) creates and returns a new immutable DFGVariable object (and didn't change a distributed factor graph object), make sure you are using the right pointers: getVariable(...).  See setTimestamp!(...) and note suggested use is at addVariable!(..., [timestamp=...]).  See DFG #315 for explanation."
    end
    return DFGVariableSummary(v.label, ts, v.tags, v.ppeDict, v.softtypename, v.bigData)
end

##------------------------------------------------------------------------------
## solvable
##------------------------------------------------------------------------------

## COMMON: solvable
# getSolvable
# setSolvable!
# isSolvable

## COMMON:

##------------------------------------------------------------------------------
## ppeDict
##------------------------------------------------------------------------------

"""
    $SIGNATURES

Get the PPE dictionary for a variable.  Recommended to use CRUD operations instead, [`getPPE`](@ref), [`addPPE!`](@ref), [`updatePPE!`](@ref), [`deletePPE!`](@ref).
"""
getPPEDict(v::VariableDataLevel1) = v.ppeDict

#TODO FIXME don't know if this should exist, should rather always update with fg object to simplify inmem vs cloud
"""
    $SIGNATURES

Get the parametric point estimate (PPE) for a variable in the factor graph.

Notes
- Defaults on keywords `solveKey` and `method`

Related

getMeanPPE, getMaxPPE, getKDEMean, getKDEFit, getPPEs, getVariablePPEs
"""
function getPPE(vari::VariableDataLevel1, solveKey::Symbol=:default)
    return  getPPEDict(vari)[solveKey]
    # return haskey(ppeDict, solveKey) ? ppeDict[solveKey] : nothing
end

# afew more aliases on PPE, brought back from deprecated DF

"""
    $SIGNATURES

Return full dictionary of PPEs in a variable, recommended to rather use CRUD: [`getPPE`](@ref),
"""
getVariablePPEDict(vari::VariableDataLevel1) = getPPEDict(vari)

getVariablePPE(args...) = getPPE(args...)

##------------------------------------------------------------------------------
## solverDataDict
##------------------------------------------------------------------------------

"""
    $SIGNATURES

Get solver data dictionary for a variable.  Advised to use graph CRUD operations instead.
"""
getSolverDataDict(v::DFGVariable) = v.solverDataDict

# TODO move to crud, don't know if this should exist, should rather always update with fg object to simplify inmem vs cloud
"""
    $SIGNATURES

Retrieve solver data structure stored in a variable.
"""
function getSolverData(v::DFGVariable, key::Symbol=:default)
    #TODO this does not fit in with some of the other error behaviour. but its used so added @error
    vnd =  haskey(v.solverDataDict, key) ? v.solverDataDict[key] : (@error "Variable does not have solver data $(key)"; nothing)
    return vnd
end

#TODO Repeated functionality? same as update
"""
    $SIGNATURES
Set solver data structure stored in a variable.
"""
setSolverData!(v::DFGVariable, data::VariableNodeData, key::Symbol=:default) = v.solverDataDict[key] = data

##------------------------------------------------------------------------------
## smallData
##------------------------------------------------------------------------------

"""
    $(SIGNATURES)
Get the small data for a variable.
"""
getSmallData(v::DFGVariable)::Dict{String, String} = v.smallData

"""
    $(SIGNATURES)
Set the small data for a variable.
This will overwrite old smallData.
"""
function setSmallData!(v::DFGVariable, smallData::Dict{String, String})::Dict{String, String}
    empty!(v.smallData)
    merge!(v.smallData, smallData)
end

##------------------------------------------------------------------------------
## bigData
##------------------------------------------------------------------------------

## see Bigdata Folder

##------------------------------------------------------------------------------
## softtypename
##------------------------------------------------------------------------------
## getter in DFGVariableSummary only
## can be utility function for others
## TODO this should return the softtype object, or try to. it should be getSofttypename for the accessor
## TODO Consider parameter N in softtype for dims, and storing constructor in softtypename
## TODO or just not having this function at all
# getSofttype(v::DFGVariableSummary) = v.softypename()
##------------------------------------------------------------------------------

"""
    $SIGNATURES
Retrieve the soft type name symbol for a DFGVariableSummary. ie :Point2, Pose2, etc.
"""
getSofttypename(v::DFGVariableSummary)::Symbol = v.softtypename

function getSofttype(v::DFGVariableSummary)::InferenceVariable
    @warn "Looking for type in `Main`. Only use if softtype has only one implementation,ie. Pose2. Otherwise use the full variable."
    return getfield(Main, v.softtypename)()
end

##==============================================================================
## Layer 2 CRUD and SET
##==============================================================================

##==============================================================================
## TAGS - See CommonAccessors
##==============================================================================

##==============================================================================
## Variable Node Data
##==============================================================================
##------------------------------------------------------------------------------
## CRUD: get, add, update, delete
##------------------------------------------------------------------------------

"""
    $(SIGNATURES)
Get variable solverdata for a given solve key.
"""
function getVariableSolverData(dfg::AbstractDFG, variablekey::Symbol, solvekey::Symbol=:default)::VariableNodeData
    v = getVariable(dfg, variablekey)
    !haskey(v.solverDataDict, solvekey) && error("Solve key '$solvekey' not found in variable '$variablekey'")
    return v.solverDataDict[solvekey]
end

"""
    $(SIGNATURES)
Add variable solver data, errors if it already exists.
"""
function addVariableSolverData!(dfg::AbstractDFG, variablekey::Symbol, vnd::VariableNodeData, solvekey::Symbol=:default)::VariableNodeData
    var = getVariable(dfg, variablekey)
    if haskey(var.solverDataDict, solvekey)
        error("VariableNodeData '$(solvekey)' already exists")
    end
    var.solverDataDict[solvekey] = vnd
    return vnd
end

"""
    $(SIGNATURES)
Add a new solver data  entry from a deepcopy of the source variable solver data.
NOTE: Copies the solver data.
"""
addVariableSolverData!(dfg::AbstractDFG, sourceVariable::DFGVariable, solvekey::Symbol=:default) =
    addVariableSolverData!(dfg, sourceVariable.label, deepcopy(getSolverData(sourceVariable, solvekey)), solvekey)

"""
    $(SIGNATURES)
Update variable solver data if it exists, otherwise add it.

Notes:
- `useCopy=true` to copy solver data and keep separate memory.
- Use `fields` to updated only a few VND.fields while adhering to `useCopy`.

Related

mergeVariableSolverData!
"""
function updateVariableSolverData!(dfg::AbstractDFG,
                                   variablekey::Symbol,
                                   vnd::VariableNodeData,
                                   solvekey::Symbol=:default,
                                   useCopy::Bool=true,
                                   fields::Vector{Symbol}=Symbol[],
                                   verbose::Bool=true  )
    #
    #This is basically just setSolverData
    var = getVariable(dfg, variablekey)
    if verbose && !haskey(var.solverDataDict, solvekey)
        @warn "VariableNodeData '$(solvekey)' does not exist, adding"
    end

    # for InMemoryDFGTypes do memory copy or repointing, for cloud this would be an different kind of update.
    usevnd = useCopy ? deepcopy(vnd) : vnd
    # should just one, or many pointers be updated?
    if haskey(var.solverDataDict, solvekey) && isa(var.solverDataDict[solvekey], VariableNodeData) && length(fields) != 0
      # change multiple pointers inside the VND var.solverDataDict[solvekey]
      for field in fields
        destField = getfield(var.solverDataDict[solvekey], field)
        srcField = getfield(usevnd, field)
        if isa(destField, Array) && size(destField) == size(srcField)
          # use broadcast (in-place operation)
          destField .= srcField
        else
          # change pointer of destination VND object member
          setfield!(var.solverDataDict[solvekey], field, srcField)
        end
      end
    else
      # change a single pointer in var.solverDataDict
      var.solverDataDict[solvekey] = usevnd
    end

    return var.solverDataDict[solvekey]
end

updateVariableSolverData!(dfg::AbstractDFG,
                          sourceVariable::DFGVariable,
                          solvekey::Symbol=:default,
                          useCopy::Bool=true,
                          fields::Vector{Symbol}=Symbol[] ) =
    updateVariableSolverData!(dfg, sourceVariable.label, getSolverData(sourceVariable, solvekey), solvekey, useCopy, fields)

function updateVariableSolverData!(dfg::AbstractDFG,
                                   sourceVariables::Vector{<:DFGVariable},
                                   solvekey::Symbol=:default,
                                   useCopy::Bool=true,
                                   fields::Vector{Symbol}=Symbol[]  )
    #I think cloud would do this in bulk for speed
    for var in sourceVariables
        updateVariableSolverData!(dfg, var.label, getSolverData(var, solvekey), solvekey, useCopy, fields)
    end
end

"""
    $(SIGNATURES)
Delete variable solver data, returns the deleted element.
"""
function deleteVariableSolverData!(dfg::AbstractDFG, variablekey::Symbol, solvekey::Symbol=:default)::VariableNodeData
    var = getVariable(dfg, variablekey)

    if !haskey(var.solverDataDict, solvekey)
        error("VariableNodeData '$(solvekey)' does not exist")
    end
    vnd = pop!(var.solverDataDict, solvekey)
    return vnd
end

"""
    $(SIGNATURES)
Delete variable solver data, returns the deleted element.
"""
deleteVariableSolverData!(dfg::AbstractDFG, sourceVariable::DFGVariable, solvekey::Symbol=:default) =
    deleteVariableSolverData!(dfg, sourceVariable.label, solvekey)

##------------------------------------------------------------------------------
## SET: list, merge
##------------------------------------------------------------------------------

"""
    $(SIGNATURES)
List all the solver data keys in the variable.
"""
function listVariableSolverData(dfg::AbstractDFG, variablekey::Symbol)::Vector{Symbol}
    v = getVariable(dfg, variablekey)
    return collect(keys(v.solverDataDict))
end

"""
    $(SIGNATURES)
Merges and updates solver and estimate data for a variable (variable can be from another graph).
If the same key is present in another collection, the value for that key will be the value it has in the last collection listed (updated).
Note: Makes a copy of the estimates and solver data so that there is no coupling between graphs.
"""
function mergeVariableSolverData!(destVariable::DFGVariable, sourceVariable::DFGVariable)::DFGVariable
    # We don't know which graph this came from, must be copied!
    merge!(destVariable.solverDataDict, deepcopy(sourceVariable.solverDataDict))
    return destVariable
end

mergeVariableSolverData!(dfg::AbstractDFG, sourceVariable::DFGVariable) =     mergeVariableSolverData!(getVariable(dfg,getLabel(sourceVariable)), sourceVariable)

##==============================================================================
## Point Parametric Estimates
##==============================================================================

##------------------------------------------------------------------------------
## CRUD: get, add, update, delete
##------------------------------------------------------------------------------

"""
    $(SIGNATURES)
Get the parametric point estimate (PPE) for a variable in the factor graph for a given solve key.

Notes
- Defaults on keywords `solveKey` and `method`

Related
getMeanPPE, getMaxPPE, getKDEMean, getKDEFit, getPPEs, getVariablePPEs
"""
function getPPE(dfg::AbstractDFG, variablekey::Symbol, ppekey::Symbol=:default)::AbstractPointParametricEst
    v = getVariable(dfg, variablekey)
    !haskey(v.ppeDict, ppekey) && error("PPE key '$ppekey' not found in variable '$variablekey'")
    return v.ppeDict[ppekey]
end

# Not the most efficient call but it at least reuses above (in memory it's probably ok)
getPPE(dfg::AbstractDFG, sourceVariable::VariableDataLevel1, ppekey::Symbol=default)::AbstractPointParametricEst = getPPE(dfg, sourceVariable.label, ppekey)

"""
    $(SIGNATURES)
Add variable PPE, errors if it already exists.
"""
function addPPE!(dfg::AbstractDFG, variablekey::Symbol, ppe::P, ppekey::Symbol=:default)::AbstractPointParametricEst where P <: AbstractPointParametricEst
    var = getVariable(dfg, variablekey)
    if haskey(var.ppeDict, ppekey)
        error("PPE '$(ppekey)' already exists")
    end
    var.ppeDict[ppekey] = ppe
    return ppe
end

"""
    $(SIGNATURES)
Add a new PPE entry from a deepcopy of the source variable PPE.
NOTE: Copies the solver data.
"""
addPPE!(dfg::AbstractDFG, sourceVariable::DFGVariable, ppekey::Symbol=:default) =
    addPPE!(dfg, sourceVariable.label, deepcopy(getPPE(sourceVariable, ppekey)), ppekey)

"""
    $(SIGNATURES)
Update PPE data if it exists, otherwise add it -- one call per `key::Symbol=:default`.
"""
function updatePPE!(dfg::AbstractDFG, variablekey::Symbol, ppe::P, ppekey::Symbol=:default)::P where P <: AbstractPointParametricEst

    var = getVariable(dfg, variablekey)
    if !haskey(var.ppeDict, ppekey)
        @warn "PPE '$(ppekey)' does not exist, adding"
    end
    #for InMemoryDFGTypes, cloud would update here
    var.ppeDict[ppekey] = ppe
    return ppe
end

"""
    $(SIGNATURES)
Update PPE data if it exists, otherwise add it.
NOTE: Copies the PPE data.
"""
updatePPE!(dfg::AbstractDFG, sourceVariable::VariableDataLevel1, ppekey::Symbol=:default) =
    updatePPE!(dfg, sourceVariable.label, deepcopy(getPPE(sourceVariable, ppekey)), ppekey)

"""
    $(SIGNATURES)
Update PPE data if it exists, otherwise add it.
"""
function updatePPE!(dfg::AbstractDFG, sourceVariables::Vector{<:VariableDataLevel1}, ppekey::Symbol=:default)
    #I think cloud would do this in bulk for speed
    for var in sourceVariables
        updatePPE!(dfg, var.label, getPPE(dfg, var, ppekey), ppekey)
    end
end

"""
    $(SIGNATURES)
Delete PPE data, returns the deleted element.
"""
function deletePPE!(dfg::AbstractDFG, variablekey::Symbol, ppekey::Symbol=:default)::AbstractPointParametricEst
    var = getVariable(dfg, variablekey)

    if !haskey(var.ppeDict, ppekey)
        error("VariableNodeData '$(ppekey)' does not exist")
    end
    vnd = pop!(var.ppeDict, ppekey)
    return vnd
end

"""
    $(SIGNATURES)
Delete PPE data, returns the deleted element.
"""
deletePPE!(dfg::AbstractDFG, sourceVariable::DFGVariable, ppekey::Symbol=:default) =
    deletePPE!(dfg, sourceVariable.label, ppekey)

##------------------------------------------------------------------------------
## SET: list, merge
##------------------------------------------------------------------------------

"""
    $(SIGNATURES)
List all the PPE data keys in the variable.
"""
function listPPEs(dfg::AbstractDFG, variablekey::Symbol)::Vector{Symbol}
    v = getVariable(dfg, variablekey)
    return collect(keys(v.ppeDict))
end

#TODO API and only correct level
"""
    $(SIGNATURES)
Merges and updates solver and estimate data for a variable (variable can be from another graph).
Note: Makes a copy of the estimates and solver data so that there is no coupling between graphs.
"""
function mergePPEs!(destVariable::AbstractDFGVariable, sourceVariable::AbstractDFGVariable)::AbstractDFGVariable
    # We don't know which graph this came from, must be copied!
    merge!(destVariable.ppeDict, deepcopy(sourceVariable.ppeDict))
    return destVariable
end


text:
module Includev
# jdad's include with verbose, like source(...,verbose=true) in R

using Printf
export includev, includeve

# Utilities : functions ends_with_semicolon, ends_with_newline

# Inspired from Base.chomp
function ends_with_newline(s::AbstractString)
    i = lastindex(s)
    (i < 1 || s[i] != '\n') && (return false)
    return true
    #(i < 1 || s[i] != '\n') && (return SubString(s, 1, i))
    #j = prevind(s,i)
    #(j < 1 || s[j] != '\r') && (return SubString(s, 1, j))
    #return SubString(s, 1, prevind(s,j))
end

# originally from julia-1.1.0/share/julia/stdlib/v1.1/REPL/src/REPL.jl
function ends_with_semicolon(line::AbstractString)
    match = findlast(isequal(';'), line)
    if match !== nothing
        # state for comment parser, assuming that the `;` isn't in a string or comment
        # so input like ";#" will still thwart this to give the wrong (anti-conservative) answer
        comment = false
        comment_start = false
        comment_close = false
        comment_multi = 0
        for c in line[(match + 1):end]
            if comment_multi > 0
                # handle nested multi-line comments
                if comment_close && c == '#'
                    comment_close = false
                    comment_multi -= 1
                elseif comment_start && c == '='
                    comment_start = false
                    comment_multi += 1
                else
                    comment_start = (c == '#')
                    comment_close = (c == '=')
                end
            elseif comment
                # handle line comments
                if c == '\r' || c == '\n'
                    comment = false
                end
            elseif comment_start
                # see what kind of comment this is
                comment_start = false
                if c == '='
                    comment_multi = 1
                else
                    comment = true
                end
            elseif c == '#'
                # start handling for a comment
                comment_start = true
            else
                # outside of a comment, encountering anything but whitespace
                # means the semi-colon was internal to the expression
                isspace(c) || return false
            end
        end
        return true
    end
    return false
end

# utility : prettytime / from BenchmarkTools
function prettytime(t)
    if t < 1e3
        value, units = t, "ns"
    elseif t < 1e6
        value, units = t / 1e3, "μs"
    elseif t < 1e9
        value, units = t / 1e6, "ms"
    else
        value, units = t / 1e9, "s"
    end
    return string(@sprintf("%.3f", value), " ", units)
end

# end utilities

"""
    includev(filetoinc::AbstractString; echo=true, elaps=false, logfile=nothing, debug=false)

"""
function includev(filetoinc::AbstractString; echo=true, elaps=false, logfile=nothing, debug=false)
  scalet=1.0e-9
  t0=time_ns()
  open(filetoinc) do f
    il=0
    ile=999
    ilx=9999
    hit_eof = false
    prompt=filetoinc * "> "
    xps1="> "
    xps2=". "
    td=(time_ns()-t0)*scalet
    tf=td;
    dt=tf-td
    while true # 1
        ps=""
        #td=time_ns()*scalet
        if elaps
            ps=@sprintf("+%9.6f =%9.6f",dt,td)
        end
        ps1= ps * xps1
        ps2= ps * xps2
        line = ""
        ast = nothing
        interrupted = false
        lined=""
        lines=""
        firstline = true
        while true # 2
            try
            debug && write(stdout,"\n$il linebefor="*line)
            oneline = readline(f, keep=true)
            line = line * oneline
            if firstline
                lined = oneline
                firstline = false
                echo && write(stdout, ps1 * oneline)
            else
                lines= lines * oneline
                echo && write(stdout, ps2 * oneline)
            end
			if !ends_with_newline(oneline)
                echo && write(stdout, '\n')
			end
            debug && write(stdout,"\n$il lineafter="*line)
            il += 1
            catch e
                if isa(e,InterruptException)
                    try # raise the debugger if present
                        ccall(:jl_raise_debugger, Int, ())
                    catch
                    end
                    line = ""
                    interrupted = true
                    break
                elseif isa(e,EOFError)
                    hit_eof = true
                    break
                else
                    rethrow()
                end
            end
            ast = Base.parse_input_line(line)
            (isa(ast,Expr) && ast.head == :incomplete) || break
            debug && @show isa(ast,Expr)
            debug && @show ast.head
            (ile>0) && (il > ile) && @info "includev : max $ilx number of lines in expr reached" && break
            (ilx>0) && (il > ilx) && @info "includev : max $ilx number of lines to process reached" && break
        end # while true 2
		# Now we have (hopefully) a group of lines with complete expression
        ast = Base.parse_input_line(line)
        if isa(ast,Expr)
            debug && write(stdout,"evaling=")
            line_to_eval = line
            td = (time_ns()-t0)*scalet
            include_string(Main,line_to_eval)
            tf = (time_ns()-t0)*scalet
            dt = tf-td
            debug && write(stdout,"...evaling complete.")
            if !ends_with_semicolon(line)
            debug && write(stdout,"...should output response?") # no ans as of now
            end
        else
            debug && write(stdout,"not_an expr-ignoring=")
        end
        ((!interrupted && isempty(line)) || hit_eof) && break
        (ilx>0) && (il > ilx) && @info "includev : max $ilx number of lines to process reached" && break
    end # while true # 1
  end # do f
  # terminate backend
  # as of now,  do not get last ans as returned result
  # nothing
end # function

includeve(filetoinc::AbstractString ;echo=true, logfile=nothing, debug=false) = includev(filetoinc; echo=echo, elaps=true, logfile=logfile, debug=debug)

end # module


text:
include("controller.jl")
include("agent.jl")

mutable struct DifferentiableAgent <: Agent
end

mutable struct DifferentiableController <: Controller
    agents::Dict{String, DifferentiableAgent}
end

function add(controller::DifferentiableController, task::Dict{String, Any})
    @info "differentiable task added" task

    # check that solution doesn't already exist
    # create
end

function update(controller::DifferentiableController, task::Dict{String, Any})
    @info "differentiable task updated" task
end

function delete(controller::DifferentiableController, task::Dict{String, Any})
    @info "differentiable task deleted" task
end

max_stars_repo_name: pbarker/kubemeta

max_issues_repo_name: pbarker/kubemeta

max_forks_repo_name: pbarker/kubemeta


text:
using Roots
using Compat.Test
import SpecialFunctions.erf

include("./test_find_zero.jl")
#include("./RootTesting.jl")
#run_benchmark_tests()

include("./test_fzero.jl")

#include("./test_fzero3.jl")
#run_robustness_test()

include("./test_newton.jl")

#include("./test_derivative_free.jl")


text:
push!(Base.LOAD_PATH,joinpath(dirname(@__FILE__),"..",".."))

using ProxSDP, MathOptInterface, Test, LinearAlgebra, Random, SparseArrays, DelimitedFiles

const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIB = MOI.Bridges
const MOIU = MOI.Utilities

const cache = MOIU.UniversalFallback(MOIU.Model{Float64}())

const optimizer = MOIU.CachingOptimizer(cache, ProxSDP.Optimizer(tol_primal = 1e-6, tol_dual = 1e-6, log_verbose = false))
const optimizer_low_acc = MOIU.CachingOptimizer(cache, ProxSDP.Optimizer(tol_primal = 1e-3, tol_dual = 1e-3, log_verbose = true, timer_verbose = true))
const optimizer_full = MOIU.CachingOptimizer(cache, ProxSDP.Optimizer(full_eig_decomp = true, tol_primal = 1e-4, tol_dual = 1e-4))
const optimizer_print = MOIU.CachingOptimizer(cache, ProxSDP.Optimizer(log_freq = 10, log_verbose = true, timer_verbose = true, tol_primal = 1e-4, tol_dual = 1e-4))
const config = MOIT.TestConfig(atol=1e-3, rtol=1e-3, infeas_certificates = false)
const config_conic = MOIT.TestConfig(atol=1e-3, rtol=1e-3, duals = false, infeas_certificates = false)

@testset "SolverName" begin
    @test MOI.get(optimizer, MOI.SolverName()) == "ProxSDP"
end

@testset "supports_default_copy_to" begin
    @test MOIU.supports_allocate_load(ProxSDP.Optimizer(), false)
    @test !MOIU.supports_allocate_load(ProxSDP.Optimizer(), true)
end

@testset "Unit" begin
    bridged = MOIB.full_bridge_optimizer(optimizer, Float64)
    MOIT.unittest(bridged, config,[
        # Quadratic functions are not supported
        "solve_qcp_edge_cases", "solve_qp_edge_cases",
        # Integer and ZeroOne sets are not supported
        "solve_integer_edge_cases", "solve_objbound_edge_cases",
        "solve_zero_one_with_bounds_1",
        "solve_zero_one_with_bounds_2",
        "solve_zero_one_with_bounds_3",
                    # `TimeLimitSec` not supported.
                    # "time_limit_sec",
                    "number_threads",
                    # ArgumentError: The number of constraints in SCSModel must be greater than 0
                    "solve_unbounded_model",
        ]
    )
end

@testset "MOI Continuous Linear" begin
    bridged = MOIB.full_bridge_optimizer(optimizer, Float64)
    MOIT.contlineartest(bridged, config, [
        # infeasible/unbounded
        # "linear8a", "linear8b", "linear8c", "linear12",
        # linear10 is poorly conditioned
        "linear10",
        # linear9 is requires precision
        "linear9",
        # primalstart not accepted
        "partial_start",
        ]
    )
    # MOIT.linear9test(MOIB.SplitInterval{Float64}(optimizer), config)
end

@testset "MOI Continuous Conic" begin
    MOIT.contconictest(MOIB.full_bridge_optimizer(optimizer, Float64), config_conic, [
        # bridge: some problem with square psd
        "rootdets",
        # exp cone
        "logdet", "exp", "dualexp",
        # pow cone
        "pow","dualpow",
        # infeasible/unbounded
        # "lin3", "lin4",
        # See https://travis-ci.com/blegat/SolverTests/jobs/268551133
        # geomean2v: Test Failed at /home/travis/.julia/dev/MathOptInterface/src/Test/contconic.jl:1328
        # Expression: MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        #  Evaluated: INFEASIBLE_OR_UNBOUNDED::TerminationStatusCode = 6 == OPTIMAL::TerminationStatusCode = 1
        "geomean2v", "geomean2f",
        "soc3", "psdt2", "normone2", "norminf2",#, "rotatedsoc2"
                # TODO - should not fail
                #"psdt1v", "soc4"
        ]
    )
end

@testset "Simple LP" begin

    bridged = MOIB.full_bridge_optimizer(optimizer, Float64)
    MOI.empty!(bridged)
    @test MOI.is_empty(bridged)

    # add 10 variables - only diagonal is relevant
    X = MOI.add_variables(bridged, 2)

    # add sdp constraints - only ensuring positivenesse of the diagonal
    vov = MOI.VectorOfVariables(X)

    c1 = MOI.add_constraint(bridged,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(2.0, X[1]),
            MOI.ScalarAffineTerm(1.0, X[2])
        ], 0.0), MOI.EqualTo(4.0))

    c2 = MOI.add_constraint(bridged,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X[1]),
            MOI.ScalarAffineTerm(2.0, X[2])
        ], 0.0), MOI.EqualTo(4.0))

    b1 = MOI.add_constraint(bridged,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X[1])
        ], 0.0), MOI.GreaterThan(0.0))

    b2 = MOI.add_constraint(bridged,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X[2])
        ], 0.0), MOI.GreaterThan(0.0))

    MOI.set(bridged,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], [X[1], X[2]]), 0.0)
        )
    MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(bridged)

    obj = MOI.get(bridged, MOI.ObjectiveValue())

    @test obj ≈ -9.33333 atol = 1e-2

    Xr = MOI.get(bridged, MOI.VariablePrimal(), X)

    @test Xr ≈ [1.3333, 1.3333] atol = 1e-2

end

@testset "Simple LP with 2 1D SDP" begin

    bridged = MOIB.full_bridge_optimizer(optimizer, Float64)
    MOI.empty!(bridged)
    @test MOI.is_empty(bridged)

    # add 10 variables - only diagonal is relevant
    X = MOI.add_variables(bridged, 2)

    # add sdp constraints - only ensuring positivenesse of the diagonal
    vov = MOI.VectorOfVariables(X)

    c1 = MOI.add_constraint(bridged,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(2.0, X[1]),
            MOI.ScalarAffineTerm(1.0, X[2])
        ], 0.0), MOI.EqualTo(4.0))

    c2 = MOI.add_constraint(bridged,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X[1]),
            MOI.ScalarAffineTerm(2.0, X[2])
        ], 0.0), MOI.EqualTo(4.0))

    b1 = MOI.add_constraint(bridged,
        MOI.VectorOfVariables([X[1]]), MOI.PositiveSemidefiniteConeTriangle(1))

    b2 = MOI.add_constraint(bridged,
        MOI.VectorOfVariables([X[2]]), MOI.PositiveSemidefiniteConeTriangle(1))

    MOI.set(bridged,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], [X[1], X[2]]), 0.0)
        )
    MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(bridged)

    obj = MOI.get(bridged, MOI.ObjectiveValue())

    @test obj ≈ -9.33333 atol = 1e-2

    Xr = MOI.get(bridged, MOI.VariablePrimal(), X)

    @test Xr ≈ [1.3333, 1.3333] atol = 1e-2

end

@testset "LP in SDP EQ form" begin

    bridged = MOIB.full_bridge_optimizer(optimizer, Float64)
    MOI.empty!(bridged)
    @test MOI.is_empty(bridged)

    # add 10 variables - only diagonal is relevant
    X = MOI.add_variables(bridged, 10)

    # add sdp constraints - only ensuring positivenesse of the diagonal
    vov = MOI.VectorOfVariables(X)
    cX = MOI.add_constraint(bridged, vov, MOI.PositiveSemidefiniteConeTriangle(4))

    c1 = MOI.add_constraint(bridged,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(2.0, X[1]),
            MOI.ScalarAffineTerm(1.0, X[3]),
            MOI.ScalarAffineTerm(1.0, X[6])
        ], 0.0), MOI.EqualTo(4.0))

    c2 = MOI.add_constraint(bridged,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X[1]),
            MOI.ScalarAffineTerm(2.0, X[3]),
            MOI.ScalarAffineTerm(1.0, X[10])
        ], 0.0), MOI.EqualTo(4.0))

    MOI.set(bridged,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], [X[1], X[3]]), 0.0)
        )
    MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(bridged)

    obj = MOI.get(bridged, MOI.ObjectiveValue())

    @test obj ≈ -9.33333 atol = 1e-2

    Xr = MOI.get(bridged, MOI.VariablePrimal(), X)

    @test Xr ≈ [1.3333, .0, 1.3333, .0, .0, .0, .0, .0, .0, .0] atol = 1e-2

end

@testset "LP in SDP INEQ form" begin

    MOI.empty!(optimizer)
    @test MOI.is_empty(optimizer)

    # add 10 variables - only diagonal is relevant
    X = MOI.add_variables(optimizer, 3)

    # add sdp constraints - only ensuring positivenesse of the diagonal
    vov = MOI.VectorOfVariables(X)
    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(2))

    c1 = MOI.add_constraint(optimizer,
        MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1,1],[
            MOI.ScalarAffineTerm(2.0, X[1]),
            MOI.ScalarAffineTerm(1.0, X[3]),
        ]), [-4.0]), MOI.Nonpositives(1))

    c2 = MOI.add_constraint(optimizer,
        MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1,1],[
            MOI.ScalarAffineTerm(1.0, X[1]),
            MOI.ScalarAffineTerm(2.0, X[3]),
        ]), [-4.0]), MOI.Nonpositives(1))

    MOI.set(optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([4.0, 3.0], [X[1], X[3]]), 0.0)
        )
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.optimize!(optimizer)

    obj = MOI.get(optimizer, MOI.ObjectiveValue())

    @test obj ≈ 9.33333 atol = 1e-2

    Xr = MOI.get(optimizer, MOI.VariablePrimal(), X)

    @test Xr ≈ [1.3333, .0, 1.3333] atol = 1e-2

    c1_d = MOI.get(optimizer, MOI.ConstraintDual(), c1)
    c2_d = MOI.get(optimizer, MOI.ConstraintDual(), c2)

    # SDP duasl error
    @test_throws ErrorException c2_d = MOI.get(optimizer, MOI.ConstraintDual(), cX)

end

@testset "SDP from MOI" begin
    # min X[1,1] + X[2,2]    max y
    #     X[2,1] = 1         [0   y/2     [ 1  0
    #                         y/2 0    <=   0  1]
    #     X >= 0              y free
    # Optimal solution:
    #
    #     ⎛ 1   1 ⎞
    # X = ⎜       ⎟           y = 2
    #     ⎝ 1   1 ⎠
    MOI.empty!(optimizer)
    @test MOI.is_empty(optimizer)

    X = MOI.add_variables(optimizer, 3)

    vov = MOI.VectorOfVariables(X)
    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(2))

    c = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[2]))], [-1.0]), MOI.Zeros(1))

    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [X[1], X[3]]), 0.0))

    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(optimizer)

    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL

    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(optimizer, MOI.DualStatus()) == MOI.FEASIBLE_POINT

    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 2 atol=1e-2

    Xv = ones(3)
    @test MOI.get(optimizer, MOI.VariablePrimal(), X) ≈ Xv atol=1e-2
    # @test MOI.get(optimizer, MOI.ConstraintPrimal(), cX) ≈ Xv atol=1e-2

    # @test MOI.get(optimizer, MOI.ConstraintDual(), c) ≈ 2 atol=1e-2
    # @show MOI.get(optimizer, MOI.ConstraintDual(), c)

end

@testset "Double SDP from MOI" begin
    # solve simultaneously two of these:
    # min X[1,1] + X[2,2]    max y
    #     X[2,1] = 1         [0   y/2     [ 1  0
    #                         y/2 0    <=   0  1]
    #     X >= 0              y free
    # Optimal solution:
    #
    #     ⎛ 1   1 ⎞
    # X = ⎜       ⎟           y = 2
    #     ⎝ 1   1 ⎠
    MOI.empty!(optimizer)
    @test MOI.is_empty(optimizer)

    X = MOI.add_variables(optimizer, 3)
    Y = MOI.add_variables(optimizer, 3)

    vov = MOI.VectorOfVariables(X)
    vov2 = MOI.VectorOfVariables(Y)
    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(2))
    cY = MOI.add_constraint(optimizer, vov2, MOI.PositiveSemidefiniteConeTriangle(2))

    c = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[2]))], [-1.0]), MOI.Zeros(1))
    c2 = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, Y[2]))], [-1.0]), MOI.Zeros(1))

    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [X[1], X[end], Y[1], Y[end]]), 0.0))

    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(optimizer)

    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL

    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(optimizer, MOI.DualStatus()) == MOI.FEASIBLE_POINT

    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 2*2 atol=1e-2

    Xv = ones(3)
    @test MOI.get(optimizer, MOI.VariablePrimal(), X) ≈ Xv atol=1e-2
    Yv = ones(3)
    @test MOI.get(optimizer, MOI.VariablePrimal(), Y) ≈ Yv atol=1e-2
    # @test MOI.get(optimizer, MOI.ConstraintPrimal(), cX) ≈ Xv atol=1e-2

    # @test MOI.get(optimizer, MOI.ConstraintDual(), c) ≈ 2 atol=1e-2
    # @show MOI.get(optimizer, MOI.ConstraintDual(), c)

end

@testset "SDP with duplicates from MOI" begin

    using MathOptInterface
    # using SCS,ProxSDP
    MOI = MathOptInterface
    MOIU = MathOptInterface.Utilities
    MOIB = MathOptInterface.Bridges

    cache = MOIU.UniversalFallback(MOIU.Model{Float64}());
    #optimizer0 = SCS.Optimizer(linear_solver=SCS.Direct, eps=1e-8);
    optimizer0 = ProxSDP.Optimizer()#linear_solver=SCS.Direct, eps=1e-8);
    MOI.empty!(cache);
    optimizer1 = MOIU.CachingOptimizer(cache, optimizer0);
    optimizer = MOIB.full_bridge_optimizer(optimizer1, Float64);

    MOI.empty!(optimizer)

    x = MOI.add_variable(optimizer)
    X = [x, x, x]

    vov = MOI.VectorOfVariables(X)
    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(2))

    c = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[2]))], [-1.0]), MOI.Zeros(1))

    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [X[1], X[3]]), 0.0))

    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(optimizer)

    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL

    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(optimizer, MOI.DualStatus()) == MOI.FEASIBLE_POINT

    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 2 atol=1e-2

    Xv = ones(3)
    @test MOI.get(optimizer, MOI.VariablePrimal(), X) ≈ Xv atol=1e-2

end

@testset "SDP from Wikipedia" begin
    # https://en.wikipedia.org/wiki/Semidefinite_programming
    MOI.empty!(optimizer)
    @test MOI.is_empty(optimizer)

    X = MOI.add_variables(optimizer, 6)

    vov = MOI.VectorOfVariables(X)
    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(3))

    cd1 = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[1]))], [-1.0]), MOI.Zeros(1))
    cd1 = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[3]))], [-1.0]), MOI.Zeros(1))
    cd1 = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[6]))], [-1.0]), MOI.Zeros(1))

    c12_ub = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[2]))], [0.1]), MOI.Nonpositives(1))   # x <= -0.1 -> x + 0.1 <= 0
    c12_lb = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(-1.0, X[2]))], [-0.2]), MOI.Nonpositives(1)) # x >= -0.2 -> -x + -0.2 <= 0

    c23_ub = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[5]))], [-0.5]), MOI.Nonpositives(1)) # x <= 0.5 ->  x - 0.5 <= 0
    c23_lb = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(-1.0, X[5]))], [0.4]), MOI.Nonpositives(1)) # x >= 0.4 -> -x + 0.4 <= 0

    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [X[4]]), 0.0))

    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(optimizer)

    obj = MOI.get(optimizer, MOI.ObjectiveValue())
    @test obj ≈ -0.978 atol=1e-2

    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    obj = MOI.get(optimizer, MOI.ObjectiveValue())
    @test obj ≈ 0.872 atol=1e-2

end

using LinearAlgebra
LinearAlgebra.symmetric_type(::Type{MathOptInterface.VariableIndex}) = MathOptInterface.VariableIndex
LinearAlgebra.symmetric(v::MathOptInterface.VariableIndex, ::Symbol) = v
LinearAlgebra.transpose(v::MathOptInterface.VariableIndex) = v

@testset "MIMO Sizes" begin
    include("base_mimo.jl")
    include("moi_mimo.jl")
    for i in 2:5
        @testset "MIMO n = $(i)" begin
            moi_mimo(optimizer, 123, i, test = true)
        end
    end
end

@testset "RANDSDP Sizes" begin
    include("base_randsdp.jl")
    include("moi_randsdp.jl")
    for n in 10:11, m in 10:11
        @testset "RANDSDP n=$n, m=$m" begin
            moi_randsdp(optimizer, 123, n, m, test = true, atol = 1e-1)
        end
    end
end

# This problems are too large for Travis
@testset "SDPLIB Sizes" begin
    datapath = joinpath(dirname(@__FILE__), "data")
    include("base_sdplib.jl")
    include("moi_sdplib.jl")
    @testset "EQPART" begin
        # badly conditioned
        # moi_sdplib(optimizer_low_acc, joinpath(datapath, "gpp124-1.dat-s"), test = true)
        moi_sdplib(optimizer_low_acc, joinpath(datapath, "gpp124-2.dat-s"), test = true)
        # moi_sdplib(optimizer, joinpath(datapath, "gpp124-3.dat-s"), test = true)
        # moi_sdplib(optimizer, joinpath(datapath, "gpp124-4.dat-s"), test = true)
    end
    @testset "MAX CUT" begin
        moi_sdplib(optimizer_low_acc, joinpath(datapath, "mcp124-1.dat-s"), test = true)
        # moi_sdplib(optimizer, joinpath(datapath, "mcp124-2.dat-s"), test = true)
        # moi_sdplib(optimizer, joinpath(datapath, "mcp124-3.dat-s"), test = true)
        # moi_sdplib(optimizer, joinpath(datapath, "mcp124-4.dat-s"), test = true)
    end
end

@testset "Sensor Localization" begin
    include("base_sensorloc.jl")
    include("moi_sensorloc.jl")
    for n in 5:5:10
        moi_sensorloc(optimizer, 0, n, test = true)
    end
end

@testset "Full eig" begin
    MOIT.psdt0vtest(
        MOIB.full_bridge_optimizer(optimizer_full, Float64),
        MOIT.TestConfig(atol=1e-3, rtol=1e-3, duals = false)
        )
end

@testset "Print" begin
    MOIT.linear15test(optimizer_print, MOIT.TestConfig(atol=1e-3, rtol=1e-3))
end

@testset "Unsupported argument" begin
    MOI.empty!(cache)
    @test_throws ErrorException  optimizer_unsupportedarg = MOIU.CachingOptimizer(cache, ProxSDP.Optimizer(unsupportedarg = 10))
    # @test_throws ErrorException MOI.optimize!(optimizer_unsupportedarg)
end

include("test_terminationstatus.jl")


text:
using EvalMetrics
using JSON

@doc """
Function to get get statistical report based on given data.

# Examples
```jldoctest
eval_report(model, Θ, data, targets)
```
where the 'Θ' is dictionary of params,
'model' is used model, 'data' is set of samples (size(dim, N) ~ (dimension of data, data count) and
targets is a Vector{Bool} of labels to given data.
""" ->
function eval_report(model, θ::Dict{Symbol, Vector},
                        data::Matrix{Float64}, targets::Vector{Bool}, show=false)
    get_probs(model, θ::Dict{Symbol, Vector},
                data::Matrix{Float64}) = tmap(x->model(θ, Array{Float64}(x)), eachcol(data))
    scores = Vector{Float64}(get_probs(model, θ, data))
    report::Dict{String, Float64} = binary_eval_report(targets, scores)
    if show JSON.print(report, 4) end
    return report["au_roccurve"]
end

@doc """
Function to get get statistical report based on given data.

# Examples
```jldoctest
auc = eval_report(model, data, targets, show=false)
```
where the 'model' is used model,
'data' is set of samples (size(dim, N) ~ (dimension of data, data count),
targets is a Vector{Bool} of labels to given data and shaw is param to print the whole
report (default false). The function returns auc::Float64 which is area under roc-curve.
""" ->
function eval_report(model, data::Matrix{Float64}, targets::Vector{Bool}, show=false)
    get_probs(model, data::Matrix{Float64}) = tmap(x->model(Array{Float64}(x)), eachcol(data))
    scores::Vector{Float64} = Vector{Float64}(get_probs(model, data))
    report::Dict{String, Float64} = binary_eval_report(targets, scores)
    if show JSON.print(report, 4) end
    return report["au_roccurve"]
end

@doc """
    Wilcoxon sing-rank non-parametric test to check if two populations have the
    same mean.

    # Examples
    ```jldoctest
    W::Int64, z::Float64, b::Bool = wilcoxon(X::Vector{Float64}, Y::Vector{Float64})
    ```
    where 'W' is W-value, 'z' is z-score and 'b' is the test result for the First type error 0.05.
    When 'b' is true, the populations 'X' and 'Y' have got the same mean.
    """ ->
function wilcoxon(X::Vector{Float64}, Y::Vector{Float64})
    xydiff::Vector{Float64} = filter(x->x!=0, X.-Y)
    df = DataFrame(diff = xydiff, absdiff = abs.(xydiff))
    sort!(df, [order(:absdiff)])
    df[!, :sgn] = sign.(df[!, :diff])
    df[!, :Rᵢ] = 1:length(df[!, :diff])
    W::Int64 = min(sum(filter(x->x.sgn==-1, df)[!, :Rᵢ]),
                    sum(filter(x->x.sgn==1, df)[!, :Rᵢ]))
    z::Float64 = mapreduce(xy->xy[1]*xy[2], +, zip(df[!,:sgn], df[!,:Rᵢ]))
    N::Int64 = size(df)[1]
    z = z/sqrt(N*(N+1)*(2*N+1)/6)
    return W, z, abs(z) <= 1.96
end

max_stars_repo_name: kozvojtex/AnomalyDetection

max_issues_repo_name: kozvojtex/AnomalyDetection

max_forks_repo_name: kozvojtex/AnomalyDetection


text:
"""
    _get_highlights_from_CSV(csv::String) -> Vector{String}

Extract the values of the `Highlight` column from the CSV file.

# Arguments
- `csv::String`: $(CSV_ARGUMENT)

# Returns
- `Vector{String}`: the highlights

# Throws
- Exceptions from: [`get_all`](@ref)

# Example
```jldoctest; output = false
using PDFHighlights
HEADER = PDFHighlights.Internal.Constants.HEADER

_file, io = mktemp()
println(io, HEADER, '\\n', "The world didn't stop spinning,,,,1")
flush(io)
file = _file * ".csv"
mv(_file, file)

PDFHighlights.Internal.CSV._get_highlights_from_CSV(file) ==
["The world didn't stop spinning"]

# output

true
```
"""
function _get_highlights_from_CSV(csv::String)::Vector{String}
    highlights, _, _, _, _ = get_all(csv)
    return highlights
end


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

"""
    DensityInterface.test_density_interface(d, x, ref_logd_at_x; kwargs...)

Test if `d` is compatible with `DensityInterface`.

Tests that [`logdensityof(d, x)`](@ref) equals `ref_logd_at_x` and
that the behavior of [`logdensityof(d)`](@ref),
[`densityof(d, x)`](@ref) and [`densityof(d)`](@ref) is consistent.

Also tests if `logfuncdensity(logdensityof(d))` returns
a density equivalent to `d` in respect to the functions above.

The results of `logdensityof(d, x)` and `densityof(d, x)` are compared to
`ref_logd_at_x` and `exp(ref_logd_at_x)` using `isapprox`. `kwargs...` are
forwarded to `isapprox`.
"""
function test_density_interface(d, x, ref_logd_at_x; kwargs...)
    @testset "test_density_interface: $d with input $x" begin
        ref_d_at_x = exp(ref_logd_at_x)

        @test hasdensity(d) == true
        @test isapprox(logdensityof(d, x), ref_logd_at_x; kwargs...)
        log_f = logdensityof(d)
        @test isapprox(log_f(x), ref_logd_at_x; kwargs...)
        @test isapprox(densityof(d,x), ref_d_at_x; kwargs...)
        @test isapprox(densityof(d)(x), ref_d_at_x; kwargs...)

        d2 = logfuncdensity(log_f)
        @test hasdensity(d2) == true
        @test isapprox(logdensityof(d2, x), ref_logd_at_x; kwargs...)
        log_f2 = logdensityof(d2)
        @test isapprox(log_f2(x), ref_logd_at_x; kwargs...)
        @test isapprox(densityof(d2,x), ref_d_at_x; kwargs...)
        @test isapprox(densityof(d2)(x), ref_d_at_x; kwargs...)
    end
end


text:
############# Controlled Gates and 2 qubit Gates#############
function CU!(C::Circuit, U::AbstractArray{ComplexF64,2},control::Int,target::Int)
    #Implentation of CU gate using Kronnecker products.
    Uf = sparse(U)

    a = length(C.qreg)

    if (target==control)
        error("Target qubit and control qubit must be Different!")
    elseif ((target>a) || (control>a)) || ((target<=0) || (control<=0))
        error("The target and control qubit must be within the size of the register!")
    end

    M1 = spzeros(ComplexF64,2,2)
    M2 = copy(M1)
    I2 = sparse(I*(1.0+0.0*im),2,2)

    conket1 = sparse([0.0+0.0*im 0.0+0.0*im;
                      0.0+0.0*im 1.0+0.0*im])
    conket0 = sparse([1.0+0.0*im 0.0+0.0*im
                      0.0+0.0*im 0.0+0.0*im])

    for i ∈ 1:a

        if (i==1)
            if (i==control)
                M1 = conket1
            elseif (i==target)
                M1 = Uf
            else
                M1 = I2
            end
        else
            if (i==control)
                M1 = kron(conket1,M1)
            elseif (i==target)
                M1 = kron(Uf,M1)
            else
                M1 = kron(I2,M1)
            end
        end
    end

    for i ∈ 1:a

        if (i==1==control)
            M2 = conket0
        elseif (i==1!==control)
            M2 = I2
        else
            if (i==control)
                M2 = kron(conket0,M2)
            else
                M2 = kron(I2,M2)
            end
        end
    end
    return push!(C.gates,M1+M2)
end

function CU!(C::Circuit,U::AbstractSparseArray{ComplexF64,Int},control::Int,target::Int)
    #Implentation of CU gate using Kronnecker products.
    a = length(C.qreg)

    if (target==control)
        error("Target qubit and control qubit must be Different!")
    elseif ((target>a) || (control>a)) || ((target<=0) || (control<=0))
        error("The target and control qubit must be within the size of the register!")
    end

    M1 = spzeros(ComplexF64,2,2)
    M2 = copy(M1)
    I2 = sparse(I*(1.0+0.0*im),2,2)

    conket1 = sparse([0.0+0.0*im 0.0+0.0*im;
                      0.0+0.0*im 1.0+0.0*im])
    conket0 = sparse([1.0+0.0*im 0.0+0.0*im
                      0.0+0.0*im 0.0+0.0*im])

    for i ∈ 1:a

        if (i==1)
            if (i==control)
                M1 = conket1
            elseif (i==target)
                M1 = U
            else
                M1 = I2
            end
        else
            if (i==control)
                M1 = kron(conket1,M1)
            elseif (i==target)
                M1 = kron(U,M1)
            else
                M1 = kron(I2,M1)
            end
        end
    end

    for i ∈ 1:a

        if (i==1==control)
            M2 = conket0
        elseif (i==1!==control)
            M2 = I2
        else
            if (i==control)
                M2 = kron(conket0,M2)
            else
                M2 = kron(I2,M2)
            end
        end
    end
    return push!(C.gates,M1+M2)
end

function CU!(C::Circuit,λ::Float64,control::Int,target::Int)
    #Implentation of CU gate using Kronnecker products.
    a = length(C.qreg)
    U = sparse([1.0  0.0
                0.0  exp(im*λ)])

    if (target==control)
        error("Target qubit and control qubit must be Different!")
    elseif ((target>a) || (control>a)) || ((target<=0) || (control<=0))
        error("The target and control qubit must be within the size of the register!")
    end

    M1 = spzeros(ComplexF64,2,2)
    M2 = copy(M1)
    I2 = sparse(I*(1.0+0.0*im),2,2)

    conket1 = sparse([0.0+0.0*im 0.0+0.0*im;
                      0.0+0.0*im 1.0+0.0*im])
    conket0 = sparse([1.0+0.0*im 0.0+0.0*im
                      0.0+0.0*im 0.0+0.0*im])

    for i ∈ 1:a

        if (i==1)
            if (i==control)
                M1 = conket1
            elseif (i==target)
                M1 = U
            else
                M1 = I2
            end
        else
            if (i==control)
                M1 = kron(conket1,M1)
            elseif (i==target)
                M1 = kron(U,M1)
            else
                M1 = kron(I2,M1)
            end
        end
    end

    for i ∈ 1:a

        if (i==1==control)
            M2 = conket0
        elseif (i==1!==control)
            M2 = I2
        else
            if (i==control)
                M2 = kron(conket0,M2)
            else
                M2 = kron(I2,M2)
            end
        end
    end
    return push!(C.gates,M1+M2)
end

###Methods for doing a number of gates in a row
function CU!(C::Circuit,U::AbstractArray{ComplexF64,2},controls::AbstractVector{Int},targets::AbstractVector{Int})
    #convenience method for doing mutliple CU operations at a time.
    a = length(C.qreg)
    b = length(controls)
    if length(targets)!==b
        error("The set of controls and targets must be the same length or the operation is ambiguous.")
    elseif ((maximum(controls)>a) || (maximum(targets)>a)) || ((minimum(controls)<=0) || (minimum(targets)<=0))
        error("The target and control qubits must all lie within the register!")
    end

    Uf = sparse(U)
    for i in 1:b
        CU!(C,Uf,controls[i],targets[i])
    end
end

function CU!(C::Circuit,U::AbstractArray{ComplexF64,2},control::Int,targets::AbstractVector{Int})
    #convenience method for doing mutliple CU operations at a time.
    a = length(C.qreg)
    if ((control>a) || (maximum(targets)>a)) || ((control<=0) || (minimum(targets)<=0))
        error("The target and control qubits must all lie within the register!")
    end

    Uf = sparse(U)
    for i in targets
        CU!(C,Uf,control,i)
    end
end

function CU!(C::Circuit,U::AbstractArray{ComplexF64,2},controls::AbstractVector{Int},target::Int)
    #convenience method for doing mutliple CU operations at a time.
    a = length(C.qreg)
    if ((maximum(controls)>a) || (target>a)) || ((minimum(controls)<=0) || (target<=0))
        error("The target and control qubits must all lie within the register!")
    end

    Uf = sparse(U)
    for i in controls
        CU!(C,Uf,i,target)
    end
end

function CU!(C::Circuit,U::AbstractSparseArray{ComplexF64,Int},controls::AbstractVector{Int},targets::AbstractVector{Int})
    #convenience method for doing mutliple CU operations at a time.
    a = length(C.qreg)
    b = length(controls)
    if length(targets)!==b
        error("The set of controls and targets must be the same length or the operation is ambiguous.")
    elseif ((maximum(controls)>a) || (maximum(targets)>a)) || ((minimum(controls)<=0) || (minimum(targets)<=0))
        error("The target and control qubits must all lie within the register!")
    end

    for i in 1:b
        CU!(C,U,controls[i],targets[i])
    end
end

function CU!(C::Circuit,U::AbstractSparseArray{ComplexF64,Int},control::Int,targets::AbstractVector{Int})
    #convenience method for doing mutliple CU operations at a time.
    a = length(C.qreg)
    if ((control>a) || (maximum(targets)>a)) || ((control<=0) || (minimum(targets)<=0))
        error("The target and control qubits must all lie within the register!")
    end

    for i in targets
        CU!(C,U,control,i)
    end
end

function CU!(C::Circuit,U::AbstractSparseArray{ComplexF64,Int},controls::AbstractVector{Int},target::Int)
    #convenience method for doing mutliple CU operations at a time.
    a = length(C.qreg)
    if ((maximum(controls)>a) || (target>a)) || ((minimum(controls)<=0) || (target<=0))
        error("The target and control qubits must all lie within the register!")
    end

    for i in controls
        CU!(C,U,i,target)
    end
end

function CU!(C::Circuit,λ::Float64,controls::AbstractVector{Int},targets::AbstractVector{Int})
    #convenience method for doing mutliple CU operations at a time.
    a = length(C.qreg)
    b = length(controls)
    if length(targets)!==b
        error("The set of controls and targets must be the same length or the operation is ambiguous.")
    elseif ((maximum(controls)>a) || (maximum(targets)>a)) || ((minimum(controls)<=0) || (minimum(targets)<=0))
        error("The target and control qubits must all lie within the register!")
    end

    for i in 1:b
        CU!(C,λ,controls[i],targets[i])
    end
end

function CU!(C::Circuit,λ::Float64,control::Int,targets::AbstractVector{Int})
    #convenience method for doing mutliple CU operations at a time.
    a = length(C.qreg)
    if ((control>a) || (maximum(targets)>a)) || ((control<=0) || (minimum(targets)<=0))
        error("The target and control qubits must all lie within the register!")
    end

    for i in targets
        CU!(C,λ,control,i)
    end
end

function CU!(C::Circuit,λ::Float64,controls::AbstractVector{Int},target::Int)
    #convenience method for doing mutliple CU operations at a time.
    a = length(C.qreg)
    if ((maximum(controls)>a) || (target>a)) || ((minimum(controls)<=0) || (target<=0))
        error("The target and control qubits must all lie within the register!")
    end

    for i in controls
        CU!(C,λ,i,target)
    end
end

function CX!(C::Circuit,control::Int,target::Int)
    #Implentation of CNOT gate using Kronnecker products.
    X = sparse([0.0+0.0*im 1.0+0.0*im;
                1.0+0.0*im 0.0+0.0*im])

    CU!(C,X,control,target)
end

function CX!(C::Circuit,controls::AbstractVector{Int},targets::AbstractVector{Int})
    #convenience method for doing mutliple CX operations at a time.
    X = sparse([0.0+0.0*im 1.0+0.0*im;
                1.0+0.0*im 0.0+0.0*im])

    CU!(C,X,controls,targets)
end

function CX!(C::Circuit,control::Int,targets::AbstractVector{Int})
    #convenience method for doing mutliple CX operations at a time.
    X = sparse([0.0+0.0*im 1.0+0.0*im;
                1.0+0.0*im 0.0+0.0*im])

    CU!(C,X,control,targets)
end

function CX!(C::Circuit,controls::AbstractVector{Int},target::Int)
    #convenience method for doing mutliple CX operations at a time.
    X = sparse([0.0+0.0*im 1.0+0.0*im;
                1.0+0.0*im 0.0+0.0*im])

    CU!(C,X,controls,target)
end

function CZ!(C::Circuit,control::Int,target::Int)
    #Implentation of CZ gate using Kronnecker products.
    Z = sparse([1.0+0.0*im  0.0+0.0*im;
                0.0+0.0*im -1.0+0.0*im])

    CU!(C,Z,control,target)
end

function CZ!(C::Circuit,controls::AbstractVector{Int},targets::AbstractVector{Int})
    #convenience method for doing mutliple CZ operations at a time.
    Z = sparse([1.0+0.0*im  0.0+0.0*im;
                0.0+0.0*im -1.0+0.0*im])

    CU!(C,Z,controls,targets)
end

function CZ!(C::Circuit,control::Int,targets::AbstractVector{Int})
    #convenience method for doing mutliple CZ operations at a time.
    Z = sparse([1.0+0.0*im  0.0+0.0*im;
                0.0+0.0*im -1.0+0.0*im])

    CU!(C,Z,control,targets)
end

function CZ!(C::Circuit,controls::AbstractVector{Int},target::Int)
    #convenience method for doing mutliple CZ operations at a time.
    Z = sparse([1.0+0.0*im  0.0+0.0*im;
                0.0+0.0*im -1.0+0.0*im])

    CU!(C,Z,controls,target)
end

function Swap!(C::Circuit,b1::Int,b2::Int)
    #quick and dirty implmentation of Swap Gate
    a = length(C.qreg)
    if b1==b2
        error("This ins't really a swap, you just tried to swap a qubit with itself!")
    elseif ((b1>a) || (b2>a)) || ((b1<=0) || (b2<=0))
        error("Both bits must be within the size of the register!")
    end
    CX!(C,b2,b1)
    CX!(C,b1,b2)
    CX!(C,b2,b1)
end

function CCX!(C::Circuit,c1::Int,c2::Int,target::Int)
    #Implementation of the Toffoli gate using Kronnecker products.
    X = sparse([0.0+0.0*im 1.0+0.0*im;
                1.0+0.0*im 0.0+0.0*im])

    a = length(C.qreg)

    if (target==c1) || (target==c2)
        error("Target qubit and control qubits must be Different!")
    elseif (c1==c2)
        error("The control qubits must be different from one another!")
    elseif ((target>a) || (c1>a) || (c2>a)) || ((target<=0) || (c1<=0) ||(c2<=0) )
        error("The target and control qubits must be within the size of the register!")
    end

    M1 = spzeros(ComplexF64,2,2)
    M2 = copy(M1)
    M3 = copy(M1)
    M4 = copy(M1)
    I2 = sparse(I*(1.0+0.0*im),2,2)

    conket1 = sparse([0.0+0.0*im 0.0+0.0*im;
                      0.0+0.0*im 1.0+0.0*im])
    conket0 = sparse([1.0+0.0*im 0.0+0.0*im
                      0.0+0.0*im 0.0+0.0*im])

    for i ∈ 1:a

        if (i==1)
            if (i==c1)||(i==c2)
                M1 = conket1
            elseif (i==target)
                M1 = X
            else
                M1 = I2
            end
        else
            if (i==c1)||(i==c2)
                M1 = kron(conket1,M1)
            elseif (i==target)
                M1 = kron(X,M1)
            else
                M1 = kron(I2,M1)
            end
        end
    end

    for i ∈ 1:a

        if (i==1)
            if (i==c1)||(i==c2)
                M2 = conket0
            else
                M2 = I2
            end
        else
            if (i==c1)||(i==c2)
                M2 = kron(conket0,M2)
            else
                M2 = kron(I2,M2)
            end
        end
    end

    for i ∈ 1:a

        if (i==1)
            if (i==c1)
                M3 = conket0
            elseif (i==c2)
                M3 = conket1
            else
                M3 = I2
            end
        else
            if (i==c1)
                M3 = kron(conket0,M3)
            elseif (i==c2)
                M3 = kron(conket1,M3)
            else
                M3 = kron(I2,M3)
            end
        end
    end

    for i ∈ 1:a

        if (i==1)
            if (i==c1)
                M4 = conket1
            elseif (i==c2)
                M4 = conket0
            else
                M4 = I2
            end
        else
            if (i==c1)
                M4 = kron(conket1,M4)
            elseif (i==c2)
                M4 = kron(conket0,M4)
            else
                M4 = kron(I2,M4)
            end
        end
    end
    return push!(C.gates,M1+M2+M3+M4)

end

function CSwap!(C::Circuit,control,bit1,bit2)
    #Controlled swap for order finding algorithm
    CX!(C,bit2,bit1)
    CCX!(C,control,bit1,bit2)
    CX!(C,bit2,bit1)
end


text:
function param_conv2d(ni,no;k=3,stride=2)
    return Dict("filter" => (k,k), "in" => ni, "out" => no, "padding" => k÷2, "stride" => stride)
end

function make_conv(d::Dict)
    return Conv(d["filter"], d["in"] => d["out"], relu, stride = d["stride"], pad = d["padding"])
end

function make_list_i_o(nfs;start=1)
    return [i==0 ? (start,nfs[1]) : (nfs[i],nfs[i+1]) for i in 0:length(nfs)-1]
end

function make_list_layers(nfs)
    conv_layers = Tuple(i==1 ? make_conv(param_conv2d(i,o,k=5)) : make_conv(param_conv2d(i,o)) for (i,o) in make_list_i_o(nfs))
    return tuple(conv_layers..., AdaptiveMeanPool((1,1)), flatten, Dense(last(nfs), 10))
end

function get_model(nfs)
    layers = make_list_layers(nfs)
    return Chain(layers...)
end


text:
using QuasiCopula, Random, Statistics, Test, LinearAlgebra, StatsFuns

### Poisson ###
@testset "Poisson(θ) * (c0 + c1 * x + c2 * x^2);" begin
Random.seed!(1234)
dist = Poisson(5)
# test if the proper c0, c1, c2 constants are stored.
γ = 1.0
d_pois = pdf_constants(γ, dist)

@test d_pois.c0 == 1 + (mean(d_pois.d)^2 * inv(var(d_pois.d)) * 0.5 * γ)
@test d_pois.c1 == 0.5 * γ * (-2 * mean(d_pois.d) * inv(var(d_pois.d)))
@test d_pois.c2 == 0.5 * γ * (inv(var(d_pois.d)))

pmf_copula!(d_pois)
reordered_k, reordered_pmf = reorder_pmf(d_pois.pmf_vec, d_pois.μ)
sum(reordered_pmf) ≈ 1

###
Random.seed!(1234)
nsample = 10_000
@info "sample $nsample points for the $dist distribution"
s = Vector{Int64}(undef, nsample)
rand!(d_pois, s) # compile
@time rand!(d_pois, s) # get time
println("sample mean = $(Statistics.mean(s)); theoretical mean = $(mean(d_pois))")
println("sample var = $(Statistics.var(s)); theoretical var = $(var(d_pois))")
end
# n = 10,0000 ::: time = 0.019234 seconds
# sample mean = 5.315; theoretical mean = 5.333333333333333
# sample var = 8.48922392239224; theoretical var = 8.555555555555557

# ### Binomial ###
@testset "Binomial(n, p) * (c0 + c1 * x + c2 * x^2);" begin
Random.seed!(123)
n = 30
p = 0.1
dist = Binomial(n, p)

# test if the proper c0, c1, c2 constants are stored.
γ = 1.0
d_binomial = pdf_constants(γ, dist)

@test d_binomial.c ≈ inv(1 + 0.5 * γ)
@test d_binomial.c0 == 1 + (mean(d_binomial.d)^2 * inv(var(d_binomial.d)) * 0.5 * γ)
@test d_binomial.c1 == 0.5 * γ * (-2 * mean(d_binomial.d) * inv(var(d_binomial.d)))
@test d_binomial.c2 == 0.5 * γ * (inv(var(d_binomial.d)))

pmf_copula!(d_binomial)
reordered_k, reordered_pmf = reorder_pmf(d_binomial.pmf_vec, d_binomial.μ)
@test sum(reordered_pmf) ≈ 1

###
Random.seed!(1234)
nsample = 10_000 #
@info "sample $nsample points for the $dist distribution using the Bisection method."
s = Vector{Int64}(undef, nsample)
rand!(d_binomial, s) # compile
@time rand!(d_binomial, s) # get time
println("sample mean = $(Statistics.mean(s)); theoretical mean = $(QuasiCopula.mean(d_binomial))")
println("sample var = $(Statistics.var(s)); theoretical var = $(QuasiCopula.var(d_binomial))")
end

# n = 10,0000 ::: time = 0.038302 seconds (40.00 k allocations: 12.817 MiB, 23.83% gc time)
# sample mean = 3.129; theoretical mean = 3.122671204329908
# sample var = 3.561715171517152; theoretical var = 3.583518347344844

# ### Geometric p = 0.2 ###
@testset "Geometric(p = 0.2) * (c0 + c1 * x + c2 * x^2);" begin
Random.seed!(123)
p = 0.2
dist = Geometric(p)

# test if the proper c0, c1, c2 constants are stored.
γ = 1.0
d_geometric = pdf_constants(γ, dist)

@test d_geometric.c ≈ inv(1 + 0.5 * γ)
@test d_geometric.c0 == 1 + (mean(d_geometric.d)^2 * inv(var(d_geometric.d)) * 0.5 * γ)
@test d_geometric.c1 == 0.5 * γ * (-2 * mean(d_geometric.d) * inv(var(d_geometric.d)))
@test d_geometric.c2 == 0.5 * γ * (inv(var(d_geometric.d)))

pmf_copula!(d_geometric)
reordered_k, reordered_pmf = reorder_pmf(d_geometric.pmf_vec, d_geometric.μ)
@test sum(reordered_pmf) ≈ 1

###
Random.seed!(1234)
nsample = 10_000 #
@info "sample $nsample points for the $dist distribution."
s = Vector{Int64}(undef, nsample)
rand!(d_geometric, s) # compile
@time rand!(d_geometric, s) # get time
println("sample mean = $(Statistics.mean(s)); theoretical mean = $(QuasiCopula.mean(d_geometric))")
println("sample var = $(Statistics.var(s)); theoretical var = $(QuasiCopula.var(d_geometric))")
end

# n = 10,0000 ::: time = 0.032615 seconds (40.00 k allocations: 34.180 MiB, 21.35% gc time)

# ### Geometric p = 0.5###
@testset "Geometric(p = 0.5) * (c0 + c1 * x + c2 * x^2);" begin
Random.seed!(123)

p = 0.5
dist = Geometric(p)

# test if the proper c0, c1, c2 constants are stored.
γ = 1.0
d_geometric = pdf_constants(γ, dist)

@test d_geometric.c ≈ inv(1 + 0.5 * γ)
@test d_geometric.c0 == 1 + (mean(d_geometric.d)^2 * inv(var(d_geometric.d)) * 0.5 * γ)
@test d_geometric.c1 == 0.5 * γ * (-2 * mean(d_geometric.d) * inv(var(d_geometric.d)))
@test d_geometric.c2 == 0.5 * γ * (inv(var(d_geometric.d)))

pmf_copula!(d_geometric)
reordered_k, reordered_pmf = reorder_pmf(d_geometric.pmf_vec, d_geometric.μ)
@test sum(reordered_pmf) ≈ 1

###
Random.seed!(1234)
nsample = 10_000 #
@info "sample $nsample points for the $dist distribution."
s = Vector{Int64}(undef, nsample)
rand!(d_geometric, s) # compile
@time rand!(d_geometric, s) # get time
println("sample mean = $(Statistics.mean(s)); theoretical mean = $(QuasiCopula.mean(d_geometric))")
println("sample var = $(Statistics.var(s)); theoretical var = $(QuasiCopula.var(d_geometric))")
end

# n = 10,0000 ::: time = 0.025624 seconds (40.00 k allocations: 34.180 MiB, 33.37% gc time)

# ### Negative Binomial r = 5, p = 0.5 ###
@testset "NegativeBinomial(r = 5, p = 0.5) * (c0 + c1 * x + c2 * x^2);" begin
Random.seed!(123)
r, p = 5, 0.5
dist = NegativeBinomial(r, p)

# test if the proper c0, c1, c2 constants are stored.
γ = 1.0
d_nb = pdf_constants(γ, dist)

@test d_nb.c ≈ inv(1 + 0.5 * γ)
@test d_nb.c0 == 1  + (mean(d_nb.d)^2 * inv(var(d_nb.d)) * 0.5 * γ)
@test d_nb.c1 == 0.5 * γ * (-2 * mean(d_nb.d) * inv(var(d_nb.d)))
@test d_nb.c2 == 0.5 * γ * (inv(var(d_nb.d)))

pmf_copula!(d_nb)
reordered_k, reordered_pmf = reorder_pmf(d_nb.pmf_vec, d_nb.μ)
@test sum(reordered_pmf) ≈ 1

###
Random.seed!(1234)
nsample = 10_000 #
@info "sample $nsample points for the $dist distribution."
s = Vector{Int64}(undef, nsample)
rand!(d_nb, s) # compile
@time rand!(d_nb, s) # get time
println("sample mean = $(Statistics.mean(s)); theoretical mean = $(QuasiCopula.mean(d_nb))")
println("sample var = $(Statistics.var(s)); theoretical var = $(QuasiCopula.var(d_nb))")
end

# n = 10,0000 ::: time = 0.142361 seconds (40.00 k allocations: 18.921 MiB, 57.42% gc time)
### Bernoulli ###
@testset "Bernoulli(θ) * (c0 + c1 * x + c2 * x^2);" begin
Random.seed!(1234)
dist = Bernoulli(0.8)
# test if the proper c0, c1, c2 constants are stored.
γ = 1.0
d_bernoulli = pdf_constants(γ, dist)

@test d_bernoulli.c0 == 1 + (mean(d_bernoulli.d)^2 * inv(var(d_bernoulli.d)) * 0.5 * γ)
@test d_bernoulli.c1 == 0.5 * γ * (-2 * mean(d_bernoulli.d) * inv(var(d_bernoulli.d)))
@test d_bernoulli.c2 == 0.5 * γ * (inv(var(d_bernoulli.d)))

pmf_copula!(d_bernoulli)
reordered_k, reordered_pmf = reorder_pmf(d_bernoulli.pmf_vec, d_bernoulli.μ)
sum(reordered_pmf) ≈ 1

###
Random.seed!(1234)
nsample = 10_000
@info "sample $nsample points for the $dist distribution"
s = Vector{Int64}(undef, nsample)
rand!(d_bernoulli, s) # compile
@time rand!(d_bernoulli, s) # get time
println("sample mean = $(Statistics.mean(s)); theoretical mean = $(mean(d_bernoulli))")
println("sample var = $(Statistics.var(s)); theoretical var = $(var(d_bernoulli))")
end
# n = 10,0000 ::: time = 0.019234 seconds
# sample mean = 5.315; theoretical mean = 5.333333333333333
# sample var = 8.48922392239224; theoretical var = 8.555555555555557


text:
using Test

using Dash, DashHtmlComponents, DashCoreComponents, DashBootstrapComponents
using PlotlyJS
using DataFrames

using Classic
#cc = Classic

function run_astar()
    maze = cc.new_maze(10, 10, (1,1), (10,10), 0.2)

    # create functions for bfs()
    is_goal(pt)     = cc.is_goal(maze.goal, pt)
    next_points(pt) = cc.next_points(maze.grid, pt)
    heuristic(pt)   = cc.manhattan_distance(pt, maze.goal)

    node = cc.a_star((1,1), is_goal, next_points, heuristic)

    if isnothing(node)
        println("Path not found")
        display(maze.grid) # print the unsovled maze
    else
        cc.mark_path!(maze.grid, node, start=maze.start, goal=maze.goal)
        display(maze.grid) # print the maze and the path
    end
end

function run_bfs()
    maze = cc.new_maze(10, 10, (1,1), (10,10), 0.2)

    # create functions for bfs()
    is_goal(pt)     = cc.is_goal(maze.goal, pt)
    next_points(pt) = cc.next_points(maze.grid, pt)

    node = cc.bfs((1,1), is_goal, next_points)

    if isnothing(node)
        println("Path not found")
        display(maze.grid) # print the unsovled maze
    else
        cc.mark_path!(maze.grid, node, start=maze.start, goal=maze.goal)
        display(maze.grid) # print the maze and the path
    end
end

function run_cell()
    x = cc.EMPTY
    y = cc.START

    println("x = $(string(x)), y = $(string(y))")
end

"""
    run_dfs()

Use depth-first search to solve a maze.
"""
function run_dfs()
    maze = cc.new_maze() #new_maze(10,10, [1,1], [10,10], 0.2)
    #@show maze
    #display(maze.grid)

    #node = ch2.Node()
    #locations = ch2.successors(maze.grid, (2,2))
    #@show locations

    # pass this to dfs(..., successors)
    successors(here) = cc.successors(maze.grid, here)
    # pass this to dfs(..., goal_test, ...)
    is_goal(here) = cc.is_goal(maze.goal, here)

    node = cc.dfs((1,1), is_goal, successors)
    if isnothing(node)
        println("Path not found")
        display(maze.grid)
    else
        cc.mark_path!(maze.grid, node, start=maze.start, goal=maze.goal)
        #@show maze
        display(maze.grid)
    end
    #@show node
end

function run_bfs_gui()
    maze = new_maze(10, 10, (1,1), (10,10), 0.2)

    # create functions for bfs()
    is_goal(pt)     = Classic.is_goal(maze.goal, pt)
    next_points(pt) = Classic.next_points(maze.grid, pt)

    node = bfs((1,1), is_goal, next_points)

    if isnothing(node)
        println("Path not found")
        display(maze.grid) # print the unsovled maze
    else
        mark_path!(maze.grid, node, start=maze.start, goal=maze.goal)
        display(maze.grid) # print the maze and the path
    end

    app = dash(external_stylesheets=[dbc_themes.SPACELAB])

    navbar = dbc_navbarsimple([
        dbc_dropdownmenu([
            dbc_dropdownmenuitem("BFS", href="#bfs", external_link=true),
            # dbc_dropdownmenuitem("", href="", external_link=true),
        ],
        in_navbar=true, label="Section", caret=true, direction="left"),
    ],
    sticky="top", expand=true, brand="WDVG", brand_href="https://www.wdvgco.com",
    )
    content =[
        dbc_container([
            dbc_badge("Line: $(@__LINE__)", color="info", className="ml-1"),
            dcc_graph(
                figure = Plot(bfs_heatmap(maze)...),
                config = Dict(),
            )
        ], className="p-3 my-2 border rounded"),
    ]
    pushfirst!(content, navbar)
    app.layout = dbc_container(content)

    run_server(
        app,
        "0.0.0.0",
        8055,
        debug=true, # enables hot reload and more
    )
end

function bfs_heatmap(maze)
    z = convert(Matrix{Int64}, maze.grid)
    # z =  [[1, 20, 30], [20, 1, 60], [30, 60, 1]]
    # z = hcat(z...)

    # Normalize color scale
    # Details at
    # https://community.plotly.com/t/colors-for-discrete-ranges-in-heatmaps/7780/10
    colorscale = [
        [0, "rgb(250,250,250)"],
        [1/4, "rgb(0,0,0)"],
        [2/4, "rgb(250, 75, 0)"],
        [3/4, "rgb(250, 75, 0)"],
        [4/4, "rgb(0, 250, 100)"],
    ]

    traces = Vector{AbstractTrace}([
        heatmap(
            z = z,
            zmin = 0,
            zmax = 4,
            colorscale = colorscale,
            showscale = false,
        )
    ])
    layout = Layout(
        title = "Solve maze with BFS",
        xaxis = Dict(
            :ticks => "",
            :showticklabels => false,
        ),
        yaxis = Dict(
            :ticks => "",
            :showticklabels => false,
        ),
        annotations = [],
    )
    annotations = layout["annotations"] # shortcut
    annotation = Dict(
        :x => 0,
        :y => 0,
        :text => "", #"S",
        #:showtext => false,
        :showarrow => false,
    )
    push!(annotations, annotation)
    annotation = Dict(
        :x => 9,
        :y => 9,
        :text => "", #"G",
        :showarrow => false,
    )
    push!(annotations, annotation)

    return traces, layout
end

#run_astar()
#run_cell()
#run_dfs()
#run_bfs()

run_bfs_gui()

nothing

#A = test_manhattan_distance()
#A # Test.Pass or not


text:
using ArgParse
using DataFrames

s=ArgParseSettings()
@add_arg_table s begin
    "--part2"
        help="Do part 2"
        action=:store_true
    "--test"
        help="Run test data"
        action=:store_true
end
parsed_args = parse_args(ARGS, s)

function parseend(endstr)
    coords = split(endstr, ",")
    (x=parse(UInt32, coords[1]), y=parse(UInt32, coords[2]))
end

function parseline(line)
    ends = split(line, " -> ")
    start = parseend(ends[1])
    stop = parseend(ends[2])
    (start=start, stop=stop)
end

function isortho(line)
    (line.start.x == line.stop.x || line.start.y == line.stop.y)
end

function getdelta(start, stop)
    if start > stop
        delta = -1
    elseif start < stop
        delta = 1
    else
        delta = 0
    end
    delta
end

function addpoints!(df, l)
    dx = getdelta(l.start.x, l.stop.x)
    dy = getdelta(l.start.y, l.stop.y)
    point = l.start
    while true
        push!(df, point)
        if point == l.stop
            break
        end
        point = (x=point.x+dx, y=point.y+dy)
    end
end

test_suffix = parsed_args["test"] ? "_test" : ""
endpoints = [parseline(line) for line in eachline("input/day_5$test_suffix.txt")]
if parsed_args["part2"]
    lines = endpoints
else
    lines = filter(isortho, endpoints)
end

df = DataFrame(x=UInt32[], y=UInt32[])
foreach(l -> addpoints!(df, l), lines)
counts = combine(groupby(df, [:x, :y]), nrow => :count)
answer = nrow(filter(:count => >=(2), counts))

println("The answer is: $answer")