text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

VOID = [[]], [[1]]

# example 1

V = [[0,0] [1,0] [2,0] [0,1] [1,1] [2,1] [0,2] [1,2] [2,2]];
FV = [[1,2,4],[2,3,5],[3,5,6],[4,5,7],[5,7,8],[6,8,9]];
pattern = repeat([1,2,-3],outer=4);
model = (V,FV);
W,FW = Lar.extrudeSimplicial(model, pattern);
GL.VIEW([ GL.GLGrid(W,FW, GL.COLORS[7], 0.5) ]);

# example 2

model = Lar.extrudeSimplicial( VOID, ones(10) )
GL.VIEW([ GL.GLGrid(model..., GL.COLORS[7], 0.5) ]); #TODO: check

model = Lar.extrudeSimplicial( model, ones(10) )
GL.VIEW([ GL.GLLines(model..., GL.COLORS[7]),
        GL.GLFrame ]);

model = Lar.extrudeSimplicial( model, ones(10) )
GL.VIEW([ GL.GLGrid(model..., GL.COLORS[7], 0.5) ]);

# example 3

model = Lar.extrudeSimplicial( VOID, repeat([1,-1],outer=10) )
GL.VIEW([ GL.GLGrid(model..., GL.COLORS[7], 0.5) ]); #TODO: check

model = Lar.extrudeSimplicial( model, repeat([1,-1],outer=10) )
GL.VIEW([ GL.GLGrid(model..., GL.COLORS[7],0.7),
        GL.GLFrame ]);

# example 4

grid_2d = Lar.simplexGrid([3,3])
GL.VIEW([ GL.GLGrid(grid_2d..., GL.COLORS[7],0.7),
        GL.GLLines(grid_2d...),    # TODO:  check lines
        GL.GLFrame ]);

grid_3d = Lar.simplexGrid([2,3,4])
GL.VIEW([ GL.GLGrid(grid_3d..., GL.COLORS[7]),
        GL.GLFrame ]);  # TODO  check


text:
# """
# Function that returns Utility for given prices and coordinates. Requires precompilation of Utility Cube. GPU version. Requires precompilation of utils cube in VRAM. Accesses values directly. Requires as many GPUs as there are statuses.
# # Arguments
# - `ind_vec::NTuple{M,<:Integer}`: coordinates of Utility.
# - `mdl::Model`: mutable struct Model.
# """

# function cube_GPU( ind_vec::Vector{<:Integer}, mdl::Model, additional_args=nothing)

#     grid = mdl.GridDef
#     index, status = ndim2lin(ind_vec, grid)
#     set_device(status-1)
#     return(mdl.Utils.qb[ status ][ index ])

# end

# """
# Function building cube. GPU version. Requires precompilation of utils cube in VRAM. Accesses values directly.Requires as many GPUs as there are statuses.
# # Arguments
# - `mdl::Model`: struct Model.
# """
# function cube_GPU_qb(mdl::Model) where M

#     n_dim = tuple( length.(mdl.GridDef.Grids)[1:(end-1)]... )
#     cube = Vector{AFArray}()

#     for status = mdl.GridDef.Grids[end]

#         set_device(status-1)
#         append!( cube, [AFArray( repeat([-Inf::Float64], n_dim... ) )] )

#     end

#     for status = mdl.GridDef.Grids[end]

#         set_device( status - 1)

#         Threads.@threads for idx = CartesianIndices(n_dim)

#             ind_vec = CartesianIndex( idx, status )

#             cube[status][idx] =  mdl.HAutility( mdl.HAconstraint( Tuple(ind_vec) , mdl )... )

#         end

#     end

#     return(cube)

# end

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(sdata::AbstractData, vars=nothing)
  # retrieve underlying domain
  sdomain = domain(sdata)

  # valid variables
  validvars = name.(variables(sdata))

  # plot all variables by default
  isnothing(vars) && (vars = validvars)
  @assert vars ⊆ validvars "invalid variable name"

  # shared plot specs
  layout --> length(vars)

  for (i, var) in enumerate(vars)
    # retrieve valid values
    vals = sdata[var]
    vals = vals isa CategoricalArray ? levelcode.(vals) : vals
    vals = replace(vals, missing => NaN)

    @series begin
      subplot := i
      title --> string(var)
      legend --> false
      sdomain, vals
    end
  end
end


text:
using Distributed
@everywhere using SetIntersectionProjection
using MAT
ENV["MPLBACKEND"]="qt5agg"
using PyPlot
using LinearAlgebra

@everywhere mutable struct compgrid
  d :: Tuple
  n :: Tuple
end

#test 2D
width=[50 100 200 400 800 1600]

test = Vector{Any}(undef,length(width))
N    = Vector{Any}(undef,length(width))

#PARSDMM options:
options              = PARSDMM_options()
options.FL           = Float32
options.evol_rel_tol = 10*eps(options.FL)

set_zero_subnormals(true)
BLAS.set_num_threads(2)

#select working precision
if options.FL==Float64
  TF = Float64
  TI = Int64
elseif options.FL==Float32
  TF = Float32
  TI = Int32
end

#define constraints
constraint = Vector{SetIntersectionProjection.set_definitions}() #initialize

#bounds:
m_min     = 1450.0
m_max     = 4750.0
set_type  = "bounds"
TD_OP     = "identity"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#vertical monotonicity
m_min     = 0.0
m_max     = 1e6
set_type  = "bounds"
TD_OP     = "D_z"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#some horizontal smoothness
m_min     = -1.0
m_max     = 1.0
set_type  = "bounds"
TD_OP     = "D_x"
app_mode  = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

#initialize arrays to save timings
log_T_serial   = Vector{Any}(undef,length(width))
T_tot_serial   = Vector{Any}(undef,length(width))
log_T_parallel = Vector{Any}(undef,length(width))
T_tot_parallel = Vector{Any}(undef,length(width))

log_T_serial_multilevel   = Vector{Any}(undef,length(width))
T_tot_serial_multilevel   = Vector{Any}(undef,length(width))
log_T_parallel_multilevel = Vector{Any}(undef,length(width))
T_tot_parallel_multilevel = Vector{Any}(undef,length(width))

for i=1:length(width)
  print(i)

  file = matopen(joinpath(dirname(pathof(SetIntersectionProjection)), "../examples/Data/compass_velocity.mat"))
  m=read(file, "Data")
  close(file)
  m = m[1:341,1:width[i]]
  m = permutedims(m,[2,1])

  comp_grid = compgrid((TF(25), TF(25)),(size(m,1), size(m,2)))
  m    = convert(Vector{TF},vec(m))
  N[i] = prod(size(m))

  #serial
  @everywhere GC.gc()
  println("")
  println("serial")
  options.parallel = false
  (P_sub,TD_OP,TD_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  (TD_OP,AtA,l,y)       = PARSDMM_precompute_distribute(TD_OP,TD_Prop,comp_grid,options)
  (x,log_PARSDMM)       = PARSDMM(m,AtA,TD_OP,TD_Prop,P_sub,comp_grid,options);
  val, t, bytes, gctime, memallocs = @timed (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,TD_Prop,P_sub,comp_grid,options);
  println(t)
  log_T_serial[i]=log_PARSDMM;
  T_tot_serial[i]=t;

  #parallel
  @everywhere GC.gc()
  println("")
  println("parallel")
  options.parallel = true
  (P_sub,TD_OP,TD_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  (TD_OP,AtA,l,y)       = PARSDMM_precompute_distribute(TD_OP,TD_Prop,comp_grid,options)
  (x,log_PARSDMM)       = PARSDMM(m,AtA,TD_OP,TD_Prop,P_sub,comp_grid,options);
  val, t, bytes, gctime, memallocs = @timed (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,TD_Prop,P_sub,comp_grid,options);
  println(t)
  log_T_parallel[i] = log_PARSDMM
  T_tot_parallel[i] = t

  #serial multilevel
  @everywhere GC.gc()
  println("")
  println("serial multilevel")
  options.parallel  = false
  n_levels          = 2
  coarsening_factor = 3
  (TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels)=setup_multi_level_PARSDMM(m,n_levels,coarsening_factor,comp_grid,constraint,options)
  (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
  val, t, bytes, gctime, memallocs = @timed (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
  println(t)
  log_T_serial_multilevel[i] = log_PARSDMM
  T_tot_serial_multilevel[i] = t

  #parallel multilevel
  @everywhere GC.gc()
  println("")
  println("parallel multilevel")
  options.parallel  = true
  n_levels          = 2
  coarsening_factor = 3
  (TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels)=setup_multi_level_PARSDMM(m,n_levels,coarsening_factor,comp_grid,constraint,options)
  (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
  val, t, bytes, gctime, memallocs = @timed (x,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
  println(t)
  log_T_parallel_multilevel[i] = log_PARSDMM
  T_tot_parallel_multilevel[i] = t

  if i==1
    figure();imshow(permutedims(reshape(x,(comp_grid.n[1],comp_grid.n[2])),[2,1]));colorbar
    savefig("projection_intersection_timings_2D_fig1.eps",bbox_inches="tight",dpi=300)
  elseif i==3
    figure();imshow(permutedims(reshape(x,(comp_grid.n[1],comp_grid.n[2])),[2,1]));colorbar
    savefig("projection_intersection_timings_2D_fig2.eps",bbox_inches="tight",dpi=300)
  elseif i==6
    figure();imshow(permutedims(reshape(x,(comp_grid.n[1],comp_grid.n[2])),[2,1]));colorbar
    savefig("projection_intersection_timings_2D_fig3.eps",bbox_inches="tight",dpi=300)
  end

end

#plot results
fig, ax = subplots()
ax[:loglog](N, T_tot_serial, marker="o", markersize=10, label="serial",linewidth=5)
ax[:loglog](N, T_tot_parallel, marker="x", markersize=10,label="parallel",linewidth=5)
ax[:loglog](N, T_tot_serial_multilevel, marker="v", markersize=10, label="serial multilevel",linewidth=5)
ax[:loglog](N, T_tot_parallel_multilevel, marker="^", markersize=10, label="parallel multilevel",linewidth=5)
ax[:legend]()
title(string("time 2D vs grid size, JuliaThreads=",Threads.nthreads(),", BLAS threads=",ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ())), fontsize=12)
xlabel("N gridpoints", fontsize=15)
ylabel("time [seconds]", fontsize=15)
savefig("projection_intersection_timings2D_1.eps",bbox_inches="tight",dpi=1200)

#######################################################################################
#######################################################################################
#
# #now with a different set of constraints:
# # transform-domain rank and bounds
# constraint=Dict()
#
# #bound constraints
# constraint["use_bounds"]=true
# constraint["m_min"]=1450
# constraint["m_max"]=4750
#
# #nuclear norm constraint on vertical derivative of the image
# constraint["use_TD_nuclear_1"]=true
# constraint["TD_nuclear_operator_1"]="D_z"
# #the max nuclear norm is adjusted to each model size, see below
#
# log_T_serial=Vector{Any}(length(width))
# T_tot_serial=Vector{Any}(length(width))
# log_T_parallel=Vector{Any}(length(width))
# T_tot_parallel=Vector{Any}(length(width))
# log_T_serial_multilevel=Vector{Any}(length(width))
# T_tot_serial_multilevel=Vector{Any}(length(width))
# log_T_parallel_multilevel=Vector{Any}(length(width))
# T_tot_parallel_multilevel=Vector{Any}(length(width))
#
# for i=1:length(width)
#   print(i)
#
#   file = matopen("compass_velocity.mat")
#   m=read(file, "Data")
#   close(file)
#   m=m[1:341,1:width[i]];
#   m=m';
#   comp_grid = compgrid((TF(25), TF(25)),(size(m,1), size(m,2)))
#
#   (Dz,dummy1,dummy2,TD_n,dummy4)=get_TD_operator(comp_grid,constraint["TD_nuclear_operator_1"],TF)
#   constraint["TD_nuclear_norm_1"]=0.25f0*norm(svdvals(reshape(Dz*vec(m),TD_n)),1)
#
#   m=convert(Vector{TF},vec(m));
#   N[i]=prod(size(m));
#
#   #serial
#   println("")
#   println("serial")
#   options.parallel=false
#   (P_sub,TD_OP,TD_Prop) = setup_constraints(constraint,comp_grid,options.FL)
#   (TD_OP,AtA,l,y) = PARSDMM_precompute_distribute(TD_OP,TD_Prop,comp_grid,options)
#   (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,TD_Prop,P_sub,comp_grid,options);
#   val, t, bytes, gctime, memallocs = @timed (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,TD_Prop,P_sub,comp_grid,options);
#   println(t)
#   log_T_serial[i]=log_PARSDMM;
#   T_tot_serial[i]=t;
#
#   #parallel
#   println("")
#   println("parallel")
#   options.parallel=true
#   (P_sub,TD_OP,TD_Prop) = setup_constraints(constraint,comp_grid,options.FL)
#   (TD_OP,AtA,l,y) = PARSDMM_precompute_distribute(TD_OP,TD_Prop,comp_grid,options)
#   (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,TD_Prop,P_sub,comp_grid,options);
#   val, t, bytes, gctime, memallocs = @timed (x,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,TD_Prop,P_sub,comp_grid,options);
#   println(t)
#   log_T_parallel[i]=log_PARSDMM;
#   T_tot_parallel[i]=t;
#
#   #serial multilevel
#   println("")
#   println("serial multilevel")
#   options.parallel=false
#   n_levels=2
#   coarsening_factor=3
#   (m_levels,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels)=setup_multi_level_PARSDMM(m,n_levels,coarsening_factor,comp_grid,constraint,options)
#   (x,log_PARSDMM) = PARSDMM_multi_level(m_levels,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
#   val, t, bytes, gctime, memallocs = @timed (x,log_PARSDMM) = PARSDMM_multi_level(m_levels,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
#   println(t)
#   log_T_serial_multilevel[i]=log_PARSDMM;
#   T_tot_serial_multilevel[i]=t;
#
#   #parallel multilevel
#   println("")
#   println("parallel multilevel")
#   options.parallel=true
#   n_levels=2
#   coarsening_factor=3
#   (m_levels,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels)=setup_multi_level_PARSDMM(m,n_levels,coarsening_factor,comp_grid,constraint,options)
#   (x,log_PARSDMM) = PARSDMM_multi_level(m_levels,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
#   val, t, bytes, gctime, memallocs = @timed (x,log_PARSDMM) = PARSDMM_multi_level(m_levels,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
#   println(t)
#   log_T_parallel_multilevel[i]=log_PARSDMM;
#   T_tot_parallel_multilevel[i]=t;
#
#   if i==1
#     figure();imshow(reshape(x,(comp_grid.n[1],comp_grid.n[2]))');colorbar
#     savefig("projection_intersection_timings_2D_fig1_b.eps",bbox_inches="tight")
#     savefig("projection_intersection_timings_2D_fig1_b.eps",bbox_inches="tight")
#   elseif i==6
#     figure();imshow(reshape(x,(comp_grid.n[1],comp_grid.n[2]))');colorbar
#     savefig("projection_intersection_timings_2D_fig3_b.eps",bbox_inches="tight")
#     savefig("projection_intersection_timings_2D_fig3_b.eps",bbox_inches="tight")
#   end
#
#
#
# end
#
# #plot results
# fig, ax = subplots()
# ax[:loglog](N, T_tot_serial, label="serial",linewidth=5)
# ax[:loglog](N, T_tot_parallel, label="parallel",linewidth=5)
# ax[:loglog](N, T_tot_serial_multilevel, label="serial multilevel",linewidth=5)
# ax[:loglog](N, T_tot_parallel_multilevel, label="parallel multilevel",linewidth=5)
# ax[:legend]()
# title(string("time 2D vs grid size, JuliaThreads=",Threads.nthreads(),", BLAS threads=",ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ())), fontsize=12)
# xlabel("N gridpoints", fontsize=15)
# ylabel("time [seconds]", fontsize=15)
# savefig("projection_intersection_timings2D_2.eps",bbox_inches="tight")
# savefig("projection_intersection_timings2D_2.eps",bbox_inches="tight")

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
# very large values
const Infinity = 100000

const OUDANSENTE =  [1200,1200,900,600, 300, -10,   0,    0,    0]::Array{Int,1}
const OUDANGOTE   = [   0,   0,  0, 10,-300,-600,-900,-1200,-1200]::Array{Int,1}

#@iprofile begin
function Eval( nextMove::Int, p::Board, gs::GameStatus)
    score::Int = 0
    i::Int = 0
    bbp::BitBoard = uint128(0)
    from::Int = 0

    bbp = p.WhitePieces
    tebanP::Int = SENTE
    while bbp > uint128(0)
        from = trailing_zeros(bbp)
        bbp $= BitSet[from+1]
        q::Int = p.square[from+1]
        koma::Int = q & 0x0f
        score += MATERIAL[koma]
    end

    bbp = p.BlackPieces
    tebanP = GOTE
    while bbp > uint128(0)
        from = trailing_zeros(bbp)
        bbp $= BitSet[from+1]
        q::Int = p.square[from+1]
        koma::Int = q & 0x0f
        score -= MATERIAL[koma]
    end

    for i = MJFU:MJHI
        score += p.WhitePiecesInHands[i] * MATERIAL[i] #MATERIAL_CAPTURE[i]
        score -= p.BlackPiecesInHands[i] * MATERIAL[i] #MATERIAL_CAPTURE[i]
    end

    if (nextMove == GOTE)
        score = -score
    end
    return score
end
#end # @iprofile

max_stars_repo_name: kimrin/JapaneseChess

max_issues_repo_name: kimrin/JapaneseChess

max_forks_repo_name: kimrin/JapaneseChess


text:
using Dates

DFG_VERSION = "0.18.1";

@enum FactorType begin
    PRIORPOSE2
    POSE2POSE2
    POSE2APRILTAG4CORNERS
end

Base.@kwdef mutable struct FactorData
    eliminated::Bool = false
    potentialused::Bool = false
    edgeIDs::Vector{String} = []
    fnc::InferenceType
    multihypo::Vector{Float64} = []
    certainhypo::Vector{Int} = []
    nullhypo::Float64 = 0.0
    solveInProgress::Int = 0
    inflation::Float64 = 3.0
end

struct Factor
    label::String
    nstime::String
    fnctype::String
    _variableOrderSymbols::Vector{String}
    data::FactorData
    solvable::Int
    tags::Vector{String}
    timestamp::String
    _version::String
end

"""
$(SIGNATURES)
Create a prior factor for a ContinuousScalar (a.k.a. Pose1) with a distribution Z representing 1D prior information,
    e.g. `Normal(0.0, 0.1)`.

Default value of Z = `Normal(0.0, 0.1)`.
"""
function PriorData(;Z::Distribution = Normal(0.0, 0.1), kwargs...)::FactorData
    data = FactorData(;fnc = ZInferenceType(Z), certainhypo = [1], kwargs...)
    return data
end

"""
$(SIGNATURES)
Create a prior factor for a Pose2 with a distribution Z representing (x,y,theta) prior information,
    e.g. `FullNormal([0.0, 0.0, 0.0], diagm([0.01, 0.01, 0.01]))`.

Default value of Z = `FullNormal([0.0, 0.0, 0.0], diagm([0.01, 0.01, 0.01]))`.
"""
function PriorPose2Data(;Z::Distribution = FullNormal([0.0, 0.0, 0.0], diagm([0.01, 0.01, 0.01])), kwargs...)::FactorData
    data = FactorData(;fnc = ZInferenceType(Z), certainhypo = [1], kwargs...)
    return data
end

"""
$(SIGNATURES)
Create a prior factor for a Point2 with a distribution Z representing (x,y) prior information,
    e.g. `FullNormal([0.0, 0.0.0], diagm([0.01, 0.01]))`.

Default value of Z = `FullNormal([0.0, 0.0], diagm([0.01, 0.01]))`.
"""
function PriorPoint2Data(;Z::Distribution = FullNormal([0.0, 0.0], diagm([0.01, 0.01])), kwargs...)::FactorData
    data = FactorData(;fnc = ZInferenceType(Z), certainhypo = [1], kwargs...)
    return data
end

"""
$(SIGNATURES)
Create a ContinousScalar->ContinousScalar (also known as Pose1->Pose1) factor with a distribution Z representing the 1D relationship
between the variables, e.g. `Normal(1.0, 0.1)`.

Default value of Z = `Normal(1.0, 0.1)`.
"""
function LinearRelativeData(;Z::Distribution = Normal(1.0, 0.1), kwargs...)::FactorData
    data = FactorData(;fnc = ZInferenceType(Z), certainhypo = [1, 2], kwargs...)
    return data
end

"""
$(SIGNATURES)
Create a Pose2->Pose2 factor with a distribution Z representing the (x,y,theta) relationship
between the variables, e.g. `FullNormal([1,0,0.3333*π], diagm([0.01,0.01,0.01]))`.

Default value of Z = `FullNormal([1,0,0.3333*π], diagm([0.01,0.01,0.01]))`.
"""
function Pose2Pose2Data(;Z::Distribution = FullNormal([1,0,0.3333*π], diagm([0.01, 0.01, 0.01])), kwargs...)::FactorData
    data = FactorData(;fnc = ZInferenceType(Z), certainhypo = [1, 2], kwargs...)
    return data
end

"""
$(SIGNATURES)
Create a Pose2->Point2 bearing+range factor with 1D distributions:
- bearing: The bearing from the pose to the point, default `Normal(0, 1)`.
- range: The range from the pose to the point, default `Normal(1, 1)`.
"""
function Pose2Point2BearingRangeData(;bearing::Distribution = Normal(0, 1), range::Distribution = Normal(1, 1), kwargs...)::FactorData
    data = FactorData(;fnc = Pose2Point2BearingRangeInferenceType(bearing, range), certainhypo = [1, 2], kwargs...)
    return data
end

"""
$(SIGNATURES)
Create a Point2->Point2 range factor with a 1D distribution:
- range: The range from the pose to the point, default `Normal(1, 1)`.
"""
function Point2Point2RangeData(;range::Distribution = Normal(1, 1), kwargs...)::FactorData
    data = FactorData(;fnc = ZInferenceType(range), certainhypo = [1, 2], kwargs...)
    return data
end

"""
$(SIGNATURES)
Create a AprilTags factor that directly relates a Pose2 to the information from an AprilTag reading.
Corners need to be provided, homography and tag length are defaulted and can be overwritten.
"""
function Pose2AprilTag4CornersData(id, corners::Vector{Float64}, homography::Vector{Float64}; K::Vector{Float64}=[300.0,0.0,0.0,0.0,300.0,0.0,180.0,120.0,1.0], taglength::Float64=0.25, kwargs...)::FactorData
    fnc = Pose2AprilTag4CornersInferenceType(;corners,homography,K,taglength,id)
    data = FactorData(; fnc, certainhypo = [1, 2], kwargs...)
    return data
end

"""
$(SIGNATURES)

Returns `<:FactorData`
"""
function ScatterAlignPose2Data(
        varType::AbstractString,
        cloud1::AbstractVector{<:AbstractVector{<:Real}},
        cloud2::AbstractVector{<:AbstractVector{<:Real}},
        bw1::AbstractVector{<:Real}=Float64[],
        bw2::AbstractVector{<:Real}=Float64[];
        mkd1 = ManifoldKernelDensity(; varType, pts=cloud1, bw=bw1 ),
        mkd2 = ManifoldKernelDensity(; varType, pts=cloud2, bw=bw2 ),
        kw_sap::NamedTuple=(;),
        kwargs...
    )
    #

    fnc = ScatterAlignPose2InferenceType(; cloud1=mkd1, cloud2=mkd2, kw_sap...)
    data = FactorData(; fnc, certainhypo = [1, 2], kwargs...)
    return data
end

"""
$(SIGNATURES)
Create a Mixture factor type with an underlying factor type, a named set of
distributions that should be mixed, the probabilities of each distribution
(the mix), and the dimensions of the underlying factor (e.g.
ContinuousScalar=1, Pose2Pose2=3, etc.).

Args:
    mechanics (Type{FactorData}): The underlying factor data type, e.g. Pose2Pose2Data. NOTE: This will change in later versions
    but for now it can be any of the FactorData classes (e,g, LinearRelative, not the object LinearRelative()).
    components (NamedTuple): The named tuple set of distributions that
    should be mixed, e.g. NamedTuple(hypo1=Normal(0, 2)), hypo2=Uniform(30, 55)).
    probabilities (List[float]): The probabilities of each distribution (the mix), e.g. [0.4, 0.6].
    dims (int): The dimensions of the underlying factor, e.g. for Pose2Pose2 it's 3.
"""
function MixtureData(
        # TODO: Need to make this type constrained to Type{<:InferenceType}
        mechanics,
        # TODO: Need to made this something cleaner and more type stable, maybe OrderedDict
        components::NamedTuple,
        probabilities::Vector{Float64},
        dims::Integer  # TODO: Confirming we can remove.
        )::FactorData #where T<:FactorData
    data = FactorData(;
        fnc = MixtureInferenceType(
            N = length(components),
            # @jim-hill-r this is why I don't like the Data suffix.
            F_ = "Packed$(replace(string(mechanics), "Data" => ""))",
            S = collect(string.(keys(components))),
            components = collect(values(components)),
            diversity = Categorical(probabilities)
        ), certainhypo = [] # This should be updated in the Factor constructor below.
    )
    return data
end

function Factor(
        label::String,
        fncType::String,
        variableOrderSymbols::Vector{String},
        data::FactorData;
        tags::Vector{String}=["FACTOR"],
        timestamp::String = string(now(Dates.UTC))*"Z",
        multihypo=nothing,
        nullhypo=nothing
    )::Factor
    #
    # TODO: Remove independent updates of this and set certainhypo here.
    data.certainhypo = Vector{Int}(1:size(variableOrderSymbols)[1])
    if multihypo !== nothing
        data.multihypo = multihypo
    end
    if nullhypo !== nothing
        data.nullhypo = float(nullhypo)
    end

    result = Factor(
        label,
        "0",
        fncType,
        variableOrderSymbols,
        data,
        1,
        tags,
        timestamp,
        DFG_VERSION
    )
    return result
end


text:
############################################################################################
"""
$(TYPEDEF)

Stores information for a subset of 'model' parameter.

# Fields
$(TYPEDFIELDS)
"""
struct Tagged{A<:NamedTuple,B<:ParameterInfo}
    "Subset of ModelWrapper parameter names."
    parameter::A
    "Information about subset of parameter distributions, transformations and constraints, see [`ParameterInfo`](@ref)."
    info::B
    function Tagged(
        model::ModelWrapper, sym::S=keys(model.val)
    ) where {S<:Union{Symbol,NTuple{k,Symbol} where k}}
        ## If sym is a single symbol, convert to tuple
        if isa(sym, Symbol)
            sym = (sym,)
        end
        ## Check if all symbols are contained in ModelWrapper
        ArgCheck.@argcheck all(haskey(model.val, param) for param in sym) "Cannot tag parameter that is not contained in ModelWrapper.val"
        ## Convert Tuple to NamedTuple with fields == value of tuple -> can subset without allocations
        sym = Tuple_to_Namedtuple(sym, true)
        ## Generate new ParameterInfo based on sym subset
        info = ParameterInfo(
            subset(model.info.constraint, sym),
            subset(model.val, sym),
            model.info.reconstruct.default,
        )
        ## Return Tagged
        return new{typeof(sym),typeof(info)}(sym, info)
    end
end

############################################################################################
# Basic functions for Tagged struct
length(tagged::Tagged) = tagged.info.reconstruct.unflatten.strict._unflatten.sz[end]
paramnames(tagged::Tagged) = keys(tagged.parameter)

#A convenient method for evaluating a prior distribution of a NamedTuple parameter
function log_prior(tagged::Tagged, θ::NamedTuple)
    return log_prior(tagged.info.constraint, subset(θ, tagged.parameter))
end

############################################################################################
# Dispatch Tagged struct for .Core functions
function fill(model::ModelWrapper, tagged::Tagged, θ::NamedTuple)
    return merge(model.val, subset(θ, tagged.parameter))
end
function fill!(model::ModelWrapper, tagged::Tagged, θ::NamedTuple)
    model.val = merge(model.val, subset(θ, tagged.parameter))
    return nothing
end
function subset(model::ModelWrapper, tagged::Tagged)
    return subset(model.val, tagged.parameter)
end

#########################################
function unconstrain(model::ModelWrapper, tagged::Tagged)
    return unconstrain(tagged.info.transform, subset(model, tagged))
end
function flatten(model::ModelWrapper, tagged::Tagged)
    return flatten(tagged.info.reconstruct, subset(model, tagged))
end
function unconstrain_flatten(model::ModelWrapper, tagged::Tagged)
    return flatten(tagged.info.reconstruct, unconstrain(model, tagged))
end

#########################################
function unflatten(
    model::ModelWrapper, tagged::Tagged, θ::AbstractVector{T}
) where {T<:Real}
    return unflatten(tagged.info.reconstruct, θ)
end
function unflatten!(
    model::ModelWrapper, tagged::Tagged, θ::AbstractVector{T}
) where {T<:Real}
    model.val = merge(model.val, unflatten(model, tagged, θ))
    return nothing
end
function unflatten_constrain(
    model::ModelWrapper, tagged::Tagged, θᵤ::AbstractVector{T}
) where {T<:Real}
    return constrain(tagged.info.transform, unflatten(model, tagged, θᵤ))
end
function unflatten_constrain!(
    model::ModelWrapper, tagged::Tagged, θᵤ::AbstractVector{T}
) where {T<:Real}
    model.val = merge(model.val, unflatten_constrain(model, tagged, θᵤ))
    return nothing
end

#########################################
function sample(_rng::Random.AbstractRNG, model::ModelWrapper, tagged::Tagged)
    return merge(model.val, sample_constraint(_rng, tagged.info.constraint, subset(model, tagged)))
end
sample(model::ModelWrapper, tagged::Tagged) = sample(Random.GLOBAL_RNG, model, tagged)

function sample!(_rng::Random.AbstractRNG, model::ModelWrapper, tagged::Tagged)
    ArgCheck.@argcheck _checkprior(subset(tagged.info.constraint, tagged.parameter)) "For inplace sample version, all constraints need to be a Distribution."
    model.val = sample(_rng, model, tagged)
    return nothing
end
sample!(model::ModelWrapper, tagged::Tagged) = sample!(Random.GLOBAL_RNG, model, tagged)

#########################################
function log_prior(model::ModelWrapper, tagged::Tagged)
    return log_prior(tagged.info.constraint, subset(model, tagged))
end
function log_prior_with_transform(model::ModelWrapper, tagged::Tagged)
    return log_prior_with_transform(tagged.info.constraint, subset(model, tagged))
end
function log_abs_det_jac(model::ModelWrapper, tagged::Tagged)
    return log_abs_det_jac(tagged.info.transform.unconstrain, subset(model, tagged))
end

############################################################################################
export Tagged,
    length,
    fill,
    fill!,
    subset,
    flatten,
    unconstrain_flatten,
    unflatten,
    unflatten!,
    unconstrain,
    unflatten_constrain,
    unflatten_constrain!,
    sample,
    sample!,
    log_prior,
    log_prior_with_transform,
    log_abs_det_jac


text:
module ProbabilisticFluxVariationGradient

    using PyPlot

    using Suppressor, Printf, ProgressMeter

    using LinearAlgebra, Distributions, Statistics, StatsFuns, Random

    using Distances

    using Optim

    using Interpolations

    using PDMats

    using ApproximateVI

    using ForwardDiff

    export bpca, noisyintersectionvi, mockline, bpcaAIS

    # Utilities

    include("UTIL/astroutil.jl")
    include("UTIL/gputil.jl")
    include("UTIL/smallutil.jl")
    include("UTIL/plotresults.jl")
    lsstwaves = setlsstwavelengths()

    # Synthetic data generation

    include("SyntheticDataset/mockline.jl")
    include("SyntheticDataset/simulatedata.jl")

    # GP approach

    # include("GP/complete_lower_bound.jl")
    # include("GP/exact_expectations.jl")
    # include("GP/complete_lower_bound_grad.jl")
    # include("GP/exact_expectations_grad.jl")
    # include("GP/gpmodel.jl")
    # include("GP/gpmodelVIexact.jl")

    # Line intersection

    include("LineIntersection/noisyintersection.jl")
    include("LineIntersection/linemarginal.jl")

    # Probabilistic principal component algorithms

    include("PPCA/bpcaAIS.jl")
    include("PPCA/bpca.jl")
    include("PPCA/ppca.jl")

end


text:
struct Molecule{M} end

Molecule(molecule::Symbol) = Molecule{molecule}()
Molecule(molecule::String) = Molecule(Symbol(molecule))
macro Molecule_str(s)
    :(Molecule{$(Expr(:quote, Symbol(s)))})
end

# Interface
molecule(::Molecule{M}) where {M} = error("`molecule` not implemented for $M")
diis(::Molecule) = true
oda(::Molecule) = true


text:
# Activation

using Base.Broadcast
using CUDA.CUDNN: cudnnActivationForward!, cudnnOpTensor!,
            CUDNN_ACTIVATION_TANH,CUDNN_ACTIVATION_SIGMOID,CUDNN_ACTIVATION_ELU,
            CUDNN_ACTIVATION_RELU,CUDNN_ACTIVATION_CLIPPED_RELU,CUDNN_OP_TENSOR_MAX

for (f, op) in [
    CUDA.tanh       => (src,dst)->cudnnActivationForward!(dst, src, mode=CUDNN_ACTIVATION_TANH),
    NNlib.σ         => (src,dst)->cudnnActivationForward!(dst, src, mode=CUDNN_ACTIVATION_SIGMOID),
    NNlib.elu       => (src,dst)->cudnnActivationForward!(dst, src, mode=CUDNN_ACTIVATION_ELU),
    NNlib.relu      => (src,dst)->cudnnActivationForward!(dst, src, mode=CUDNN_ACTIVATION_RELU),
    # NNlib.relu6     => (src,dst)->cudnnActivationForward!(dst, src, mode=CUDNN_ACTIVATION_CLIPPED_RELU, coef=6.0),
    # NNlib.leakyrelu => (src,dst)->cudnnOpTensor!(dst, src, src; op=CUDNN_OP_TENSOR_MAX, alpha1=0.01),
    ]

    @eval begin
        # in-place
        function Base.materialize!(dst::DenseCuArray{<:CUDNNFloat},
                                   bc::Broadcast.Broadcasted{<:Any,<:Any,typeof($f),<:Tuple{DenseCuArray}})
            $op(bc.args[1], dst)
            return dst
        end

        # out of place
        function Base.materialize(bc::Broadcast.Broadcasted{<:Any,<:Any,typeof($f),<:Tuple{DenseCuArray}})
            ElType = Broadcast.combine_eltypes(bc.f, bc.args)
            dst = similar(bc, ElType)
            $op(bc.args[1], dst)
            return dst
        end
    end
end

# CUDNN_ACTIVATION_IDENTITY does not work with cudnnActivationForward
# FIXME: put this optimization in GPUArrays' `copyto!` (like Base.Broadcast's `copyto!`)
Base.broadcasted(::typeof(identity), x::DenseCuArray{T}) where {T<:CUDNNFloat} = x


text:
abstract type AbstractSpatial end
abstract type AbstractSpatialPoints <: AbstractSpatial end

mutable struct Spatial{B<:bbox, P} <: AbstractSpatial
    bbox::B
    proj::P
end

mutable struct SpatialPoints{C<:AbstractArray{<:Real, 2}, B<:bbox, P} <: AbstractSpatialPoints
    coords::C
    bbox::B
    proj::P
    # spatial::Spatial
end

"""
    SpatialPixelsDataFrame{
        T<:AbstractArray{<:Real, 2},
        C<:AbstractArray{<:Real, 2}, B<:bbox, P}  <: AbstractSpatialPoints

    SpatialPixelsDataFrame(data::AbstractArray{<:Real, 2}, coords::AbstractArray{<:Real, 2}, b::bbox, proj = 4326)

# Usage
`SpatialPixelsDataFrame(data, coords, bbox, proj)`
"""
Base.@kwdef mutable struct SpatialPixelsDataFrame{
    T<:AbstractArray{<:Real, 2},
    C<:AbstractArray{<:Real, 2},
    B<:bbox, P}  <: AbstractSpatialPoints

    data::T
    coords::C
    bbox::B
    proj::P = 4326
end

# SpatialPixelsDataFrame(data::AbstractArray{<:Real, 2},
#     coords::AbstractArray{<:Real, 2}, b::bbox, proj = 4326) =
#     SpatialPixelsDataFrame(data, coords, b, proj)

st_bbox(sp::AbstractSpatial) = sp.bbox
st_proj(sp::AbstractSpatial) = sp.proj
st_coords(sp::AbstractSpatialPoints) = sp.coords

export Spatial, SpatialPoints, SpatialPixelsDataFrame,
    st_bbox, st_proj, st_coords


text:
include("../src/LarSurf.jl")
# include("arr_fcn.jl")

import SparseArrays.spzeros
using Plasm
using LinearAlgebraicRepresentation

Lar = LinearAlgebraicRepresentation

# 3D 3x3x3

ibf = LarSurf.inner_boundary_filter([3, 3, 3])


text:
@doc raw"""
    Circle2

An object `c` of type [`Circle2`](@ref) is a circle in the two-dimensional
Euclidean plane ``Ε^2``.

The circle is oriented, i.e. its boundary has clockwise or counterclockwise
orientation. The boundary splits ``Ε^2`` into a positive and a negative
side, where the positive side is to the left of the boundary. The boundary also
splits ``Ε^2`` into a bounded and an unbounded side. Note that the circle
can be degenerated, i.e. the squared radius may be zero.
"""
Circle2

@doc raw"""
    Circle2(center::Point2, squared_radius::Real, ori::Orientation = COUNTERCLOCKWISE)

Introduces a variable `c` of type [`Circle2`](@ref).

It is initialized to the circle with center `center`, squared radius
`squared_radius`, and orientation `ori`.

!!! info "Precondition"

    `ori` ``≠`` `COLLINEAR`, and further, `squared_radius` ``≥ 0``.
"""
Circle2(center::Point2, squared_radius::Real, ori::Orientation = COUNTERCLOCKWISE)

@cxxdereference Circle2(c::Point2, r²::Real, ori::Union{Orientation,Ref{Orientation}} = COUNTERCLOCKWISE) =
    Circle2(c, convert(FT, r²), ori)

@doc raw"""
    Circle2(p::Point2, q::Point2, r::Point2)

Introduces a variable `c` of type [`Circle2`](@ref).

It is initialized to the unique circle which passes through the points `p`, `q`
and `r`. The orientation of the circle is the orientation of the point triple
`p`, `q`, `r`.

!!! info "Precondition"

    `p`, `q`, and `r` are not collinear.
"""
Circle2(p::Point2, q::Point2, r::Point2)

@doc raw"""
    Circle2(p::Point2, q::Point2, ori::Orientation = COUNTERCLOCKWISE)

Introduces a variable `c` of type [`Circle2`](@ref).

It is initialized to the circle with diameter ``\overline{pq}`` and orientation
`ori`.

!!! info "Precondition"

    `ori` ``≠`` `COLLINEAR`.
"""
Circle2(p::Point2, q::Point2, ori::Orientation = COUNTERCLOCKWISE)

@doc raw"""
    Circle2(p::Point2, ori::Orientation = COUNTERCLOCKWISE)

Introduces a variable `c` of type [`Circle2`](@ref).

It is initialized to the circle with center `center`, squared radius zero, and
orientation `ori`.

!!! info "Precondition"

    `ori` ``≠`` `COLLINEAR`.

!!! info "Postcondition"

    `is_degenerate(c) == true`.
"""
Circle2(p::Point2, ori::Orientation = COUNTERCLOCKWISE)

"""
    center(c::Circle2)

Returns the center of `c`.
"""
center(c::Circle2)

"""
    squared_radius(c::Circle2)

Returns the squared radius of `c`.
"""
squared_radius(c::Circle2)

"""
    orientation(c::Circle2)

Returns the orientation of `c`.
"""
orientation(c::Circle2)

"""
    ==(c::Circle2, circle2::Circle2)

Returns `true`, iff `c` and `circle2` are equal, i.e. if they have the same
center, same squared radius and same orientation.
"""
==(c::Circle2, circle2::Circle2)

"""
    is_degenerate(c::Circle2)

Returns `true`, iff `c` is degenerate, i.e. if `c` has squared radius zero.
"""
is_degenerate(c::Circle2)

"""
    oriented_side(c::Circle2, p::Point2)

Returns either the constant [`ON_ORIENTED_BOUNDARY`](@ref),
[`ON_POSITIVE_SIDE`](@ref), or [`ON_NEGATIVE_SIDE`](@ref), iff `p` lies on the
boundary, properly on the positive side, or properly on the negative side of
`c`, resp.
"""
oriented_side(c::Circle2, p::Point2)

"""
    bounded_side(c::Circle2, p::Point2)

Returns [`ON_BOUNDED_SIDE`](@ref), [`ON_BOUNDARY`](@ref), or
[`ON_UNBOUNDED_SIDE`](@ref) iff p lies properly inside, on the boundary, or
properly outside of `c`, resp.
"""
bounded_side(c::Circle2, p::Point2)

"""
    has_on_positive_side(c::Circle2, p::Point2)

Returns `true`, iff `p` lies on the positive side of `c`.
"""
has_on_positive_side(c::Circle2, p::Point2)

"""
    has_on_negative_side(c::Circle2, p::Point2)

Returns `true`, iff `p` lies on the negative side of `c`.
"""
has_on_negative_side(c::Circle2, p::Point2)

"""
    has_on_boundary(c::Circle2, p::Point2)

Returns `true`, iff `p` lies on the boundary of `c`.
"""
has_on_boundary(c::Circle2, p::Point2)

"""
    has_on_bounded_side(c::Circle2, p::Point2)

Returns `true`, iff `p` lies properly inside of `c`.
"""
has_on_bounded_side(c::Circle2, p::Point2)

"""
    has_on_unbounded_side(c::Circle2, p::Point2)

Returns `true`, iff `p` lies properly outside of `c`.
"""
has_on_unbounded_side(c::Circle2, p::Point2)

"""
    opposite(c::Circle2)

Returns the circle with the same center and squared radius as `c`, but with
opposite orientation.
"""
opposite(c::Circle2)

"""
    orthogonal_transform(c::Circle2, at::AffTransformation2)

Returns the circle obtained by applying `at` on `c`.

!!! info "Precondition"

    `at` is an orthogonal transformation.
"""
orthogonal_transform(c::Circle2, at::AffTransformation2)

"""
    bbox(c::Circle2)

Returns a bounding box containing `c`.
"""
bbox(c::Circle2)


text:
module Statistics

import Distributions
using Distributions
using StatsBase
# import Debug
# using Debug

import Base.mean, Base.std, Base.var

export mean, var, std,
		gaussian_regret, gaussian_percentile, gaussian_top_percentile,  gaussian_shortfall, gaussian_average_shortfall,
		gaussian_expected_shortfall, gaussian_downside_variance, gaussian_downside_deviation, gaussian_value_at_risk,
		semi_variance, semi_deviation, percentile, top_percentile, downside_variance, downside_deviation,
		shortfall, average_shortfall, expected_shortfall, value_at_risk, regret

#if !isdefined(Main.StatsBase,:weights)
#   weights(x) = x
#   mean(v::AbstractVector, w::AbstractVector) = sum(v.*w)/sum(w)
#end

function var(v::AbstractVector, w::AbstractVector, m)
	n=length(v)
	@assert (n>1)
	s=expectation_value(v, w, (x)->(x-m)^2, (x)->true) [1]

	return s*n/(n-1)
end
std(v::AbstractVector, w::AbstractVector) = sqrt(var(v, w))

var(v::AbstractVector, w::AbstractVector) = var(v, w, mean(v,weights(w)))

gaussian_downside_variance(v::AbstractVector, w::AbstractVector) = gaussian_regret(v,w,0)
gassian_downside_deviation(v::AbstractVector, w::AbstractVector) = sqrt(gaussian_downside_variance(v,w))

#Dembo, Freeman "The Rules Of Risk", Wiley (2001)
function gaussian_regret(v::AbstractVector, w::AbstractVector, target::Real)
	m=mean(v,weights(w))
	s = std(v, w)
	variance = var(v, w)
	f1 = variance + m*m -2*target*m + target*target
	d=Normal(m, s)
	alpha = cdf(d, target)
	f2 = m-target
	beta = variance*pdf(d, target)
	(alpha*f1 - beta*f2) / alpha
end

function gaussian_percentile(v::AbstractVector, w::AbstractVector, percentile::Real)
	@assert percentile > 0 && percentile < 1.0
	m=mean(v,weights(w))
	s=stdm(v,m)
	quantile(Normal(m,s), percentile)
end

gaussian_top_percentile(v::AbstractVector, w::AbstractVector, percentile::Real)=gaussian_percentile(v,w,1-percentile)

function gaussian_potential_upside(v::AbstractVector, w::AbstractVector, percentile::Real)
	@assert percentile <1.0 && percentile >= 0.9
	max(gaussian_percentile(v,w,percentile),0.0)
end

function gaussian_value_at_risk(v::AbstractVector, w::AbstractVector, percentile::Real)
	@assert percentile <1.0 && percentile >= 0.9
	-1*min(gaussian_percentile(v,w,1-percentile),0.0)
end

#Artzner, Delbaen, Eber and Heath, "Coherent measures of risk", Mathematical Finance 9 (1999)
function gaussian_expected_shortfall(v::AbstractVector, w::AbstractVector, percentile::Real)
	@assert percentile <1.0 && percentile >= 0.9
	m=mean(v,weights(w))
	s=stdm(v,m)
	d=Normal(m,s)
	var = quantile(1-percentile)
	g=pdf(d, var)
	r=m-s*s*g/(1-percentile)
	return -1*min(r,0)
end

function gaussian_shortfall(v::AbstractVector, w::AbstractVector, target::Real)
	m=mean(v,weights(w))
	s=stdm(v,m)
	cdf(Normal(m,s), target)
end

function gaussian_average_shortfall(v::AbstractVector, w::AbstractVector, target::Real)
	m=mean(v,weights(w))
	s=stdm(v,m)
	d=Normal(m,s)
	g=pdf(d,target)
	gi=cdf(d, target)

	(target - m) + s*s*g/gi
end

error_estimate(v::AbstractVector, w::AbstractVector) = sqrt(var(v,w)/length(v))

semi_variance(v::AbstractVector, w::AbstractVector) = regret(v, mean(v,weights(w)))
semi_deviation(v::AbstractVector, w::AbstractVector) = sqrt(semi_variance(v,weights(w)))

downside_variance(v) = regret(v, 0)
downside_deviation(v) = sqrt(downside_variance(v))

#Dembo and Freeman, "The Rules Of Risk", Wiley (2001)
function regret(v::AbstractVector, w::AbstractVector, target::Real)
	n=length(v)
	(r,_)=expectation_value(v, w, (x)->(x-target)^2, (x)->(x<target))
	r*n/(n-1)
end

potential_upside(v::AbstractVector, w::AbstractVector, p::Real) = max(percentile(v,w,p),0)
value_at_risk(v::AbstractVector, w::AbstractVector, p::Real) = -1*min(percentile(v,w,1-p),0)

#Artzner, Delbaen, Eber and Heath, "Coherent measures of risk", Mathematical Finance 9 (1999)
function expected_shortfall(v::AbstractVector, w::AbstractVector, p::Real)
	@assert p>= 0.9 && p <1
	target = -1* value_at_risk(v, w, p)

	(r,n)=expectation_value(v, w, (x)->(x), (x)->(x<target))
	return -1 * min(r, 0.0)
end

function average_shortfall(v::AbstractVector, w::AbstractVector, target::Real)
	return expectation_value(v, w, (x)->(target-x), (x)->(x<target)) [1]
end

function shortfall(v::AbstractVector, w::AbstractVector, target::Real)
	return expectation_value(v, w, (x)->(1), (x)->(x<target)) [1]
end

function expectation_value(v::AbstractVector, w::AbstractVector, mapfun::Function, filterfun::Function)
	@assert length(v) == length(w)
	@assert length(v) >0

	s=0
	sw=0
	n=0
	for i in 1:length(v)
		if filterfun(v[i])
			s = s + mapfun(v[i])*w[i]
			sw = sw + w[i]
			n=n+1
		end
	end

	if n==0
		 return NaN, 0
	end
	return s/sw , n
end

#Weighted percentile.
#Warning! Both value and weight arrays are copied and sorted
function percentile(v::AbstractVector, w::AbstractVector, p::Real)
	@assert p>0 && p<1
	wts = sum(w)
	@assert wts > 0
	vs, ord = Base.sortperm(v)

	k=1; l=length(v)
	i = w[ord[k]]
	target = p*wts

	while (i<target && k!=l)
		k=k+1
		i = i + w[ord[k]]
	end
	return vs[k]
end

function top_percentile(v::AbstractVector, w::AbstractVector, p::Real)
	@assert p>0 && p<1
	wts = sum(w)
	@assert wts > 0
	vs, ord = Base.sortperm(v)

	l=1; k=length(v)
	i = w[ord[k]]
	target = p*wts

	while (i<target && k!=l)
		k=k-1
		i = i + w[ord[k]]
	end
	return vs[k]
end

function kurtosis(v::AbstractVector, w::AbstractVector, m, corrected::Bool)
	n=length(v)
	vv=v-m
	y,_=expectation_value(v, w, (x)->(x-m)^3, (x)->true)
	if corrected
		@assert n>3
		c1=(n/(n-1.0)) * (n/(n-2.0)) * ((n+1.0)/(n-3.0))
		c2=3.0 * ((n-1.0)/(n-2.0)) * ((n-1.0)/(n-3.0));
		#(n*(n+1)*k)/((n-1)*(n-2)*(n-3))  - 3*(n-1)*(n-1)/( (n-2)*(n-3))
		c1*y/(var(v)^2) -c2
	else
		k - 3
	end
end

end


text:
function agm(x::T, y::T) where {T<:Real}
    a, b = x, y
    epsilon = eps(min(a,b))

    signbit(a) != signbit(b) && throw(DomainError("$a * $b must be nonnegative"))

    while abs(a-b) > epsilon
     c = (a + b) * 0.5
     b = sqrt(a * b)
     a = c
    end

    return a
end

agm1(x::T) where {T<:Real} = agm(one(T), x)


text:
#=

Define all the scientific constants needed in this module here

=#

# Avogadro's number:
Nₐ = 6.0221415e23;


text:
include("utils.jl")

function _test_d_proj(cone, tol)
    n = MOI.dimension(cone)
    x = randn(n)
    dx = 1e-6 * randn(n)
    proj_x = ConeProgramDiff._proj(x, cone)
    proj_xdx = ConeProgramDiff._proj(x+dx, cone)
    dproj_finite_diff = proj_xdx - proj_x

    dproj_test = ConeProgramDiff._d_proj(x, cone) * dx
    if !isapprox(dproj_finite_diff, dproj_test, atol=tol)
        println("Error!")
        error("Error for x = $x, dx = $dx, $cone\nnorm=$(norm(dproj_finite_diff -dproj_test))")
    end
end

function test_d_proj(cone::MOI.AbstractVectorSet; tol=1e-8)
    Random.seed!(0)
    for _ in 1:100
        _test_d_proj(cone, tol)
        _test_d_proj(MOI.dual_set(cone), tol)
    end
end

function test_d_proj(n::Int; tol=2e-6)
    cones = [
        MOI.Zeros(n),
        MOI.Nonnegatives(n),
        MOI.SecondOrderCone(n),
        MOI.PositiveSemidefiniteConeTriangle(n)
    ]
    for cone in cones
        # println("Testing $cone")
        test_d_proj(cone; tol=tol)
    end
    return true
end

function test_d_proj_exp(tol)
    function _helper(x, dx, tol; dual)
        proj_x = ConeProgramDiff._proj_exp_cone(x,  dual=dual)
        proj_xdx = ConeProgramDiff._proj_exp_cone(x+dx,  dual=dual)
        dproj_finite_diff = proj_xdx - proj_x
        dproj_test = ConeProgramDiff._d_proj_exp_cone(x; dual=dual) * dx
        # println(dproj_finite_diff)
        # println(dproj_test)
        @assert isapprox(dproj_finite_diff, dproj_test, atol=tol)
    end

    Random.seed!(0)
    case_p = zeros(4)
    case_d = zeros(4)
    for _ in 1:100
        x = randn(3)
        dx = 1e-6 * randn(3)
        case_p[det_case_exp_cone(x; dual=false)] += 1
        # println(det_case_exp_cone(x; dual=false))
        _helper(x, dx, tol; dual=false)

        case_d[det_case_exp_cone(x; dual=true)] += 1
        # println(det_case_exp_cone(x; dual=true))
        # println(x)
        _helper(x, dx, tol; dual=true)
    end
    @assert all(case_p .> 0) && all(case_d .> 0)
    return true
end

function test_d_proj_pow(tol)
    function _helper(x, α, dx, tol; dual)
        # println("$x \t $dx \t $(x+dx)")
        proj_x = ConeProgramDiff._proj_pow_cone(x, α; dual=dual)
        proj_xdx = ConeProgramDiff._proj_pow_cone(x+dx, α; dual=dual)
        dproj_finite_diff = proj_xdx - proj_x
        dproj_test = ConeProgramDiff._d_proj_pow_cone(x, α; dual=dual) * dx
        # println(dproj_finite_diff)
        # println(dproj_test)
        @assert isapprox(dproj_finite_diff, dproj_test, atol=tol)
    end

    Random.seed!(0)
    case_p = zeros(4)
    case_d = zeros(4)
    for _ in 1:100
        x = randn(3)
        dx = 1e-6 * randn(3)
        α = rand()*0.6 + 0.2

        # Need to get some into case 3
        if rand(1:10) == 1
            x[3] = 0
        end

        case_p[det_case_pow_cone(x, α; dual=false)] += 1
        _helper(x, α, dx, tol; dual=false)

        case_d[det_case_pow_cone(x, α; dual=true)] += 1
        _helper(x, α, dx, tol; dual=true)
    end
    @assert all(case_p .> 0) && all(case_d .> 0)
    return true
end

function test_d_project_onto_cone()
    function _test_d_proj_all_cones(x, cone_prod; tol=1e-6, dual=false)
        cones = dual ? [MOI.dual_set(c) for c in cone_prod] : cone_prod
        dx = 1e-7 * randn(length(x))
        Dpi = d_project_onto_cone(x, cones)
        proj_x = ConeProgramDiff.project_onto_cone(x, cones)
        proj_xdx = ConeProgramDiff.project_onto_cone(x+dx, cones)
        dproj_finite_diff = proj_xdx - proj_x

        dproj_test = Dpi * dx
        @assert isapprox(dproj_finite_diff, dproj_test, atol=tol)
        # if ~isapprox(dproj_finite_diff, dproj_test, atol=tol)
        #     println("Test failed:")
        #     println("determinant: $(det(Dpi))")
        #     inds = abs.(dproj_finite_diff - dproj_test) .>= tol
        #     println((dproj_finite_diff - dproj_test)[inds])
        #     println()
        # end
    end

    for _ in 1:10
        cone_prod = get_random_product_cone(10)
        N = sum([MOI.dimension(c) for c in cone_prod])
        local x = randn(N)
        _test_d_proj_all_cones(x, cone_prod; dual=false)
        _test_d_proj_all_cones(x, cone_prod; dual=true)
    end
    return true
end

@testset "D_Projections Test (Zero, Pos, SOC, PSD)" begin
    @test test_d_proj(10)
end

@testset "D_Projection Test: Exp Cone" begin
    @test test_d_proj_exp(1e-5)
end

@testset "D_Projection Test: Pow Cone" begin
    @test test_d_proj_pow(1e-6)
end

@testset "D_Projection: d_project_onto_cone()" begin
    test_d_project_onto_cone()
end


text:
# groupby aggregators

#
# Aggregators are structs with fit! and value methods
#
# fit!(aggregator, row) updates aggregator with row
# value(aggregator) returns current value
#

export fit!,
    value,
    Sum,
    Count,
    Avg

abstract type AbstractAggregator end

mutable struct Sum{F} <: AbstractAggregator
    valfunc::F
    val::Float64
    Sum(valfunc::F, val = 0.0) where {F <: Function} = new{F}(valfunc, val)
end

Base.eltype(::Type{Sum{F}}) where {F} = Float64

"""
    Sum(column::Symbol, val = 0.0)

Aggregator for sum of `column`.
"""
function Sum(column::Symbol, val = 0.0)
    valfunc = s -> s[column]
    Sum(valfunc, val)
end

function fit!(s::Sum, x)
    s.val += s.valfunc(x)
end

value(s::Sum) = s.val

"""
    Count()

Aggregator for count.
"""
mutable struct Count <: AbstractAggregator
    val::Int64
    Count() = new(0)
end

Base.eltype(::Type{Count}) = Int64

function fit!(c::Count, x)
    c.val += 1
end

value(c::Count) = c.val

mutable struct Avg{F} <: AbstractAggregator
    keyfunc::F
    total::Float64
    n::Int64
    Avg(keyfunc::F) where {F <: Function} = new{F}(keyfunc, 0.0, 0)
end

Base.eltype(::Type{Avg{F}}) where {F} = Float64

"""
    Avg(column::Symbol)

Aggregator for average of `column`.
"""
function Avg(column::Symbol)
    valfunc = s -> s[column]
    Avg(valfunc)
end

function fit!(a::Avg, x)
    a.total += a.keyfunc(x)
    a.n += 1
end

value(a::Avg) = a.total / a.n

# d is collection of key -> Aggregator pairs
# fit!(d, row) fit!s all aggregators with row
function fit!(d, x)
    for (k, v) in d
        fit!(v, x)
    end
end


text:
abstract(::AType{typeof(getfield)}, m::Const{Module}, f::Const{Symbol}) =
  Const(getfield(m.value, f.value))

abstract(::AType{typeof(getfield)}, m::AType{T}, f::Const{Symbol}) where T =
  fieldtype(T, f.value)

abstract(::AType{typeof(Core.apply_type)}, Ts::Const...) =
  Const(Core.apply_type(map(T -> T.value, Ts)...))

abstract(::AType{typeof(typeof)}, x::Const) = Const(widen(x))
abstract(::AType{typeof(typeof)}, x::AType{T}) where T =
  isconcretetype(T) ? Const(T) : Type

abstract(::AType{typeof(fieldtype)}, T::Const{<:Type}, f::Const{Symbol}) =
  Const(fieldtype(T.value, f.value))

for op in :[+, -, *, /, &, |, ^].args
  @eval abstract(::AType{typeof($op)}, ::AType{S}, ::AType{T}) where {S<:Number,T<:Number} =
          promote_type(S, T)
  @eval abstract(::AType{typeof($op)}, a::Const{<:Number}, b::Const{<:Number}) =
          Const($op(a.value, b.value))
end

for op in :[>, >=, <, <=, ==, !=].args
  @eval abstract(::AType{typeof($op)}, ::AType{S}, ::AType{T}) where {S<:Number,T<:Number} =
          Bool
  @eval abstract(::AType{typeof($op)}, a::Const{<:Number}, b::Const{<:Number}) =
          Const($op(a.value, b.value))
end

abstract(::AType{typeof(===)}, a, b) = Bool
abstract(::AType{typeof(===)}, a::Const, b::Const) = Const(a.value == b.value)

abstract(::AType{typeof(==)}, a, b) = Bool
abstract(::AType{typeof(==)}, a::Const, b::Const) = Const(a.value == b.value)

abstract(::AType{typeof(float)}, x::Const{<:Number}) = Const(float(x.value))

abstract(::AType{typeof(convert)}, ::Const{Type{T}}, x::Const{<:Number}) where T<:Number =
  Const(convert(T, x.value))

abstract(::AType{typeof(!)}, x::AType{Bool}) = Bool
abstract(::AType{typeof(!)}, x::Const{Bool}) = Const(!(x.value))

abstract(::AType{typeof(repr)}, x) = String
abstract(::AType{typeof(repr)}, x::Const) = Const(repr(x.value))
abstract(::AType{typeof(println)}, xs...) = Nothing
abstract(::AType{typeof(print)}, xs...) = Nothing


text:
# # Consumption Savings Problem
# In this case study, a household endowed with ``B_0`` dollars of wealth
# must decide how much to consume and save
# to maximize its utility over its finite lifecycle.

# ## Formulation

# The corresponding dynamic optimization problem is expressed:
# ```math
# \begin{aligned}
# 	&V(B_0,0) = &&\underset{c(t), B(t)}{\text{max}} \int_{t = 0}^{t=T} e^{-\rho t} u(c(t)) dt \\
# 	&\text{s.t.} \\
# 	&&& \frac{dB}{dt} = r \times B(t) - c(t), && t \in [0,T] \\
#   &&& B(0) = B_0 \\
#   &&& B(T) = 0
# \end{aligned}
# ```
# where
# the household lives during time ``t \in [0,T]``,
# the state variable ``B(t)`` is the household's stock of wealth at time ``t``,
# the choice variable ``c(t)`` is the household's consumption at time ``t``,
# ``r`` is the interest rate,
# ``B_0`` is the household's wealth endowment (initial condition),
# ``B(T) = 0`` is the terminal condition,
# and ``T`` is the time
# horizon.

# ## Model Definition

# Let's implement this in `InfiniteOpt` and first import the packages we need:
using InfiniteOpt, Ipopt

# We set the preference and constraint parameters:
ρ = 0.025  # discount rate
k = 100.0  # utility bliss point
T = 10.0   # life horizon
r = 0.05   # interest rate
B0 = 100.0 # endowment
u(c; k=k) = -(c - k)^2       # utility function
discount(t; ρ=ρ) = exp(-ρ*t) # discount function
BC(B, c; r=r) = r*B - c      # budget constraint

# We set the hyperparameters:
opt = Ipopt.Optimizer   # desired solver
ns = 1_000;             # number of points in the time grid

# We initialize the infinite model and choose the Ipopt solver:
m = InfiniteModel(opt)

# Let's specify our infinite parameter which is time ``t \in [0, T]``:
@infinite_parameter(m, t in [0, T], num_supports = ns)

# Now let's specify the variables:
@variable(m, B, Infinite(t)) ## state variables
@variable(m, c, Infinite(t)) ## control variables

# Specify the objective:
@objective(m, Max, integral(u(c), t, weight_func = discount))

# Set the initial/terminal conditions:
@constraint(m, B == B0, DomainRestrictions(t => 0))
@constraint(m, B == 0, DomainRestrictions(t => T))

# Set the budget constraint:
@constraint(m, c1, deriv(B, t) == BC(B, c; r=r))

# ## Problem Solution

# Optimize the model:
optimize!(m)
termination_status(m)

# Extract the results:
c_opt = value(c)
B_opt = value(B)
ts = supports(t)
opt_obj = objective_value(m) # V(B0, 0)

# Plot the results:
using Plots
ix = 2:(length(ts)-1) # index for plotting
plot(ts[ix],   B_opt[ix], lab = "B: wealth balance")
plot!(ts[ix],  c_opt[ix], lab = "c: consumption")

# That's it, now we have our optimal trajectory!

# This very simple problem has a closed form solution:
λ1 = exp((r)T)
λ2 = exp(-(r-ρ)T)
den = (λ1-λ2)r
Ω1 = (k + (r*B0-k)λ2)/den
Ω2 = (k + (r*B0-k)λ1)/den
c0 = r*B0 + (r)Ω1 + (r-ρ)Ω2
BB(t; k=k,r=r,ρ=ρ,Ω1=Ω1,Ω2=Ω2) = (k/r) - Ω1*exp((r)t) + Ω2*exp(-(r-ρ)t)
cc(t; k=k,r=r,ρ=ρ,c0=c0)       = k + (c0-k)*exp(-(r-ρ)t)

# Compare the solution given by `InfiniteOpt` with the closed form:
plot(legend=:topright);
plot!(ts[ix], c_opt[ix], color = 1, lab = "c: consumption, InfiniteOpt");
plot!(ts[ix], cc, color = 1, linestyle=:dash, lab = "c: consumption, closed form");
plot!(ts[ix], B_opt[ix], color = 4, lab = "B: wealth balance, InfiniteOpt");
plot!(ts[ix], BB, color = 4, linestyle=:dash, lab = "B: wealth balance, closed form")

# Not bad!

# ### Maintenance Tests
# These are here to ensure this example stays up to date.
using Test
@test termination_status(m) == MOI.LOCALLY_SOLVED
@test has_values(m)
@test B_opt isa Vector{<:Real}
@test c_opt isa Vector{<:Real}


text:
function Insertionsort(x::AbstractVector)
    for i = 2:length(x)
        index = x[i]
        j = i

        while j > 0 && x[j-1] > index
            x[j] = x[j-1]
            j -= 1
        end

        x[j] = index

    end

    return x
end


text:
module CancerData

export getcancerdata

using HTTP, CSV

"""
    getcancerdata()

Gets the cancer dataset for binary classification. Returns `features` and
`binary_response`. The feature matrix is standardized, the first column is a
vector of ones as the intercept. The binary responses have true denote Malignant
and false Benign. Internet connection is needed for this function.
"""
function getcancerdata()
    io = HTTP.get("https://raw.githubusercontent.com/MichielStock/STMO/master/chapters/04.Unconstrained/Data/BreastCancer.csv")
    cancer_data = CSV.read(io.body)
    binary_response = cancer_data.y .== "M"
    # extract feature matrix X
    features = Matrix(cancer_data[:,3:end])
    # standarizing features
    # this is needed for gradient descent to run faster
    features .-= sum(features, dims=1) / size(features, 1)
    features ./= sum(features.^2, dims=1).^0.5 / (size(features, 1).^0.5)
    # add intercept using dummy variable
    features = [ones(size(features, 1)) features]
    return features, binary_response
end

end #Words


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    AbstractProgressiveHedgingAttribute

Abstract supertype for attribute objects specific to the progressive-hedging algorithm.
"""
abstract type AbstractProgressiveHedgingAttribute <: AbstractStructuredOptimizerAttribute end
"""
    PrimalTolerance

An optimizer attribute for specifying the primal tolerance in the progressive-hedging algorithm.
"""
struct PrimalTolerance <: AbstractProgressiveHedgingAttribute end
"""
    DualTolerance

An optimizer attribute for specifying the dual tolerance in the progressive-hedging algorithm.
"""
struct DualTolerance <: AbstractProgressiveHedgingAttribute end
"""
    Regularizer

An optimizer attribute for specifying a penalization procedure to be used in the progressive-hedging algorithm. Options are:

- [`Fixed`](@ref):  Fixed penalty (default) ?Fixed for parameter descriptions.
- [`Adaptive`](@ref): Adaptive penalty update ?Adaptive for parameter descriptions.
"""
struct Penalizer <: AbstractProgressiveHedgingAttribute end
"""
    PenaltyTerm

An optimizer attribute for specifying what proximal term should be used in subproblemsof the progressive-hedging algorithm. Options are:

- [`Quadratic`](@ref) (default)
- [`Linearized`](@ref)
- [`InfNorm`](@ref)
- [`ManhattanNorm`](@ref)
"""
struct PenaltyTerm <: AbstractProgressiveHedgingAttribute end
"""
    PenalizationParameter

Abstract supertype for penalization-specific attributes.
"""
abstract type PenalizationParameter <: AbstractProgressiveHedgingAttribute end


text:
module CondNormalMMSE

include("FastMMSE.jl")
include("DiscreteMMSE.jl")
include("StructuredMMSE.jl")
include("MLEstimators.jl")

end

max_stars_repo_name: tum-msv/learning-mmse-est


text:
module test_26_TooManyLegends

using ModiaResult
using OrderedCollections
using Unitful
ModiaResult.@usingModiaPlot

t = range(0.0, stop=10.0, length=100)

result = OrderedDict{String,Any}()

result["time"] = t*u"s"
result["phi"]  = sin.(t)u"rad"
result["phi2"] = 0.5 * sin.(t)u"rad"
result["w"]    = cos.(t)u"rad/s"
result["w2"]   = 0.6 * cos.(t)u"rad/s"
result["r"]    = [[0.4 * cos(t[i]),
                   0.5 * sin(t[i]),
                   0.3 * cos(t[i])] for i in eachindex(t)]*u"m"

println("\n... test_26_TooManyLegends:")
ModiaResult.printResultInfo(result)

plot(result, ("phi", "r", "w", "w2"),
     maxLegend = 5,
     heading   = "Too many legends")

end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
@testset "TR_003.tif" begin
filename = joinpath(testbase, "TR_003.tif")
    # $ ScanImageTiffReader image shape TR_003.tif
    # Shape: 512 x 512 x 10 @ i16
    tsize = ScanImageTiffReader.open(filename) do io
        size(io)
    end
    @test tsize == (512, 512, 10)
    ttype = ScanImageTiffReader.open(filename) do io
        pxtype(io)
    end
    @test ttype == Int16

    # $ ScanImageTiffReader image bytes TR_003.tif
    # 5 MB
    dat = ScanImageTiffReader.open(filename) do io
        data(io)
    end
    @test abs(sizeof(dat)/2^20 - 5.0) < 1e-5
    # @test ScanImageTiffReader.open(filename, length) == size(dat, 3) # TODO: check this assumption

    # $ ScanImageTiffReader descriptions --frame 0 TR_003.tif
    # "Frame Tag = 00000001"
    # "scanimage.SI4.acqFrameBufferLength = 11"
    # "scanimage.SI4.acqFrameBufferLengthMin = 2"
    # "scanimage.SI4.acqFramesDone = 0"
    # "scanimage.SI4.acqNumAveragedFrames = 1"
    # "scanimage.SI4.acqNumFrames = 10"
    # "scanimage.SI4.acqState = 'idle'"
    # "scanimage.SI4.beamDirectMode = false"
    # "scanimage.SI4.beamFillFracAdjust = 14"
    # "scanimage.SI4.beamFlybackBlanking = true"
    # [...]
    desc = ScanImageTiffReader.open(filename) do io
        split(description(io, 1), "\n")
    end
    @test desc[1] == "Frame Tag = 00000001"
    @test desc[2] == "scanimage.SI4.acqFrameBufferLength = 11"
    @test desc[3] == "scanimage.SI4.acqFrameBufferLengthMin = 2"
    @test desc[4] == "scanimage.SI4.acqFramesDone = 0"
    @test desc[5] == "scanimage.SI4.acqNumAveragedFrames = 1"
    @test desc[6] == "scanimage.SI4.acqNumFrames = 10"
    @test desc[7] == "scanimage.SI4.acqState = 'idle'"
    @test desc[8] == "scanimage.SI4.beamDirectMode = false"
    @test desc[9] == "scanimage.SI4.beamFillFracAdjust = 14"
    @test desc[10] == "scanimage.SI4.beamFlybackBlanking = true"

    # $ ScanImageTiffReader metadata TR_003.tif
    #
    @test ScanImageTiffReader.open(metadata, filename) == ""
end


text:
# fitting

"""
    struct FittedCylinder{A<:AbstractArray, R<:Real} <: FittedShape

Cylinder primitive, defined by its axis direction,
a point that lies on its axis,
and its radius.
Also stored, if the normals point outwards of the shape.
"""
struct FittedCylinder{A<:AbstractArray, R<:Real} <: FittedShape
    axis::A
    center::A
    radius::R
    outwards::Bool
end

Base.show(io::IO, x::FittedCylinder) =
    print(io, """cylinder, R: $(x.radius)""")

Base.show(io::IO, ::MIME"text/plain", x::FittedCylinder{A, R}) where {A, R} =
    print(io, """FittedCylinder{$A, $R}\ncenter: $(x.center), axis: $(x.axis), R: $(x.radius), $(x.outwards ? "outwards" : "inwards" )""")

strt(x::FittedCylinder) = "cylinder"

function defaultshapeparameters(::Type{FittedCylinder})
    return (cylinder=(ϵ=0.3, α=deg2rad(5)),)
end

function setcylinderOuterity(fc, b)
    FittedCylinder(fc.axis, fc.center, fc.radius, b)
end

function fit2pointcylinder(p, n, params)
    #@unpack α_cylinder, ϵ_cylinder, parallelthrdeg = params
    @extract params : params_sphere=sphere
    @extract params_sphere : ϵ_cylinder=ϵ α_cylinder=α
    @extract params.common : parallelthrdeg
    # the normals are too parallel so cannot fit cylinder to it
    if abs(dot(n[1], n[2])) > cosd(parallelthrdeg)
        return nothing
    end

    an = normalize(n[1] × n[2])

    function lineplaneintersect(n, u, w)
        # http://geomalgorithms.com/a05-_intersect-1.html
        # n=plane normal
        # u= direction of the line
        # w=

        return dot(-n, w)/dot(n,u)
    end

    function project2plane(n, w)
        #n: normal
        #w: helyvektor
        return w+n*lineplaneintersect(n,n,w)
    end

    function projectto2d(xaxis, yaxis, zaxis, p1)
        # xax: pl a ponban muataó vektor
        xx = xaxis[1];
        xy = xaxis[2];
        xz = xaxis[3];

        yx = yaxis[1];
        yy = yaxis[2];
        yz = yaxis[3];

        zx = zaxis[1];
        zy = zaxis[2];
        zz = zaxis[3];

        px = p1[1];
        py = p1[2];
        pz = p1[3];

        r1 = -((-(pz*yy*zx) + py*yz*zx + pz*yx*zy - px*yz*zy - py*yx*zz + px*yy*zz)/(xz*yy*zx - xy*yz*zx - xz*yx*zy + xx*yz*zy + xy*yx*zz - xx*yy*zz))
        r2 = -((pz*xy*zx - py*xz*zx - pz*xx*zy + px*xz*zy + py*xx*zz - px*xy*zz)/(xz*yy*zx - xy*yz*zx - xz*yx*zy + xx*yz*zy + xy*yx*zz - xx*yy*zz))
        r3 = -((pz*xy*yx - py*xz*yx - pz*xx*yy + px*xz*yy + py*xx*yz - px*xy*yz)/(-(xz*yy*zx) + xy*yz*zx + xz*yx*zy - xx*yz*zy - xy*yx*zz + xx*yy*zz))

        return SVector(r1, r2)
    end

    function lineintersectionpoint(l1, l2)
        # l1 = [[p1x, p1y],[p2x,p2y]]
        a = l1[1]
        b = l1[2]
        c = l2[1]
        d = l2[2]

        amb_ = a-b
        cmd_ = c-d

        d1 = det(vcat(a',b'))
        d2 = det(vcat(c',d'))
        d3 = det(vcat(amb_',cmd_'))
        return (d1*cmd_-d2*amb_)/d3
    end

    xax = normalize(project2plane(an, p[1]))
    yax = normalize(an × xax)

    p11proj = projectto2d(xax, yax, an, project2plane(an, p[1]))
    p12proj = projectto2d(xax, yax, an, project2plane(an, p[1]+n[1]))

    p21proj = projectto2d(xax, yax, an, project2plane(an, p[2]))
    p22proj = projectto2d(xax, yax, an, project2plane(an, p[2]+n[2]))

    interc = lineintersectionpoint([p11proj,p12proj], [p21proj,p22proj])

    c = interc[1]*xax+interc[2]*yax

    nnormies = [norm(pt - c - an*dot( an, pt-c )) for pt in p[1:2]]
    R = (nnormies[1] + nnormies[2])/2
    # R = (norm(interc-p11proj) + norm(interc-p21proj))/2

    #p[1] - an*dot( an, p[1]-c )

    outw = dot(p12proj-p11proj, p11proj-interc) > 0

    return FittedCylinder(an, c, R, outw)
end

"""
    fit(::Type{FittedCylinder}, p, n, pc, params)

Fit a cylinder to 2 points.
Additional points and their normals are used to validate the fit.
Normals are expected to be normalized.
Return `nothing` if points do not fit to a cylinder.
"""
function fit(::Type{FittedCylinder}, p, n, pc, params)
    #@unpack ϵ_cylinder, α_cylinder, parallelthrdeg = params
    @extract params : params_sphere=sphere
    @extract params_sphere : ϵ_cylinder=ϵ α_cylinder=α
    @extract params.common : parallelthrdeg
    pl = length(p)
    @assert pl == length(n) "Size must be the same."
    @assert pl > 2 "Size must be at least 3."

    # "forcefit" a cylinder
    fc = fit2pointcylinder(p, n, params)

    fc === nothing && return nothing

    thr = cos(α_cylinder)
    vert_ok = falses(pl)
    norm_ok = falses(pl)
    invnorm_ok = falses(pl)
    for i in 1:pl
        # current normal
        curr_norm = p[i] - fc.axis*dot( fc.axis, p[i]-fc.center ) - fc.center
        # vertice check
        vert_ok[i] = abs(norm(curr_norm)-fc.radius) < ϵ_cylinder
        # normal check
        dotp = dot( normalize(curr_norm), n[i] )
        norm_ok[i] = dotp > thr
        invnorm_ok[i] = dotp < -thr
    end
    vert_ok == trues(pl) || return nothing
    norm_ok == trues(pl) && return setcylinderOuterity(fc, true)
    invnorm_ok == trues(pl) && return setcylinderOuterity(fc, false)
    return nothing
end

# bitmapping

function scorecandidate(pc, candidate::FittedCylinder, subsetID, params)
    ps = @view pc.vertices[pc.subsets[subsetID]]
    ns = @view pc.normals[pc.subsets[subsetID]]
    ens = @view pc.isenabled[pc.subsets[subsetID]]

    cp = compatiblesCylinder(candidate, ps, ns, params)
    inder = cp.&ens
    inpoints = (pc.subsets[subsetID])[inder]
    #inpoints = ((pc.subsets[1])[ens])[cp]
    score = estimatescore(length(pc.subsets[subsetID]), pc.size, length(inpoints))
    return (score, inpoints)
end

"""
    compatiblesCylinder(cylinder, points, normals, eps, alpharad)

Create a bool-indexer array for those points that are compatible to the cylinder.
Give back the projected points too for parameter space magic.

Compatibility is measured with an `eps` distance to the cylinder
and an `alpharad` angle to it's normal.
"""
function compatiblesCylinder(cylinder, points, normals, params)
    #@unpack ϵ_cylinder, α_cylinder = params
    @extract params : params_sphere=sphere
    @extract params_sphere : ϵ_cylinder=ϵ α_cylinder=α
    @assert length(points) == length(normals) "Size must be the same."

    c = cylinder.center
    R = cylinder.radius
    a = cylinder.axis

    comp = falses(length(points))

    for i in 1:length(points)
        curr_norm = points[i] - a*dot( a, points[i]-c ) - c
        # if the radius is correct
        if abs(norm(curr_norm)-R) < ϵ_cylinder
            if cylinder.outwards
                comp[i] = isparallel(normalize(curr_norm), normals[i], α_cylinder)
            else
                comp[i] = isparallel(-normalize(curr_norm), normals[i], α_cylinder)
            end
        end

        # playing with parameters:
        #TODO: implement it
    end
    return comp
end

"""
    refit(s::T, pc, params) where {T<:FittedCylinder}

Refit cylinder.
"""
function refit(s::T, pc, params) where {T<:FittedCylinder}
    # TODO: use octree for that
    pcv = @view pc.vertices[pc.isenabled]
    pcn = @view pc.normals[pc.isenabled]
    cp = compatiblesCylinder(s, pcv, pcn, params)
    return ExtractedShape(s, ((1:pc.size)[pc.isenabled])[cp])
end


text:
third32  = Float32(1)/Float32(3)
thirdR32 = Robust32(1)/Robust32(3)
two32  = Float32(2)
twoR32 = Robust32(2)
fourthirds = Float32(4)/Float32(3)

@testset "add/sub" begin
  @test two32 + two32 == twoR32 + twoR32
  @test two32 - third32 == twoR32 - thirdR32

  # @test fourthirds - (two32 - 2*third32) > 2.75 * (fourthirds - (twoR32 - 2*thirdR32))
end

@testset "mul/div" begin
  @test abs(2/9 -  (third32 * 2*third32)) > 9 *( abs(2/9 - (thirdR32 * 2*thirdR32)))
  @test two32 - third32 == twoR32 - thirdR32

  #@test fourthirds - (two32 - 2*third32) > 2.75 * (fourthirds - (twoR32 - 2*thirdR32))
end

@testset "root/power" begin
  @test abs(two32 - sqrt(sqrt(two32))^4) > abs(twoR32 - sqrt(sqrt(twoR32))^4)
  @test abs(two32 - sqrt(cbrt(two32))^6) > abs(twoR32 - sqrt(cbrt(twoR32))^6)
end


text:
#
# v.0.7 (with errors, to be fixed in Pccf2)
#
module Roc

function makeBinsOneSide(changes, maxDelay)
    # returns an array of tuples
    [(c,b) for (c,b) in zip(changes, [c + maxDelay for c in changes])]
end

function makeBinsTwoSide(changes, maxDelay)
    # returns an array of tuples
    right=[c+maxDelay for c in changes]
    left=[c-maxDelay for c in changes]
    [(l,r) for (l,r) in zip(left,right)]
end

function makeBins(changes,maxDelay)
    makeBinsTwoSide(changes,maxDelay)
end

function fallsInBins(val,bins::Array{Tuple{Int64,Int64}})
    # returns bins as an array of tuples: [(left,right)]
    bins[[e[1]<=val<=e[2] for e in bins]]
end
# println(fallsInBins(18, makeBins([1,20,30],3)))

function ifCdeFallsIntoBins(cde::Int64, bins::Array{Tuple{Int64,Int64}})
    length(fallsInBins(cde, bins))>0 ? true : false
end
# println(ifCdeFallsIntoBins(18, makeBins([1,20,30],3)))

function tpNum(detections,changes,maxDelay)
    # Calculate number of true positives
    bins = makeBins(changes, maxDelay)
    cdesInBins = Set(filter(cde->ifCdeFallsIntoBins(cde, bins), detections))
    length(cdesInBins)
end
#
function tpNum(detections,changes)
    intersect(Set(detections), Set(changes))
end

function faNum(detections,changes,maxDelay)
    length(detections)-tpNum(detections,changes,maxDelay)
end

function fnNum(detections,changes,maxDelay)
    # Calculate number of false negatives
    bins=makeBins(changes.maxDelay)
    length(bins)-tpNum(detections, changes, maxDelay)
end

function tpRate(detections,changes,maxDelay)
    # True positive TP rate = TP/(TP + FN); In case if there are
    # several true changes within the same bin we take min(1, tprate)
    tpn=tpNum(detections,changes,maxDelay)
    min(1, tpn/length(changes))
end

function faRate(detections,changes,maxDelay,totalNum)
    # False alarm/positive rate FPR = FP / (FP + TN) ; == specificity
    fan = faNum(detections,changes,maxDelay)
    tpn = tpNum(detections,changes,maxDelay)
    trueNegNum = totalNum - length(changes) * maxDelay
    fan/(fan+trueNegNum)
end

end

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
using Plots, Dates

func_evals=10000
n_vars=10
lb=-10ones(n_vars)
ub=15ones(n_vars)
x=lb+(ub-lb).*rand(n_vars)
x_evals=zeros(func_evals,n_vars)
x_evals[1,:]=lb
x_evals[2,:]=x
x_evals[3,:]=ub
function func(x)
    # sum(x.^2)

    n = length(x)
    1 + (1/4000)*sum(abs2, (x.+7)) - prod(cos.((x.+7) ./ sqrt.(1:n)))
end
opti_evals=zeros(func_evals)
opti_evals[1]=eval(func)(lb)
opti_evals[2]=eval(func)(x)
opti_evals[3]=eval(func)(ub)
iopt=sortperm(opti_evals[1:3])
opti_f=opti_evals[iopt[1]]
opti_x=x_evals[iopt[1],:]
x=copy(opti_x)
iter=3

while iter<func_evals
    global opti_f,opti_x,x_evals,opti_evals,iter
    for j=rand(1:n_vars)
        iter+=1; if iter>func_evals @goto escape_label end
        if iter>func_evals/3 alpha1=Int64(floor((func_evals/3)*(1-iter/func_evals)))+2 else alpha1=iter end
        # println(Dates.format(now(), "HH:MM:SS")," iter= ",iter," alpha1= ",alpha1)
        indsBEST=sortperm(opti_evals[1:iter-1])[1:alpha1-1]
        inp_x=x_evals[indsBEST,j];
        rr=minimum(inp_x)+rand()*(maximum(inp_x)-minimum(inp_x));
        x[j]=rr; fu=eval(func)(x); opti_evals[iter]=fu; x_evals[iter,:].=copy(x)

        if fu<opti_f
            opti_x.=copy(x);
            opti_f=copy(fu);
            println(Dates.format(now(), "HH:MM:SS")," iter= ",iter," otpi_f= ",opti_f," alpha1= ",alpha1)
        else
            x.=copy(vec(opti_x))
        end
    end
    @label escape_label
end
opti_evals.-=minimum(opti_evals); opti_evals./=maximum(opti_evals)
plot(sort(opti_evals,rev=true).^(1/10))

# using Optim
# res = Optim.optimize(func, fill(lb[1], n_vars), fill(ub[1], n_vars), x_evals[2,:], SAMIN(verbosity=0), Optim.Options(f_calls_limit=func_evals,store_trace =true));
# opti_x=res.minimizer; println(func(opti_x))
# res = Optim.optimize(func, fill(lb[1], n_vars), fill(ub[1], n_vars), x_evals[2,:], NelderMead(), Optim.Options(f_calls_limit=func_evals,store_trace =true))
# opti_x=res.minimizer; println(func(opti_x))
# res = Optim.optimize(func, fill(lb[1], n_vars), fill(ub[1], n_vars), x_evals[2,:], ParticleSwarm(), Optim.Options(f_calls_limit=func_evals,store_trace =true))
# opti_x=res.minimizer; println(func(opti_x))


text:
function chess!(r::Robot)
    move_till_border!(r)
    paint_vertival!(r)
    move!(r,Nord)
    paint_Horizon!(r)
    move_till_border!(r)
end

function move_till_border!(r::Robot)
    while isborder(r, Sud) == false
        move!(r, Sud)
    end
    while isborder(r, West) == false
        move!(r, West)
    end
end

function paint_vertival!(r::Robot)
    for i in 1:10
        for j in 1:5
            if isborder(r, Nord) == false
                putmarker!(r)
                move!(r, Nord)
                if isborder(r,Nord) == false
                    move!(r, Nord)
                end
            end
        end
        putmarker!(r)
        for k in 1:11
            if isborder(r, Sud) == false
                move!(r, Sud)
            end
        end
        if isborder(r, Ost) == false
            move!(r, Ost)
            if isborder(r,Ost) == false
                move!(r, Ost)
            end
        end
        if isborder(r,Ost) == true
            break
        end
    end
end

function paint_Horizon!(r::Robot)
    for i in 1:10
        for j in 1:5
            if isborder(r, West) == false
                putmarker!(r)
                move!(r, West)
                if isborder(r,West) == false
                    move!(r, West)
                end
            end
        end
        putmarker!(r)
        for k in 1:11
            if isborder(r, Ost) == false
                move!(r, Ost)
            end
        end
        if isborder(r, Nord) == false
            move!(r, Nord)
            if isborder(r,Nord) == false
                move!(r, Nord)
            end
        end
        if isborder(r,Nord) == true
            break
        end
    end
end


text:
function allocate(datakind::Symbol, dtype::DataType, dims...)
    if datakind == :local
        return zeros(dtype, dims...)
    elseif datakind == :shared
        return SharedArray{dtype}(dims...)
    else
        ErrorException("Unknown kind: " * string(datakind)) |> throw
    end
end

#=
"""
This function checks if CFL criteria is satisfied which is required by Euler Forward Scheme. Explicitly,

 K Δt       1
------  <= ---
(Δz)^2      2

for every layer. This function returns true every layer is satisfied, returns false if any of the layers is not.

"""
function checkAllDiffusionStability(;
    Δzs:: AbstractArray{Float64, 1},
    K  :: Float64,
    Δt :: Float64,
)

    return all( Δzs .>= √(2.0 * K * Δt) )
end

function checkDiffusionStability(;
    Δz :: Float64,
    K  :: Float64,
    Δt :: Float64,
)

    return Δz >= √(2.0 * K * Δt)
end

function checkDiffusionStability(oc::OceanColumn; Δt)
    return checkAllDiffusionStability(Δzs=oc.Δzs, K=oc.K, Δt=Δt)
end

function minΔz(;
    K :: Float64,
    Δt:: Float64,
)
    return √(2.0 * K * Δt)
end

=#

function boundMLD(h_ML::Float64; h_ML_min::Float64, h_ML_max::Float64)
    return max(min(h_ML, h_ML_max), h_ML_min)
end

"""
    getTKE(fric_u)

# Description
This function returns the TKE (turbulent kinetic energy) `k = 0.5 * (v'^2)` of ML. This parameterization is given by Kim 1976: "A Generalized Bulk Model of the Oceanic Mixed Layer" in its equation (11)

"""
function getTKE(;
    fric_u :: Float64
)
    cm = max(3e-2, 3.0 * fric_u)
    return  0.5 * cm^2.0
end

function updateFLDO!(
        ocn :: Ocean,
    )

    @loop_hor ocn i j let
        OC_updateFLDO!(ocn, i, j)
    end
end

function OC_updateFLDO!(
        ocn :: Ocean,
        i   :: Integer,
        j   :: Integer,
    )
    ocn.FLDO[i, j] = getFLDO(zs=ocn.cols.zs[i, j], h_ML=ocn.h_ML[i, j], Nz=ocn.Nz[i, j])
end

"""

    Returns the FLDO. If mixed-layer depth is equal to the total depth
    of ocean column, -1 will be returned.

"""
function getFLDO(;
    zs   :: AbstractArray{Float64,1},
    h_ML :: Float64,
    Nz   :: Integer,
)
    for i = 1:Nz
        if h_ML < - zs[i+1]  # I don't use equality in order to avoid Δb = 0 during some initialization
            return i
        end
    end

    return -1
end

function getLayerFromDepth(;
    zs   :: AbstractArray{Float64,1},
    z    :: Float64,
    Nz   :: Integer,
)
    for i = 1:Nz
        if zs[i+1] < z  # I don't use equality in order to avoid Δb = 0 during some initialization
            return i
        end
    end

    return -1
end

function getWindStress(;
    u10::Float64
)

    return u10 * 1e-3 * ( (u10 < 25.0)
                    ? 2.7 + 0.142 * u10 + 0.0764 * u10^2.0
                    : u10 * (2.16 + 0.5406 * (1.0 - exp(- (u10 - 25.0) / 7.5)))
    )

end

function getFricU(;
    ua::Float64
)
    return √(getWindStress(u10=ua) / ρ)
end

function OC_updateB!(
    ocn :: Ocean,
    i   :: Integer,
    j   :: Integer,
)

    ocn.b_ML[i, j] = TS2b(ocn.T_ML[i, j], ocn.S_ML[i, j])
    for k=1:ocn.Nz[i, j]
        ocn.bs[k, i, j] = TS2b(ocn.Ts[k, i, j], ocn.Ss[k, i, j])
    end

end

function updateB!(ocn::Ocean)

    @loop_hor ocn i j let
        OC_updateB!(ocn, i, j)
    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# Julia wrapper for header: /usr/include/scip/cons_disjunction.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrDisjunction(scip)
    ccall((:SCIPincludeConshdlrDisjunction, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPcreateConsDisjunction(scip, cons, name, nconss, conss, relaxcons, initial, enforce, check, _local, modifiable, dynamic)
    ccall((:SCIPcreateConsDisjunction, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_CONS}}, Ptr{SCIP_CONS}, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, nconss, conss, relaxcons, initial, enforce, check, _local, modifiable, dynamic)
end

function SCIPcreateConsBasicDisjunction(scip, cons, name, nconss, conss, relaxcons)
    ccall((:SCIPcreateConsBasicDisjunction, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_CONS}}, Ptr{SCIP_CONS}), scip, cons, name, nconss, conss, relaxcons)
end

function SCIPaddConsElemDisjunction(scip, cons, addcons)
    ccall((:SCIPaddConsElemDisjunction, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_CONS}), scip, cons, addcons)
end


text:
using LazySets, MathematicalPredicates, SparseArrays

function mna1_specification()
    # initial set: xᵢ ∈ [0.001, 0.0015] if i ≤ 2 and xᵢ = 0 otherwise
    X0 = Hyperrectangle(low=[fill(0.001, 2); zeros(576)],
                        high=[fill(0.0015, 2); zeros(576)])

    # safety property: x1 ≤ 0.5
    property = is_contained_in(HalfSpace(sparsevec([1], [1.0], 578), 0.5))

    # time horizon: 20 time units
    time_horizon = 20.0

    # specification
    O = Dict(:T => time_horizon, :property => property)

    return X0, O
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
function setup_primal_infeasibility()
    options = Dict(:verbose => false,
                   :eps_abs => 1e-05,
                   :eps_rel => 1e-05,
                   :eps_dual_inf => 1e-18,
                   :scaling => true)
    return options
end

tol = 1e-5

@testset "primal_infeasibility" begin

    @testset "primal_infeasible_problem" begin
        srand(1)

        n = 50
        m = 500
        P = sprandn(n, n, 0.6)
        P = P' * P
        q = randn(n)
        A = sprandn(m, n, 0.6)
        u = 3 .+ randn(m)
        l = -3 .+ randn(m)

        # Make problem infeasible
        A[Int(n/2), :] = A[Int(n/2) + 1, :]
        l[Int(n/2)] = u[Int(n/2) + 1] + 10 * rand()
        u[Int(n/2)] = l[Int(n/2)] + 0.5

        options = setup_primal_infeasibility()

        model = OSQP.Model()
        OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)
        results = OSQP.solve!(model)

        @test results.info.status == :Primal_infeasible

    end

    @testset "primal_dual_infeasible_problem" begin
        srand(1)

        n = 2
        m = 4
        P = spzeros(n, n)
        q = [-1., -1]
        A = sparse([1. -1; -1. 1.; 1. 0.; 0. 1.])
        l = [1., 1., 0., 0.]
        u = Inf * ones(m)

        options = setup_primal_infeasibility()

        model = OSQP.Model()
        OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)
        results = OSQP.solve!(model)

        @test results.info.status == :Primal_infeasible

    end
end


text:
# This file was generated, do not modify it. # hide
G3 = SearchGraph(; dist, db=MatrixDatabase(X))
callbacks3 = SearchGraphCallbacks(hyperparameters=OptimizeParameters(kind=MinRecall(0.95)))
@elapsed index!(G3; callbacks=callbacks3)

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
include("ClientAI.jl")
using ClientAI
using Base.Test

# Test basic coordinate type stuff
a = HexGridCoordinate(1.5, 2.5)
@test a.x == 1.5
@test a.y == 2.5

b = Position(1, 2)
d = Dict("x"=>1, "y"=>2)
@test b.x == 1 && b.y == 2
@test b == HexGridCoordinate(1, 2)
@test b == Position(d)
@test to_dict(b) == d

@test zero(Position) == Position(0,0)

# cube coordinates
@test HexGridCoordinate(CubeGridCoordinate(b)) == b

# test distance function
@test distance(Position(0,0), Position(1,0)) == 1
@test distance(Position(0,0), Position(0,1)) == 1
@test distance(Position(0,0), Position(-1,1)) == 1

# grid iteration
area = GeneralCircularArea(Position(1,1), 0)
@test collect(area) == [Position(1,1)]
area = GeneralCircularArea(Position(0,0), 1)
@test collect(area) == [Position(-1,0), Position(-1,1), Position(0,-1), Position(0,0), Position(0,1), Position(1,-1), Position(1,0)]
area = GeneralCircularArea(Position(0,0), 3)
@test length(area) == 37
for p in area
  @test distance(p, Position(0,0)) <= radius(area)
end

i = intersect(GeneralCircularArea(Position(1,0), 1), OriginCircularArea(1))
@test collect(i) == [Position(0,0), Position(0,1), Position(1, -1), Position(1, 0)]

@test length(collect(intersect(circle(Position(0,0), 2), circle(10)))) == length(collect(circle(Position(0,0), 2)))

max_stars_repo_name: Condimentum/spacetyckiting

max_issues_repo_name: Condimentum/spacetyckiting

max_forks_repo_name: Condimentum/spacetyckiting


text:
"""
  Automatic sweeping of function evaluations over a sample space.

  # AUTHORSHIP
    * Author    : Eduardo J. Alvarez
    * Email     : Edo.AlvarezR@gmail.com
    * Created   : Nov 2018
    * License   : MIT
"""
module SweepFunction

# ------------ EXPOSED FUNCTIONS------------------------------------------------
export ArgSweep, get_nvals
export Sweep, runsweep, printsweep, get_nevals

# ------------ GENERIC MODULES -------------------------------------------------
# None required

# ------------ FLOW LAB MODULES ------------------------------------------------
# None required

# ------------ GLOBAL VARIABLES AND DATA STRUCTURES ----------------------------
const module_path = splitdir(@__FILE__)[1]      # Path to this module
                                                # Default path to data files

# ------------ HEADERS ---------------------------------------------------------
# None required
for header_name in []
  include("SweepFunction_"*header_name*".jl")
end

################################################################################
# ARGSWEEP TYPE
################################################################################
"""
  `ArgSweep{T}(arg::Symbol, vals::Array{T, 1})`

Defines an argument `arg` to sweep and the set of values `vals` that the
argument will take.

NOTE: This is only for optional (keyword) arguments.
"""
type ArgSweep{T}
  arg::Symbol
  vals::Array{T, 1}
end

"""
  `get_nvals(self::ArgSweep)`

Returns the number of argument values
"""
function get_nvals(self::ArgSweep)
  return size(self.vals, 1)
end

"""
  `get_optarg(self::ArgSweep, i::Int64)`

Returns the i-th value of the argument
"""
function get_optarg(self::ArgSweep, i::Int64)
  if i>get_nvals(self)
    error("Invalid index $i; maximum is $(get_nvals(self))")
  end

  return (self.arg, self.vals[i])
end
##### END OF ARGSWEEP ##########################################################

################################################################################
# SWEEP TYPE
################################################################################

"""
  `Sweep(fun::Function, argsweep::Array{ArgSweep, 1},
cons_args::Array{Any, 1}, cons_optargs::Array{Tuple{Symbol, Any}, 1})`

Defines a sweep of function `fun` over multiple arguments `argsweep` with
constant required and optional arguments `cons_args` and `cons_optargs` (if
none, `cons_args` and `cons_optargs` can be omitted)
"""
type Sweep
  fun::Function                  # Function to sweep
  argsweep::Array{ArgSweep, 1}   # Variables to sweep

  # Constant required arguments
  cons_args::Array{Any, 1}

  # Constat optional arguments
  cons_optargs::Array{Tuple{Symbol, Any}, 1}

  Sweep(
          fun, argsweep,
          cons_args=Any[],
          cons_optargs=Tuple{Symbol, Any}[]
       ) = new(
          fun, argsweep,
          cons_args,
          cons_optargs
      )
end

"""
  `runsweep(self::Sweep)`

Runs the sweep evaluation of the function.

```jldoctest
julia> function fun(cons1, cons2; arg1=1.0, arg2=1.0, arg3=1.0)
           return arg1*arg2*arg3 + cons1 + cons2
       end;

julia> arg_sweep1 = ArgSweep{Int64}(:arg1, [1, 2, 3]);

julia> arg_sweep2 = ArgSweep{Float64}(:arg2, [4.0, 5.0]);

julia> arg_sweeps = [arg_sweep1, arg_sweep2];

julia> cons_args = [10, 20];

julia> sweep = Sweep(fun, arg_sweeps, cons_args);

julia> res = runsweep(sweep);

julia> printsweep(res)

"f(10, 20; arg1=1, arg2=4.0) = 34.0\n
f(10, 20; arg1=2, arg2=4.0) = 38.0\n
f(10, 20; arg1=3, arg2=4.0) = 42.0\n
f(10, 20; arg1=1, arg2=5.0) = 35.0\n
f(10, 20; arg1=2, arg2=5.0) = 40.0\n
f(10, 20; arg1=3, arg2=5.0) = 45.0\n"
```
"""
function runsweep(self::Sweep; threaded_ites=0, _out=[])

  # Base case: no variables left to sweep
  if size(self.argsweep, 1)==0

    # Evaluates the function
    val = self.fun(self.cons_args...; self.cons_optargs...)

    # Saves this evaluation
    this_run = Dict(
                    "arguments"              => self.cons_args,
                    "optional_arguments"     => self.cons_optargs,
                    "result"                 => val
                    )
    push!(_out, this_run)

  # Recursive case
  else

    # Sweeps last argument
    argsweep = self.argsweep[end]

    if get_nvals(argsweep)==0
      error("LOGIC ERROR: Received an argument sweep with 0 values! $argsweep")
    end

    # Parallelized runs
    if threaded_ites >= 1
        Threads.@threads for i in 1:get_nvals(argsweep) # Iterates over argument values

          # Get the optional argument pointer and its value
          this_optarg = get_optarg(argsweep, i)

          # Creates a new sweep with this argument being constant
          sweep = Sweep(self.fun,
                          self.argsweep[1:end-1],
                          self.cons_args,
                          vcat(this_optarg, self.cons_optargs)
                        )

          # Recursive call
          runsweep(sweep; threaded_ites=threaded_ites-1, _out=_out)
        end

    # Serial runs
    else
        for i in 1:get_nvals(argsweep) # Iterates over argument values

          # Get the optional argument pointer and its value
          this_optarg = get_optarg(argsweep, i)

          # Creates a new sweep with this argument being constant
          sweep = Sweep(self.fun,
                          self.argsweep[1:end-1],
                          self.cons_args,
                          vcat(this_optarg, self.cons_optargs)
                        )

          # Recursive call
          runsweep(sweep; threaded_ites=threaded_ites-1, _out=_out)
        end
    end
  end

  return _out
end

"""
  `printsweep(runsweep_output; verbose=true, lvl=0)`

Given the output of `runsweep()`, it will format it and print it as a string.
"""
function printsweep(runsweep_output; verbose=true, lvl=0, newln=false)

  out = ""

  for this_run in runsweep_output # Iterates over function calls

    # Builds the string showing the function
    str = "\t"^lvl*"f("

    # Required arguments
    for (i, val) in enumerate(this_run["arguments"])
      str *= string(val)
      str *= i!=size(this_run["arguments"], 1) ? ", " : ""
    end

    # Optional arguments
    for (i, (arg, val)) in enumerate(this_run["optional_arguments"])

      if i==1
        str *= "; "
      end

      str *= string(arg)*"="string(val)
      str *= i!=size(this_run["optional_arguments"], 1) ? ", " : ""
    end

    str *= ")"*("\n"*"\t"^(lvl+1))^newln

    # Function output
    str *= " = " * string(this_run["result"])

    out *= str*"\n"
  end

  if verbose
    println(out)
  end

  return out
end

"""
  `get_nevals(self::Sweep)`

Returns the number of function evaluations this sweep will take.
"""
function get_nevals(self::Sweep)
  return prod( [ i for i in
                [get_nvals(argsweep) for argsweep in self.argsweep
                ] if i!=0])
end
##### END OF SWEEP #############################################################

end # END OF MODULE


text:
const data = UInt[]
function create_lidict(target, data)
    udata = unique(data)
    lidict = Dict{UInt,Base.Profile.LineInfo}()
    for addr in udata
        if addr == 0
            continue
        end
        icxx"""
            lldb_private::Address addr;
            $target->ResolveLoadAddress($addr,addr);
            auto module = addr.GetModule();
            auto *function = addr.CalculateSymbolContextFunction();
            lldb_private::LineEntry line_entry;
            addr.CalculateSymbolContextLineEntry(line_entry);
            $:(begin
                funcnameptr = icxx"return function;" != C_NULL ?
                    icxx"return function->GetName().GetCString();" : C_NULL
                funcname = funcnameptr == C_NULL ? "???" :
                    bytestring(funcnameptr)
                filenameptr =  icxx"return !line_entry.file;" ?
                    C_NULL : icxx"return line_entry.file.GetFilename();"
                filename = filenameptr == C_NULL ? "???" :
                    bytestring(filenameptr)
                isjl = icxx"return module.get();" != C_NULL ?
                    !Gallium.isJuliaModule(icxx"return module;") :
                    true;
                line = icxx"return line_entry.line;"
                lidict[addr] = Base.Profile.LineInfo(
                funcname, filename, line,
                "", # For now, easy to do correctly
                0,
                isjl,
                reinterpret(Int64,addr))
                nothing
            end);
        """
    end
    lidict
end
function callback(env, ctx, id, loc)
    ctx = pcpp"lldb_private::StoppointCallbackContext"(ctx)
    exe_ctx = icxx"&$ctx->exe_ctx_ref;"
    process = icxx"$exe_ctx->GetProcessSP();"
    target = icxx"$exe_ctx->GetTargetSP().get();"
    profiler_thread = icxx"$exe_ctx->GetThreadSP();"
    tlist = icxx"$process->GetThreadList();"
    append!(data,map(frame->icxx"$frame->GetFrameCodeAddress().GetLoadAddress($target);",tlist[0]))
    push!(data,0)
    # Skip the rest of the data collection on the profiler thread
    Gallium.JumpToLine(profiler_thread,"signals-apple.c",295)
    return false
end
SetBreakpointAtLoc(callback,targets(dbg)[0],"signals-apple.c",238)
#=
using ProfileView
lidict = create_lidict(targets(dbg)[0],data)
ProfileView.view(data,C=true,lidict=lidict)
=#