text:
using CPLEX, JuMP
function cplex_raw_bigM(thePortfolio::SparsePortfolioData, ΔT_max::Float64=3600, gap::Float64=1e-4, debugLevel::Int=0)

  master_stock_problem_cplex=Model(solver=CplexSolver(CPX_PARAM_MIQCPSTRAT=1, CPX_PARAM_TILIM=ΔT_max, CPX_PARAM_EPGAP=gap, CPX_PARAM_THREADS=1))

  @variable(master_stock_problem_cplex, z[1:thePortfolio.n], Bin)
  @variable(master_stock_problem_cplex, x[1:thePortfolio.n]>=0)
  @constraint(master_stock_problem_cplex, sum(z[i] for i=1:thePortfolio.n)<=thePortfolio.k)
  @constraint(master_stock_problem_cplex, sparsity[i=1:thePortfolio.n], x[i]<=z[i])
  @constraint(master_stock_problem_cplex, minInvestment[i=1:thePortfolio.n], x[i]>=z[i]*thePortfolio.min_investment[i])
  @constraint(master_stock_problem_cplex, sumsToOne, sum(x[i] for i=1:thePortfolio.n)==1.0)
  @constraint(master_stock_problem_cplex, thePortfolio.A*x.>=thePortfolio.l)
  @constraint(master_stock_problem_cplex, thePortfolio.A*x.<=thePortfolio.u)
  @objective(master_stock_problem_cplex, Min, (0.5*sum((1.0/(thePortfolio.γ[i])*x[i]^2 for i=1:thePortfolio.n)) +0.5*x'*thePortfolio.X'*thePortfolio.X*x-thePortfolio.μ'*x+0.5*thePortfolio.Y'*thePortfolio.Y)')

  solve(master_stock_problem_cplex)

  optimalStocks=getvalue(x)

  ofv = getobjectivevalue(master_stock_problem_cplex)

  optimalIndicies=findall(a->a>0.5, getvalue(z))

  return CardinalityConstrainedPortfolio(optimalIndicies, optimalStocks[optimalIndicies], 0.0, zeros(thePortfolio.n), zeros(thePortfolio.m), zeros(thePortfolio.m), zeros(thePortfolio.m), getobjgap(master_stock_problem_cplex), getobjgap(master_stock_problem_cplex), 0, thePortfolio.μ'*optimalStocks, optimalStocks'*thePortfolio.X'*thePortfolio.X*optimalStocks, getsolvetime(master_stock_problem_cplex), getnodecount(master_stock_problem_cplex))

end

function cplex_raw_MISOCP(thePortfolio::SparsePortfolioData, ΔT_max::Float64=3600, gap::Float64=1e-4, debugLevel::Int=0)

  misocp_cplex=Model(solver=CplexSolver(CPX_PARAM_MIQCPSTRAT=1, CPX_PARAM_TILIM=ΔT_max, CPX_PARAM_EPGAP=gap, CPX_PARAM_THREADS=1))

  @variable(misocp_cplex, z[1:thePortfolio.n], Bin) # defines whether a stock is selected (big-M)
  @variable(misocp_cplex, x[1:thePortfolio.n]>=0.0)    # defines the quantity selected for each stock (short selling is permitted)
  @variable(misocp_cplex, θ[1:thePortfolio.n]>=0.0)
  @variable(misocp_cplex, τ>=0.0)

  @constraint(misocp_cplex, socpconstraint[i=1:thePortfolio.n], θ[i]+z[i]>=norm([2.0*x[i]; θ[i]-z[i]]))
  @constraint(misocp_cplex, sum(z[i] for i=1:thePortfolio.n)<=thePortfolio.k)
  @constraint(misocp_cplex, thePortfolio.A*x.>=thePortfolio.l)
  @constraint(misocp_cplex, thePortfolio.A*x.<=thePortfolio.u)
  @constraint(misocp_cplex, minInvestment[i=1:thePortfolio.n], x[i]>=z[i]*thePortfolio.min_investment[i])
  @constraint(misocp_cplex, sumsToOne, sum(x[i] for i=1:thePortfolio.n)==1.0)
  @constraint(misocp_cplex, τ+1.0>=norm([2.0*thePortfolio.X*x;τ-1.0]))
  @objective(misocp_cplex, Min, (0.5*sum(θ[i]/thePortfolio.γ[i] for i=1:thePortfolio.n)+0.5*τ-thePortfolio.μ'*x+0.5*thePortfolio.Y'*thePortfolio.Y)')

  solve(misocp_cplex)

  optimalStocks=getvalue(x)

  ofv = getobjectivevalue(misocp_cplex)
  optimalIndicies=findall(a->a>0.5, getvalue(z))
  numIters::Int128=0

  return CardinalityConstrainedPortfolio(optimalIndicies, optimalStocks[optimalIndicies], 0.0, zeros(thePortfolio.n), zeros(thePortfolio.m), zeros(thePortfolio.m),  zeros(thePortfolio.m), getobjgap(misocp_cplex), getobjgap(misocp_cplex), numIters, thePortfolio.μ'*optimalStocks, optimalStocks'*thePortfolio.X'*thePortfolio.X*optimalStocks, getsolvetime(misocp_cplex), getnodecount(misocp_cplex))

end

function cplex_raw_MISOCP_inout(thePortfolio::SparsePortfolioData, ΔT_max::Float64=3600, gap::Float64=1e-4, debugLevel::Int=0)

  misocp_cplex_inout=Model(solver=CplexSolver(CPX_PARAM_MIQCPSTRAT=1, CPX_PARAM_TILIM=ΔT_max, CPX_PARAM_EPGAP=gap, CPX_PARAM_THREADS=1))

  @variable(misocp_cplex_inout, z[1:thePortfolio.n], Bin) # defines whether a stock is selected (big-M)
  @variable(misocp_cplex_inout, x[1:thePortfolio.n]>=0.0)    # defines the quantity selected for each stock (short selling is permitted)
  @variable(misocp_cplex_inout, θ[1:thePortfolio.n]>=0.0)
  @variable(misocp_cplex_inout, τ>=0.0)
  @variable(misocp_cplex_inout, ofv)

  @constraint(misocp_cplex_inout, socpconstraint[i=1:thePortfolio.n], θ[i]+z[i]>=norm([2.0*x[i]; θ[i]-z[i]]))
  @constraint(misocp_cplex_inout, sum(z[i] for i=1:thePortfolio.n)<=thePortfolio.k)
  @constraint(misocp_cplex_inout, thePortfolio.A*x.>=thePortfolio.l)
  @constraint(misocp_cplex_inout, thePortfolio.A*x.<=thePortfolio.u)
  @constraint(misocp_cplex_inout, minInvestment[i=1:thePortfolio.n], x[i]>=z[i]*thePortfolio.min_investment[i])
  @constraint(misocp_cplex_inout, sumsToOne, sum(x[i] for i=1:thePortfolio.n)==1.0)
  @constraint(misocp_cplex_inout, τ+1.0>=norm([2.0*thePortfolio.X*x;τ-1.0]))
  @constraint(misocp_cplex_inout, ofv>=(0.5*sum(θ[i]/thePortfolio.γ[i] for i=1:thePortfolio.n)+0.5*τ-thePortfolio.μ'*x+0.5*thePortfolio.Y'*thePortfolio.Y)')
  @objective(misocp_cplex_inout, Min, ofv)

  # Add cuts from in-out method, to match what we do with the outer-approximation method...
  kelleyCutPool=Cut[]
  zSOCP=cplex_MISOCP_relaxation(thePortfolio, ΔT_max)
  stabilizationPoint=zSOCP
  kelleyCutPool=getKelleyPrimalCuts(thePortfolio, true, true, stabilizationPoint, 200)
  if length(kelleyCutPool) > 0
    for c in kelleyCutPool
       @constraint(misocp_cplex_inout, ofv >= c.p + dot(c.∇s, z))
    end
  end
  solve(misocp_cplex_inout)

  optimalStocks=getvalue(x)

  ofv = getobjectivevalue(misocp_cplex_inout)
  optimalIndicies=findall(a->a>0.5, getvalue(z))
  numIters::Int128=0

  return CardinalityConstrainedPortfolio(optimalIndicies, optimalStocks[optimalIndicies], 0.0, zeros(thePortfolio.n), zeros(thePortfolio.m), zeros(thePortfolio.m),  zeros(thePortfolio.m), getobjgap(misocp_cplex_inout), getobjgap(misocp_cplex_inout), numIters, thePortfolio.μ'*optimalStocks, optimalStocks'*thePortfolio.X'*thePortfolio.X*optimalStocks, getsolvetime(misocp_cplex_inout), getnodecount(misocp_cplex_inout))

end

# This function is currently used for diagnostic purposes only (i.e. not part of the main code base)
function cplex_raw_MISOCP_subset(thePortfolio::SparsePortfolioData, z::Array{Float64})

  socp_mosek=Model(solver=MosekSolver(MSK_DPAR_INTPNT_QO_TOL_PFEAS=1e-6, MSK_DPAR_INTPNT_QO_TOL_DFEAS=1e-6, MSK_IPAR_LOG=0, MSK_IPAR_MAX_NUM_WARNINGS=0))

  @variable(socp_mosek, x[1:thePortfolio.n]>=0.0)    # defines the quantity selected for each stock (short selling is permitted)
  @variable(socp_mosek, θ[1:thePortfolio.n]>=0.0)
  @variable(socp_mosek, τ>=0.0)

  @constraint(socp_mosek, socpconstraint[i=1:thePortfolio.n], θ[i]+z[i]>=norm([2.0*x[i]; θ[i]-z[i]]))
  @constraint(socp_mosek, sum(z[i] for i=1:thePortfolio.n)<=thePortfolio.k)
  @constraint(socp_mosek, thePortfolio.A*x.>=thePortfolio.l)
  @constraint(socp_mosek, thePortfolio.A*x.<=thePortfolio.u)
  @constraint(socp_mosek, minInvestment[i=1:thePortfolio.n], x[i]>=z[i]*thePortfolio.min_investment[i])
  @constraint(socp_mosek, sumsToOne, sum(x[i] for i=1:thePortfolio.n)==1.0)
  @constraint(socp_mosek, τ+1.0>=norm([2.0*thePortfolio.X*x;τ-1.0]))
  @objective(socp_mosek, Min, (0.5*sum(θ[i]/thePortfolio.γ[i] for i=1:thePortfolio.n)+0.5*τ-thePortfolio.μ'*x+0.5*thePortfolio.Y'*thePortfolio.Y)')

  solve(socp_mosek)

  ofv = getobjectivevalue(socp_mosek)
end

function cplex_raw_MISOCP_subset2(thePortfolio::SparsePortfolioData, z::Array{Float64})

  socp_mosek=Model(solver=MosekSolver(MSK_DPAR_INTPNT_QO_TOL_PFEAS=1e-6, MSK_DPAR_INTPNT_QO_TOL_DFEAS=1e-6, MSK_IPAR_LOG=0, MSK_IPAR_MAX_NUM_WARNINGS=0))

  @variable(socp_mosek, x[1:thePortfolio.n]>=0.0)    # defines the quantity selected for each stock (short selling is permitted)
  @variable(socp_mosek, θ[1:thePortfolio.n]>=0.0)
  @variable(socp_mosek, τ>=0.0)

  @constraint(socp_mosek, socpconstraint[i=1:thePortfolio.n], θ[i]+z[i]>=norm([2.0*x[i]; θ[i]-z[i]]))
  @constraint(socp_mosek, sum(z[i] for i=1:thePortfolio.n)<=thePortfolio.k)
  @constraint(socp_mosek, thePortfolio.A*x.>=thePortfolio.l)
  @constraint(socp_mosek, thePortfolio.A*x.<=thePortfolio.u)
  @constraint(socp_mosek, minInvestment[i=1:thePortfolio.n], x[i]>=z[i]*thePortfolio.min_investment[i])
  @constraint(socp_mosek, sumsToOne, sum(x[i] for i=1:thePortfolio.n)==1.0)
  @constraint(socp_mosek, τ+1.0>=norm([2.0*thePortfolio.X*x;τ-1.0]))
  @objective(socp_mosek, Min, (0.5*sum(θ[i]/thePortfolio.γ[i] for i=1:thePortfolio.n)+0.5*τ-thePortfolio.μ'*x+0.5*thePortfolio.Y'*thePortfolio.Y)')

  theStatus=solve(socp_mosek)

  ofv = getobjectivevalue(socp_mosek)
  return theStatus, getvalue(x)
end

function cplex_MISOCP_relaxation(thePortfolio::SparsePortfolioData, ΔT_max::Float64=3600)

  misocp_cplex=Model(solver=CplexSolver(CPX_PARAM_SCRIND=0, CPX_PARAM_TILIM=ΔT_max))

  @variable(misocp_cplex, z[1:thePortfolio.n]>=0.0) # defines whether a stock is selected (big-M)
  @constraint(misocp_cplex, z.<=ones(thePortfolio.n))
  @variable(misocp_cplex, x[1:thePortfolio.n]>=0.0)    # defines the quantity selected for each stock (short selling is permitted)
  @variable(misocp_cplex, θ[1:thePortfolio.n]>=0.0)
  @variable(misocp_cplex, τ>=0.0)

  @constraint(misocp_cplex, socpconstraint[i=1:thePortfolio.n], θ[i]+z[i]>=norm([2.0*x[i]; θ[i]-z[i]]))
  @constraint(misocp_cplex, sum(z[i] for i=1:thePortfolio.n)<=thePortfolio.k)
  @constraint(misocp_cplex, thePortfolio.A*x.>=thePortfolio.l)
  @constraint(misocp_cplex, thePortfolio.A*x.<=thePortfolio.u)
  @constraint(misocp_cplex, minInvestment[i=1:thePortfolio.n], x[i]>=z[i]*thePortfolio.min_investment[i])
  @constraint(misocp_cplex, sumsToOne, sum(x[i] for i=1:thePortfolio.n)==1.0)
  @constraint(misocp_cplex, τ+1.0>=norm([2.0*thePortfolio.X*x;τ-1.0]))
  @objective(misocp_cplex, Min, (0.5*sum(θ[i]/thePortfolio.γ[i] for i=1:thePortfolio.n)+0.5*τ-thePortfolio.μ'*x+0.5*thePortfolio.Y'*thePortfolio.Y)')

  solve(misocp_cplex)

  optimalStocks=getvalue(x)

  ofv = getobjectivevalue(misocp_cplex)
  optimalIndicies=findall(a->a>0.5, getvalue(z))
  numIters::Int128=0

  return getvalue(z)

end


text:
using IndexFunArrays, Documenter

 # set seed fixed for documentation
DocMeta.setdocmeta!(IndexFunArrays, :DocTestSetup, :(using IndexFunArrays); recursive=true)
makedocs(modules = [IndexFunArrays],
         sitename = "IndexFunArrays.jl",
         pages = Any[
            "IndexFunArrays.jl" => "index.md",
            "Misc" => "misc.md",
            "Distance Functions" => "distance.md",
            "Window Functions" => "window.md",
         ]
        )

deploydocs(repo = "github.com/bionanoimaging/IndexFunArrays.jl.git", devbranch="master")

ext: jl


text:
module TestGraphviz

using Test
import LightGraphs, MetaGraphs
using LightGraphs: add_vertex!, add_vertices!, add_edge!
using MetaGraphs: MetaGraph, MetaDiGraph

using Catlab.Graphics.Graphviz

# MetaGraphs
############

# Undirected simple graph
g = MetaGraph()
add_vertex!(g, :label, "v")
add_vertices!(g, 2)
add_edge!(g, 1, 2, :xlabel, "e")
gv = to_graphviz(g)::Graph
@test !gv.directed
nodes = filter(s -> isa(s,Node), gv.stmts)
edges = filter(s -> isa(s,Edge), gv.stmts)
@test length(nodes) == 3
@test length(edges) == 1
@test nodes[1].attrs[:label] == "v"
@test edges[1].attrs[:xlabel] == "e"

# Directed simple graph
g = MetaDiGraph()
add_vertices!(g, 3); add_edge!(g, 1, 2); add_edge!(g, 2, 3)
gv = to_graphviz(g)::Graph
@test gv.directed
@test length(filter(s -> isa(s,Node), gv.stmts)) == 3
@test length(filter(s -> isa(s,Edge), gv.stmts)) == 2

# Directed multigraph
g = MetaDiGraph()
add_vertices!(g, 2)
add_edge!(g, 1, 2, :edges, [Dict(:label => "e1"), Dict(:label => "e2")])
gv = to_graphviz(g; multigraph=true)::Graph
@test gv.directed
nodes = filter(s -> isa(s,Node), gv.stmts)
edges = filter(s -> isa(s,Edge), gv.stmts)
@test length(nodes) == 2
@test length(edges) == 2
@test edges[1].attrs[:label] == "e1"
@test edges[2].attrs[:label] == "e2"

# Pretty-print
##############

spprint(expr::Expression) = sprint(pprint, expr)

# Node statement
@test spprint(Node("n")) == "n;"
@test spprint(Node("n"; label="foo")) == "n [label=\"foo\"];"
@test spprint(Node("n"; shape="box", style="filled")) ==
  "n [shape=\"box\",style=\"filled\"];"

# Edge statement
@test spprint(Edge("n1","n2")) == "n1 -- n2;"
@test spprint(Edge("n1","n2"; label="foo")) == "n1 -- n2 [label=\"foo\"];"
@test spprint(Edge("n1","n2"; style="dotted", weight="10")) ==
  "n1 -- n2 [style=\"dotted\",weight=\"10\"];"
@test spprint(Edge(NodeID("n1","p1"), NodeID("n2","p2"))) == "n1:p1 -- n2:p2;"
@test spprint(Edge(NodeID("n1","p1"), NodeID("n2","p2"); label="bar")) ==
  "n1:p1 -- n2:p2 [label=\"bar\"];"
@test spprint(Edge(NodeID("n1","p1","w"), NodeID("n2", "p2", "e"))) ==
  "n1:p1:w -- n2:p2:e;"
@test spprint(Edge("n1","n2","n3")) == "n1 -- n2 -- n3;"
@test spprint(Edge(NodeID("n1"), NodeID("n2"), NodeID("n3"))) ==
  "n1 -- n2 -- n3;"

# Graph statement
graph = Graph("G",
  Node("n1"),
  Node("n2"),
  Edge("n1","n2")
)
@test spprint(graph) == """
graph G {
  n1;
  n2;
  n1 -- n2;
}
"""

graph = Digraph("G",
  Node("n1"),
  Node("n2"),
  Edge("n1","n2")
)
@test spprint(graph) == """
digraph G {
  n1;
  n2;
  n1 -> n2;
}
"""

graph = Digraph("G",
  Node("n1"; label="foo"),
  Node("n2"; label="bar"),
  Edge("n1","n2");
  graph_attrs = Attributes(:rankdir => "LR"),
  node_attrs = Attributes(:shape => "box", :style => "filled"),
  edge_attrs = Attributes(:style => "dotted")
)
@test spprint(graph) == """
digraph G {
  graph [rankdir="LR"];
  node [shape="box",style="filled"];
  edge [style="dotted"];
  n1 [label="foo"];
  n2 [label="bar"];
  n1 -> n2;
}
"""

# Subgraph statement
subgraph = Subgraph("sub",
  Node("n1"),
  Node("n2"),
  Edge("n1","n2")
)
@test spprint(subgraph) == """
subgraph sub {
  n1;
  n2;
  n1 -- n2;
}"""

subgraph = Subgraph(
  Node("n1"),
  Node("n2"),
  graph_attrs = Attributes(:rank => "same")
)
@test spprint(subgraph) == """
{
  graph [rank="same"];
  n1;
  n2;
}"""

end


text:
# Fréchet Distribution

immutable Frechet <: ContinuousUnivariateDistribution
    shape::Float64
    scale::Float64
    function Frechet(sh::Real, sc::Real)
    	zero(sh) < sh && zero(sc) < sc || error("Both shape and scale must be positive")
    	new(float64(sh), float64(sc))
    end
end

Frechet(sh::Real) = Frechet(sh, 1.0)

## Support
@distr_support Frechet 0.0 Inf

## Properties
mean(d::Frechet) = d.shape > 1.0 ? d.scale * gamma(1.0 - 1.0 / d.shape) : Inf
median(d::Frechet) = d.scale * log(2)^(-1.0 / d.shape)

mode(d::Frechet) = (ik = -1.0/d.shape; d.scale * (1.0-ik)^ik)

function var(d::Frechet)
    ia = 1.0 / d.shape
    d.shape > 2.0 ? abs2(d.scale) * (gamma(1.0 - 2.0 * ia) - abs2(gamma(1.0 - ia))) : Inf
end

function skewness(d::Frechet)
    d.shape <= 3.0 && return NaN
    tmp_mean = mean(d)
    tmp_var = var(d)
    tmp = gamma(1.0 - 3.0 / d.shape) * d.scale^3
    tmp -= 3.0 * tmp_mean * tmp_var
    tmp -= tmp_mean^3
    return tmp / tmp_var / sqrt(tmp_var)
end

function kurtosis(d::Frechet)
    d.shape <= 4.0 && return NaN
    λ, k = d.scale, d.shape
    μ = mean(d)
    σ = std(d)
    γ = skewness(d)
    den = λ^4 * gamma(1.0 - 4.0 / k) -
          4.0 * γ * σ^3 * μ -
          6.0 * μ^2 * σ^2 - μ^4
    num = σ^4
    return den / num - 3.0
end

function entropy(d::Frechet)
    const γ = 0.57721566490153286060
    1.0 + γ / d.shape + γ + log(d.scale / d.shape)
end

## Functions
function pdf(d::Frechet, x::Float64)
    x < zero(x) && return 0.0
    a = d.scale/x
    d.shape/d.scale * a^(d.shape+1.0) * exp(-a^d.shape)
end
function logpdf(d::Frechet, x::Float64)
    x < zero(x) && return -Inf
    a = d.scale/x
    log(d.shape/d.scale) + (d.shape+1.0)*log(a) - a^d.shape
end

cdf(d::Frechet, x::Float64) = x <= zero(x) ? 0.0 : exp(-((d.scale / x)^d.shape))
ccdf(d::Frechet, x::Float64) = x <= zero(x) ? 1.0 : -expm1(-((d.scale / x)^d.shape))
logcdf(d::Frechet, x::Float64) = x <= zero(x) ? -Inf : -(d.scale / x)^d.shape
logccdf(d::Frechet, x::Float64) = x <= zero(x) ? 0.0 :  log1mexp(-((d.scale / x)^d.shape))

quantile(d::Frechet, p::Float64) = @checkquantile p d.scale*(-log(p))^(-1/d.shape)
cquantile(d::Frechet, p::Float64) = @checkquantile p d.scale*(-log1p(-p))^(-1/d.shape)
invlogcdf(d::Frechet, lp::Float64) = lp > zero(lp) ? NaN : d.scale*(-lp)^(-1/d.shape)
invlogccdf(d::Frechet, lp::Float64) = lp > zero(lp) ? NaN : d.scale*(-log1mexp(lp))^(-1/d.shape)

function gradlogpdf(d::Frechet, x::Float64)
  insupport(Frechet, x) ? -(d.shape + 1.0) / x + d.shape * (d.scale^d.shape) * x^(-d.shape - 1.0)  : 0.0
end

## Sampling
rand(d::Frechet) = d.scale*Base.Random.randmtzig_exprnd()^(-1/d.shape)


text:
using Compose
mutable struct Frame{A, S, T, SA}
    add::A
    show::S
    tree::T
    save::SA
    tag::String
    objects::Array
function Frame(width::Int64, height::Int64, lm, bm)
    base = string("compose(context(units=UnitBox(0,0,",
    width, ",",
    height, ",",
    lm, ",", 0mm,",",
    0mm, ",", bm,
        ")),")
    objects = []
    tag = base
    composition = nothing
    add(object) = composition,objects,tag = _frameup(base,objects,object)
    show() = composition
    tree() = introspect(composition)
    save(name) = draw(SVG(name), composition)
    A, S, T, SA = typeof(add), typeof(show), typeof(tree), typeof(save)
    new{A, S, T, SA}(add, show, tree, save, tag, objects)
end
function _frameup(tag,objects,object)
    push!(objects, object.tag)
    objecttags = ""
    for o in objects
       objecttags = string(objecttags,"(context(),",o, "),")
    end
    tag = string(tag,objecttags,")")
    exp = Meta.parse(tag)
    println(exp)
    composition = eval(exp)
    return(composition,objects,tag)
end
end
function Container(width, height, lm, rm)
    tag = string("context(units=UnitBox(0,0,",
    width, ",",
    height, ",",
    0px, ",", 0px,",",
    0px, ",", 0px,
        ")),")
    objects = []
    composition = nothing
    add(object) = composition,objects,tag = _containerup(tag,objects,object)
    show() = composition
    tree() = introspect(composition)
    save(name) = draw(SVG(name), composition)
    (var) -> (add;show;base;objects;composition;save;tree;width;height;lm;rm;tag)
end
function _containerup(tag,objects,object)
    if typeof(object.tag) == Core.Box
        tagger = Core.getfield(object.tag, :contents)::Any
    else
        tagger = object.tag
    end
    push!(objects, object.tag)
    objecttags = ""
    for o in objects
       objecttags = string(objecttags,"(context(),",o, "),")
    end
    tag = string(tag,objecttags,"")
    println(tag)
    return(tag,objects,tag)
end


text:
using LinearAlgebra

# Vectors are created with '[]'
x = [ -1.1, 0.0, 3.6, -7.2 ]    # using commas
length(x)
z = [ -1.1; 0.0; 3.6; -7.2 ]    # using semicolons
length(z)

# Vector Operations
x[3]                            # 1-based indexing
x[3] = 4.0;                     # assignment
y = x;                          # y references x
y = copy(x);                    # y makes a new copy of x
y == x;                         # check vector equality

a = [ 1.3]
b = 1.3
a == b                          # false, scalars != 1-vector arrays
a[1] == based                   # true, first element same

# Stack Vectors
x = [1, -2]; y = [1,1,0];       # concat stacked vectors with ';'
z = [x;y]
g = [1, -2]; h = [1,1,0];       # using vertical concat
z = vcat(g,h)

#Watch out!
z = (x,y)                       # LIST of 2 vectors
z = [x,y]                       # ARRAY of 2 vectors

# Vector of first differences
m = [ 1, 0, 0, -2, 2 ];
d = m[2:end] - m[1:end-1]       # produces [ -1, 0, -2, 4 ]

max_stars_repo_name: lynnlangit/julia-linear-algebra

max_issues_repo_name: lynnlangit/julia-linear-algebra

max_forks_repo_name: lynnlangit/julia-linear-algebra


text:
"""
theta_w(p[Pa],Temp[K],qv[kg/kg]) from eqn 3.8 of Davies-Jones 2008 and
eqn. 43 of Bolton, 1980, MWR, 108, 1046-1053.
"""
theta_w(p::Number,T::Number,qv::Number) = theta_w(uconvert(hPa,p),uconvert(K,T),NoUnits(uconvert(u"kg/kg",qv))) # wrapper handles diverse units
theta_w(p,T,qv)=theta_w.(p,T,qv) # array broadcasting
function theta_w(p::Pressure,T::Temperature,qv::NondimensionalQuantity) # specific function proper
  # Davies-Jones rational function coefficients
  a0=  7.101574
  a1=-20.68208
  a2= 16.11182
  a3=  2.574631
  a4= -5.205688
  b1= -3.552497
  b2=  3.781782
  b3= -0.6899655
  b4= -0.5929340

  # compute Theta_e from Bolton as first guess approx of thw
  thw = theta_e(p,T,qv);

  # evaluates Davies-Jones 2008 rational function fit for thw>=173.15 [here in K]
  if thw>=173.15
    X=thw/273.15
    X2=X*X
    X3=X2*X
    X4=X2*X2
    thw=thw-exp((a0+a1*X+a2*X2+a3*X3+a4*X4)/(1.0+b1*X+b2*X2+b3*X3+b4*X4));
  end

end

#=
# constants
global Rd Cp
A = 2675;  %K
C = 273.15; %K
lam = Cp/Rd;
ep=0.622
es0=6.112 % hPa
a=17.67;
b=243.5;

% Davies-Jones 2008 method resembles Bolton
if the<=257
   dlnesdlnthe=a*b/((Temp-C+b).*(Temp-C+b));
   thw=the-C-A*rs/(1+A*rs.*dlnesdlnthe));
elseif the
elseif 377<=the & the<674
end

% doesn't work at all
% initialize thw
thw=max(min(-35,Twet(th-273.15,qv,1e5,5)),85)+273.15;
thw=th;

for iter=1:niter
   % es(T) in hPa from es.m Buck
   esat=6.1121*(1.0007 + 3.46e-8*1000).*exp((17.502*(thw-273.15))./(240.97 + thw-273.15));
   %hPa                          %hPa
   rs=622.*esat./(1000-esat); % saturation mixing ratio in hPa, Bolton eqn. 41
   the=thw.*exp( (3.376./thw-0.00254).*rs.*(1+0.81e-3*rs) );
   rath=thetae./the; % ratio of Bolton's thetae to guess thetae
   thw=thw.*rath.^0.287;
end
=#


text:
module Intrinsics

using YaoAPI: AbstractRegister
using YaoHIR
using YaoLocations
using CompilerPluginTools

export main, apply, measure, barrier, expect

# true intrinsics
@intrinsic_stub device apply(r::AbstractRegister, gate::Routine)
@intrinsic_stub device apply(r::AbstractRegister, gate::Routine, ::Locations)
@intrinsic_stub device apply(r::AbstractRegister, gate::Routine, ::Locations, ::CtrlLocations)
@intrinsic_stub device measure(r::AbstractRegister, ::Locations)
@intrinsic_stub device barrier(r::AbstractRegister, ::Locations)
@intrinsic_stub device expect(r::AbstractRegister, ::Locations, nshots::Int)

# this is just for doing overlay since we don't
# have explicit register semantic inside @device
@intrinsic_stub device apply(gate::Routine)
@intrinsic_stub device apply(gate::Routine, ::Locations)
@intrinsic_stub device apply(gate::Routine, ::Locations, ::CtrlLocations)
@intrinsic_stub device measure(::Locations)
@intrinsic_stub device barrier(::Locations)
@intrinsic_stub device expect(::Locations, nshots::Int)

using YaoHIR.IntrinsicOperation

export X, Y, Z, H, S, T, SWAP, shift, Rx, Ry, Rz, UGate

end


text:
# Autogenerated wrapper script for FFMPEG_jll for x86_64-linux-gnu
export libavdevice, libavresample, libswresample, libavutil, libavcodec, ffprobe, libavformat, ffmpeg, libpostproc, libswscale, libavfilter

using libass_jll
using libfdk_aac_jll
using FriBidi_jll
using FreeType2_jll
using LAME_jll
using libvorbis_jll
using Ogg_jll
using LibVPX_jll
using x264_jll
using x265_jll
using Bzip2_jll
using Zlib_jll
using OpenSSL_jll
using Opus_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `libavdevice`
const libavdevice_splitpath = ["lib", "libavdevice.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libavdevice_path = ""

# libavdevice-specific global declaration
# This will be filled out by __init__()
libavdevice_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libavdevice = "libavdevice.so.58"

# Relative path to `libavresample`
const libavresample_splitpath = ["lib", "libavresample.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libavresample_path = ""

# libavresample-specific global declaration
# This will be filled out by __init__()
libavresample_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libavresample = "libavresample.so.4"

# Relative path to `libswresample`
const libswresample_splitpath = ["lib", "libswresample.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libswresample_path = ""

# libswresample-specific global declaration
# This will be filled out by __init__()
libswresample_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libswresample = "libswresample.so.3"

# Relative path to `libavutil`
const libavutil_splitpath = ["lib", "libavutil.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libavutil_path = ""

# libavutil-specific global declaration
# This will be filled out by __init__()
libavutil_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libavutil = "libavutil.so.56"

# Relative path to `libavcodec`
const libavcodec_splitpath = ["lib", "libavcodec.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libavcodec_path = ""

# libavcodec-specific global declaration
# This will be filled out by __init__()
libavcodec_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libavcodec = "libavcodec.so.58"

# Relative path to `ffprobe`
const ffprobe_splitpath = ["bin", "ffprobe"]

# This will be filled out by __init__() for all products, as it must be done at runtime
ffprobe_path = ""

# ffprobe-specific global declaration
function ffprobe(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(ffprobe_path)
    end
end

# Relative path to `libavformat`
const libavformat_splitpath = ["lib", "libavformat.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libavformat_path = ""

# libavformat-specific global declaration
# This will be filled out by __init__()
libavformat_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libavformat = "libavformat.so.58"

# Relative path to `ffmpeg`
const ffmpeg_splitpath = ["bin", "ffmpeg"]

# This will be filled out by __init__() for all products, as it must be done at runtime
ffmpeg_path = ""

# ffmpeg-specific global declaration
function ffmpeg(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(ffmpeg_path)
    end
end

# Relative path to `libpostproc`
const libpostproc_splitpath = ["lib", "libpostproc.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libpostproc_path = ""

# libpostproc-specific global declaration
# This will be filled out by __init__()
libpostproc_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libpostproc = "libpostproc.so.55"

# Relative path to `libswscale`
const libswscale_splitpath = ["lib", "libswscale.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libswscale_path = ""

# libswscale-specific global declaration
# This will be filled out by __init__()
libswscale_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libswscale = "libswscale.so.5"

# Relative path to `libavfilter`
const libavfilter_splitpath = ["lib", "libavfilter.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libavfilter_path = ""

# libavfilter-specific global declaration
# This will be filled out by __init__()
libavfilter_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libavfilter = "libavfilter.so.7"

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"FFMPEG")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (libass_jll.PATH_list, libfdk_aac_jll.PATH_list, FriBidi_jll.PATH_list, FreeType2_jll.PATH_list, LAME_jll.PATH_list, libvorbis_jll.PATH_list, Ogg_jll.PATH_list, LibVPX_jll.PATH_list, x264_jll.PATH_list, x265_jll.PATH_list, Bzip2_jll.PATH_list, Zlib_jll.PATH_list, OpenSSL_jll.PATH_list, Opus_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (libass_jll.LIBPATH_list, libfdk_aac_jll.LIBPATH_list, FriBidi_jll.LIBPATH_list, FreeType2_jll.LIBPATH_list, LAME_jll.LIBPATH_list, libvorbis_jll.LIBPATH_list, Ogg_jll.LIBPATH_list, LibVPX_jll.LIBPATH_list, x264_jll.LIBPATH_list, x265_jll.LIBPATH_list, Bzip2_jll.LIBPATH_list, Zlib_jll.LIBPATH_list, OpenSSL_jll.LIBPATH_list, Opus_jll.LIBPATH_list,))

    global libavdevice_path = normpath(joinpath(artifact_dir, libavdevice_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libavdevice_handle = dlopen(libavdevice_path)
    push!(LIBPATH_list, dirname(libavdevice_path))

    global libavresample_path = normpath(joinpath(artifact_dir, libavresample_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libavresample_handle = dlopen(libavresample_path)
    push!(LIBPATH_list, dirname(libavresample_path))

    global libswresample_path = normpath(joinpath(artifact_dir, libswresample_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libswresample_handle = dlopen(libswresample_path)
    push!(LIBPATH_list, dirname(libswresample_path))

    global libavutil_path = normpath(joinpath(artifact_dir, libavutil_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libavutil_handle = dlopen(libavutil_path)
    push!(LIBPATH_list, dirname(libavutil_path))

    global libavcodec_path = normpath(joinpath(artifact_dir, libavcodec_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libavcodec_handle = dlopen(libavcodec_path)
    push!(LIBPATH_list, dirname(libavcodec_path))

    global ffprobe_path = normpath(joinpath(artifact_dir, ffprobe_splitpath...))

    push!(PATH_list, dirname(ffprobe_path))
    global libavformat_path = normpath(joinpath(artifact_dir, libavformat_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libavformat_handle = dlopen(libavformat_path)
    push!(LIBPATH_list, dirname(libavformat_path))

    global ffmpeg_path = normpath(joinpath(artifact_dir, ffmpeg_splitpath...))

    push!(PATH_list, dirname(ffmpeg_path))
    global libpostproc_path = normpath(joinpath(artifact_dir, libpostproc_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libpostproc_handle = dlopen(libpostproc_path)
    push!(LIBPATH_list, dirname(libpostproc_path))

    global libswscale_path = normpath(joinpath(artifact_dir, libswscale_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libswscale_handle = dlopen(libswscale_path)
    push!(LIBPATH_list, dirname(libswscale_path))

    global libavfilter_path = normpath(joinpath(artifact_dir, libavfilter_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libavfilter_handle = dlopen(libavfilter_path)
    push!(LIBPATH_list, dirname(libavfilter_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
const live_mode = true

import ..debug_header

@static if live_mode

    struct IODebug{T <: IO} <: IO
        io::T
    end

    logwrite(iod::IODebug, x) = show_io_debug(compat_stdout(), "➡ ", x)
    logread(iod::IODebug, x) = show_io_debug(compat_stdout(), "⬅ ", x)
    logunread(iod::IODebug, x) = show_io_debug(compat_stdout(), "♻ ", x)

else

    struct IODebug{T <: IO} <: IO
        io::T
        log::Vector{Tuple{String,String}}
    end

    IODebug(io::T) where T <: IO = IODebug{T}(io, [])

    logwrite(iod::IODebug, x) = push!(iod.log, ("➡ ", x))
    logread(iod::IODebug, x) = push!(iod.log, ("⬅ ", x))
    logunread(iod::IODebug, x) = push!(iod.log, ("♻ ", x))

end

Base.write(iod::IODebug, a...) = (logwrite(iod, join(a)); write(iod.io, a...))

Base.write(iod::IODebug, x::String) = (logwrite(iod, x); write(iod.io, x))

Base.unsafe_write(iod::IODebug, x::Ptr{UInt8}, n::UInt) =
    (logwrite(iod, unsafe_string(x,n));
     unsafe_write(iod.io, x, n))

Base.read(iod::IODebug, n::Integer) =
    (r = read(iod.io, n);
     logread(iod, String(r)); r)

Base.readavailable(iod::IODebug) =
    (r = readavailable(iod.io);
     logread(iod, String(r)); r)

IOExtras.unread!(iod::IODebug, bytes) =
    (logunread(iod, String(bytes));
     unread!(iod.io, bytes))

Base.eof(iod::IODebug) = eof(iod.io)
Base.close(iod::IODebug) = close(iod.io)
Base.isopen(iod::IODebug) = isopen(iod.io)
Base.iswritable(iod::IODebug) = iswritable(iod.io)
Base.isreadable(iod::IODebug) = isreadable(iod.io)
IOExtras.startread(iod::IODebug) = startread(iod.io)
IOExtras.startwrite(iod::IODebug) = startwrite(iod.io)
IOExtras.closeread(iod::IODebug) = closeread(iod.io)
IOExtras.closewrite(iod::IODebug) = closewrite(iod.io)

@static if isdefined(Base, :bytesavailable)
    Base.bytesavailable(iod::IODebug) = bytesavailable(iod.io)
else
    Base.nb_available(iod::IODebug) = nb_available(iod.io)
end

Base.show(io::IO, iod::IODebug) = show(io, iod.io)

function show_log(io::IO, iod::IODebug)
    lock(io)
    println(io, "$(typeof(iod)):\nio:     $(iod.io)")
    prevop = ""
    for (operation, bytes) in iod.log
        if prevop != "" && prevop != operation
            println(io)
        end
        show_io_debug(io, operation, bytes)
        prevop = operation
    end
    println(io)
    unlock(io)
end

function show_io_debug(io::IO, operation, bytes)
    prefix = string(debug_header(), rpad(operation, 4))
    i = j = 1
    while i < length(bytes)
        j = findnext(bytes, '\n', i)
        if j == nothing || j == 0
            j = prevind(bytes, length(bytes)+1)
        end
        println(io, prefix, "\"", escape_string(bytes[i:j]), "\"")
        if i == 1
            prefix = rpad("DEBUG:", length(prefix) - 1)
        end
        i = nextind(bytes, j)
    end
end


text:
include("./sequencingModel.jl")
include("./dpModel.jl")


text:
const ORDERED = Dict(
    :concept => issorted
)


text:
# modifier for converting to gray-scale
const CCV_IO_GRAY      = 0x100
# modifier for converting to color
const CCV_IO_RGB_COLOR = 0x300

# modifier for not copy the data over when read raw in-memory data
const CCV_IO_NO_COPY = 0x10000

# read self-describe in-memory data
const CCV_IO_ANY_STREAM     = 0x010
const CCV_IO_BMP_STREAM     = 0x011
const CCV_IO_JPEG_STREAM    = 0x012
const CCV_IO_PNG_STREAM     = 0x013
const CCV_IO_PLAIN_STREAM   = 0x014
const CCV_IO_DEFLATE_STREAM = 0x015
# read self-describe on-disk data
const CCV_IO_ANY_FILE       = 0x020
const CCV_IO_BMP_FILE       = 0x021
const CCV_IO_JPEG_FILE      = 0x022
const CCV_IO_PNG_FILE       = 0x023
const CCV_IO_BINARY_FILE    = 0x024
# read not-self-describe in-memory data (a.k.a. raw data)
# you need to specify rows cols or scanline for these data
const CCV_IO_ANY_RAW        = 0x040
const CCV_IO_RGB_RAW        = 0x041
const CCV_IO_RGBA_RAW       = 0x042
const CCV_IO_ARGB_RAW       = 0x043
const CCV_IO_BGR_RAW        = 0x044
const CCV_IO_BGRA_RAW       = 0x045
const CCV_IO_ABGR_RAW       = 0x046
const CCV_IO_GRAY_RAW       = 0x047

const CCV_IO_FINAL = 0x00
const CCV_IO_CONTINUE = 0x01
const CCV_IO_ERROR = 0x02
const CCV_IO_ATTEMPTED = 0x03
const CCV_IO_UNKNOWN = 0x04

const CCV_8U  = 0x01000
const CCV_32S = 0x02000
const CCV_32F = 0x04000
const CCV_64S = 0x08000
const CCV_64F = 0x10000

const CCV2JULIA_TYPE =
    Dict(CCV_8U => UInt8,
         CCV_32S => Int32,
         CCV_32F => Float32,
         CCV_64S => Int64,
         CCV_64F => Float64)

const JULIA2CCV_TYPE = Dict(zip(values(CCV2JULIA_TYPE), keys(CCV2JULIA_TYPE)))


text:
abstract type MonteCarloSampler <: Sampler end
abstract type FullSpaceSampler <: Sampler end

parallel_type(s::Sampler) = true

abstract type SamplerCache{T}  end

"""
    caches(sampler, v, net, [parallel_type=NotParallel()])

Creates a `SamplerCache` for this sampler and state.
"""
cache(s::Sampler, hilb::AbstractHilbert, net, par=NotParallel()) =
    _sampler_cache(s, state(hilb, net), hilb, net, par)

cache(s::Sampler, σ, hilb::AbstractHilbert, net, par=NotParallel()) =
    _sampler_cache(s, σ, hilb, net, par)

"""
    init_sampler!(sampler, net, σ, [c=cache(sampler, σ, net)]) -> c

Initializes the sampler `sampler` and state `σ`. If no `SamplerCache` is
provided, one will be initialized and returned. The state σ is the first in the
list of sampled states.
"""
init_sampler!(s::Sampler, net, hilb::AbstractHilbert, σ, par_cache::ParallelType=NotParallel()) =
    init_sampler!(s, net, σ, cache(s, σ, hilb, net, par_cache))

"""
    chain_length(sampler, sampler_cache) -> Int

Returns the estimated length of the chain.
"""
function chain_length end

"""
    sampling_function(s::SamplerCache) -> function

Returns the sampling function. By default it's 2.0*real().
"""
function sampling_function end

samplenext!(σ, sampl::Sampler, net, sampler_cache) =
    samplenext!(σ, σ, sampl, net, sampler_cache)


text:
#Generate GLM Formula
#function genFmula(iv::Array{Symbol},isLastOffset::Bool=false)  # Y=1st, if isLastOffset:: offset last
#    isLastOffset ?  genFmula(iv[1],iv,last(iv)) : genFmula(iv[1],iv)
#end
#function genFmula(y::Symbol, iv::Array{Symbol},offset::Symbol=:Missing)
#    xvars=setdiff(iv,[y,offset])
#    eval(parse( string(y)*" ~ 1"* reduce(*, [ " + "*  string(c) for c in xvars ] )  ))
#end
#function genFmula(y::Symbol, iv::Array{Symbol})
#    xvars=setdiff(iv,[y])
#    eval(parse( string(y)*" ~ 1"* reduce(*, [ " + "*  string(c) for c in xvars ] )  ))
#end
function genF(y::Symbol, iv::Array{Symbol},ranef::Array{Symbol}=Symbol[])
        vars=setdiff(iv,vcat([y],ranef))
        f = string(y)*" ~ 1"
        f = length(vars) > 0 ?  f * reduce(*, [ " + "*  string(c) for c in unique(vars) ] ) : f
        f = length(ranef) > 0 ?  f * reduce(*, [ " + "*  "(1 | "*string(c)*")" for c in unique(ranef) ] )  : f
        eval(parse(f))
end

# remove single level vars
function FS_singleLevel(dfd::DataFrame,vars::Array{Symbol}=Symbol[])
    if length(vars)==0  vars=names(dfd)  end
    return [c for c in filter(x -> length(unique( dfd[x] ))==1, vars)]
end

function raneffect(m::MixedModels.GeneralizedLinearMixedModel)
    DataFrame(parameter=levels(m.LMM.trms[1].f), coef=vec(ranef(m, named=true)'[1]), stderr = sqrt.(vec(condVar(m)[1])) )  # stderrX =vec(condVar(m)[1])
end

isDict(::Dict) = true
#isDict(::OrderedDict) = true
isDict(::Any) = false
isArray(::AbstractArray) = true
isArray(::Any) = false
isNum{T<:Number}(::AbstractArray{T}) = true
isNum(::Any) = false
isFloat{T<:Float64}(::AbstractArray{T}) = true
isFloat(::Any) = false
isInt{T<:Int64}(::AbstractArray{T}) = true
isInt(::Any) = false
isString{T<:String}(::AbstractArray{T}) = true
isString(::String) = true
isString(::Any) = false
isBool{T<:Bool}(::AbstractArray{T}) = true
isBool(::Any) = false
isSymbol{T<:Symbol}(::AbstractArray{T}) = true
isSymbol(::Any) = false

function dfcor(dfr::DataFrame, cols::Vector{Symbol})
    nms = Symbol[]
    arr = Float64[]
    for (n, v) in eachcol(dfr)
        if n ∈ cols && isNum(v)
            push!(nms, n)
            append!(arr, v)
        end
    end
    result = NamedArray(cor(reshape(arr, (size(dfr, 1), length(nms)))))
    NamedArrays.setnames!(result, string.(nms), 1)
    NamedArrays.setnames!(result, string.(nms), 2)
    result
end
dfcor(dfr::DataFrame) = dfcor(dfr, names(dfr))

function dfcor1(dfr::DataFrame, cols::Vector{Symbol})
    m = size(dfr, 1)
    nms = Symbol[]
    arr = ones(m)
    for (n, v) in eachcol(dfr)
        if n ∈ cols && isNum(v)
            push!(nms, n)
            append!(arr, v)
        end
    end
    n = length(nms) + 1
    R = full(UpperTriangular(view(qrfact!(reshape(arr, (m, n)))[:R], 2:n, 2:n)))
    for j in 1:size(R, 2)
        colj = view(R, 1:j, j)
        colj ./= norm(colj)
    end
    Rt = UpperTriangular(R)
    result = NamedArray(Rt'Rt)
    NamedArrays.setnames!(result, string.(nms), 1)
    NamedArrays.setnames!(result, string.(nms), 2)
    result
end
dfcor1(dfr::DataFrame) = dfcor1(dfr, names(dfr))

function blank(dfx::DataFrame,d::Dict=Dict())
    arr = Any[]
    for (n, v) in eachcol(dfx)
        if haskey(d,n)
            arr = vcat(arr, [d[n]])
        else
            if isFloat(v)
                arr = vcat(arr, [0.0])
            elseif isInt(v)
                arr = vcat(arr, [0])
            elseif isString(v)
                arr = vcat(arr, [""])
            elseif isBool(v)
                        arr = vcat(arr, [false])
            elseif isSymbol(v)
                        arr = vcat(arr, [:empty])
            else
                println("ERROR : datatype Not Found!!!!!!!!")
            end
        end
    end
    return  arr
end

Base.lowercase(df::DataFrame) = names!(df, map(x->Symbol(lowercase(string(x))),names(df)))
Base.lowercase(da::DataArray{String}) = [lowercase(x) for x in da]

function filler(dfx::DataFrame, rowIn::Array{Any})
    row=filler(dfx)
    for i in 1:length(rowIn)
        row[i] = rowIn[i]
    end
    return row
end

function filler(dfx::DataFrame, cols::Vector{Symbol}=Symbol[])
    a = Any[]
    if length(cols)==0 cols=names(dfx) end
    for (n, v) in eachcol(dfx)
        if n ∈ cols && isFloat(v)
            push!(a, 0.0)
        elseif n ∈ cols && isInt(v)
            push!(a, 0)
        elseif n ∈ cols && isString(v)
            push!(a, "")
        elseif n ∈ cols && isBool(v)
            push!(a, false)
        end
    end
    return a
end

function cov2cor(m)
    d = Diagonal([inv(sqrt(m[i, i])) for i in 1:size(m, 2)])
    A_mul_B!(d, m * d)
end

function vif2(m)
    v = vcov(m)
    assign = m.mm.assign
    nms = coefnames(m.mf)
    if (ipos = findfirst(nms, "(Intercept)")) > 0
        inds = deleteat!(collect(eachindex(nms)), ipos)
        v = view(v, inds, inds)
        assign = view(assign, inds)
    else
        warning("No intercept: vifs may not be sensible.")
    end
    terms = convert(Vector{Symbol}, m.mf.terms.terms)
    n_terms = length(terms)
    if n_terms < 2
        error("model contains fewer than 2 terms")
    end
    R = cov2cor(v)
    detR = det(Symmetric(R))
    result = NamedArray(Array(eltype(R), (length(terms), 3)))
    setnames!(result, string.(terms), 1)
    setnames!(result, ["GVIF", "Df", "GVIF^(1/(2*Df))"], 2)
    inds = collect(eachindex(terms))
    ainds = collect(eachindex(assign))
    for i in inds
        subs = find(x -> x == i, assign)
        notsubs = setdiff(ainds, subs)
        result[i, 1] = det(Symmetric(R[notsubs, notsubs])) / detR
        if (result[i, 2] = length(subs)) > 1
            result[i, 1] *= det(Symmetric(R[subs, subs]))
        end
    end
    if all(result[:, 2] .== 1)
        return result[:, 1:1]
    end
    result
end

function vifDF2(m::RegressionModel)
    sdf=coefDF(m,true)
    v=vif2(m)
    vdf=DataFrame(parameter=convert(Array{Symbol},names(v)[1]), vif=convert(Array{Float64}, v)[:,1]   )
    #vdf[:parameter] = convert(Array{Symbol}, vdf[:variable])
    join(sdf,vdf[[:parameter,:vif]], on = :parameter, kind=:outer)
end

""" Doug
    profileθ(m::GeneralizedLinearMixedModel)

Profile the deviance as a function of θ for models with a single, scalar r.e. term
"""
function profileθ{T}(m::GeneralizedLinearMixedModel{T})
    θ = getθ(m)
    isa(θ, Vector{T}) && length(θ) == 1 ||
        throw(ArgumentError("m must have a single scalar r.e. term"))
    th = zeros(T, 1)
    δ = zeros(T, 1)
    fill!(θ, 0)
    decreasing = false
    basedev = pirls!(setθ!(m, θ), true)
    for t in 0.001:0.001:2.0
        θ[1] = t
        dd = pirls!(setθ!(m, θ), true) - basedev
        sgn = sign(dd - δ[end])
        push!(δ, dd)
        push!(th, t)
        dd ≤ 4.0 || break
        decreasing && sgn > 0 && break
        decreasing = sgn < 0
    end
    th, δ
end

function gethostworkers()
    d = Dict()
    @sync @async for (idx, pid) in enumerate(workers())
        d[pid] = remotecall_fetch(getipaddr,pid)
    end
    iout=Dict()
    for ip in [string(ip) for ip in unique(values(d))]
        a=Int64[]
        for (key, value) in d
            if string(value) == ip
                push!(a,key)
            end
        end
        iout[ip] = sort(a)
    end
    l = get(iout,string(getipaddr()),[])
    iout[string(getipaddr())] = convert(Array{Int64},vcat([1],l))
    return iout
end

function gethostModelDistribution(raneff::Array{Symbol})
    function allocateCore(a::Array{Int64})
        proc = a[1]
        ao = length(setdiff(a,proc)) > 0 ?  setdiff(a,proc) : a
        return proc, ao
    end
    nodehosts = gethostworkers()
    hosts = [k for k in keys(nodehosts)]
    hostcnt = length(hosts)
    #raneff = [:creative, :publisher1, :placement]
    h=1
    m=:pen
    penprocs = Dict()
    for r in raneff
        proc, nodehosts[hosts[h]] = allocateCore(nodehosts[hosts[h]])
        penprocs[r] = Dict(:proc=>proc, :host=>hosts[h])
        println(r,"  :  ", hosts[h], "  ",proc)
        h = h == hostcnt ? h=1 : h+1
    end

    static_nodehosts = deepcopy(nodehosts)

    println("Pen : ",penprocs )
    println(nodehosts)
end

function checksingularity(form::Formula, data::DataFrame, tolerance = 1.e-8)
    mf = ModelFrame(form, data)
    mm = ModelMatrix(mf)
    qrf = qrfact!(mm.m, Val{true})
    vals = abs.(diag(qrf[:R]))
    firstbad = findfirst(x -> x < min(tolerance, 0.5) * vals[1], vals)
    if firstbad == 0
        return Symbol[]
    end
    mf.terms.terms[view(mm.assign[qrf[:p]], firstbad:length(vals))]
end

function coefDF(m::RegressionModel, paramsAsSymbols::Bool=false)
    vout = DataFrame( #vars=vcat([:intercept],m.mf.terms.terms)  #g.model.mm.assign # m.mm.assign
               #vars=vcat([:intercept],[m.mf.terms.terms[x] for x in m.mm.assign[2:end]])
               parameter=DataFrames.coefnames(m.mf)
             , coef=GLM.coef(m)
             , stderr=stderr(m)
             , zval=GLM.coef(m)./stderr(m)
             , pval= ccdf(FDist(1, dof_residual(m)), abs2(GLM.coef(m)./stderr(m))  )
             )
    if paramsAsSymbols == true
        vout[:parameter] = map(x->Symbol(replace(replace(x,"group: 1","group"),"(Intercept)","Intercept")), vout[:parameter])
    end
    return vout
end

"""
Relevel DF for GLM
"""
function relevel{T,I}(v::PooledDataVector{T,I}, lev::T)
    pool = v.pool
    i = convert(I, findfirst(pool, lev))
    i ≠ 0 || throw(ArgumentError("\"$lev\" is not in v.pool"))
    if i == 1
        return copy(v)
    end
    orig = convert(Vector{I}, 1:length(pool))
    perm = append!([i], deleteat!(orig, i))
    PooledDataArray(DataArrays.RefArray(invperm(perm)[v.refs]), v.pool[perm])
end
#dfd[ranef] = relevel(dfd[ranef], "none")

#GLM REsiduals
function xResiduals(g::DataFrames.DataFrameRegressionModel)
    resp = g.model.rr
    #sign(resp.y - resp.mu) .* sqrt(resp.devresid)    # Original
    #sign(resp.y - resp.mu) .* sqrt(complex(resp.devresid))
    sign(resp.y - resp.mu) .* sqrt(((resp.devresid)+1)-1)
    #sign(resp.y - resp.mu) .* sqrt(abs(resp.devresid))
end

function unpool(dfd::DataFrame)
    v_factors = Symbol[]
    for c in names(dfd)
        if typeof(dfd[c]) in [ PooledDataArray{String,UInt8,1} ]
            println("         Converting : ", c )
            push!(v_factors,c)
            dfd[c] = Array(dfd[c])
        end
    end
    return v_factors
end

function poolit!(dfd::DataFrame, vars::Array{Symbol}=Symbol[])
    if length(vars) == 0   vars=names(dfd) end
    for (i, c) in enumerate(dfd.columns)
        if typeof(c) == DataVector{String}
            if dfd.colindex.names[i] in vars
                #println("[",i,"]  : ",dfd.colindex.names[i])
                #unique(dfd[i])
                dfd[i] = pool(dfd[i])
            end
        end
        #typeof(c) == DataVector{String} && println("[",i,"]  : ",dfd.colindex.names[i])
    end
end

function getColswithType(dtype::String, dfd::DataFrame, vars::Array{Symbol}=Symbol[])
    clst=Symbol[]
    typLst = Dict("num" => [DataVector{Float64},DataVector{Int64}], "str" => [DataVector{String}] )
    if length(vars) == 0   vars=names(dfd) end
    if dtype in ["num","str"]
        for (i, c) in enumerate(dfd.columns)
            if dfd.colindex.names[i] in vars
                    if typeof(c) in typLst[dtype]
                #if typeof(c) in [ DataVector{Float64}, DataVector{Int64} ] # = = DataVector{String}
                    #println("[",i,"]  : ",dfd.colindex.names[i])
                    push!(clst,dfd.colindex.names[i])
                end
             end
         end
    else
        println("ERROR: getColswithType : Invalid dtype!")
    end
    return clst
end

function corDFD(dfd::DataFrame, vars::Array{Symbol}=Symbol[])
    if length(vars) == 0
        vars=names(dfd)
    #else
    #    vars=convert(Array{Symbol},vars)
    end
    cols=getColswithType("num", dfd, convert(Array{Symbol},vars) )
    carr= cor(Array(dfd[cols]))
    tcarr = triu(carr,1)
    d = DataFrame(tcarr)
    names!(d,cols)
    d[:vars] = cols
    return stack(d,cols)
end
#function corDFD(dfd::DataFrame, vars::Array{Any}=[])
#    varsS=Symbol[]
#    if length(vars) > 0 varsS=convert(Array{Symbol},vars) end
#    corrDF(dfd,varsS)
#end

function cmparr(a1::Array{Symbol},a2::Array{Symbol},ignore::Array{Symbol}=Symbol[])
    a1=setdiff(a1,ignore)
    a2=setdiff(a2,ignore)
    a1_only=setdiff(a1,a2)
    a2_only=setdiff(a2,a1)
    l1=length(a1_only)
    l2=length(a2_only)
    println(l1,"--",l2)
    l1>l2 ? DataFrame(a1_only=a1_only,a2_only=vcat(a2_only,fill(:x, l1-l2))) : DataFrame(DataFrame(a1_only=vcat(a1_only,fill(:x, l2-l1)),a2_only=a2_only))
end

function cmparr(a1::Array{Any},a2::Array{Symbol},ignore::Array{Symbol}=Symbol[])
    cmparr(convert(Array{Symbol},a1), a2)
end

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
##### Beginning of file

function _print_welcome_message(a::AbstractDict = ENV)::Nothing
    predictmdfull_version::VersionNumber = version()
    predictmdfull_pkgdir::String = package_directory()
    if !_suppress_welcome_message(a)
        @info(string("This is PredictMDFull, version ",predictmdfull_version,),)
        @info(string("For help, please visit https://predictmd.net",),)
        @debug(string("PredictMDFull package directory: ",predictmdfull_pkgdir,),)
    end
    return nothing
end

function _suppress_welcome_message(a::AbstractDict = ENV)::Bool
    return get(a, "SUPPRESS_PREDICTMDFULL_WELCOME_MESSAGE", "false") == "true"
end

##### End of file


text:
const VertexBuffer = Allocated{Buffer,DeviceMemory}
const IndexBuffer = Allocated{Buffer,DeviceMemory}

"""
Describes data that an object needs to be drawn, but without having a pipeline created yet.
"""
struct ShaderDependencies
    vertex_buffer::VertexBuffer
    index_buffer::Optional{IndexBuffer}
    descriptor_sets::Vector{DescriptorSet}
end

function Vulkan.update_descriptor_sets(device, shader_dependencies::ShaderDependencies, resources)
    update_descriptor_sets(
        device,
        map(Base.Fix1(WriteDescriptorSet, shader_dependencies), resources),
        [],
    )
end

struct DescriptorInfo
    type::DescriptorType
    "1-based index into a descriptor set."
    index::Int
    binding::Int
end

struct Descriptor
    set::DescriptorSet
    info::DescriptorInfo
end


text:
# Communicator

const NCCL_UNIQUE_ID_BYTES = 128
const ncclUniqueId_t = NTuple{NCCL_UNIQUE_ID_BYTES, Cchar}

struct UniqueID
    internal::ncclUniqueId_t

    function UniqueID()
        buf = zeros(Cchar, NCCL_UNIQUE_ID_BYTES)
        @apicall(:ncclGetUniqueId, (Ptr{Cchar},), buf)
        new(Tuple(buf))
    end
end

Base.convert(::Type{ncclUniqueId_t}, id::UniqueID) = id.internal

const ncclComm_t = Ptr{Cvoid}

struct Communicator
    handle::ncclComm_t
end

# creates a new communicator (multi thread/process version)
function Communicator(nranks, comm_id, rank)
    handle_ref = Ref{ncclComm_t}()
    @apicall(:ncclCommInitRank, (Ptr{ncclComm_t}, Cint, ncclUniqueId_t, Cint),
             handle_ref,nranks, comm_id, rank)

    Communicator(handle_ref[])
end

# creates a clique of communicators (single process version)
function Communicator(devices)
    ndev = length(devices)
    comms = Vector{ncclComm_t}(undef, ndev)
    devlist = Cint[dev.ordinal for dev in devices]
    @apicall(:ncclCommInitAll, (Ptr{ncclComm_t}, Cint, Ptr{Cint}), comms, ndev, devlist)
    Communicator.(comms)
end


text:
"""
    unix_socket_server(socket::AbstractString, io_port::Integer, start::Condition)

Starts a bi-directional unix socket server that uses a UNIX-socket `socket`
and communicates with the search server through the TCP port `io_port`.
The server is started once the condition `start` is triggered.
"""
function unix_socket_server(socket::AbstractString,
                            io_port::Integer,
                            start::Condition)
    # Checks
    if issocket(socket)
        rm(socket)
    elseif isempty(socket)
        @error "No socket file specified, cannot create Garamond socket."
    elseif isfile(socket)
        @error "$socket already exists, cannot create Garamond socket."
    elseif isdir(socket)
        @error "$socket is a directory, cannot create Garamond socket."
    else
        socket = abspath(socket)
        _path = strip.(split(socket, "/"))
        directory = join(_path[1:end-1], "/")
        !isdir(directory) && mkpath(directory)
    end

    # Wait for search server to be ready
    wait(start)
    server = listen(socket)
    @info "UNIX-Socket server online @$socket..."

    # Start serving requests
    while true
        connection = accept(server)
        @async while isopen(connection)
            request = readline(connection, keep=true)
            if !isempty(request)
                # Send request to search server and get response
                ssconn = connect(Sockets.localhost, io_port)
                print(ssconn, request)
                response = ifelse(isopen(ssconn), readline(ssconn), "")  # expects a "\n" terminator
                close(ssconn)
                # Return response
                println(connection, response)
            end
        end
    end
    return nothing
end


text:
#=
import Pkg; Pkg.activate(joinpath(@__DIR__, "../../../"))
using CytofDensityEstimation
const CDE = CytofDensityEstimation
=#

function scenarios(n)
  if n == 1
    return CDE.Model.generate_samples(NC=1000, NT=1000, gammaC=0.3, gammaT=0.2,
                                      etaC=[.5, .5, 0], etaT=[.5, .2, .3],
                                      loc=[-1, 1, 3.], scale=[.7, .7, .7],
                                      df=[15, 30, 10], skew=[-20, -5, 0.])
  elseif n == 2
    return CDE.Model.generate_samples(NC=1000, NT=1000, gammaC=0.3, gammaT=0.3,
                                      etaC=[.5, .5, 0], etaT=[.5, .2, .3],
                                      loc=[-1, 1, 3.], scale=[.7, .7, .7],
                                      df=[15, 30, 10], skew=[-20, -5, 0.])
  elseif n == 3
    return CDE.Model.generate_samples(NC=1000, NT=1000, gammaC=0.3, gammaT=0.3,
                                      etaC=[.5, .5, 0], etaT=[.5, .5, .0],
                                      loc=[-1, 1, 3.], scale=[.7, .7, .7],
                                      df=[15, 30, 10], skew=[-20, -5, 0.])
  end
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
export start, next, done, length, haskey, keys, values, getindex

## Overloads for transparent iteration on `Exponent`, `Polyomial`, `Point`

## Exponent iterator
start(expo::Exponent) = start(expo.expo)
next(expo::Exponent, state) = next(expo.expo, state)
done(expo::Exponent, state) = done(expo.expo, state)
length(expo::Exponent) = length(expo.expo)
haskey(expo::Exponent, key) = haskey(expo.expo, key)
keys(expo::Exponent) = keys(expo.expo)
values(expo::Exponent) = values(expo.expo)
getindex(expo::Exponent, var::Variable) = expo.expo[var]

## Polynomial
start(poly::Polynomial) = start(poly.poly)
next(poly::Polynomial, state) = next(poly.poly, state)
done(poly::Polynomial, state) = done(poly.poly, state)
length(poly::Polynomial) = length(poly.poly)
haskey(poly::Polynomial, key) = haskey(poly.poly, key)
keys(poly::Polynomial) = keys(poly.poly)
values(poly::Polynomial) = values(poly.poly)
getindex(poly::Polynomial, expo::Exponent) = poly.poly[expo]

## Point
start(pt::Point) = start(pt.coords)
next(pt::Point, state) = next(pt.coords, state)
done(pt::Point, state) = done(pt.coords, state)
length(pt::Point) = length(pt.coords)
haskey(pt::Point, key) = haskey(pt.coords, key)
keys(pt::Point) = keys(pt.coords)
values(pt::Point) = values(pt.coords)
getindex(pt::Point, var::Variable) = pt.coords[var]


text:
#using Revise, IndividualDisplacements
#include(joinpath(dirname(pathof(IndividualDisplacements)),"examples.jl"))

using MeshArrays, DifferentialEquations, Plots

"""
    ex_1()

Global ocean case -- just reading from file for now.

```
df=IndividualDisplacements.ex_1()

p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"PlotIndDisp.jl"))
PyPlot.figure(); PlotMapProj(df,300); gcf()
```
"""
function ex_1()
   dirIn="run_offflt/"
   prec=Float32
   df=IndividualDisplacements.ReadDisplacements(dirIn,prec)
end

"""
    ex_2()

Reproducing `MITgcm/verification/flt_example/` case. This is based on an
extended and modified configuration of the standard MITgcm test case.

```
(df,pl)=IndividualDisplacements.ex_2(); display(pl)

p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"PlotIndDisp.jl"))
PyPlot.figure(); PlotBasic(df,300,100000.0); gcf()
```
"""
function ex_2()
   dirIn="flt_example/"
   prec=Float32
   df=IndividualDisplacements.ReadDisplacements(dirIn,prec)
   uvetc=IndividualDisplacements.ReadGriddedFields()
   #
   tmp=df[df.ID .== 200, :]
   nSteps=Int32(tmp[end,:time]/3600)-2
   ref=transpose([tmp[1:nSteps,:lon] tmp[1:nSteps,:lat]])
   maxLon=80*5.e3
   maxLat=42*5.e3
   for i=1:nSteps-1
       ref[1,i+1]-ref[1,i]>maxLon/2 ? ref[1,i+1:end]-=fill(maxLon,(nSteps-i)) : nothing
       ref[1,i+1]-ref[1,i]<-maxLon/2 ? ref[1,i+1:end]+=fill(maxLon,(nSteps-i)) : nothing
       ref[2,i+1]-ref[2,i]>maxLat/2 ? ref[2,i+1:end]-=fill(maxLat,(nSteps-i)) : nothing
       ref[2,i+1]-ref[2,i]<-maxLat/2 ? ref[2,i+1:end]+=fill(maxLat,(nSteps-i)) : nothing
   end
   ref=ref./uvetc["dx"]
   #
   comp_vel=IndividualDisplacements.VelComp
   get_vel=IndividualDisplacements.VelCopy
   uInit=[tmp[1,:lon];tmp[1,:lat]]./uvetc["dx"]
   du=fill(0.0,2)
   #
   tspan = (0.0,nSteps*3600.0)
   #prob = ODEProblem(get_vel,uInit,tspan,tmp)
   prob = ODEProblem(comp_vel,uInit,tspan,uvetc)
   sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
   #
   Plots.plot(sol[1,:],sol[2,:],linewidth=5,lc=:black, title="One Trajectory Example",
   xaxis="x",yaxis="y",label="Julia Solution") # legend=false
   pl=Plots.plot!(ref[1,:],ref[2,:],lw=3,ls=:dash,lc=:red,label="MITgcm Solution")
   #
   return df,pl
end


text:
s = open("input.txt") do file
    readlines(file)
end

function GetID(seatCode)
    charArr = seatCode |> collect;
    binString = String(map(x->ifelse(x=='F','0','1'),charArr[1:7]));

    row = parse(Int, binString, base=2);

    binString = String(map(x->ifelse(x=='L','0','1'),charArr[8:10]));
    col = parse(Int, binString, base=2);

    return row * 8 + col;
end

foundSeats = zeros(930);

for seat in s
    foundSeats[GetID(seat)] = 1;
end

for i = 2:929
    if foundSeats[i] == 0 && foundSeats[i-1] == 1 && foundSeats[i+1] == 1
        println(i)
    end
end


text:
cmd(command) = ngSpice_Command(command)

function init()
    pvoid = convert(Ptr{Nothing}, 0)
    ngSpice_Init(gen_psendchar(), gen_psendstat(),
        gen_pcontrolledexit(),
        gen_psenddata(),
        gen_psendinitdata(),
        gen_pbgthread(), pvoid)
    ngSpice_Init_JLExtensions(
        @cfunction(path_resolve, Cstring, (Cstring, Cstring))
    )

end

isrunning() = ngSpice_running() # always returns 0

run()     = (println("Running the simulator"); cmd("run"))
bgrun()   = (println("Running the simulator in a background thread");
                cmd("bg_run"))

set_breakpoint(bkpt::Float64) = ngSpice_SetBkpt(bkpt)

stop()    = (println("Stopping the simulator"); cmd("stop"))
bghalt()  = (println("Halting the simulator in a background thread");
                cmd("bg_halt"))

reset()   = (println("Resetting the simulator"); cmd("reset"))
resume()  = (println("Resuming the simulator"); cmd("resume"))

quit()    = cmd("quit")
exit()    = (println("Quitting immediately"); cmd("unset askquit"); cmd("quit"))

max_stars_repo_name: JuliaComputing/Ngspice

max_issues_repo_name: JuliaComputing/Ngspice

max_forks_repo_name: JuliaComputing/Ngspice


text:
"""
Phone data

# Components
- `year::Integer`: years from 1950 to 1973.
- `calls::Float64`: phone calls (in millions).

# Reference
P. J. Rousseeuw and A. M. Leroy (1987) _Robust Regression &
     Outlier Detection._ Wiley.
"""
const phones = DataFrame(
 year=collect(50:73),
 calls=[4.4, 4.7, 4.7, 5.9, 6.6, 7.3, 8.1, 8.8, 10.6, 12.0, 13.5, 14.9, 16.1, 21.2, 119.0, 124.0,
          142.0, 159.0, 182.0, 212.0, 43.0, 24.0, 27.0, 29.0]
)

"""
Hawkins & Bradu & Kass data

# Components
- `x1::Float64`: first independent variable.
- `x2::Float64`: second independent variable.
- `x3::Float64`: third independent variable.
- `y::Float64`: dependent (response) variable.

# Reference
Hawkins, D.M., Bradu, D., and Kass, G.V. (1984) Location of several outliers in multiple regression data using elemental sets. Technometrics 26, 197–208.
"""
const hbk = DataFrame(
    x1=[10.1, 9.5, 10.7, 9.9, 10.3, 10.8, 10.5, 9.9, 9.7, 9.3, 11, 12, 12, 11, 3.4, 3.1, 0, 2.3, 0.8, 3.1, 2.6, 0.4, 2, 1.3, 1, 0.9, 3.3, 1.8, 1.2, 1.2, 3.1, 0.5, 1.5, 0.4, 3.1, 1.1, 0.1, 1.5, 2.1, 0.5, 3.4, 0.3, 0.1, 1.8, 1.9, 1.8, 3, 3.1, 3.1, 2.1, 2.3, 3.3, 0.3, 1.1, 0.5, 1.8, 1.8, 2.4, 1.6, 0.3, 0.4, 0.9, 1.1, 2.8, 2, 0.2, 1.6, 0.1, 2, 1, 2.2, 0.6, 0.3, 0, 0.3],
    x2=[19.6, 20.5, 20.2, 21.5, 21.1, 20.4, 20.9, 19.6, 20.7, 19.7, 24, 23, 26, 34, 2.9, 2.2, 1.6, 1.6, 2.9, 3.4, 2.2, 3.2, 2.3, 2.3, 0, 3.3, 2.5, 0.8, 0.9, 0.7, 1.4, 2.4, 3.1, 0, 2.4, 2.2, 3, 1.2, 0, 2, 1.6, 1, 3.3, 0.5, 0.1, 0.5, 0.1, 1.6, 2.5, 2.8, 1.5, 0.6, 0.4, 3, 2.4, 3.2, 0.7, 3.4, 2.1, 1.5, 3.4, 0.1, 2.7, 3, 0.7, 1.8, 2, 0, 0.6, 2.2, 2.5, 2, 1.7, 2.2, 0.4],
    x3=[28.3, 28.9, 31, 31.7, 31.1, 29.2, 29.1, 28.8, 31, 30.3, 35, 37, 34, 34, 2.1, 0.3, 0.2, 2, 1.6, 2.2, 1.9, 1.9, 0.8, 0.5, 0.4, 2.5, 2.9, 2, 0.8, 3.4, 1, 0.3, 1.5, 0.7, 3, 2.7, 2.6, 0.2, 1.2, 1.2, 2.9, 2.7, 0.9, 3.2, 0.6, 3, 0.8, 3, 1.9, 2.9, 0.4, 1.2, 3.3, 0.3, 0.9, 0.9, 0.7, 1.5, 3, 3.3, 3, 0.3, 0.2, 2.9, 2.7, 0.8, 1.2, 1.1, 0.3, 2.9, 2.3, 1.5, 2.2, 1.6, 2.6],
    y=[9.7, 10.1, 10.3, 9.5, 10, 10, 10.8, 10.3, 9.6, 9.9, -0.2, -0.4, 0.7, 0.1, -0.4, 0.6, -0.2, 0, 0.1, 0.4, 0.9, 0.3, -0.8, 0.7, -0.3, -0.8, -0.7, 0.3, 0.3, -0.3, 0, -0.4, -0.6, -0.7, 0.3, -1, -0.6, 0.9, -0.7, -0.5, -0.1, -0.7, 0.6, -0.7, -0.5, -0.4, -0.9, 0.1, 0.9, -0.4, 0.7, -0.5, 0.7, 0.7, 0, 0.1, 0.7, -0.1, -0.3, -0.9, -0.3, 0.6, -0.3, -0.5, 0.6, -0.9, -0.7, 0.6, 0.2, 0.7, 0.2, -0.2, 0.4, -0.9, 0.2]
)

"""
Animals data

# Components
- `names::AbstractString`: names of animals.
- `body::Float64`: body weight in kg.
- `brain::Float64`: brain weight in g.

# References
     Venables, W. N. and Ripley, B. D. (1999) _Modern Applied
     Statistics with S-PLUS._ Third Edition. Springer.

     P. J. Rousseeuw and A. M. Leroy (1987) _Robust Regression and
     Outlier Detection._ Wiley, p. 57.
"""
const animals = DataFrame(
    names=["Mountain beaver", "Cow", "Grey wolf", "Goat", "Guinea pig", "Dipliodocus", "Asian elephant", "Donkey", "Horse", "Potar monkey", "Cat", "Giraffe", "Gorilla", "Human", "African elephant", "Triceratops", "Rhesus monkey", "Kangaroo", "Golden hamster", "Mouse", "Rabbit", "Sheep", "Jaguar", "Chimpanzee", "Rat", "Brachiosaurus", "Mole", "Pig"],
    body=[1.35, 465, 36.33, 27.66, 1.04, 11700, 2547, 187.1, 521, 10, 3.3, 529, 207, 62, 6654, 9400, 6.8, 35, 0.12, 0.023, 2.5, 55.5, 100, 52.16, 0.28, 87000, 0.122, 192],
    brain=[8.1, 423, 119.5, 115, 5.5, 50, 4603, 419, 655, 115, 25.6, 680, 406, 1320, 5712, 70, 179, 56, 1, 0.4, 12.1, 175, 157, 440, 1.9, 154.5, 3, 180]
)

"""
Weight loss data

# Components
- `days::Integer`: time in days since the start of the diet program.
- `weight::Float64`: weight in kg.

# Reference
     Venables, W. N. and Ripley, B. D. (1999) _Modern Applied
     Statistics with S-PLUS._ Third Edition. Springer.
"""
const weightloss = DataFrame(
    days=[0, 4, 7, 7, 11, 18, 24, 30, 32, 43, 46, 60, 64, 70, 71, 71, 73, 74, 84, 88, 95, 102, 106, 109, 115, 122, 133, 137, 140, 143, 147, 148, 149, 150, 153, 156, 161, 164, 165, 165, 170, 176, 179, 198, 214, 218, 221, 225, 233, 238, 241, 246],
    weight=[184.35, 182.51, 180.45, 179.91, 177.91, 175.81, 173.11, 170.06, 169.31, 165.1, 163.11, 158.3, 155.8, 154.31, 153.86, 154.2, 152.2, 152.8, 150.3, 147.8, 146.1, 145.6, 142.5, 142.3, 139.4, 137.9, 133.7, 133.7, 133.3, 131.2, 133, 132.2, 130.8, 131.3, 129, 127.9, 126.9, 127.7, 129.5, 128.4, 125.4, 124.9, 124.9, 118.2, 118.2, 115.3, 115.7, 116, 115.5, 112.6, 114, 112.6]
)

"""
Stack loss data

# Components
- `airflow::Float64`: flow of cooling air (independent variable).
- `watertemp::Float64`: cooling water inlet temperature (independent variable).
- `acidcond::Float64`: concentration of acid (independent variable).
- `stackloss::Float64`: stack loss (dependent variable).

# Outliers
    Observations 1, 3, 4, and 21 are outliers.

# References
    Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) _The New S Language_.  Wadsworth & Brooks/Cole.

    Dodge, Y. (1996) The guinea pig of multiple regression. In: _Robust Statistics, Data Analysis, and Computer Intensive Methods;
    In Honor of Peter Huber's 60th Birthday_, 1996, _Lecture Notes in Statistics_ *109*, Springer-Verlag, New York.
"""
const stackloss = DataFrame(
    airflow=[80.0, 80, 75, 62, 62, 62, 62, 62, 58, 58, 58, 58, 58, 58, 50, 50, 50, 50, 50, 56, 70],
    watertemp=[27.0, 27, 25, 24, 22, 23, 24, 24, 23, 18, 18, 17, 18, 19, 18, 18, 19, 19, 20, 20, 20],
    acidcond=[89.0, 88, 90, 87, 87, 87, 93, 93, 87, 80, 89, 88, 82, 93, 89, 86, 72, 79, 80, 82, 91],
    stackloss=[42.0, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15]
)

"""
Hadi & Simonoff (1993) Random data

# Components
- `x1::Float64`: Random values.
- `x2::Float64`: Random values.
- `y::Float64`: Random values (independent variable).

# Outliers
    Observations 1, 2, and 3 are outliers.

# References
Hadi, Ali S., and Jeffrey S. Simonoff. "Procedures for the identification of
multiple outliers in linear models." Journal of the American Statistical
Association 88.424 (1993): 1264-1272.
"""
const hs93randomdata = DataFrame(
    x1=[15.0, 14.95, 14.9, 13.7909, 3.60424, 6.3905, 5.13756, 6.98242, 8.54968, 9.41802, 8.58942, 2.32221, 14.3409, 8.08558, 9.38206, 4.32664, 12.7987, 7.92987, 2.26946, 12.8386, 3.66188, 0.254591, 3.91701, 13.1687, 0.953868],
    x2=[15.0, 14.95, 14.9, 17.0999, 3.00793, 11.0249, 9.11747, 6.70008, 6.58946, 7.14941, 3.39223, -0.115934, 10.935, 5.28895, 2.28385, 3.50448, 15.7696, 13.8133, 6.30972, 5.66232, 3.55028, 0.230256, 11.5185, 9.82886, 4.29971],
    y=[34.0, 33.9, 33.8, 32.6016, 6.38553, 16.1119, 14.8969, 14.0065, 15.1297, 17.9406, 10.22806, 1.6473, 25.1874, 13.764, 11.9967, 8.2127, 28.0132, 22.8697, 7.15197, 18.2589, 7.09652, 0.17605, 15.8921, 22.4458, 4.65196]
)

"""
Modified Wood Gravity Data

# Components
- `x1::Float64`: Random values.
- `x2::Float64`: Random values.
- `x3::Float64`: Random values.
- `x4::Float64`: Random values.
- `x5::Float64`: Random values.
- `y::Float64`: Random values (independent variable).

# References
P. J. Rousseeuw and A. M. Leroy (1987) Robust Regression and Outlier Detection. Wiley, p.243, table 8.
"""
const woodgravity = DataFrame(
    x1=[0.573, 0.651, 0.606, 0.437, 0.547, 0.444, 0.489, 0.413, 0.536, 0.685, 0.664, 0.703, 0.653, 0.586, 0.534, 0.523, 0.58, 0.448, 0.417, 0.528],
    x2=[0.1059, 0.1356, 0.1273, 0.1591, 0.1135, 0.1628, 0.1231, 0.1673, 0.1182, 0.1564, 0.1588, 0.1335, 0.1395, 0.1114, 0.1143, 0.132, 0.1249, 0.1028, 0.1687, 0.1057],
    x3=[0.465, 0.527, 0.494, 0.446, 0.531, 0.429, 0.562, 0.418, 0.592, 0.631, 0.506, 0.519, 0.625, 0.505, 0.521, 0.505, 0.546, 0.522, 0.405, 0.424],
    x4=[0.538, 0.545, 0.521, 0.423, 0.519, 0.411, 0.455, 0.43, 0.464, 0.564, 0.481, 0.484, 0.519, 0.565, 0.57, 0.612, 0.608, 0.534, 0.415, 0.566],
    x5=[0.841, 0.887, 0.92, 0.992, 0.915, 0.984, 0.824, 0.978, 0.854, 0.914, 0.867, 0.812, 0.892, 0.889, 0.889, 0.919, 0.954, 0.918, 0.981, 0.909],
    y=[0.534, 0.535, 0.57, 0.45, 0.548, 0.431, 0.481, 0.423, 0.475, 0.486, 0.554, 0.519, 0.492, 0.517, 0.502, 0.508, 0.52, 0.506, 0.401, 0.568]
)

"""
Scottish Hill Races Data

# Components
- `dist::Array{Float64, 1}`: Distance in miles (Independent).
- `climb::Array{Float64, 1}`: Heights in feet (Independent).
- `time::Array{Float64, 1}`: Record times in hours (Dependent).

# Model
time ~ dist + climb

# References
A.C. Atkinson (1986) Comment: Aspects of diagnostic regression analysis. Statistical Science 1, 397-402.
"""
const hills = DataFrame(
    dist=[2.4, 6.0, 6.0, 7.5, 8.0, 8.0, 16.0, 6.0, 5.0, 6.0, 28.0, 5.0, 9.5, 6.0, 4.5, 10.0, 14.0, 3.0, 4.5, 5.5, 3.0, 3.5, 6.0, 2.0, 3.0, 4.0, 6.0, 5.0, 6.5, 5.0, 10.0, 6.0, 18.0, 4.5, 20.0],
    climb=[650.0, 2500.0, 900.0, 800.0, 3070.0, 2866.0, 7500.0, 800.0, 800.0, 650.0, 2100.0, 2000.0, 2200.0, 500.0, 1500.0, 3000.0, 2200.0, 350.0, 1000.0, 600.0, 300.0, 1500.0, 2200.0, 900.0, 600.0, 2000.0, 800.0, 950.0, 1750.0, 500.0, 4400.0, 600.0, 5200.0, 850.0, 5000.0],
    time=[0.268055555555556, 0.805833333333333, 0.560833333333333, 0.76, 1.03777777777778, 1.22027777777778, 3.41027777777778, 0.606111111111111, 0.495833333333333, 0.6625, 3.21111111111111, 0.7175, 1.08333333333333, 0.735555555555556, 0.448888888888889, 1.20416666666667, 1.64027777777778, 1.31083333333333, 0.290277777777778, 0.542777777777778, 0.265833333333333, 0.465, 0.794166666666667, 0.298888888888889, 0.311388888888889, 0.436944444444444, 0.573888888888889, 0.476111111111111, 0.841666666666667, 0.349166666666667, 1.42638888888889, 0.539722222222222, 2.8375, 0.468333333333333, 2.66388888888889]
)

"""
    Soft Drink Delivery Data

# Components
- `cases::Array{Float64, 1}`: Independent variable.
- `distance::Array{Float64, 1}`: Independent variable.
- `time::Array{Float64, 1}`: Dependent variable.

# Model
time ~ distance + cases

# Reference
D. C. Montgomery and E. A. Peck (1992) Introduction to Regression Analysis. Wiley, New York.
"""
const softdrinkdelivery = DataFrame(
    time=[16.68, 11.5, 12.03, 14.88, 13.75, 18.11, 8.0, 17.83, 79.24, 21.5, 40.33, 21.0, 13.5, 19.75, 24.0, 29.0, 15.35, 19.0, 9.5, 35.1, 17.9, 52.32, 18.75, 19.83, 10.75],
    distance=[560.0, 220.0, 340.0, 80.0, 150.0, 330.0, 110.0, 210.0, 1460.0, 605.0, 688.0, 215.0, 255.0, 462.0, 448.0, 776.0, 200.0, 132.0, 36.0, 770.0, 140.0, 810.0, 450.0, 635.0, 150.0],
    cases=[7.0, 3.0, 3.0, 4.0, 6.0, 7.0, 2.0, 7.0, 30.0, 5.0, 16.0, 10.0, 4.0, 6.0, 9.0, 10.0, 6.0, 7.0, 3.0, 17.0, 10.0, 26.0, 9.0, 8.0, 4.0]
)

max_stars_repo_name: angeris/LinRegOutliers

max_issues_repo_name: angeris/LinRegOutliers

max_forks_repo_name: angeris/LinRegOutliers


text:
@testset "Index" begin
    buf = UInt8[]
    writer = BSONWriter(buf)
    writer["x"] = Int64(1)
    writer["y"] = w -> begin
        w["x"] = Int64(2)
    end
    close(writer)
    reader = BSONReader(buf)
    index = BSONIndex(128)
    index[] = reader
    entry = index[LightBSON.BSONIndexKey("x", 0)]
    @test entry !== nothing
    @test entry.type == BSON_TYPE_INT64
    @test entry.offset == 7
    entry = index[LightBSON.BSONIndexKey("y", 0)]
    @test entry !== nothing
    @test entry.type == BSON_TYPE_DOCUMENT
    @test entry.offset == 18
    entry = index[LightBSON.BSONIndexKey("x", 18)]
    @test entry !== nothing
    @test entry.type == BSON_TYPE_INT64
    @test entry.offset == 25
end


text:
using StochasticDiffEq, Zygote
using DiffEqSensitivity, Test

function dt!(du, u, p, t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end

function dW!(du, u, p, t)
  du[1] = 0.1u[1]
  du[2] = 0.1u[2]
end

u0 = [1.0,1.0]
tspan = (0.0, 10.0)
p = [2.2, 1.0, 2.0, 0.4]
prob_sde = SDEProblem(dt!, dW!, u0, tspan,p)

condition(u,t,integrator) = integrator.t >9.0 #some condition
function affect!(integrator)
	 #println("Callback")  #some callback
end
cb = DiscreteCallback(condition,affect!,save_positions=(false,false))

function predict_sde(p)
  return Array(solve(prob_sde, EM(), p=p, saveat = 0.1,sensealg = ForwardDiffSensitivity(), dt=0.001, callback=cb))
end

loss_sde(p)= sum(abs2, x-1 for x in predict_sde(p))

loss_sde(p)
@time dp = gradient(p) do p
	loss_sde(p)
end

@test !iszero(dp[1])


text:
module MacroTest

using Test
using ExprTools: splitdef, combinedef
using Pretend: Pretend, apply, @mockable

# Define a macro that also does splitdef/combinedef
macro before(ex)
    def = splitdef(ex)
    pushfirst!(def[:body].args, :(print("Entry: ")))
    return esc(combinedef(def))
end

function test()
    Pretend.activate()
    @testset "Inner Macro" begin
        # Set up mock for function that is defined with another macro
        @mockable @before hello() = 1

        # Should be patched properly
        apply(hello => () -> 2) do
            @test hello() == 2
        end
    end
end

end # module

using .MacroTest: test
MacroTest.test()


text:
include("PRM.jl")
#using Distributions
pfont = Plots.Font("sans-serif", 8, :left, :left, 0.0, RGB{U8}(0.0, 0.0, 0.0))
glvisualize(titlefont = pfont, guidefont = pfont)
#Plots.resetfontsizes()
#Plots.scalefontsizes(0.9)

####################################
# Randomly Generate Clutter (obstacles)
####################################

# implement function to create clutter
# Let's say on average we want to create 4 obstacles... (otherwise, *most* of
# the time we will product one large rectangle that is falling out of the room

#targetNumObs = 2
function genClutter(roomperimeter, param, targetNumObs, clutterPercentage, roomWidth, roomHeight)
    obstacles = Vector{HyperRectangle}()
    roomArea = roomWidth * roomHeight
    targetSumObsArea= roomArea * clutterPercentage
    sumObsArea = 0
    while sumObsArea < targetSumObsArea
        x,y = rand(Uniform(0, roomWidth),2)

        randWidth, randHeight = rand(Uniform(1, roomWidth/targetNumObs),2)
        #randWidth, randHeight = rand(1:roomWidth/targetNumObs,2)
        protoObstacle = HyperRectangle(Vec(x, y), Vec(randWidth, randHeight)) #Todo
        if contains(roomperimeter, protoObstacle)
            push!(obstacles, protoObstacle)
            sumObsArea += randWidth*randHeight
        end
    end
    actualSumObsArea = sumObsArea
    #print("obstacles generated")
    return actualSumObsArea, obstacles
end
####

function testGenClutter()

    ####################################
    ## PARAMETERS
    ####################################
    numPts = 100
    connectRadius = 5
    param = algT.AlgParameters(numPts, connectRadius)

    targetNumObs = 3
    clutterPercentage = 0.1
    roomWidth,roomHeight  = 20,20

    goalstate = Point(19.,19)

    ####################################
    ## INIT
    ####################################
    obstacles = Vector{HyperRectangle}()
    perimeter = HyperRectangle(Vec(0.,0), Vec(roomWidth,roomHeight))
    wallsPerimeter = algfxn.decompRect(perimeter)

    walls = Vector{LineSegment}()
    for l in wallsPerimeter
        push!(walls, l)
    end
    #plot!(walls, color =:black)
    #plot!(obstacles, fillalpha=0.5)

####################################
# Run PRM once and display config space plot
####################################

    obsArea, obstacles = genClutter(perimeter, param, targetNumObs, clutterPercentage, roomWidth, roomHeight)
    r = algT.Room(roomWidth,roomHeight,walls,obstacles)
    roomPlot = plotfxn.plotRoom(r)

    ## Run preprocessing
    nodeslist, edgeslist = preprocessPRM(r, param)
    print("\n ---- pre-processed --- \n")

    ## Query created RM

    pathcost, isPathFound, solPath = queryPRM(startstate, goalstate, nodeslist, edgeslist, obstacles)
    print("\n ---- queried ---- \n")

    ## Plot path found
    timestamp = Base.Dates.now()
    title = "PRM with # Pts =$numPts, maxDist=$connectRadius, \npathcost = $pathcost,
            timestamp=$timestamp)\n\n"
    roadmap = algT.roadmap(startstate, goalstate, nodeslist, edgeslist)

    #prmPlot= plotfxn.plotPRM(roomPlot, roadmap, solPath, title::String, afont)
    prmPlot= plotfxn.plotPRM(roomPlot, roadmap, solPath, title::String)

    print("\n --- Time --- \n")
    @show timestamp
    #print("\n --- Obstacles --- \n")
    #@show obstacles
    #print("\n --- Nodes --- \n")
    #@show nodeslist
    #print("\n --- Edges --- \n")
    #@show edgeslist
    # print("\n -------- \n")

    #f = Plots.Font(05,"Liberation Serif")
    #plot!(prmPlot,titlefont = f)

    gui(prmPlot)
        print("\n --- Press ENTER to continue --- \n")
        #readline(STDIN)
        # todo: save one "representative" figure from each run
        # print("\n --- continue --- \n")
    @show  isPathFound

end

function clutterExp()
    ####################################
    # Run multiple trials and scatterplot all cost vs area for all runs
    # due to varying absolute area
    ####################################
    ####################################
    ## PARAMETERS - RRT
    ####################################
    connectRadius = 8.
    numPts = 100.
    param = algT.AlgParameters(numPts, connectRadius)
    targetNumObs = 5
    #clutterPercentage = 0.15
    roomWidth,roomHeight  = 20,20

    startstate = Point(1.,1)
    goalstate = Point(19.,19)

    ####################################
    ## PARAMETERS - plot clutter
    ####################################
    N = 50

    clutterPercentages = [0 0.01 0.05 0.1 0.12 0.15 0.2 0.25 0.3 0.4 0.5]
    clutterPercentages = [0 0.1 0.2 0.3 0.5]
    #clutterPercentages = 0 : 0.05 :0.3

    ####################################
    ## Init
    ####################################
    listCosts = Vector{Float32}()
    listpSucc= Vector{Float32}()

    obstacles = Vector{HyperRectangle}()
    perimeter = HyperRectangle(Vec(0.,0), Vec(roomWidth,roomHeight))
    wallsPerimeter = algfxn.decompRect(perimeter)

    walls = Vector{LineSegment}()
    for l in wallsPerimeter
        push!(walls, l)
    end

    ####################################
    # Run Trials
    ####################################

    plot()
    tic()

    stddevs = Vector{Float64}()
    avgCosts = Vector{Float32}() #one per %clutter

    print("\n ------ \n")
    for pClutter in clutterPercentages
        idx, nSuccess, pathcost = 0,0,0.
        pathcosts = Vector{Float32}() #one per %clutter
        stddev_n = 0.

        param = algT.AlgParameters(numPts, connectRadius)

        while idx < N
            idx += 1

            obstacles = Vector{HyperRectangle}()
            obsArea, obstacles = genClutter(perimeter, param, targetNumObs, pClutter, roomWidth, roomHeight)
            r = algT.Room(roomWidth,roomHeight,walls,obstacles)
            nodeslist, edgeslist = preprocessPRM(r, param)
            pathcost, isPathFound, solPath = queryPRM(startstate, goalstate, nodeslist, edgeslist, obstacles)

            if isPathFound
                nSuccess += 1
            end
            if !(pathcost == Void)
                push!(pathcosts, pathcost)
            end
        end

        totalcost = sum(pathcosts)
        avgCost = totalcost /  nSuccess # average cost across the nTrials
        #@show nSuccess
        pSucc = nSuccess / N

        sumSqdError = sum([cost - avgCost for cost in pathcosts].^2)
        #@show sumSqdError
        #@show pathcosts
        #@show avgCost
        stddev_n = sqrt( sumSqdError / (nSuccess -1))
        #@show stddev_n

        #@printf("For the iter of %d the avg cost was %d across %d trials", maxIter, avgCost, nTrials)
        push!(avgCosts, avgCost) #todo combine pt1
        push!(stddevs, stddev_n) #todo combine pt2
        push!(listpSucc, pSucc)
        print("\n ------ \n")
    end

    timeExperiment = toc();
    print("Time to run experiment: $timeExperiment\n")
    timestamp = Base.Dates.now()

    ####################################
    # Plot Results
    ####################################
    # clutter (area) on X
    # pathcost on Y
    sizeplot = (1000,1000)

    #supTitle= ""
    # why does glvisualize() have terrible title() issues? oh well comment out for now
    supTitle="prm with maxdist=$connectRadius, "* "numPts = $numPts,"* "Avgd across #samples = $N" *
            "\n(time to run:$timeExperiment, "*
            "\ntimestamp=$timestamp)"
    costTitle= "\nclutter % vs pathcost\n"

    yerrCost = 10
    print("type of stddevs, $(typeof(stddevs))")
    stddevs = stddevs'
    yerrCost = [0 1.0 2.1 5.1 10.0]
    print("type of yerrCost , $(typeof(yerrCost))")

    pPRMcost = scatter(clutterPercentages, avgCosts',
        markercolor = :black,
        title = supTitle * costTitle, ylabel = "euclidean path cost", xlabel = "clutter %",
        yaxis=((0,50), 0:5:40), yerr = stddevs) #this automatically centers the 1 stddev, I think. So each side is +- 0.5 stddev?
         # yerr = yerrCost)

    successTitle = "clutter % vs pSuccess"
    pPRMsuccess = scatter(clutterPercentages, listpSucc',
        color = :orange, markersize= 6,
        title = successTitle, title_location=:center, ylabel = ("P(success)=numSucc/$N trials"), xlabel = "clutter %",
        yaxis=((0, 1.2), 0:0.1:1))

    plot(pPRMcost, pPRMsuccess, layout=(2,1), legend=false,
        xaxis=((-0.05, 0.6), 0:0.05:0.6),
        size = sizeplot)
end
#testGenClutter()

clutterExp()


text:
# # Matrix Completion

# We present another example that is about matrix completion. The idea is, given a partially observed matrix ``Y\in\mathbb{R}^{m\times n}``, to find
# ``X\in\mathbb{R}^{m\times n}`` to minimize the sum of squared errors from the observed entries while 'completing' the matrix ``Y``, i.e. filling the unobserved
# entries to match ``Y`` as good as possible. A detailed explanation can be found in section 4.2 of
# [the paper](https://arxiv.org/pdf/2104.06675.pdf).
# We will try to solve
# ```math
# \min_{||X||_*\le \tau} \sum_{(i,j)\in\mathcal{I}} (X_{i,j}iY_{i,j})^2,
# ```
# where ``\tau>0`` and ``\mathcal{I}`` denote the indices of the observed entries. We will use [`FrankWolfe.NuclearNormLMO`](@ref) and compare our
# Frank-Wolfe implementation with a Projected Gradient Descent (PGD) algorithm which, after each gradient descent step, projects the iterates back onto the nuclear
# norm ball. We use a movielens dataset for comparison.

using FrankWolfe
using ZipFile, DataFrames, CSV

using Random
using Plots

using Profile

import Arpack
using SparseArrays, LinearAlgebra

using LaTeXStrings

temp_zipfile = download("http://files.grouplens.org/datasets/movielens/ml-latest-small.zip")

zarchive = ZipFile.Reader(temp_zipfile)

movies_file = zarchive.files[findfirst(f -> occursin("movies", f.name), zarchive.files)]
movies_frame = CSV.read(movies_file, DataFrame)

ratings_file = zarchive.files[findfirst(f -> occursin("ratings", f.name), zarchive.files)]
ratings_frame = CSV.read(ratings_file, DataFrame)

users = unique(ratings_frame[:, :userId])
movies = unique(ratings_frame[:, :movieId])

@assert users == eachindex(users)
movies_revert = zeros(Int, maximum(movies))
for (idx, m) in enumerate(movies)
    movies_revert[m] = idx
end
movies_indices = [movies_revert[idx] for idx in ratings_frame[:, :movieId]]

const rating_matrix = sparse(
    ratings_frame[:, :userId],
    movies_indices,
    ratings_frame[:, :rating],
    length(users),
    length(movies),
)

missing_rate = 0.05

Random.seed!(42)

const missing_ratings = Tuple{Int,Int}[]
const present_ratings = Tuple{Int,Int}[]
let
    (I, J, V) = SparseArrays.findnz(rating_matrix)
    for idx in eachindex(I)
        if V[idx] > 0
            if rand() <= missing_rate
                push!(missing_ratings, (I[idx], J[idx]))
            else
                push!(present_ratings, (I[idx], J[idx]))
            end
        end
    end
end

function f(X)
    r = 0.0
    for (i, j) in present_ratings
        r += 0.5 * (X[i, j] - rating_matrix[i, j])^2
    end
    return r
end

function grad!(storage, X)
    storage .= 0
    for (i, j) in present_ratings
        storage[i, j] = X[i, j] - rating_matrix[i, j]
    end
    return nothing
end

function test_loss(X)
    r = 0.0
    for (i, j) in missing_ratings
        r += 0.5 * (X[i, j] - rating_matrix[i, j])^2
    end
    return r
end

function project_nuclear_norm_ball(X; radius=1.0)
    U, sing_val, Vt = svd(X)
    if (sum(sing_val) <= radius)
        return X, -norm_estimation * U[:, 1] * Vt[:, 1]'
    end
    sing_val = FrankWolfe.projection_simplex_sort(sing_val, s=radius)
    return U * Diagonal(sing_val) * Vt', -norm_estimation * U[:, 1] * Vt[:, 1]'
end

norm_estimation = 10 * Arpack.svds(rating_matrix, nsv=1, ritzvec=false)[1].S[1]

const lmo = FrankWolfe.NuclearNormLMO(norm_estimation)
const x0 = FrankWolfe.compute_extreme_point(lmo, ones(size(rating_matrix)))
const k = 10

FrankWolfe.benchmark_oracles(
    f,
    (str, x) -> grad!(str, x),
    () -> randn(size(rating_matrix)),
    lmo;
    k=10,
)

gradient = spzeros(size(x0)...)
gradient_aux = spzeros(size(x0)...)

function build_callback(trajectory_arr)
    return function callback(state)
        return push!(trajectory_arr, (Tuple(state)[1:5]..., test_loss(state.x)))
    end
end

# The smoothness constant is estimated:

num_pairs = 100
L_estimate = -Inf
for i in 1:num_pairs
    global L_estimate
    u1 = rand(size(x0, 1))
    u1 ./= sum(u1)
    u1 .*= norm_estimation
    v1 = rand(size(x0, 2))
    v1 ./= sum(v1)
    x = FrankWolfe.RankOneMatrix(u1, v1)
    u2 = rand(size(x0, 1))
    u2 ./= sum(u2)
    u2 .*= norm_estimation
    v2 = rand(size(x0, 2))
    v2 ./= sum(v2)
    y = FrankWolfe.RankOneMatrix(u2, v2)
    grad!(gradient, x)
    grad!(gradient_aux, y)
    new_L = norm(gradient - gradient_aux) / norm(x - y)
    if new_L > L_estimate
        L_estimate = new_L
    end
end

# We can now perform projected gradient descent:

xgd = Matrix(x0)
function_values = Float64[]
timing_values = Float64[]
function_test_values = Float64[]

time_start = time_ns()
for _ in 1:k
    f_val = f(xgd)
    push!(function_values, f_val)
    push!(function_test_values, test_loss(xgd))
    push!(timing_values, (time_ns() - time_start) / 1e9)
    @info f_val
    grad!(gradient, xgd)
    xgd_new, vertex = project_nuclear_norm_ball(xgd - gradient / L_estimate, radius=norm_estimation)
    gamma, _ = FrankWolfe.backtrackingLS(f, gradient, xgd, xgd - xgd_new, 1.0)
    @. xgd -= gamma * (xgd - xgd_new)
end

trajectory_arr_fw = Vector{Tuple{Int64,Float64,Float64,Float64,Float64,Float64}}()
callback = build_callback(trajectory_arr_fw)
xfin, _, _, _, traj_data = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0;
    epsilon=1e-9,
    max_iteration=10 * k,
    print_iter=k / 10,
    verbose=false,
    linesearch_tol=1e-8,
    line_search=FrankWolfe.Adaptive(),
    emphasis=FrankWolfe.memory,
    gradient=gradient,
    callback=callback,
)

trajectory_arr_lazy = Vector{Tuple{Int64,Float64,Float64,Float64,Float64,Float64}}()
callback = build_callback(trajectory_arr_lazy)
xlazy, _, _, _, _ = FrankWolfe.lazified_conditional_gradient(
    f,
    grad!,
    lmo,
    x0;
    epsilon=1e-9,
    max_iteration=10 * k,
    print_iter=k / 10,
    verbose=false,
    linesearch_tol=1e-8,
    line_search=FrankWolfe.Adaptive(),
    emphasis=FrankWolfe.memory,
    gradient=gradient,
    callback=callback,
)

trajectory_arr_lazy_ref = Vector{Tuple{Int64,Float64,Float64,Float64,Float64,Float64}}()
callback = build_callback(trajectory_arr_lazy_ref)
xlazy, _, _, _, _ = FrankWolfe.lazified_conditional_gradient(
    f,
    grad!,
    lmo,
    x0;
    epsilon=1e-9,
    max_iteration=50 * k,
    print_iter=k / 10,
    verbose=false,
    linesearch_tol=1e-8,
    line_search=FrankWolfe.Adaptive(),
    emphasis=FrankWolfe.memory,
    gradient=gradient,
    callback=callback,
)

fw_test_values = getindex.(trajectory_arr_fw, 6)
lazy_test_values = getindex.(trajectory_arr_lazy, 6)

results = Dict("svals_gd"=>svdvals(xgd),
"svals_fw"=>svdvals(xfin),
"svals_lcg"=>svdvals(xlazy),
"fw_test_values"=>fw_test_values,
"lazy_test_values"=>lazy_test_values,
"trajectory_arr_fw"=>trajectory_arr_fw,
"trajectory_arr_lazy"=>trajectory_arr_lazy,
"function_values_gd"=>function_values,
"function_values_test_gd"=>function_test_values,
"timing_values_gd"=>timing_values,
"trajectory_arr_lazy_ref"=>trajectory_arr_lazy_ref)

ref_optimum = results["trajectory_arr_lazy_ref"][end][2]

iteration_list = [
    [x[1] + 1 for x in results["trajectory_arr_fw"]],
    [x[1] + 1 for x in results["trajectory_arr_lazy"]],
    collect(1:1:length(results["function_values_gd"])),
]
time_list = [
    [x[5] for x in results["trajectory_arr_fw"]],
    [x[5] for x in results["trajectory_arr_lazy"]],
    results["timing_values_gd"],
]
primal_gap_list = [
    [x[2] - ref_optimum for x in results["trajectory_arr_fw"]],
    [x[2] - ref_optimum for x in results["trajectory_arr_lazy"]],
    [x - ref_optimum for x in results["function_values_gd"]],
]
test_list =
    [results["fw_test_values"], results["lazy_test_values"], results["function_values_test_gd"]]

label = [L"\textrm{FW}", L"\textrm{L-CG}", L"\textrm{GD}"]

FrankWolfe.plot_results(
    [primal_gap_list, primal_gap_list, test_list, test_list],
    [iteration_list, time_list, iteration_list, time_list],
    label,
    [L"\textrm{Iteration}", L"\textrm{Time}", L"\textrm{Iteration}", L"\textrm{Time}"],
    [
        L"\textrm{Primal Gap}",
        L"\textrm{Primal Gap}",
        L"\textrm{Test Error}",
        L"\textrm{Test Error}",
    ],
    xscalelog=[:log, :identity, :log, :identity],
    legend_position=[:bottomleft, nothing, nothing, nothing]
)