text:
module CurryRTS

export Node, runCurry, makeChoice, setChoice, makeFree, setRHS, setConst,
       setForwardNode, bindVarNode, voidNode, rew_nf

##############################################################################
# Definition of graph nodes
mutable struct Node
  tag    :: UInt8       # 0: operation
                        # 1: forward
                        # 2: choice
                        # 3: fail
                        # 4: free
                        # 5: integer
                        # 6: float
                        # 7: char
                        # 8: partial function call
                        # 9: partial constructor call
                        # >=10: constructor
  value  :: Int64       # demanded argument (0: nothing demanded) for functions,
                        # choice identifier for choices,
                        # index of constructors or integer/char value
                        # missing arguments for partial calls
  otsk   :: Int64       # identifier of the owner task
  taskns :: Union{Nothing,Dict{Int64,Node}} # forward nodes for different tasks
  fcode                 # function code
  args   :: Array{Node} # arguments or forward node
  symbol :: Union{String,Float64} # string representation of node symbol
                                  # or value of a float
end

##############################################################################
# Definition of tasks
mutable struct Task
  control :: Node              # the root node to be evaluated by this task
  fprint  :: Dict{Int64,UInt8} # the fingerprint for this task
  tasknum :: Int64             # task number (for debugging)
  parents :: Array{Int32,1}    # the numbers of the parent tasks
end

##############################################################################

# Evaluate all (non-deterministic) results of a main operation.
function runCurry(mainname::String, maincode, withnf::Bool, interactive::Bool,
                  withbfs::Bool, onevalue::Bool, quiet::Bool)
  tasks::Vector{Task} = []
  fp = Dict{Int64,UInt8}()
  taskcounter = 0 # counter to enumerate new tasks
  nd = Node(0, 0, taskcounter, nothing, maincode, [], mainname)
  mainnd = withnf ? makeNormalForm(nd) : nd
  currtask = Task(mainnd, fp, taskcounter, [])
  createdtasks = 1
  finished = false
  if !quiet ; print("EVALUATE: "); printNode(currtask, mainnd); println() ; end
  while ! finished
    #println("EVALUATING TASK " * string(currtask.tasknum) * "...")
    result = hnfNode(currtask, taskcounter)
    if result.tag >= 3 # fail, free, constructor, integer,...
      if result.tag > 3 # no failing computation
        if !quiet ; print("RESULT: ") end
        printNode(currtask, result); println()
        finished = onevalue
        if !finished && interactive && ! isempty(tasks)
          answer = "?"
          while answer != "" && answer != "y" && answer != "n" && answer != "a"
            print("More values? [Y(es)/n(o)/a(ll)] ")
            answer = chomp(readline())
            if answer == "a"
              interactive = false
            elseif answer == "n"
              finished = true
            end
          end
        end
      end
      if isempty(tasks)
        finished = true
      else
        currtask = withbfs ? popfirst!(tasks) : pop!(tasks)
      end
    elseif result.tag == 2 # choice
      fp = currtask.fprint
      if haskey(fp, result.value)  # selection exists for this choice
        currtask.control = result.args[fp[result.value]]
      else
        #println("PULL TAB AT ROOT, NEW TASK:" * string(taskcounter+1))
        fp2 = copy(currtask.fprint)
        fp2[result.value] = 2
        newparents = vcat([currtask.tasknum], currtask.parents)
        currtask.fprint[result.value] = 1
        currtask.control = result.args[1]
        currtask.tasknum = taskcounter+1
        currtask.parents = newparents
        if withbfs
          push!(tasks, currtask)
          currtask = popfirst!(tasks)
        end
        push!(tasks, Task(result.args[2], fp2, taskcounter+2, newparents))
        taskcounter += 2
        createdtasks += 1
      end
    end
  end
  if !quiet ; println("NUMBER OF CREATED TASKS: " * string(createdtasks)) ; end
end

# The current task during the evaluation of hnfNode.
# It is used for makeFree and makeChoice if a free variable is instantiated
# by a narrowing step. It is initialized with a dummy value.
currentHnfTask = Task(Node(0,0,0,nothing,nothing,[],""),
                      Dict{Int64,UInt8}(), 0, [])

# Compute the head-normal form of the node contained in the control
# of the given task.
# The node is updated (maybe as a forward node) and the result node
# (which might be constructor-rooted or a choice or failure node) is returned.
function hnfNode(task :: Task, tskcnt :: Int64) :: Node
  global currentHnfTask
  currentHnfTask = task
  stack ::Vector{Node} = []
  cn = task.control
  while !isempty(stack) || cn.tag == 0 || cn.tag == 1
    #traceTask("ROOT", task, isempty(stack) ? cn : stack[1])
    if cn.tag == 4 && cn.otsk < task.tasknum && cn.taskns != nothing
      # follow the possible instantiation of a free variable:
      cn1 = followTaskForwardNode(task, cn)
      if cn1 != cn && !isempty(stack)
        # binding found, replace corresponding argument in top stack function:
        topfun = last(stack)
        if task.tasknum > topfun.otsk
          # insert forward node in topfun:
          if topfun.taskns == nothing
            #println("CREATE TASKNS")
            topfun.taskns = Dict{Int64,Node}()
          end
          fwdnode = Node(0,topfun.value,task.tasknum,nothing,topfun.fcode,
                         copy(topfun.args),topfun.symbol)
          topfun.taskns[task.tasknum] = fwdnode
          stack[end] = fwdnode
        end
        last(stack).args[topfun.value] = cn1
      end
      cn = cn1
    end
    if cn.tag >= 4      # free, constructor, integer,...
      cn = pop!(stack)

      if !isempty(stack)
        topfun = last(stack)
        if cn.otsk > topfun.otsk
          # insert forward node in topfun:
          if topfun.taskns == nothing
            #println("CREATE TASKNS")
            topfun.taskns = Dict{Int64,Node}()
          end
          fwdnode = Node(0,topfun.value,cn.otsk,nothing,topfun.fcode,
                         copy(topfun.args),topfun.symbol)
          topfun.taskns[cn.otsk] = fwdnode
          stack[end] = fwdnode
        end
        # set argument of top function (just to be safe for forward nodes)
        topfun = last(stack)
        topfun.args[topfun.value] = cn
      end

      cn.fcode(cn) # apply rewrite rule

    elseif cn.tag == 1  # forward node?
      #println("FORWARD")
      cn = cn.args[1]
      if !isempty(stack)
        topfun = last(stack)
        topfun.args[topfun.value] = cn
      end
    elseif cn.tag == 2  # choice node? follow choice or pull-tab step:
      fp = task.fprint
      if haskey(fp, cn.value)  # selection exists for this choice
        #println("FOLLOW CHOICE SELECTION")
        cn = cn.args[fp[cn.value]]
        if !isempty(stack)
          topfun = last(stack)
          if task.tasknum > topfun.otsk
            # insert forward node in topfun:
            if topfun.taskns == nothing
              #println("CREATE TASKNS")
              topfun.taskns = Dict{Int64,Node}()
            end
            fwdnode = Node(0,topfun.value,task.tasknum,nothing,topfun.fcode,
                           copy(topfun.args),topfun.symbol)
            topfun.taskns[task.tasknum] = fwdnode
            stack[end] = fwdnode
          end
          last(stack).args[topfun.value] = cn
        end
      else # pull tab step:
        fn = pop!(stack)
        args1 = copy(fn.args)
        args1[fn.value] = cn.args[1]
        args2 = copy(fn.args)
        args2[fn.value] = cn.args[2]
        c1 = Node(0,fn.value,tskcnt+1,nothing,fn.fcode,args1,fn.symbol)
        c2 = Node(0,fn.value,tskcnt+2,nothing,fn.fcode,args2,fn.symbol)
        fn.tag    = 2
        fn.value  = cn.value
        fn.taskns = nothing
        fn.fcode  = nothing
        empty!(fn.args); push!(fn.args,c1); push!(fn.args,c2)
        fn.symbol = "?"
        cn = fn
      end
    elseif cn.tag == 3  # failure node?
      empty!(stack)
    # final alternative: cn.tag == 0, i.e., operation
    else
      hasfwdnode = false
      if cn.otsk < task.tasknum && cn.taskns != nothing
        # check whether there is a forward node in this or some parent task:
        if haskey(cn.taskns, task.tasknum)
          hasfwdnode = true
          cn = cn.taskns[task.tasknum]
        else
          for i in 1:length(task.parents)
            if haskey(cn.taskns, task.parents[i])
              hasfwdnode = true
              cn = cn.taskns[task.parents[i]]
              break
            end
          end
        end
      end
      if hasfwdnode
        #println("FOLLOW RESULT FORWARD NODE")
        if !isempty(stack)
          topfun = last(stack)
          topfun.args[topfun.value] = cn
        end
      elseif cn.value == 0 # no argument demanded?
        cn.fcode(cn) # apply rewrite rule
      else
        push!(stack, cn)   # push current function on stack
        cn = cn.args[cn.value] # proceed with demanded argument
      end
    end
  end
  #traceTask("HNF", task, cn)
  return cn
end

##############################################################################
# Auxiliary operations for dealing with choice nodes

# The global counter for fresh choice identifiers
freshChoiceId = 0

# Create a choice node with a fresh choice identifier for a given task
# identifier.
function makeChoice(otsk, l, r) :: Node
  global freshChoiceId
  freshChoiceId += 1
  return Node(2,freshChoiceId,otsk,nothing,nothing,[l,r],"?")
end

# replace argument node by a choice node with a fresh choice identifier
function setChoice(root :: Node, l :: Node, r :: Node)
  global freshChoiceId
  freshChoiceId += 1
  root.symbol = "?"
  root.tag    = 2
  root.value  = freshChoiceId
  empty!(root.args); push!(root.args,l); push!(root.args,r)
end

##############################################################################
# Auxiliary operations

# Node representing a void value. Used to initialize Node variables.
function voidNode(otsk)
  return Node(10, 0, otsk, nothing, nothing, [], "()")
end

# Create a free variable node
function makeFree(otsk) :: Node
  return Node(4, 0, otsk, nothing, nothing, [], "FREE")
end

# Set the first argument node (which is the root of a redex)
# to the right-hand side provided as the remaining arguments.
function setRHS(root :: Node, ntag, nvalue, nfcode, nargs, nsymbol)
  root.tag    = ntag
  root.value  = nvalue
  root.fcode  = nfcode
  root.args   = nargs
  root.symbol = nsymbol
end

# Set the first argument node (which is the root of a redex)
# to a constant value provided as the remaining arguments.
function setConst(root :: Node, ntag, nvalue)
  root.tag     = ntag
  root.value   = nvalue
end

# Set the first argument node to a forward node.
function setForwardNode(root :: Node, fnode :: Node)
  root.tag  = 1
  empty!(root.args); push!(root.args,fnode)
end

# Bind the first argument node (a free variable) to another node.
# Note that this binding is made only for the current task.
function bindVarNode(varnode :: Node, n :: Node)
  global currentHnfTask
  if currentHnfTask.tasknum > varnode.otsk
    # insert forward node in varnode:
    if varnode.taskns == nothing
      #println("CREATE TASKNS FOR FREE VARIABLE")
      varnode.taskns = Dict{Int64,Node}()
    end
    varnode.taskns[currentHnfTask.tasknum] = n
  else
    setForwardNode(varnode, n)
  end
end

##############################################################################
# ADDITIONAL STUFF TO SUPPORT THE REPL:
##############################################################################
# Normal form computation.

# Wrap a node with the `normalForm` operation.
function makeNormalForm(root :: Node)
  return Node(0,1,root.otsk,nothing,rew_nf,[root],"normalForm")
end

# Operation to compute the normal form of its argument.
# Implemented by wrapping all arguments with `normalForm` and
# evaluating successively all arguments (by operation `rew_nfArgs`).
function rew_nf(root :: Node)
  x = root.args[1]
  if x.tag == 4 # free variable: create forward node
    root.tag  = 1
    empty!(root.args); push!(root.args,x)
  elseif x.tag < 10 || length(x.args) == 0 # no argument evaluation necessary
    root.tag    = x.tag
    root.symbol = x.symbol
    root.value  = x.value
    root.otsk   = x.otsk
    root.taskns = x.taskns
    root.fcode  = nothing
    root.args   = x.args
  else
    root.value = 1
    root.fcode = rew_nfArgs
    root.args  = vcat(map(makeNormalForm,x.args), [x]) # remember constructor
  end
end

function rew_nfArgs(root :: Node)
  darg = root.value + 1
  if darg == length(root.args) # all constructor arguments evaluated
    cterm       = pop!(root.args)
    root.tag    = cterm.tag
    root.symbol = cterm.symbol
    root.value  = cterm.value
    root.otsk   = cterm.otsk
    root.taskns = cterm.taskns
    root.fcode  = nothing
  else
    root.value = darg
  end
end

##############################################################################
# Auxiliaries for printing results and tracing

# Print a graph node as a term in a task, i.e., fingerprints are considered
# when printing choice nodes.
function printNode(tsk:: Task, n :: Node)
  hasfwdnode = false
  if (n.tag == 0 || n.tag == 4) && n.otsk < tsk.tasknum && n.taskns != nothing
    n1 = followTaskForwardNode(tsk, n)
    if n1 != n
      n = n1
      hasfwdnode = true
    end
  end
  if hasfwdnode
    # follow result forward node:
    print("%") ; printNode(tsk, n)
  elseif n.tag == 1 # forward node
    print("^") ; printNode(tsk, n.args[1])
  elseif n.tag == 2 # choice node
    fp = tsk.fprint
    if haskey(fp, n.value)
      printNode(tsk, n.args[fp[n.value]])
    else
      print("(")
      printNode(tsk, n.args[1])
      print(" ?" * string(n.value) * " ")
      printNode(tsk, n.args[2])
      print(")")
    end
  elseif n.tag == 3 # fail node
    print("!")
  elseif n.tag == 4 # free node
    #print("FREE")
    #print("_x" * string(objectid(n)))
    print("_x" * string(objectid(n))[end-3:end])
    #printNodeSymbol(n)
  elseif n.tag == 5 # integer node
    print(string(n.value))
  elseif n.tag == 7 # character node
    if n.value > 31
      print("'" * Char(n.value) * "'")
    else
      print("(char " * string(n.value) * ")")
    end
  elseif length(n.args) == 0
    printNodeSymbol(n)
  elseif length(n.args) == 2 && !isletter(n.symbol[1]) # infix operator
    print("(")
    printNode(tsk,n.args[1])
    print(" ") ; printNodeSymbol(n) ; print(" ")
    printNode(tsk,n.args[2])
    print(")")
  else
    print("(")
    printNodeSymbol(n)
    for i = 1 : length(n.args)
      print(" ") ; printNode(tsk,n.args[i])
    end
    print(")")
  end
end

function followTaskForwardNode(tsk:: Task, n :: Node)
  if haskey(n.taskns, tsk.tasknum)
    n = n.taskns[tsk.tasknum]
  else
    for i in 1:length(tsk.parents)
      if haskey(n.taskns, tsk.parents[i])
        n = n.taskns[tsk.parents[i]]
        break
      end
    end
  end
  return n
end

# Print a node symbol
function printNodeSymbol(n :: Node)
  #print(n.symbol * (n.tag == 0 ? "[" * string(objectid(n)) * "]" : ""))
  #print(n.symbol * (n.tag == 0 ? "[" * string(n.otsk) * "]" : ""))
  print(n.symbol)
end

# Show a finger print
function showFPrint(fp :: Dict{Int64,UInt8}) :: String
  s = "["
  for (index, value) in pairs(fp)
    s = s * string(index) * "=>" * (value==1 ? "L" : "R") * " "
  end
  return rstrip(s) * "]"
end

function parentsString(p) :: String
  s = ""
  for i in 1:length(p)
    s = s * "," * string(p[i])
  end
  return s
end

function traceTask(ctxt :: String, task :: Task, cn :: Node)
  for i = 1:length(task.parents)
    print("  ")
  end
  print(ctxt * "[" * string(task.tasknum) * parentsString(task.parents) * "]")
  print(showFPrint(task.fprint))
  print(": ")
  printNode(task,cn)
  println()
end

##############################################################################

end

max_stars_repo_name: cau-placc/julia-curry

max_issues_repo_name: cau-placc/julia-curry

max_forks_repo_name: cau-placc/julia-curry


text:
function _filterdown(f, root; kwargs...)
    founds = String[]
    walkdown(root; kwargs..., thfrec = 0.0) do path
        f(path) && let
            push!(founds, path)
        end
        false
    end
    founds
end

function _filterdown_th(f, root; thfrec = 1.0, kwargs...)
    founds_pool = [String[] for _ in 1:nthreads()]
    walkdown(root; thfrec, kwargs...) do path
        f(path) && push!(founds_pool[threadid()], path)
        false
    end
    vcat(founds_pool...)
end

function filterdown(f::Function, root;
        thfrec::Float64 = 0.0, kwargs...
    )
    thfrec > 0.0 ?
        _filterdown_th(f, root; thfrec, kwargs...) :
        _filterdown(f, root; kwargs...)
end

ext: jl


text:
using Documenter
using Glimmer, Glimmer.FlexUI

makedocs(
    sitename = "Glimmer.jl",
    format = Documenter.HTML(
        # prettyurls = get(ENV, "CI", nothing) == "true",
        assets = [asset("assets/GlimmerJulia.png", class = :ico, islocal = true)],
    ),
    modules = [Glimmer],
    pages = [
        "Home" => "index.md",
        "Overview" => "overview.md",
        "Examples" => [
            "Meshing Example" => "examples/Meshing.md",
        ],
        "Reference" => "ref.md",
    ],
    expandfirst = [])

# deploydocs(
#     repo = "github.com/galran/Glimmer.jl.git",
#     devbranch = "main",
#     push_preview = true,
#     target = "build",
#     deps = nothing,
#     make = nothing,
# )

deploydocs(
    repo = "github.com/galran/Glimmer.jl.git",
    branch = "gh-pages",
    devbranch = "main",
    push_preview = true,
)


text:
using Clustering
using Random

function ind_to_indexes(c)
    a = similar(c, Int)
    count = 1
    @inbounds for i in eachindex(c)
        a[count] = i
        count += (c[i] != zero(eltype(c)))
    end
    return resize!(a, count-1)
end

function sample_ind_wo_replacement(c, k)
    a = falses(size(c))
    idxs = ind_to_indexes(c)
    perm = randperm(length(idxs))
    a[idxs[perm[1:k]]] .= 1
    return a
end

function kmeans_filter1(reps, eps)
    n = size(reps)[2]
    kmass = kmeans(reps, 2).assignments
    c1, c2 = kmass .== 1, kmass .== 2
    c1c, c1p, c2c, c2p = (
        sum(c1[1:end-eps]),
        sum(c1[end-eps+1:end]),
        sum(c2[1:end-eps]),
        sum(c2[end-eps+1:end])
    )
    if c1p/(c1p + c1c) > c2p/(c2p + c2c)
        good, bad = c1, c2
    else
        good, bad = c2, c1
    end
    return sample(ind_to_indexes(good))
end

function kmeans_filter2(reps, eps, k=100, limit=1.5)
    reps_pca = pca(reps, 100)[1]
    to_remove = Set{Int}()
    n = 0
    while length(to_remove) < round(limit*eps)
        i = kmeans_filter1(reps_pca, eps)
        push!(to_remove, i)
        n += 1
        if n > 10000
            break
        end
    end
    ind = trues(size(reps)[2])
    for i in to_remove
        ind[i] = 0
    end
    return ind
end

max_stars_repo_name: TDteach/spectre-defense

max_issues_repo_name: TDteach/spectre-defense

max_forks_repo_name: TDteach/spectre-defense


text:
export SuperHalley

function SuperHalley(nlp_stop,
				     PData :: PDataFact,
				     dₙ :: Vector,
				     gt :: Vector)

	nlp_at_x = nlp_stop.current_state
	x = nlp_at_x.x

	n = length(x)

	e₁ = zeros(n)
    e₁[1] = 1.0
    ∇³fdₙ = ∇f³xuv(nlp_stop.pb, x, dₙ, e₁)
    for i = 2:n
        eᵢ = zeros(n)
        eᵢ[i] = 1.0
        ∇³fdₙ = hcat(∇³fdₙ, ∇f³xuv(nlp_stop.pb, x, dₙ, eᵢ))
    end

	# hess_plus_tenseur = (nlp_at_x.Hx + ∇³fdₙ)
	if PData.λ .== 0.0
		hess_plus_tenseur = (nlp_at_x.Hx .+  ∇³fdₙ)
	else
		hess_plus_tenseur = (nlp_at_x.Hx .+ ∇³fdₙ) .+ (PData.λ .* Matrix(1.0I, n, n))
	end

	if isposdef(hess_plus_tenseur)
        (L, D, pp, ρ, ncomp) = ldlt_symm(hess_plus_tenseur, 'r')
   	else
		# println("on est ici")
	   	return dₙ
   	end

	# try
    #    (L, D, pp, ρ, ncomp) = ldlt_symm(hess_plus_tenseur, 'r')
	# catch
 	# 	println("*******   Problem in LDLt Halley")
    #    	#res = PDataLDLt()
    #    	#res.OK = false
    #    	#return res
    #    	d = NaN * zeros(n)
    #    	return d, xdemi
   	# end

	if true in isnan.(D)
   	   println("*******   Problem in D from LDLt: NaN")
	   println(" cond (H) = $(cond(H))")
	   #res = PDataLDLt()
	   #res.OK = false
	   #return res
	   d = NaN * zeros(n)
	   return d, xdemi
   	end

	D = Hermitian(D)
   	X = eigen(D)
  	Δ = X.values
  	Q =  X.vectors

	ϵ =  sqrt(eps(eltype(gt)))
   	Γ = max.(abs.(Δ), ϵ)

	d̃ = L \ nlp_at_x.gx[pp]
    d̂ = L' \ (Q * (Q' * d̃ ./ Γ))
    dV = -d̂[invperm(pp)]

	grad_plus_tuv = -(nlp_at_x.gx + 0.5 .* ∇f³xuv(nlp_stop.pb, x, dₙ, dV))
	grad_plus_tuv = -grad_plus_tuv

	ϵ2 =  sqrt(eps(eltype(gt)))
   	Γ2 = max.(abs.(PData.Δ), ϵ2)

	d̃ = PData.L \ grad_plus_tuv[pp]
    d̂ = PData.L' \ (PData.Q * (PData.Q' * d̃ ./ Γ2))
    dSH = -d̂[invperm(PData.pp)]

	return dSH
end


text:
include("cheb.jl")
using .cheb
(a,b) = projmat(5)

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
# This file was generated, do not modify it.

using HTTP, CSV, MLJ, StatsBase, PyPlot
req = HTTP.get("https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data")
data = CSV.read(req.body,
                header=["Class", "Alcool", "Malic acid",
                        "Ash", "Alcalinity of ash", "Magnesium",
                        "Total phenols", "Flavanoids",
                        "Nonflavanoid phenols", "Proanthcyanins",
                        "Color intensity", "Hue",
                        "OD280/OD315 of diluted wines", "Proline"])
# the target is the Class column, everything else is a feature
y, X = unpack(data, ==(:Class), colname->true);

scitype(y)

yc = coerce(y, OrderedFactor);

scitype(X)

sch = schema(X)
println(rpad(" Name", 28), "| Scitype")
println("-"^45)
for (name, scitype) in zip(sch.names, sch.scitypes)
    println(rpad("$name", 30), scitype)
end

X[1:5, :Proline]

Xc = coerce(X, :Proline=>Continuous, :Magnesium=>Continuous);

describe(Xc, :mean, :std)

@load KNNClassifier pkg="NearestNeighbors"
@load MultinomialClassifier pkg="MLJLinearModels";

@pipeline KnnPipe(std=Standardizer(), clf=KNNClassifier()) is_probabilistic=true
@pipeline MnPipe(std=Standardizer(), clf=MultinomialClassifier()) is_probabilistic=true

train, test = partition(eachindex(yc), 0.8, shuffle=true, rng=111)
Xtrain = selectrows(Xc, train)
Xtest = selectrows(Xc, test)
ytrain = selectrows(yc, train)
ytest = selectrows(yc, test);

knn = machine(KnnPipe(), Xtrain, ytrain)
multi = machine(MnPipe(), Xtrain, ytrain)

opts = (resampling=Holdout(fraction_train=0.9), measure=cross_entropy)
res = evaluate!(knn; opts...)
round(res.measurement[1], sigdigits=3)

res = evaluate!(multi; opts...)
round(res.measurement[1], sigdigits=3)

mcr_k = misclassification_rate(predict_mode(knn, Xtrain), ytrain)
mcr_m = misclassification_rate(predict_mode(multi, Xtrain), ytrain)
println(rpad("KNN mcr:", 10), round(mcr_k, sigdigits=3))
println(rpad("MNC mcr:", 10), round(mcr_m, sigdigits=3))

# @pipeline PCAPipe(std=Standardizer(), t=PCA(maxoutdim=2))
# pca = machine(PCAPipe(), Xtrain)
# fit!(pca, Xtrain)
# W = transform(pca, Xtrain)

@load PCA

pca = Xc |> Standardizer() |> PCA(maxoutdim=2)
fit!(pca)
W = pca(rows=train);

x1 = W.x1
x2 = W.x2

mask_1 = ytrain .== 1
mask_2 = ytrain .== 2
mask_3 = ytrain .== 3

figure(figsize=(8, 6))
plot(x1[mask_1], x2[mask_1], linestyle="none", marker="o", color="red")
plot(x1[mask_2], x2[mask_2], linestyle="none", marker="o", color="blue")
plot(x1[mask_3], x2[mask_3], linestyle="none", marker="o", color="magenta")

xlabel("PCA dimension 1", fontsize=14)
ylabel("PCA dimension 2", fontsize=14)
legend(["Class 1", "Class 2", "Class 3"], fontsize=12)
xticks(fontsize=12)
yticks(fontsize=12)

savefig("assets/literate/EX-wine-pca.svg") # hide

perf_k = misclassification_rate(predict_mode(knn, Xtest), ytest)
perf_m = misclassification_rate(predict_mode(multi, Xtest), ytest)
println(rpad("KNN mcr:", 10), round(perf_k, sigdigits=3))
println(rpad("MNC mcr:", 10), round(perf_m, sigdigits=3))


text:
abstract type AbstractTreeViewNode end
@Base.kwdef struct TreeViewNode <: AbstractTreeViewNode
    id::String = string(UUIDs.uuid4())
    label::String = "nothing"
    children::Vector{TreeViewNode} = TreeViewNode[]
    data = nothing # insert arbitrary data here
end
label(t::TreeViewNode) = t.label
tree_id(t::TreeViewNode) = t.id
children(t::TreeViewNode) = t.children

"""
    TreeViewNode(::String)

Create a TreeViewNode
"""
function TreeViewNode(text::String)
    TreeViewNode(string(UUIDs.uuid4()), text, TreeViewNode[], nothing)
end

Base.push!(parent::TreeViewNode, child::TreeViewNode) = push!(children(parent), child)

"""
TreeViewRoot(::String)

Create a TreeViewRoot for use as content in a TreeView widget

# Examples

Create a TreeViewRoot
```
    tree = TreeViewRoot()
    beverages = TreeViewNode("Beverages")
    push!(tree, beverages)
    push!(tree, TreeViewNode("Food"))
```
And push some children into a node if you want
```
    push!(beverages, TreeViewNode("Water"))
    push!(beverages, TreeViewNode("Tea"))
    push!(beverages, TreeViewNode("Coffee"))
```

"""
struct TreeViewRoot
    nodes::AbstractArray{<:AbstractTreeViewNode}
end
TreeViewRoot() = TreeViewRoot(TreeViewNode[])

Base.push!(tree::TreeViewRoot, node::TreeViewNode) = push!(tree.nodes, node)

function example_treeview_root()
    tree = TreeViewRoot()

    beverages = TreeViewNode("Beverages")
    push!(tree, beverages)
    push!(beverages, TreeViewNode("Water"))
    push!(beverages, TreeViewNode("Tea"))

    coffee = TreeViewNode("Coffee")
    push!(beverages, coffee)
    push!(coffee, TreeViewNode("Espresso"))
    push!(coffee, TreeViewNode("Cappuccino"))

    meat = TreeViewNode("Meat")
    push!(tree, meat)
    push!(meat, TreeViewNode("Chicken"))
    push!(meat, TreeViewNode("Steak"))

    vegetables = TreeViewNode("Vegetables")
    push!(tree, vegetables)

    return tree
end

struct TreeView
    content::TreeViewRoot
    selection::Observable
    scope::Scope
end

"""
TreeView(::TreeViewRoot)

Create a TreeView widget from a TreeViewRoot structure

# Examples

Create a TreeView
```
    tree = TreeViewRoot()
    beverages = TreeViewNode("Beverages")
    push!(beverages, TreeViewNode("Water"))
    push!(beverages, TreeViewNode("Tea"))
    push!(tree, beverages)
    push!(tree, TreeViewNode("Food"))

    treeview = TreeView(tree)
```
To get the Observable selection of elements:
```
    treeview[]
```
"""
function TreeView(content::TreeViewRoot)
    treeview_scope = create_scope(tohtml(content))
    obs = Observable(treeview_scope, SELECTION_OBSERVABLE_NAME, Dict{String, String}[])
    return TreeView(content, obs, treeview_scope)
end

# standard behavior of widgets: get the observable via the getindex: treeview[]
Base.getindex(tv::TreeView) = tv.selection

Base.show(io::IO, x::TreeView) = show(io, x.scope)

# just forward all MIME to WebIO.Scope
Base.show(io::IO, m::MIME, x::TreeView) = show(io, m, x.scope)


text:
include(joinpath(pwd(), "src/models/hopper.jl"))
include(joinpath(pwd(), "src/constraints/contact.jl"))
include(joinpath(pwd(), "src/constraints/loop.jl"))
include(joinpath(pwd(), "src/constraints/free_time.jl"))

# Free-time model
model_ft = free_time_model(model)

function fd(model::Hopper, x⁺, x, u, w, h, t)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2⁺ = view(x⁺, 1:model.nq)
	q2⁻ = view(x, model.nq .+ (1:model.nq))
	q1 = view(x, 1:model.nq)
	u_ctrl = view(u, model.idx_u)
	λ = view(u, model.idx_λ)
	b = view(u, model.idx_b)
	h = u[end]

	[q2⁺ - q2⁻;
	((1.0 / h) * (M_func(model, q1) * (SVector{4}(q2⁺) - SVector{4}(q1))
	- M_func(model, q2⁺) * (SVector{4}(q3) - SVector{4}(q2⁺)))
	+ transpose(B_func(model, q3)) * SVector{2}(u_ctrl)
	+ transpose(N_func(model, q3)) * SVector{1}(λ)
	+ transpose(P_func(model, q3)) * SVector{2}(b)
	- h * G_func(model, q2⁺))]
end

function maximum_dissipation(model::Hopper, x⁺, u, h)
	q3 = x⁺[model.nq .+ (1:model.nq)]
	q2 = x⁺[1:model.nq]
	ψ = u[model.idx_ψ]
	ψ_stack = ψ[1] * ones(model.nb)
	η = u[model.idx_η]
	h = u[end]
	return P_func(model, q3) * (q3 - q2) / h + ψ_stack - η
end

# Horizon
T = 21

# Time step
tf = 1.0
h = tf / (T - 1)

# Bounds
_uu = Inf * ones(model_ft.m)
_uu[end] = 5.0 * h
_ul = zeros(model_ft.m)
_ul[model_ft.idx_u] .= -Inf
_ul[end] = 0.5 * h
ul, uu = control_bounds(model_ft, T, _ul, _uu)

# Initial and final states
q1 = [0.0, 0.25 , 0.0, 0.25]
q_right = [0.0, 0.5 + 0.25, pi / 2.0, 0.25]
q_top = [0.0, 0.5 + 0.5, pi, 0.25]
q_left = [0.0, 0.5 + 0.25, 3.0 * pi / 2.0, 0.25]
qT = [0.0, 0.5,  2.0 * pi, 0.5]

xl, xu = state_bounds(model_ft, T,
		[model_ft.qL; model_ft.qL],
		[model_ft.qU; model_ft.qU],
        x1 = [q1; Inf * ones(model_ft.nq)],
		xT = [Inf * ones(model_ft.nq); qT])

q_ref = [linear_interp(q1, q_right, 6)...,
         linear_interp(q_right, q_top, 6)[2:end]...,
         linear_interp(q_top, q_left, 6)[2:end]...,
         linear_interp(q_left, qT, 6)[2:end]...]

x_ref = configuration_to_state(q_ref)

# Objective
Qq = 1.0 * Diagonal([1.0, 1.0, 1.0, 1.0])
Q = cat(0.5 * Qq, 0.5 * Qq, dims = (1, 2))
# QT = cat(0.5 * Qq, 100.0 * Diagonal(ones(model_ft.nq)), dims = (1, 2))
R = Diagonal([1.0e-1, 1.0e-3, zeros(model_ft.m - model_ft.nu)...])

obj_tracking = quadratic_time_tracking_objective(
    [Q for t = 1:T],
    [R for t = 1:T-1],
    [x_ref[t] for t = 1:T],
    [zeros(model_ft.m) for t = 1:T],
    1.0)

obj_contact_penalty = PenaltyObjective(100.0, model_ft.m - 1)

obj = MultiObjective([obj_tracking, obj_contact_penalty])

# Constraints
con_free_time = free_time_constraints(T)
con_contact = contact_constraints(model_ft, T)
con_loop = loop_constraints(model_ft, 1, T)
con = multiple_constraints([con_free_time, con_contact])

# Problem
prob = trajectory_optimization_problem(model_ft,
               obj,
               T,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con
               )

# Trajectory initialization
X0 = deepcopy(x_ref) # linear interpolation on state
U0 = [[1.0e-5 * rand(model_ft.m-1); h] for t = 1:T-1] # random controls

# Pack trajectories into vector
Z0 = pack(X0, U0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
include("/home/taylor/.julia/dev/SNOPT7/src/SNOPT7.jl")

@time Z̄ = solve(prob, copy(Z0),
	nlp = :ipopt,
	tol = 1.0e-3, c_tol = 1.0e-3, mapl = 5)

X̄, Ū = unpack(Z̄, prob)

@show Ū[4][end]
@show check_slack(Z̄, prob)

using Plots
plot(hcat(Ū...)[end,:], linetype=:steppost)

include(joinpath(pwd(), "src/models/visualize.jl"))
vis = Visualizer()
open(vis)
visualize!(vis, model_ft, state_to_configuration(X̄), Δt = Ū[1][end])

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
# sample to show robot module's movement
module RobotMoveSample
  # include external modules
  include(joinpath(split(@__FILE__, "src")[1], "src/model/world/puddle_world.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/agent/agent.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/robot/differential_wheeled_robot/differential_wheeled_robot.jl"))

  # main function
  # show world's graph
  function main(delta_time=0.1, end_time=20; is_test=false)
    # save path of output
    path = "src/model/robot/robot_move_sample.gif"

    # world instance
    world = PuddleWorld(-5.0, 5.0, -5.0, 5.0,
                        delta_time=delta_time,
                        end_time=end_time,
                        is_test=is_test,
                        save_path=path)

    # robot instance
    # straight movement
    straight = Agent(0.2, 0.0)
    robot1 = DifferentialWheeledRobot([2.0, 3.0, pi/6], 0.2, "black",
                                      straight, delta_time)
    append(world, robot1) # add in world
    # circling movement
    circling = Agent(0.2, 10.0/180*pi)
    robot2 = DifferentialWheeledRobot([-2.0, -1.0, pi/5*6], 0.2, "red",
                                      circling, delta_time)
    append(world, robot2) # add in world
    # staying
    staying = Agent(0.0, 0.0)
    robot3 = DifferentialWheeledRobot([0.0, 0.0, 0.0], 0.2, "blue",
                                      staying, delta_time)
    append(world, robot3) # add in world

    # draw graph
    draw(world)
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
module problem109

description = """
In the game of darts a player throws three darts at a target board which is split into twenty equal sized sections numbered one to twenty.
The score of a dart is determined by the number of the region that the dart lands in. A dart landing outside the red/green outer ring scores zero. The black and cream regions inside this ring represent single scores. However, the red/green outer ring and middle ring score double and treble scores respectively.
At the centre of the board are two concentric circles called the bull region, or bulls-eye. The outer bull is worth 25 points and the inner bull is a double, worth 50 points.
There are many variations of rules but in the most popular game the players will begin with a score 301 or 501 and the first player to reduce their running total to zero is a winner. However, it is normal to play a "doubles out" system, which means that the player must land a double (including the double bulls-eye at the centre of the board) on their final dart to win; any other dart that would reduce their running total to one or lower means the score for that set of three darts is "bust".
When a player is able to finish on their current score it is called a "checkout" and the highest checkout is 170: T20 T20 D25 (two treble 20s and double bull).

There are exactly eleven distinct ways to checkout on a score of 6:

D3
D1 	D2
S2 	D2
D2 	D1
S4 	D1
S1 	S1 	D2
S1 	T1 	D1
S1 	S3 	D1
D1 	D1 	D1
D1 	S2 	D1
S2 	S2 	D1

Note that D1 D2 is considered different to D2 D1 as they finish on different doubles. However, the combination S1 T1 D1 is considered the same as T1 S1 D1.
In addition we shall not include misses in considering combinations; for example, D3 is the same as 0 D3 and 0 0 D3.
Incredibly there are 42336 distinct ways of checking out in total.

How many distinct ways can a player checkout with a score less than 100?
"""

# Use a triple nested loop to count all the solutions
# The outer loop is all the doubles that represent the final dart
# The two inner loops represent the preceding two throws.
# The count is incremented at each nested loop, as the game can be won in 1,2 or 3 throws.
function solve()
  count = 0
  singles = [1:20;25] # 25 for outer bull
  doubles = 2 * singles
  triples = 3 * singles[1:end-1] # no triple bull
  throws = [singles;doubles;triples]

  for i = 1:length(doubles)
    count += 1
    for j = 1:length(throws)
      if doubles[i] + throws[j] < 100
        count += 1
        for k = j:length(throws)
          if doubles[i] + throws[j] + throws[k] < 100
            count += 1
          end
        end
      end
    end
  end
  println("$count ways to go out on less than 100")
end

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
module PiGPIOC

const pigpio = "libpigpio"
include("pigpio_common.jl")
include("pigpio.jl")

end # module


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCellDG
include("../useful_routines.jl")
include("plane_incoherent_interface_transformation_solver.jl")
APS = AnalyticalPlaneSolver

L = 1.0
W = 1.0
lambda1, mu1 = 200.0, 80.0
lambda2, mu2 = 150.0, 40.0
penaltyfactor = 1e2
eta = 1
nelmts = 4
dx = L / nelmts
dy = W / nelmts
penalty = penaltyfactor / dx * (lambda1 + mu1 + lambda2 + mu2) * 0.5

theta0 = -0.01

stiffness = CutCellDG.HookeStiffness(lambda1, mu1, lambda2, mu2)
transfstress =
    CutCellDG.plane_strain_transformation_stress(lambda1, mu1, theta0)

polyorder = 1
numqp = 2

elasticitybasis = LagrangeTensorProductBasis(2, polyorder)
levelsetbasis = HermiteTensorProductBasis(2)

dim, nf = size(interpolation_points(levelsetbasis))
refpoints = interpolation_points(elasticitybasis)

dgmesh = CutCellDG.DGMesh([0.0, 0.0], [L, W], [nelmts, nelmts], refpoints)
cgmesh = CutCellDG.CGMesh([0.0, 0.0], [L, W], [nelmts, nelmts], nf)

R = 0.51
x0 = [R, 0.0]
normal = [1.0, 0.0]
levelset = CutCellDG.LevelSet(
    x -> plane_distance_function(x, normal, x0),
    cgmesh,
    levelsetbasis,
)

cutmesh = CutCellDG.CutMesh(dgmesh, levelset)
cellquads = CutCellDG.CellQuadratures(cutmesh, levelset, numqp)
interfacequads = CutCellDG.InterfaceQuadratures(cutmesh, levelset, numqp)
facequads = CutCellDG.FaceQuadratures(cutmesh, levelset, numqp)
mergedmesh = CutCellDG.MergedMesh!(
    cutmesh,
    cellquads,
    facequads,
    interfacequads,
    tinyratio = 0.3,
)

@assert CutCellDG.has_merged_cells(mergedmesh)

sysmatrix = CutCellDG.SystemMatrix()
sysrhs = CutCellDG.SystemRHS()

CutCellDG.assemble_displacement_bilinear_forms!(
    sysmatrix,
    elasticitybasis,
    cellquads,
    stiffness,
    mergedmesh,
)
CutCellDG.assemble_bulk_transformation_linear_form!(
    sysrhs,
    transfstress,
    elasticitybasis,
    cellquads,
    mergedmesh,
)
CutCellDG.assemble_interelement_condition!(
    sysmatrix,
    elasticitybasis,
    facequads,
    stiffness,
    mergedmesh,
    penalty,
    eta,
)
CutCellDG.assemble_interelement_transformation_linear_form!(
    sysrhs,
    transfstress,
    elasticitybasis,
    facequads,
    mergedmesh,
)

CutCellDG.assemble_incoherent_interface_condition!(
    sysmatrix,
    elasticitybasis,
    interfacequads,
    stiffness,
    mergedmesh,
    penalty,
    eta,
)
CutCellDG.assemble_incoherent_interface_transformation_linear_form!(
    sysrhs,
    transfstress,
    elasticitybasis,
    interfacequads,
    mergedmesh,
)

CutCellDG.assemble_penalty_displacement_component_bc!(
    sysmatrix,
    sysrhs,
    x -> 0.0,
    elasticitybasis,
    facequads,
    stiffness,
    mergedmesh,
    x -> x[1] ≈ 0.0,
    [1.0, 0.0],
    penalty,
)
CutCellDG.assemble_penalty_displacement_component_transformation_linear_form!(
    sysrhs,
    transfstress,
    elasticitybasis,
    facequads,
    mergedmesh,
    x -> x[1] ≈ 0.0,
    [1.0, 0.0],
)

CutCellDG.assemble_penalty_displacement_component_bc!(
    sysmatrix,
    sysrhs,
    x -> 0.0,
    elasticitybasis,
    facequads,
    stiffness,
    mergedmesh,
    x -> x[2] ≈ 0.0,
    [0.0, 1.0],
    penalty,
)
CutCellDG.assemble_penalty_displacement_component_transformation_linear_form!(
    sysrhs,
    transfstress,
    elasticitybasis,
    facequads,
    mergedmesh,
    x -> x[2] ≈ 0.0,
    [0.0, 1.0],
)

CutCellDG.assemble_penalty_displacement_component_bc!(
    sysmatrix,
    sysrhs,
    x -> 0.0,
    elasticitybasis,
    facequads,
    stiffness,
    mergedmesh,
    x -> x[1] ≈ L,
    [1.0, 0.0],
    penalty,
)
CutCellDG.assemble_penalty_displacement_component_transformation_linear_form!(
    sysrhs,
    transfstress,
    elasticitybasis,
    facequads,
    mergedmesh,
    x -> x[1] ≈ L,
    [1.0, 0.0],
)

matrix = CutCellDG.sparse_displacement_operator(sysmatrix, mergedmesh)
rhs = CutCellDG.displacement_rhs_vector(sysrhs, mergedmesh)

sol = matrix \ rhs
disp = reshape(sol, 2, :)

solver = APS.PlaneSolver(L, W, R, lambda1, mu1, lambda2, mu2, theta0)

err = mesh_L2_error(
    disp,
    x -> APS.displacement_field(solver, x),
    elasticitybasis,
    cellquads,
    mergedmesh,
)

@test maximum(err) < 1e2eps()


text:
using Visualize: CairoWindow, AbstractCairoWindow, Canvas, Scene, ColorBuffer, DepthBuffer, clip2pixel_space, CairoContext
using Visualize.JLRasterization: JLRasterizer, edge_function, interpolate,standard_transparency, broadcastmin, broadcastmax
using FileIO, GeometryTypes, Colors
using Visualize: normalmesh, Light, Shading, MeshUniforms, Camera

x = CairoWindow()

function (r::JLRasterizer{Vert, Args, FragN}){Vert, Args, FragN}(
        window::AbstractCairoWindow, vertex_array::AbstractArray{Vert}, uniforms::Args
    )
    canvas = window[Canvas]
    framebuffers = canvas[ColorBuffer]; depthbuffer = canvas[DepthBuffer]
    resolution = Vec2f0(size(framebuffers[1]))
    # hoisting out functions... Seems to help inference a bit. Or not?
    vshader = r.vertexshader
    gshader = r.geometryshader
    fshader = r.fragmentshader
    FragNVal = Val{FragN}()
    fragments_drawn = 0
    clip_triangles = NTuple{3, Vec2f0}[]
    for face in vertex_array
        vertex_stage = map(reverse(face)) do f
            vshader(f, uniforms...)
        end
        geom_stage = if isa(r.geometryshader, Void)
            (vertex_stage,)
        else
            reset!(r.geometry_view)
            gshader(r.emit, vertex_stage, uniforms...)
            r.geometry_view.view
        end
        for geom_face in geom_stage
            fdepth = map(geom_face) do vert
                fv = first(vert)
                p = clip2pixel_space(fv, resolution)
                p[Vec(1, 2)], p[3]
            end
            f = map(first, fdepth)

            depths = map(last, fdepth)
            vertex_out = map(last, geom_face)
            # Bounding rectangle
            mini = floor(Int, max.(reduce(broadcastmin, f), 1f0))
            maxi = ceil(Int, min.(reduce(broadcastmax, f), resolution))
            area = edge_function(f[1], f[2], f[3])
            fragment_count = fragments_drawn
            for y = mini[2]:maxi[2]
                for x = mini[1]:maxi[1]
                    p = Vec(Float64(x), Float64(y))
                    w = Vec(
                        edge_function(f[2], f[3], p),
                        edge_function(f[3], f[1], p),
                        edge_function(f[1], f[2], p)
                    )
                    yi, xi = (y), (x)
                    if all(w-> w <= 0f0, w) && checkbounds(Bool, framebuffers[1], yi, xi)
                        bary = w / area
                        depth = bary[1] * depths[1] + bary[2] * depths[2] + bary[3] * depths[3]

                        if depth <= depthbuffer[yi, xi]
                            depthbuffer[yi, xi] = depth
                            fragment_in = interpolate(bary, vertex_out, FragNVal)
                            fragment_out = fshader(fragment_in, uniforms...)
                            for i = eachindex(fragment_out)
                                src_color = framebuffers[i][yi, xi]
                                dest_color = fragment_out[i]
                                fragments_drawn += 1
                                color = standard_transparency(
                                    RGBA{Float32}(src_color),
                                    RGBA{Float32}(dest_color[1], dest_color[2], dest_color[3], dest_color[4])
                                )
                                framebuffers[i][yi, xi] = mapc(x-> clamp(x, 0, 1), color)
                            end
                        end
                    end
                end
            end
            if fragments_drawn != fragment_count # only add triangle if it appeared on screen!
                push!(clip_triangles, f)
            end
        end
    end
    cr = canvas[CairoContext]
    Cairo.save(cr)
    Cairo.set_source_rgba(cr, 0.5, 0, 0, 1.0)
    for face in clip_triangles
        for (i, point) in enumerate(face)
            if i == 1
                Cairo.move_to(cr, point[1], point[2])
            else
                Cairo.line_to(cr, point[1], point[2])
            end
        end
        # Cairo.close_path(cr)
    end
    # Cairo.fill(cr)
    Cairo.clip(cr)
    # Cairo.reset_clip(cr)
    Cairo.set_source_surface(cr, Cairo.CairoImageSurface(rotl90(framebuffers[1])), 0, 0)
    Cairo.paint(cr)
    Cairo.restore(cr)
    println("fragments drawn: ", fragments_drawn)
    return
end

catmesh = normalmesh(load(Pkg.dir("GLVisualize", "assets", "cat.obj")))
catmesh = Base.view(
    reinterpret(Visualize.VertexN, catmesh.parent), catmesh.indexes[1]
)

args = MeshUniforms()

light = Light(
    Vec3f0(10),
    Vec3f0(0.0),
    Vec3f0(0.9),
    1.0f0,
    Vec3f0(0.2),
    1.0f0
)
shading = Shading(
    Vec3f0(0.0),
    Vec3f0(0.3),
    8.0f0
)
scene = Visualize.SceneUniforms(x[Camera])

uniforms = (scene, light, shading, args)
raster = JLRasterizer(
    catmesh,
    uniforms,
    Visualize.vert_mesh, Visualize.frag_mesh
)

sf = canvas[Visualize.CairoSurface]
canvas[Visualize.CairoContext] = Cairo.CairoContext(sf)
canvas[DepthBuffer] = zeros(Float32, 500, 500)

Cairo.save(cr);
Cairo.set_source_rgb(cr,1.0,1.0,1.0);    # light gray
Cairo.rectangle(cr,0.0,0.0, get(x, Visualize.Resolution)...); # background
Cairo.fill(cr);
Cairo.restore(cr);
Cairo.save(cr);
Cairo.reset_clip(cr)
raster(x, catmesh, uniforms)
Cairo.write_to_png(canvas[Visualize.CairoSurface], "test.png")
# save("test.png", canvas[Visualize.ColorBuffer][1])

#
# c = CairoPDFSurface("test.pdf", resolution...);
# cr = CairoContext(c);
# Cairo.save(cr);
# set_source_rgb(cr, 1.0, 1.0, 1.0);    # light gray
# rectangle(cr, 0.0, 0.0, resolution...); # background
# fill(cr);
# restore(cr);


text:
# test Pose2DPoint2D constraint evaluation function

using LinearAlgebra, Statistics
using Test
using RoME
# , IncrementalInference, Distributions

##

@testset "test pose and point combinations..." begin

##

N = 100
fg = initfg()

initCov = Matrix(Diagonal([0.03;0.03;0.001]))
odoCov = Matrix(Diagonal([3.0;3.0;0.01]))

# Some starting position
v1 = addVariable!(fg, :x0, Pose2, N=N)
initPosePrior = PriorPose2(MvNormal(zeros(3), initCov))
f1  = addFactor!(fg,[v1], initPosePrior)

@test Pose2Pose2(MvNormal(randn(2), Matrix{Float64}(LinearAlgebra.I, 2,2))) !== nothing
@test Pose2Pose2(MvNormal(randn(2), Matrix{Float64}(LinearAlgebra.I, 2,2))) !== nothing
@test Pose2Pose2(MvNormal(randn(2), Matrix{Float64}(LinearAlgebra.I, 2,2))) !== nothing

# and a second pose
v2 = addVariable!(fg, :x1, Pose2, N=N)
ppc = Pose2Pose2(MvNormal([50.0;0.0;pi/2], odoCov))
f2 = addFactor!(fg, [:x0;:x1], ppc)

# test evaluation of pose pose constraint
pts = approxConv(fg, :x0x1f1, :x1)

pts[3,:] .= TU.wrapRad.(pts[3,:])
@show mv = Statistics.mean(pts,dims=2)
@test norm(mv[1:2] - [50.0;0.0]) < 10.0
@test abs(mv[3] - pi/2) < 0.5

##

# @show ls(fg)
tree, smt, hist = solveTree!(fg)
# ensureAllInitialized!(fg)
# tree = wipeBuildNewTree!(fg)
# inferOverTreeR!(fg, tree,N=N)
# inferOverTree!(fg, tree, N=N)

# test post evaluation values are correct
pts = getVal(fg, :x0)
pts[3,:] .= TU.wrapRad.(pts[3,:])

@test norm(Statistics.mean(pts, dims=2)[1:2] - [0.0;0.0]) < 10.0
@test abs(Statistics.mean(pts, dims=2)[3]) < 0.5

pts = getVal(fg, :x1)
pts[3,:] .= TU.wrapRad.(pts[3,:])

@test norm(Statistics.mean(pts, dims=2)[1:2]-[50.0;0.0]) < 10.0
@test abs(Statistics.mean(pts, dims=2)[3]-pi/2) < 0.5

# check that yaw is working
v3 = addVariable!(fg, :x2, Pose2, N=N)
ppc = Pose2Pose2(MvNormal([50.0;0.0;0.0], odoCov))
f3 = addFactor!(fg, [v2;v3], ppc)

solveTree!(fg)
solveTree!(fg)

# test post evaluation values are correct
pts = getVal(fg, :x0)
@test norm(Statistics.mean(pts, dims=2)[1:2]-[0.0;0.0]) < 20.0
pts[3,:] .= TU.wrapRad.(pts[3,:])
@test abs(Statistics.mean(pts, dims=2)[3]) < 0.5

pts = getVal(fg, :x1)
@test norm(Statistics.mean(pts, dims=2)[1:2]-[50.0;0.0]) < 20.0
pts[3,:] .= TU.wrapRad.(pts[3,:])
@test abs(Statistics.mean(pts, dims=2)[3] - pi/2) < 0.5

pts = getVal(fg, :x2)
@test norm(Statistics.mean(pts, dims=2)[1:2]-[50.0;50.0]) < 20.0
pts[3,:] .= TU.wrapRad.(pts[3,:])
@test abs(Statistics.mean(pts, dims=2)[3]-pi/2) < 0.5

println("test bearing range evaluations")

# new landmark
l1 = addVariable!(fg, :l1, Point2, N=N)
# and pose to landmark constraint
rhoZ1 = norm([10.0;0.0])
ppr = Pose2Point2BearingRange(Uniform(-pi,pi),Normal(rhoZ1,1.0))
f4 = addFactor!(fg, [:x0;:l1], ppr)

pts = approxConv(fg, :x0l1f1, :l1, N=N)
# pts = evalFactor(fg, f4, l1.index, N=N)
## res = zeros(2)
## meas = getSample(ppr)
## ppr(res,nothing,1,meas,xi,lm )

# all points should lie in a ring around 0,0
@test sum(sqrt.(sum(pts.^2, dims=1 )) .< 5.0) == 0
@test sum(sqrt.(sum(pts.^2, dims=1 )) .< 15.0) == N

pts = approxConv(fg, :x0l1f1, :x0)
# pts = evalFactor(fg, f4, v1.index, N=200)
# @show sum(sqrt(sum(pts.^2, 1 )) .< 5.0)
# @test sum(sqrt(sum(pts.^2, 1 )) .< 5.0) == 0

# add a prior to landmark
pp2 = PriorPoint2(MvNormal([10.0;0.0], Matrix(Diagonal([1.0;1.0]))))

f5 = addFactor!(fg,[l1], pp2)
pts = approxConv(fg, f5.label, l1.label)
# pts = evalFactor(fg, f5, l1.label)

@test norm(Statistics.mean(pts, dims=2)[:]-[10.0;0.0]) < 5.0

# println("test Pose2D plotting")

# drawPoses(fg);
# drawPosesLandms(fg);

# using KernelDensityEstimate
# using Gadfly
#
# pts = getVal(fg, :l1)
#
# p1= kde!(pts)
# plotKDE( p1 , dimLbls=["x";"y";"z"])

# plotKDE( [marginal(p1c,[1;2]);marginal(p1,[1;2])] , dimLbls=["x";"y";"z"],c=["red";"black"],levels=3)
# p1c = deepcopy(p1)

# plotKDE( marginal(getBelief(fg, :x2),[1;2]) , dimLbls=["x";"y";"z"])
#
# axis = [[1.5;3.5]';[-1.25;1.25]';[-1.0;1.0]']
# draw( PDF("/home/dehann/Desktop/test.pdf",30cm,20cm),
#       plotKDE( p1, dimLbls=["x";"y";"z"], axis=axis)  )
#

end

#


text:
using Revise
using PyPlot
using PyCall
using Statistics
using LinearAlgebra
using Infiltrator
using Bem2d

"""
    plotbcsUeff(els, titlestring)

Plot boundary conditions and Ueff
"""
function plotbcsUeff(els, bcs, Ueff, suptitlestring)
    figure(figsize=(20, 20))
    # subplot(2, 2, 1)
    # quiver(els.xcenter[1:els.endidx], els.ycenter[1:els.endidx],
    #        bcs[1:2:end], bcs[2:2:end])
    # gca().set_aspect("equal")
    # title("bcs")

    subplot(2, 2, 1)
    plot(bcs[1:2:end], "rx", label="x")
    plot(bcs[2:2:end], "b+", label="y")
    legend()
    title("bcs")

    subplot(2, 2, 3)
    quiver(els.xcenter[1:els.endidx], els.ycenter[1:els.endidx],
           Ueff[1:2:end], Ueff[2:2:end])
    gca().set_aspect("equal")
    title("Ueff")

    subplot(2, 2, 4)
    plot(Ueff[1:2:end], "rx", label="x")
    plot(Ueff[2:2:end], "b+", label="y")
    legend()
    title("Ueff")

    suptitle(suptitlestring)
end

"""
    plotgeometry(els, titlestring)

els structure geometry plotting for diagnostics
"""
function plotgeometry(els, titlestring)
    scale = 500
    figure(figsize=(15,10))
    for i in 1:els.endidx-1
        plot([els.x1[i], els.x2[i]], [els.y1[i], els.y2[i]], "-g")
        quiver(els.xcenter[i], els.ycenter[i],
               els.xnormal[i], els.ynormal[i],
               width=1e-3, scale=1e2)
        text(els.xcenter[i] + scale * els.xnormal[i],
             els.ycenter[i] + scale * els.ynormal[i],
             string(i), fontsize=10)
    end
    title(titlestring)
    gca().set_aspect("equal")
end

"""
    flipnormals!(els, idx)

Flip the normals of the elements idx in els
"""
function flipnormals!(els, idx)
    els.xnormal[idx] *= -1
    els.ynormal[idx] *= -1
end

"""
    cloneels!(els, idx, label)

Create a copy of the elements with indices idx adding them to the
else structure and giving them the new name label.
"""
function cloneels!(els, idx, label)
    addelsez!(els, els.x1[idx], els.y1[idx], els.x2[idx], els.y2[idx], label)
end

"""
    layeredboxfault()

Comparing half-space and dislocation in a box solutions
"""
function layeredboxfault()
    close("all")
    PLOTGEOMETRY = false
    mu1 = 3e10
    nu1 = 0.25
    mu2 = mu1
    nu2 = nu1
    npts = 100
    offset = 100 # meters

    ###
    ### Single domain box BEM problem
    ###
    elsbox = Elements(Int(1e5))
    boxB = -30e3
    boxR = 30e3
    boxT = 0
    boxL = -30e3
    nside = 20
    nfault = 1
    x1, y1, x2, y2 = discretizedline(boxL, boxB, boxR, boxB, nside) # Bottom
    addelsez!(elsbox, x1, y1, x2, y2, "B")
    x1, y1, x2, y2 = discretizedline(boxR, boxB, boxR, boxT, nside) # Right hand side
    addelsez!(elsbox, x1, y1, x2, y2, "R")
    x1, y1, x2, y2 = discretizedline(boxR, boxT, boxL, boxT, nside) # Top
    addelsez!(elsbox, x1, y1, x2, y2, "T")
    x1, y1, x2, y2 = discretizedline(boxL, boxT, boxL, boxB, nside) # Left hand side
    addelsez!(elsbox, x1, y1, x2, y2, "L")
    x1, y1, x2, y2 = discretizedline(-10e3, -10e3, 0, 0, nfault)
    addelsez!(elsbox, x1, y1, x2, y2, "F")

    idxbox = getidxdict(elsbox)
    PLOTGEOMETRY && plotgeometry(elsbox, "Box boundaries and normals")

    # Fault contribution
    T_TB_F, H_TB_F = PUTC(slip2dispstress, elsbox,
                          [idxbox["B"] ; idxbox["R"] ; idxbox["T"] ; idxbox["L"]],
                          idxbox["F"], mu1, nu1)
    Fslip = [1; 1]; # y-direction slip only
    Uslip = T_TB_F * Fslip
    Tslip = H_TB_F * Fslip

    # Kernels and solve
    T_B_BRTL, _ = PUTC(slip2dispstress, elsbox, idxbox["B"],
                       [idxbox["B"] ; idxbox["R"] ; idxbox["T"] ; idxbox["L"]],
                       mu1, nu1)
    T_R_BRTL, _ = PUTC(slip2dispstress, elsbox, idxbox["R"],
                       [idxbox["B"] ; idxbox["R"] ; idxbox["T"] ; idxbox["L"]],
                       mu1, nu1)
    T_L_BRTL, _ = PUTC(slip2dispstress, elsbox, idxbox["L"],
                       [idxbox["B"] ; idxbox["R"] ; idxbox["T"] ; idxbox["L"]],
                       mu1, nu1)
    _, H_T_BRTL = PUTC(slip2dispstress, elsbox, idxbox["T"],
                       [idxbox["B"] ; idxbox["R"] ; idxbox["T"] ; idxbox["L"]],
                       mu1, nu1)

    bcsbox = zeros(8 * nside)
    bcsbox[1:2*nside] = -Uslip[1:2*nside] # Bottom
    bcsbox[2*nside+1:4*nside] = -Uslip[2*nside+1:4*nside] # Right
    bcsbox[4*nside+1:6*nside] = -Tslip[4*nside+1:6*nside] # Top
    bcsbox[6*nside+1:8*nside] = -Uslip[6*nside+1:8*nside] # Left

    THbox = [T_B_BRTL ; T_R_BRTL; H_T_BRTL ; T_L_BRTL]
    Ueffbox =THbox \ bcsbox
    # plotbcsUeff(elsbox, bcsbox, Ueffbox, "Box")
    xgrid, ygrid = obsgrid(boxL+offset, boxB+offset, boxR-offset, boxT-offset, npts)
    UTB, STB = constdispstress(slip2dispstress, xgrid, ygrid, elsbox,
                               [idxbox["B"] ; idxbox["R"] ; idxbox["T"] ; idxbox["L"]],
                               Ueffbox[1:2:end], Ueffbox[2:2:end], mu1, nu1)
    UF, SF = constdispstress(slip2dispstress, xgrid, ygrid, elsbox, idxbox["F"],
                             Fslip[1:2:end], Fslip[2:2:end], mu1, nu1)
    U1 = UTB .+ UF
    S1 = STB .+ SF
    plotfields(elsbox, reshape(xgrid, npts, npts), reshape(ygrid, npts, npts),
               U1, S1, "single layer")

    ###
    ### Layered box
    ###
    elslayer = Elements(Int(1e5))
    l1B = -15e3
    l1R = 30e3
    l1T = 0
    l1L = -30e3
    l2B = -30e3
    l2R = 30e3
    l2T = -15e3
    l2L = -30e3
    # nside = 100
    x1, y1, x2, y2 = discretizedline(l1L, l1B, l1R, l1B, nside) # Layer 1 bottom
    addelsez!(elslayer, x1, y1, x2, y2, "B1")
    x1, y1, x2, y2 = discretizedline(l1R, l1B, l1R, l1T, nside) # Layer 1 right
    addelsez!(elslayer, x1, y1, x2, y2, "R1")
    x1, y1, x2, y2 = discretizedline(l1R, l1T, l1L, l1T, nside) # Layer 1 top
    addelsez!(elslayer, x1, y1, x2, y2, "T1")
    x1, y1, x2, y2 = discretizedline(l1L, l1T, l1L, l1B, nside) # Layer 1 left
    addelsez!(elslayer, x1, y1, x2, y2, "L1")
    x1, y1, x2, y2 = discretizedline(l2L, l2B, l2R, l2B, nside) # Layer 2 bottom
    addelsez!(elslayer, x1, y1, x2, y2, "B2")
    x1, y1, x2, y2 = discretizedline(l2R, l2B, l2R, l2T, nside) # Layer 2 right
    addelsez!(elslayer, x1, y1, x2, y2, "R2")
    x1, y1, x2, y2 = discretizedline(l2L, l2T, l2R, l2T, nside) # Layer 2 top
    addelsez!(elslayer, x2, y2, x1, y1, "T2") # Note ordering!!!
    x1, y1, x2, y2 = discretizedline(l2L, l2T, l2L, l2B, nside) # Layer 2 left
    addelsez!(elslayer, x1, y1, x2, y2, "L2")
    x1, y1, x2, y2 = discretizedline(-10e3, -10e3, 0, 0, nfault) # Fault
    addelsez!(elslayer, x1, y1, x2, y2, "F")
    idxlayer = getidxdict(elslayer)
    plotgeometry(elslayer, "Layer boundaries and normals")

    # Build design matrix
    C1idx = [idxlayer["B1"] ; idxlayer["R1"] ; idxlayer["T1"] ; idxlayer["L1"]]
    C2idx = [idxlayer["B2"] ; idxlayer["R2"] ; idxlayer["T2"] ; idxlayer["L2"]]

    # Fault contribution
    T_C1_F, H_C1_F = PUTC(slip2dispstress, elslayer, C1idx, idxlayer["F"], mu1, nu1)
    UC1slip = T_C1_F * Fslip
    TC1slip = H_C1_F * Fslip

    # Boundaries shared by upper (C1) and lower (C2) layers
    T_B1_C1, H_B1_C1 = PUTC(slip2dispstress, elslayer, idxlayer["B1"], C1idx, mu1, nu1)
    T_T2_C2, H_T2_C2 = PUTC(slip2dispstress, elslayer, idxlayer["T2"], C2idx, mu2, nu2)

    # Upper layer (C1) only boundaries
    T_R1_C1, H_R1_C1 = PUTC(slip2dispstress, elslayer, idxlayer["R1"], C1idx, mu1, nu1)
    T_T1_C1, H_T1_C1 = PUTC(slip2dispstress, elslayer, idxlayer["T1"], C1idx, mu1, nu1)
    T_L1_C1, H_L1_C1 = PUTC(slip2dispstress, elslayer, idxlayer["L1"], C1idx, mu1, nu1)

    # Lower layer (C2) only boundaries
    T_B2_C2, H_B2_C2 = PUTC(slip2dispstress, elslayer, idxlayer["B2"], C2idx, mu2, nu2)
    T_R2_C2, H_R2_C2 = PUTC(slip2dispstress, elslayer, idxlayer["R2"], C2idx, mu2, nu2)
    T_L2_C2, H_L2_C2 = PUTC(slip2dispstress, elslayer, idxlayer["L2"], C2idx, mu2, nu2)

    # Place submatrices into larger matrix (zero slip sides)
    TH = zeros(16 * nside, 16 * nside)
    TH[1:2*nside, 1:8*nside] = T_B1_C1 # Displacements at B1 due to C1
    TH[1:2*nside, 8*nside+1:16*nside] = -T_T2_C2 # Displacements at T2 due to C2
    TH[2*nside+1:4*nside, 1:8*nside] = H_B1_C1 # Tractions at B1 due to C1
    TH[2*nside+1:4*nside, 8*nside+1:16*nside] = H_T2_C2 # Tractions at T2 due to C2
    TH[4*nside+1:6*nside, 1:8*nside] = T_R1_C1 # Displacements at R1 due to C1
    TH[6*nside+1:8*nside, 1:8*nside] = H_T1_C1 # Tractions at T1 due to C1
    TH[8*nside+1:10*nside, 1:8*nside] = T_L1_C1 # Displacements at L1 due to C1
    TH[10*nside+1:12*nside, 8*nside+1:16*nside] = T_B2_C2 # Displacments at B2 due to C2
    TH[12*nside+1:14*nside, 8*nside+1:16*nside] = T_R2_C2 # Tractions at R2 due to C2
    TH[14*nside+1:16*nside, 8*nside+1:16*nside] = T_L2_C2 # Tractions at L2 due to C2

    # BCS (from fault) and solve
    bcslayer = zeros(16 * nside)
    bcslayer[1:2*nside] = -UC1slip[1:2*nside] # Displacments at B1 due to F
    bcslayer[2*nside+1:4*nside] = -TC1slip[1:2*nside] # Tractions at B1 due to F
    bcslayer[4*nside+1:6*nside] = -UC1slip[2*nside+1:4*nside] # Displacments at R1 due to F
    bcslayer[6*nside+1:8*nside] = -TC1slip[4*nside+1:6*nside] # Tractions at T1 due to F
    bcslayer[8*nside+1:10*nside] = -UC1slip[6*nside+1:8*nside] # Displacements at L1 due to F
    Uefflayer = TH \ bcslayer
    plotbcsUeff(elslayer, log10.(abs.(bcslayer)), Uefflayer, "Two-layer")
    UeffT2 = Uefflayer[1:8*nside]
    UeffB2 = Uefflayer[8*nside+1:16*nside]

    # Volume visualization for two domain case
    Tidx = findall(x -> x > -15e3, ygrid)
    Bidx = findall(x -> x < -15e3, ygrid)
    Tx = xgrid[Tidx]
    Ty = ygrid[Tidx]
    Bx = xgrid[Bidx]
    By = ygrid[Bidx]
    UT2, ST2 = constdispstress(slip2dispstress, Tx, Ty, elslayer, C1idx,
                               UeffT2[1:2:end], UeffT2[2:2:end], mu1, nu1)
    UB2, SB2 = constdispstress(slip2dispstress, Bx, By, elslayer, C2idx,
                               UeffB2[1:2:end], UeffB2[2:2:end], mu2, nu2)
    UT2F, ST2F = constdispstress(slip2dispstress, Tx, Ty, elslayer, idxlayer["F"],
                                 Fslip[1:2:end], Fslip[2:2:end], mu1, nu1)
    UT2 = UT2 .+ UT2F
    ST2 = ST2 .+ ST2F
    U2 = [UB2 ; UT2] # Note B, T ordering
    S2 = [SB2 ; ST2] # Note B, T ordering
    plotfields(elslayer, reshape(xgrid, npts, npts), reshape(ygrid, npts, npts),
               U2, S2, "two-layer")
    plotfields(elslayer, reshape(xgrid, npts, npts), reshape(ygrid, npts, npts),
               U2 .- U1, S2 .- S1, "residuals")
end
layeredboxfault()


text:
# sse(x) gives the L2 norm of x
sse(x) = (x'*x)[1]

function levenberg_marquardt(f::Function, g::Function, x0; tolX=1e-8, tolG=1e-12, maxIter=100, lambda=100.0, show_trace=false)
	# finds argmin sum(f(x).^2) using the Levenberg-Marquardt algorithm
	#          x
	# The function f should take an input vector of length n and return an output vector of length m
	# The function g is the Jacobian of f, and should be an m x n matrix
	# x0 is an initial guess for the solution
	# fargs is a tuple of additional arguments to pass to f
	# available options:
	#   tolX - search tolerance in x
	#   tolG - search tolerance in gradient
	#   maxIter - maximum number of iterations
	#   lambda - (inverse of) initial trust region radius
	#   show_trace - print a status summary on each iteration if true
	# returns: x, J
	#   x - least squares solution for x
	#   J - estimate of the Jacobian of f at x

	# other constants
	const MAX_LAMBDA = 1e16 # minimum trust region radius
	const MIN_LAMBDA = 1e-16 # maximum trust region radius
	const MIN_STEP_QUALITY = 1e-3
	const MIN_DIAGONAL = 1e-6 # lower bound on values of diagonal matrix used to regularize the trust region step

	converged = false
	iterCt = 0
	x = x0
	delta_x = copy(x0)
	f_calls = 0
	g_calls = 0

	fcur = f(x)
	f_calls += 1
	residual = sse(fcur)

	# Maintain a trace of the system.
	tr = OptimizationTrace()
	if show_trace
		d = {"lambda" => lambda}
		os = OptimizationState(x, sse(fcur), iterCt, d)
		push!(tr, os)
		println(os)
	end

	while ( ~converged && iterCt < maxIter )
		J = g(x)
		g_calls += 1
		# we want to solve:
		#    argmin 0.5*||J(x)*delta_x + f(x)||^2 + lambda*||diagm(J'*J)*delta_x||^2
		# Solving for the minimum gives:
		#    (J'*J + lambda*DtD) * delta_x == -J^T * f(x), where DtD = diagm(sum(J.^2,1))
		# Where we have used the equivalence: diagm(J'*J) = diagm(sum(J.^2, 1))
		# It is additionally useful to bound the elements of DtD below to help
		# prevent "parameter evaporation".
		DtD = diagm(Float64[max(x, MIN_DIAGONAL) for x in sum(J.^2,1)])
		delta_x = ( J'*J + sqrt(lambda)*DtD ) \ -J'*fcur
		# if the linear assumption is valid, our new residual should be:
		predicted_residual = sse(J*delta_x + fcur)
		# check for numerical problems in solving for delta_x by ensuring that the predicted residual is smaller
		# than the current residual
		if predicted_residual > residual
			error("Error solving for delta_x: predicted residual increase.")
		end
		# try the step and compute its quality
		trial_f = f(x + delta_x)
		f_calls += 1
		trial_residual = sse(trial_f)
		# step quality = residual change / predicted residual change
		rho = (trial_residual - residual) / (predicted_residual - residual)

		if rho > MIN_STEP_QUALITY
			x += delta_x
			fcur = trial_f
			residual = trial_residual
			# increase trust region radius
			lambda = max(0.1*lambda, MIN_LAMBDA)
		else
			# decrease trust region radius
			lambda = min(10*lambda, MAX_LAMBDA)
		end
		iterCt += 1

		# show state
		if show_trace
			d = {"g(x)" => norm(J'*fcur, Inf), "dx" => delta_x, "lambda" => lambda}
			os = OptimizationState(x, sse(fcur), iterCt, d)
			push!(tr, os)
			println(os)
		end

		# check convergence criteria:
		# 1. Small gradient: norm(J^T * fcur, Inf) < tolG
		# 2. Small step size: norm(delta_x) < tolX
		if norm(J' * fcur, Inf) < tolG
			converged = true
		elseif norm(delta_x) < tolX*(tolX + norm(x))
			converged = true
		end
	end

	# give the user info about the stopping condition
	if ~converged
		println("Exceeded maximum number of iterations")
	end

	OptimizationResults("Levenberg-Marquardt", x0, x, sse(fcur), iterCt, converged, tr, f_calls, g_calls)
end


text:
include("../src/components/neteconomy_component.jl")

@testset "neteconomy" begin

    Precision = 1.0e-11
    T = 60
    f = readxlsx(joinpath(@__DIR__, "..", "data", "DICE2010_082710d.xlsx"))

    m = Model()

    set_dimension!(m, :time, MimiDICE2010.model_years)

    add_comp!(m, neteconomy, :neteconomy)

# Set the parameters that would normally be internal connection from their Excel values
    set_param!(m, :neteconomy, :YGROSS, read_params(f, "B92:BI92", T))
    set_param!(m, :neteconomy, :DAMFRAC, read_params(f, "B93:BI93", T))

# Load the rest of the external parameters
    p = dice2010_excel_parameters(joinpath(@__DIR__, "..", "data", "DICE2010_082710d.xlsx"))
    set_param!(m, :neteconomy, :cost1, p[:cost1])
    set_param!(m, :neteconomy, :MIU, p[:MIU])
    set_param!(m, :neteconomy, :expcost2, p[:expcost2])
    set_param!(m, :neteconomy, :partfract, p[:partfract])
    set_param!(m, :neteconomy, :pbacktime, p[:pbacktime])
    set_param!(m, :neteconomy, :S, p[:S])
    set_param!(m, :neteconomy, :l, p[:l])

# Run the one-component model
    run(m)

# Extract the generated variables
    ABATECOST   = m[:neteconomy, :ABATECOST]
    C           = m[:neteconomy, :C]
    CPC         = m[:neteconomy, :CPC]
    CPRICE      = m[:neteconomy, :CPRICE]
    I           = m[:neteconomy, :I]
    Y           = m[:neteconomy, :Y]
    YNET        = m[:neteconomy, :YNET]

# Extract the true values
    True_ABATECOST  = read_params(f, "B97:BI97", T)
    True_C          = read_params(f, "B125:BI125", T)
    True_CPC        = read_params(f, "B126:BI126", T)
    True_CPRICE     = read_params(f, "B134:BI134", T)
    True_I          = read_params(f, "B101:BI101", T)
    True_Y          = read_params(f, "B98:BI98", T)
    True_YNET       = read_params(f, "B95:BI95", T)

# Test that the values are the same
    @test maximum(abs, ABATECOST .- True_ABATECOST) ≈ 0. atol = Precision
    @test maximum(abs, C .- True_C) ≈ 0. atol = Precision
    @test maximum(abs, CPC .- True_CPC) ≈ 0. atol = Precision
    @test maximum(abs, CPRICE .- True_CPRICE) ≈ 0. atol = Precision
    @test maximum(abs, I .- True_I) ≈ 0. atol = Precision
    @test maximum(abs, Y .- True_Y) ≈ 0. atol = Precision
    @test maximum(abs, YNET .- True_YNET) ≈ 0. atol = Precision

end


text:
# This file is a part of Julia, but is derived from
# https://github.com/google/double-conversion which has the following license
#
# Copyright 2006-2014, the V8 project authors. All rights reserved.
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
#     * Redistributions of source code must retain the above copyright
#       notice, this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above
#       copyright notice, this list of conditions and the following
#       disclaimer in the documentation and/or other materials provided
#       with the distribution.
#     * Neither the name of Google Inc. nor the names of its
#       contributors may be used to endorse or promote products derived
#       from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
# OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

const kDoubleSignificandSize = 53

function filldigits32fixedlength(n1,requested_len,buffer,len)
    for i = (requested_len-1):-1:0
        buffer[len+i] = 0x30 + n1 % 10
        n1 = div(n1,10)
    end
    return len + requested_len
end

function filldigits32(n,buffer,len)
    n_len = 0
    while n != 0
        digit = n % 10
        n = div(n,10)
        buffer[len+n_len] = 0x30 + digit
        n_len += 1
    end
    i,j = len, len + n_len - 1
    while i < j
        buffer[i], buffer[j] = buffer[j], buffer[i]
        i += 1
        j -= 1
    end
    return len + n_len
end

function filldigits64fixedlength(n2,buffer,len)
    kTen7 = 10000000
    part2 = n2 % kTen7
    n2 = div(n2,kTen7)
    part0, part1 = divrem(n2,kTen7)
    len = filldigits32fixedlength(part0, 3, buffer, len)
    len = filldigits32fixedlength(part1, 7, buffer, len)
    len = filldigits32fixedlength(part2, 7, buffer, len)
    return len
end

function filldigits64(n3,buffer,len)
    kTen7 = 10000000
    part2 = n3 % kTen7
    n3 = div(n3,kTen7)
    part0, part1 = divrem(n3,kTen7)
    if part0 != 0
        len = filldigits32(part0, buffer, len)
        len = filldigits32fixedlength(part1, 7, buffer, len)
        len = filldigits32fixedlength(part2, 7, buffer, len)
    elseif part1 != 0
        len = filldigits32(part1, buffer, len)
        len = filldigits32fixedlength(part2, 7, buffer, len)
    else
        len = filldigits32(part2, buffer, len)
    end
    return len
end

function roundup(buffer, len, decimal_point)
    if len == 1
        buffer[1] = 0x31
        decimal_point = 1
        len = 2
        return len, decimal_point
    end
    buffer[len - 1] += 1
    for i = (len-1):-1:2
        buffer[i] != 0x30 + 10 && return len, decimal_point
        buffer[i] = 0x30
        buffer[i - 1] += 1
    end
    if buffer[1] == 0x30 + 10
        buffer[1] = 0x31
        decimal_point += 1
    end
    return len, decimal_point
end

function fillfractionals(fractionals, exponent,
                         fractional_count, buffer,
                         len, decimal_point)
    if -exponent <= 64
        point = -exponent
        for i = 1:fractional_count
            fractionals == 0 && break
            fractionals *= 5
            point -= 1
            digit = fractionals >> point
            buffer[len] = 0x30 + digit
            len += 1
            fractionals -= UInt64(digit) << point
        end
        if ((fractionals >> (point - 1)) & 1) == 1
            len, decimal_point = roundup(buffer, len, decimal_point)
        end
    else
        fract128 = UInt128(fractionals) << 64
        fract128 = shift(fract128,-exponent - 64)
        point = 128
        for i = 1:fractional_count
            fract128 == 0 && break
            fract128 *= 5
            point -= 1
            digit, fract128 = divrem2(fract128,point)
            buffer[len] = 0x30 + digit
            len += 1
        end
        if bitat(fract128,point - 1) == 1
            len, decimal_point = roundup(buffer, len, decimal_point)
        end
    end
    return len, decimal_point
end

low(x) = UInt64(x&0xffffffffffffffff)
high(x) = UInt64(x >>> 64)
bitat(x::UInt128,y) = y >= 64 ? (Int32(high(x) >> (y-64)) & 1) : (Int32(low(x) >> y) & 1)
function divrem2(x,power)
    h = high(x)
    l = low(x)
    if power >= 64
        result = Int32(h >> (power - 64))
        h -= UInt64(result) << (power - 64)
        return result, (UInt128(h) << 64) + l
    else
        part_low::UInt64 = l >> power
        part_high::UInt64 = h << (64 - power)
        result = Int32(part_low + part_high)
        return result, UInt128(l - (part_low << power))
    end
end
function shift(x::UInt128,amt)
    if amt == 0
      return x
    elseif amt == -64
        return x << 64
    elseif amt == 64
        return x >> 64
    elseif amt <= 0
        h = high(x); l = low(x)
        h <<= -amt
        h += l >> (64 + amt)
        l <<= -amt
        return (UInt128(h) << 64) + l
    else
        h = high(x); l = low(x)
        l >>= amt
        l += h << (64 - amt)
        h >>= amt
        return (UInt128(h) << 64) + l
    end
end

function trimzeros(buffer, len, decimal_point)
    while len > 1 && buffer[len - 1] == 0x30
        len -= 1
    end
    first_non_zero::Int32 = 1
    while first_non_zero < len && buffer[first_non_zero] == 0x30
        first_non_zero += 1
    end
    if first_non_zero != 1
        for i = first_non_zero:(len-1)
            buffer[i - first_non_zero + 1] = buffer[i]
        end
        len -= first_non_zero-1
        decimal_point -= first_non_zero-1
    end
    return len, decimal_point
end

function fastfixedtoa(v,mode,fractional_count,buffer)
    v = Float64(v)
    significand::UInt64 = _significand(v)
    exponent = _exponent(v)
    exponent > 20 && return false, 0, 0
    fractional_count > 20 && return false, 0, 0
    len = 1
    if exponent + kDoubleSignificandSize > 64
        kFive17 = divisor = Int64(5)^17
        divisor_power = 17
        dividend = significand
        if exponent > divisor_power
            dividend <<= exponent - divisor_power
            quotient = div(dividend,divisor)
            remainder = (dividend % divisor) << divisor_power
        else
            divisor <<= divisor_power - exponent
            quotient = div(dividend,divisor)
            remainder = (dividend % divisor) << exponent
        end
        len = filldigits32(quotient, buffer, len)
        len = filldigits64fixedlength(remainder, buffer, len)
        decimal_point = len-1
    elseif exponent >= 0
        significand <<= exponent
        len = filldigits64(significand, buffer, len)
        decimal_point = len-1
    elseif exponent > -kDoubleSignificandSize
        integrals = significand >> -exponent
        fractionals = significand - (integrals << -exponent)
        if integrals > 0xFFFFFFFF
            len = filldigits64(integrals,buffer,len)
        else
            len = filldigits32(integrals%UInt32,buffer,len)
        end
        decimal_point = len-1
        len, decimal_point = fillfractionals(fractionals,exponent,fractional_count,
                                             buffer,len, decimal_point)
    elseif exponent < -128
        len = 1
        decimal_point = -fractional_count
    else
        decimal_point = 0
        len, decimal_point = fillfractionals(significand,exponent,fractional_count,
                                             buffer,len, decimal_point)
    end
    len, decimal_point = trimzeros(buffer,len,decimal_point)
    buffer[len] = 0
    if (len-1) == 0
        decimal_point = -fractional_count
    end
    return true, len, decimal_point
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
# Tests async mode operation.
# Both server and client run in the same process.
# In this mode the transport can be in-memory.
include("srvr.jl")
include("clnt.jl")

function test_asyncsrvr()
    run_srvr(JuliaWebAPI.JSONMsgFormat(), JuliaWebAPI.ZMQTransport(SRVR_ADDR, ZMQ.REP, true), true)
    run_clnt(JuliaWebAPI.JSONMsgFormat(), JuliaWebAPI.ZMQTransport("127.0.0.1", 9999, ZMQ.REQ, false))
end

# run tests if invoked with run flag
!isempty(ARGS) && (ARGS[1] == "--runasyncsrvr") && test_asyncsrvr()


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, Distributions, Plots
	Random.seed!(0)
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing5.12"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	alpha = 0.05
	L(obs) = obs - (1-sqrt(alpha))
	U(obs) = obs + (1-sqrt(alpha))
end

# ╔═╡ ae6c7528-3af7-11eb-14a4-4bb1f292614b
begin
	mu = 5.57
	observation = rand(TriangularDist(mu-1,mu+1,mu))
	Text("Lower bound L: $(L(observation)))")
end

# ╔═╡ ae6cb47a-3af7-11eb-0264-bb161b7731a1
Text("Upper bound U: $(U(observation)))")

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing5.12"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═ae6c7528-3af7-11eb-14a4-4bb1f292614b
# ╠═ae6cb47a-3af7-11eb-0264-bb161b7731a1
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
#using Colors: Iterators
using Luxor, IterTools
using Colors
using Test

#import Luxor: text

## http://juliagraphics.github.io/Luxor.jl/stable/examples/#The-obligatory-%22Hello-World%22
function learn_hello_world()
    Drawing(1000, 1000, "hello-world.png")
    origin()
    background("black")
    sethue("red")
    fontsize(50)
    text("hello world")
    finish()
    preview()
end

# circle with hello
function learn_hello_circle()
    plt = @png begin
        fontsize(50)
        Luxor.circle(Point(0,0), 150, :stroke)
        text("hello world", halign=:center, valign=:middle)
    end

    return plt
end

# The arrow thing
function learn_hello_arrow()
    plt = @svg begin
        sethue("red")
        randpoint = Point(rand(-200:200), rand(-200:200))
        Luxor.circle(randpoint, 2, :fill)
        sethue("black")
        foreach(
            f -> arrow(f, between(f, randpoint, 0.1), arrowheadlength=6),
            first.(collect(Table(fill(20, 15), fill(20, 15))))
        )
    end

    return plt
end

# @draw, @drawsvg
# Display in Code, Jupyter, Pluto

function learn_hello_julia_sector()
    plt = @png begin
        setopacity(0.85)
        steps = 20
        gap = 2
        for (n, θ) in enumerate(range(0, step=2π/steps, length=steps))
            sethue([Luxor.julia_green,
                Luxor.julia_red,
                Luxor.julia_purple,
                Luxor.julia_blue][mod1(n, 4)])
            sector(Point(0, 0), 50, 250 + 2n, θ, θ + 2π/steps - deg2rad(gap), :fill)
        end
    end

    return plt
end

function learn_hello_pluto()
    S = 15
    plt = @svg begin
        background("slateblue4")
        # Notice O = Point(0,0)

        points = between.(O + (0,S), O - (0,S), range(0,1,length=3)) # broadcast
        hues = [Luxor.julia_red, Luxor.julia_purple, Luxor.julia_green]
        # `enumerate` can be replaced with `pairs`
        for (i, point) in enumerate(points)
            sethue(hues[i])
            Luxor.circle(point, S, :fill)
        end
    end 100 100

    return plt
end

function learn_hello_julia_logos()
    plt = @svg begin
        origin()
        background("white")

        for theta in range(0, step=pi/8, length=16)
            gsave()
            scale(0.2)
            rotate(theta)
            translate(350,0)
            julialogo(action = :fill, bodycolor = randomhue())
            grestore()
        end

        gsave()
        scale(0.3)
        juliacircles()
        grestore()

        translate(150,-150)
        scale(0.3)
        julialogo()

    end 600 400

    return plt
end

function triangle(points, degree, color)
    sethue(color)
    poly(points, :fill)
end

function sierpinski(points, degree, colors)
    triangle(points, degree, colors[degree])
    if degree > 1
        p1, p2, p3 = points
        sierpinski(
            [p1, midpoint(p1,p2), midpoint(p1,p3)],
            degree - 1,
            colors
        )
        sierpinski(
            [p2, midpoint(p1,p2), midpoint(p2,p3)],
            degree - 1,
            colors
        )
        sierpinski(
            [p3, midpoint(p3,p2), midpoint(p1,p3)],
            degree - 1,
            colors
        )

    end
end

function learn_hello_sierpinski_triangle()
    depth = 5
    colors = distinguishable_colors(depth)

    plt = @svg begin
        background("white")
        origin()

        Luxor.circle(Point(0,0), 75, :clip)
        points = ngon(Point(0,0), 150, 3, -pi/2, vertices=true)
        sierpinski(points, depth, colors)

    end

    return plt
end

function learn_hello_jupyter(r, m)
    plt = @svg begin
        origin()
        background("white")
        setcolor"antiquewhite"

        squircle(O, 124, 124, :fill)
        squircle(O, 124, 124, :stroke)
        ls = range(1, stop=120, length=18)
        @show ls
        setline(0.5)
        for (n,pt) in enumerate(ngon(O,100,3,vertices=true))
            @layer begin
                rotate(pi/6 + n * 2/3*pi)
                c = m
                for (s,e) in partition(ls, 2, 1)
                    @show (s,e)
                    sethue([Luxor.darker_red, Luxor.darker_green, Luxor.darker_purple][n])
                    sector(O, s, e, -pi/8 + r/10, pi * c, 15, :fill)
                    setcolor"white"
                    sector(O, s, e, -pi/8 + r/10, pi * c, 15, :stroke)
                    c *= m
                end
            end
        end
    end 300 300

    return plt
end

## http://juliagraphics.github.io/Luxor.jl/stable/tutorial/
function learn_tutorial_first_steps_1()
    @svg begin
        text("Hello world")
        circle(O, 200, :stroke)
    end
end

function learn_tutorial_first_steps_2()
    plt = @svg begin
        text("Hello again, world!", Point(0,250))
        circle(O, 200, :fill)
        rulers()
    end

    return plt
end

function dot(P; r = 2)
    circle(P, r, :fill)
end

## http://juliagraphics.github.io/Luxor.jl/stable/tutorial/#Euclidean-eggs
function learn_euclidean_eggs()
    plt = @svg begin
        radius = 80
        setdash("dot")
        sethue("gray30") # does not affect the current opacity setting
        A, B = [Point(x,0) for x in [-radius, radius]]
        line(A, B, :stroke)

        # Circle O
        # O == Point(0,0)
        circle(O, radius, :stroke)

        label("A", :NW, A)
        label("O", :N,  O)
        label("B", :NE, B)

        dot.([A, O, B]) # black dots
        circle.([A, B], 2*radius, :stroke) # big circle

        # Intersection of circle A, B
        count, C, D = intersectionlinecircle(Point(0, -2*radius), Point(0, 2*radius), A, 2*radius)

        if 2 == count
            #circle.([C,D], 2, :fill)
            dot.([C,D])
            label.(["C", "D"], :N, [C, D])
        end

        ## Upper circle
        # C1, the tutorial use intersection but it is just `radius` above O
        C1 = O + Point(0,-radius)
        dot(C1,r=3)
        label("C1", :N, C1)

        # ip1, ip2

        src = B; dst = C1
        count, I, J = intersectionlinecircle(src, dst, src, 2*radius)
        dot.([I,J])
        ip = distance(C1,I) < distance(C1,J) ? I : J
        label("ip1", :N, ip)
        ip1 = ip

        src = A; dst = C1
        count, I, J = intersectionlinecircle(src, dst, src, 2*radius)
        dot.([I,J])
        ip = distance(C1,I) < distance(C1,J) ? I : J
        label("ip2", :N, ip)
        ip2 = ip

        circle(C1, distance(C1,ip2), :stroke)

        ## Eggs at the ready
        setline(5)
        setdash("solid")

        arc2r( B,   A, ip1, :path)
        arc2r(C1, ip1, ip2, :path)
        arc2r( A, ip2,   B, :path)
        arc2r( O,   B,   A, :path)

        strokepreserve() # stokepath()
        setopacity(0.6)
        sethue("ivory")
        fillpath() # fillpreserve(), preserve the path after fill
    end 700 700

    return plt
end

function euclidean_egg(radius, action=:none)
    A, B = [Point(x,0) for x in [-radius, radius]]
    count, C, D = intersectionlinecircle(Point(0, -2*radius), Point(0, 2*radius), A, 2*radius)

    C1 = O + Point(0,-radius)

    src = B; dst = C1
    count, I, J = intersectionlinecircle(src, dst, src, 2*radius)
    ip = distance(C1,I) < distance(C1,J) ? I : J
    ip1 = ip

    src = A; dst = C1
    count, I, J = intersectionlinecircle(src, dst, src, 2*radius)
    ip = distance(C1,I) < distance(C1,J) ? I : J
    ip2 = ip

    newpath()
    arc2r( B,   A, ip1, :path)
    arc2r(C1, ip1, ip2, :path)
    arc2r( A, ip2,   B, :path)
    arc2r( O,   B,   A, :path)
    closepath()

    do_action(action)
end

function learn_euclidean_eggs_2()
    plt = @svg begin
        setopacity(0.7)
        for theta in range(0, step=pi/6, length=12)
            @layer begin # gsave() ... grestore()
                rotate(theta)
                translate(0,-150)
                euclidean_egg(50, :path)
                setline(10)
                randomhue()
                fillpreserve()

                randomhue()
                strokepath()
            end
        end
    end 800 800
end

function learn_euclidean_eggs_spiky()
    plt = @svg begin
        setlinejoin("round")
        euclidean_egg(160, :path)

        path = pathtopoly()

        # pathtopoly() -> Vector{Vector{Point}}
        # the first item is all points but the last one is not connected to
        # the first point
        # the second item is 1 point which is the starting point
        # by combining the 2, we crate a closed loop

        pgon = collect(Iterators.flatten(path))
        # pgon = vcat(path[1], path[2]) # same result as above

        pc = polycentroid(pgon)
        circle(pc, 5, :fill)

        for pt in 1:2:length(pgon)
            pgon[pt] = between(pc, pgon[pt], 0.5)
        end

        poly(pgon, :stroke)
    end

    return plt
end

function learn_euclidean_eggs_offsetpoly()
    plt = @svg begin
        euclidean_egg(80, :path)
        path = pathtopoly()
        ## offsetpoly() does not like endpoint connects to first point
        #pgon = collect(Iterators.flatten(path))
        pgon = first(path)
        pc = polycentroid(pgon)

        for pt in 1:2:length(pgon)
            pgon[pt] = between(pc, pgon[pt], 0.8)
        end
        poly(pgon, :stroke)

        for i in 30:-3:1
            randomhue()
            op = offsetpoly(pgon, i)
            poly(op, :stroke, close=true)
        end
    end

    return plt
end

function learn_euclidean_eggs_clipping()
    plt = @svg begin
        setopacity(0.5)

        egg(a) = euclidean_egg(150, a)
        sethue("gold")
        egg(:fill)
        egg(:clip)
        @layer begin
            for i in 360:-4:1
                sethue(Colors.HSV(i, 1.0, 0.8))
                rotate(pi/90)
                ngon(O, i, 5, 0, :stroke)
            end
        end
        clipreset()
        setopacity(1)
        setline(5)
        sethue("red")
        egg(:stroke)
    end

    return plt
end

## http://juliagraphics.github.io/Luxor.jl/stable/basics/
function learn_basics()
    P = Point(12,13)
    @test 12 == P.x
    @test 13 == P.y
    @test P.x isa Float64
    Q = Point(4,5)
    @test P.x + Q.x == (P + Q).x
    @test P.y + Q.y == (P + Q).y

    @test 10 * P.x == (10*P).x
    @test 10 * P.y == (10*P).y
end

## https://discourse.julialang.org/t/luxor-questions/41946/9
ytext(t; kwargs...) = ytext(t, O, yflip; kwargs...)
ytext(t, xpos, ypos; kwargs...) = ytext(t, Point(xpos, ypos); kwargs...)
function ytext(t, pt::Point; kwargs...)
    @layer begin
        translate(pt) # must translate first
        scale(1,-1)
        text(t; kwargs...)
    end
end

function learn_basics_triangle()
    plt = @svg begin
        background("white")
        origin()

        # @layer begin
        #     #rulers()
        #     points = ngon(Point(0,0), 250, 3, 90/360*2pi, :fill, vertices=true)
        #     poly(points, :stroke, close=true)
        # end
        @layer begin
            scale(1,-1)
            rulers()
            sethue("gold")
            # positive angle from +x to +y
            points = ngon(Point(0,0), 250, 5, 90/360*2pi, :fill, vertices=true)
            poly(points, :stroke, close=true)
        end

    end 500 500

    return plt
end
function learn_basics_box(;text=Luxor.text)
    plt = @svg begin
        scale(1,-1)
        rulers()
        # circle(O, 20, :stroke)
        l = 50
        theta = 2pi/l
        for dist in range(0, stop=100, length=l )
            #@show theta

            #rotate(theta)
            p = O + (dist,0)
            box(p, 20, 20, :stroke)

            ## https://discourse.julialang.org/t/luxor-questions/41946/9
            # @layer begin
            #     translate(p)
            #     scale(1,-1)
            #     text("$dist", valign=:middle, halign=:left)
            # end

            #text("$dist", p, valign=:middle, halign=:left)
            rotate(theta)
        end

        #text("TEST", Point(100,100))
    end 500 500

    return plt
end

## A few examples
# plt = learn_hello_world();               display(plt)
# plt = learn_hello_circle();              display(plt)
# plt = learn_hello_arrow();               display(plt)
# plt = learn_hello_julia_sector();        display(plt)
# plt = learn_hello_pluto();               display(plt)
# plt = learn_hello_julia_logos();         display(plt)
# plt = learn_hello_sierpinski_triangle(); display(plt)
# plt = learn_hello_jupyter(2, 0.95);      display(plt)
# TODO: Not finished yet

## Tutorial
# plt = learn_tutorial_first_steps_1();    display(plt)
# plt = learn_tutorial_first_steps_2();    display(plt)
# plt = learn_euclidean_eggs();            display(plt)
# plt = learn_euclidean_eggs_2();          display(plt)
# plt = learn_euclidean_eggs_spiky();      display(plt)
# plt = learn_euclidean_eggs_offsetpoly(); display(plt)
# plt = learn_euclidean_eggs_clipping();   display(plt)

## http://juliagraphics.github.io/Luxor.jl/stable/basics/

@testset "Base" begin
    learn_basics()
end

#plt = learn_basics_box(text=Luxor.text); display(plt)
plt = learn_basics_box(text=ytext); display(plt)
#plt = learn_basics_triangle(); display(plt)

nothing


text:
import accrue.cryptoerase.runtime.Condition;

class C {
    static void foo(Object c) {}

    public static void main(String[] args) {
	final Condition c = new Condition();
	foo(c);
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
using Memento
export m_debug, m_info, m_notice, m_warn, m_error

remove_handler(get_logger(Mocha), "console")
add_handler(get_logger(Mocha),
            DefaultHandler(STDOUT,
                           DefaultFormatter("[{date} | {level} | {name}]: {msg}")))
set_level(get_logger(Mocha), "info")

function m_debug(msg :: AbstractString...)
  debug(get_logger(Mocha), prod(msg))
end

function m_info(msg :: AbstractString...)
  info(get_logger(Mocha), prod(msg))
end

function m_notice(msg :: AbstractString...)
  notice(get_logger(Mocha), prod(msg))
end

function m_warn(msg :: AbstractString...)
  warn(get_logger(Mocha), prod(msg))
end

function m_error(msg :: AbstractString...)
  error(get_logger(Mocha), prod(msg))
end


text:
# Source this script as e.g.
#
#     include("PATH/TO/devrepl.jl")
#
# from *any* Julia REPL or run it as e.g.
#
#     julia -i --banner=no PATH/TO/devrepl.jl
#
# from anywhere. This will change the current working directory and
# activate/initialize the correct Julia environment for you.
#
# You may also run this in vscode to initialize a development REPL
#
using Pkg
using Downloads: download

cd(@__DIR__)
Pkg.activate("test")

function _instantiate()
    installorg_script = joinpath("..", "scripts", "installorg.jl")
    if !isfile(installorg_script)
        @warn "$(@__DIR__) should be inside the JuliaQuantumControl development environment. See https://github.com/JuliaQuantumControl/JuliaQuantumControl#readme"
        installorg_script = download(
            "https://raw.githubusercontent.com/JuliaQuantumControl/JuliaQuantumControl/master/scripts/installorg.jl",
        )
    end
    if !isfile(joinpath("..", ".JuliaFormatter.toml"))
        download(
            "https://raw.githubusercontent.com/JuliaQuantumControl/JuliaQuantumControl/master/.JuliaFormatter.toml",
            ".JuliaFormatter.toml"
        )
    end
    include(installorg_script)
    eval(:(installorg()))
end

if !isfile(joinpath("test", "Manifest.toml"))
    _instantiate()
end
include("test/init.jl")

if abspath(PROGRAM_FILE) == @__FILE__
    help()
end


text:
######################################################################
# policies.jl
# Does basic policy stuff.
# Policy's don't need to remember vehicle's max_step,
#  as the action function ensures the step is normalized to max_step.
######################################################################

abstract type Policy end

"""
`action(m, x, o, f, p)`

Returns an action given SearchDomain `m`, vehicle `x`, observation `o`, filter `f`, and policy `p`.
"""
function action(m::SearchDomain, x::Vehicle, o, f::AbstractFilter,p::Policy)
	return error(typeof(p), " does not yet implement action.")
end

# TODO: should ensure the heading is not too much either
function normalize(a::Action, x::Vehicle)
	ax = a[1]; ay = a[2]
	den = sqrt(ax*ax + ay*ay)
	return x.max_step * ax / den, x.max_step * ay / den, a[3]
end

function reset!(p::Policy) end

# to make action lists
include("action_list.jl")

######################################################################
# greedy policies
######################################################################
include("greedy.jl")
include("greedyp.jl")
include("greedyp2.jl")
include("greedy2.jl")

######################################################################
# GaussianMPC
######################################################################
include("gaussian_mpc.jl")

######################################################################
# Waypoint policy
######################################################################
include("waypoint.jl")

######################################################################
# Random policy
######################################################################
type RandomPolicy <: Policy end

function action(m::SearchDomain, x::Vehicle, o, f::AbstractFilter, p::RandomPolicy)
	ax = rand() - 0.5
	ay = rand() - 0.5
	ah = 10.0
	if rand() < 0.5
		ah = -10.0
	end
	return normalize((ax,ay,ah), x)
end

# random policy for rotate for bearing method
include("random2.jl")

######################################################################
# Sit policy: Just sit and don't move
######################################################################
type SitPolicy <: Policy end

action(m::SearchDomain, x::Vehicle, o, f::AbstractFilter, p::SitPolicy) = (0.0,0.0,0.0)

######################################################################
# Spin Policy: Just sit there and spin
######################################################################
type SpinPolicy <: Policy end

function action(m::SearchDomain, x::Vehicle, o, f::AbstractFilter, p::SpinPolicy)
	return (0.0, 0.0, 10.0)
end

######################################################################
# CirclePolicy: Move orthogonally to last measurement

# This ends up tracing a circle around the jammer
######################################################################
type CirclePolicy <: Policy
	last::Action

	CirclePolicy() = new( (0.0, 0.0, 0.0) )
end

# Remembers last action to ensure we follow same direction around circle
# Otherwise, it will "chatter" back and forth
# TODO: take into account distance from edge
function action(m::SearchDomain, x::Vehicle, o, f::AbstractFilter, p::CirclePolicy)
	#ax = -1.0 / sind(o)
	#ay = 1.0 / cosd(o)

	ax = -cosd(o);
	ay = sind(o);

	if ax*p.last[1] + ay*p.last[2] < 0.
		ax = -ax
		ay = -ay
	end

	temp_a = normalize((ax,ay,0.0), x)
	temp_x = x.x + temp_a[1]
	temp_y = x.y + temp_a[2]
	ax = temp_a[1]
	ay = temp_a[2]
	# If this moves you out of bounds, go the other way
	if temp_x < 0. || temp_x > m.length || temp_y < 0. || temp_y > m.length
		ax *= -1.
		ay *= -1.
	end

	p.last = (ax, ay, 0.0)
	return (ax, ay, 0.0)
	#return normalize((ax,ay,0.0), x)
end

######################################################################
# SpiralPolicy: Like CirclePolicy, but also move towards jammer
#
# TODO: actually implement this
######################################################################
# Move orthogonally to last measurement
type SpiralPolicy <: Policy
	last::Action

	SpiralPolicy() = new( (0.0, 0.0) )
end

# Remembers last action to ensure we follow same direction around circle
# Otherwise, it will "chatter" back and forth
# TODO: take into account distance from edge
function action(m::SearchDomain, x::Vehicle, o, f::DF, p::SpiralPolicy)
	ax = -1.0 / sind(o)
	ay = 1.0 / cosd(o)

	if ax*p.last[1] + ay*p.last[2] < 0.
		ax = -ax
		ay = -ay
	end

	p.last = (ax, ay)

	return normalize((ax,ay), x)
end


text:
# Set up MATLAB COM server if not in CI environment
if !parse(Bool, lowercase(get(ENV, "CI", "false")))
    Sys.iswindows() && run(`matlabCOM.bat`)
end


text:
export AbstractAffinePatch,
    nequations,
    state,
    onpatch!,
    setup!,
    changepatch!

"""
    AbstractAffinePatch

An affine patch is a hyperplane defined by ``v⋅x-1=0``.
"""
abstract type AbstractAffinePatch end

"""
    AbstractAffinePatchState{N}

This holds the actual patch information. `N` indicates the number of products of projective spaces.
"""
abstract type AbstractAffinePatchState{N} end

"""
    nequations(::AbstractAffinePatchState)

Number of equations an affine patch adds.
"""
nequations(::AbstractAffinePatchState{N}) where N = N

"""
    state(::AbstractAffinePatch, x)::AbstractAffinePatchState

Construct the state of the path from `x`.
"""
state(p::AbstractAffinePatch, x) = throw(MethodError(state, (p,x)))

"""
    onpatch!(x::AbstractVector, ::AbstractAffinePatchState)

Scale a vector `x` such that it is on the affine patch.
"""
function onpatch! end

"""
    setup!(::AbstractAffinePatchState, x::AbstractVector)

Setup the affine patch depending on `x` and modify `x` if necessary.
This is only called once at the beginning of a tracked path.
"""
setup!(state::AbstractAffinePatchState, x::AbstractVector) = onpatch!(x, state)

"""
    changepatch!(::AbstractAffinePatch, x::AbstractVector)

The same as [`setup!`](@ref) but only called during the path tracking.
"""
changepatch!(::AbstractAffinePatchState, x::AbstractVector) = nothing

include("affine_patches/common.jl")
include("affine_patches/orthogonal_patch.jl")
include("affine_patches/embedding_patch.jl")
include("affine_patches/random_patch.jl")
include("affine_patches/fixed_patch.jl")