text:
import Base.isequal
import Base.delete!, Base.union
import Base.join

export add!, delete!, contract!, induce, add_edges!
export line_graph, complement, complement!, adjoint
export cartesian, lex, relabel, trim
export disjoint_union, union, join, subdivide

# check for graph equality
function SimpleGraphs.delete!(G::SimpleGraph, H::SimpleGraph)
    if G.V != H.V || NE(G) != NE(H)
        return false
    end
    try  # check vertex sets for sortability
        sort(collect(G.V))
        # sort(collect(H.V))
    catch  #if not sortable, resort to a slow method
        for e in G.E
            if !has(H, e[1], e[2])
                return false
            end
        end
        return true
    end
    # otherwise faster to check edge sets this way
    return G.E == H.E
end

function SimpleGraphs.isequal(G::SimpleGraph, H::SimpleGraph)
    return G.V == H.V && G.E == H.E
end

function ==(G::SimpleGraph, H::SimpleGraph)
    return isequal(G, H)
end

# adding vertices
"""
`add!(G,v)` adds a new vertex `v` to the graph.

`add!(G,v,w)` adds a new edge `(v,w)` to the graph. If one (or both)
of those vertices is not already in the graph, it is added to the
vertex set.
"""
function add!(G::SimpleGraph{T}, v) where {T}
    if has(G, v)
        return false
    end
    cache_clear(G)
    push!(G.V, v)
    if G.Nflag
        G.N[v] = Set{T}()
    end
    return true
end

# adding edges
function add!(G::SimpleGraph{T}, v, w) where {T}
    if v == w
        return false
    end

    try
        if v > w
            v, w = w, v
        end
    catch
        if hash(v) > hash(w)
            v, w = w, v
        end
    end

    if ~has(G, v)
        add!(G, v)
    end
    if ~has(G, w)
        add!(G, w)
    end
    if has(G, v, w)
        return false
    end
    cache_clear(G)
    push!(G.E, (v, w))
    if G.Nflag
        push!(G.N[v], w)
        push!(G.N[w], v)
    end
    return true
end

"""
`add_edges!(G,edge_table)` adds edges to the graph `G`. The argument
`edge_table` is an iterable list (array or set) of edges to be
added (two-tuples). Returns the number of edges added to the graph.

This works both when `G` is a `SimpleGraph` and when `G` is
a `SimpleDigraph`.
"""
function add_edges!(G::AbstractSimpleGraph, edge_list)::Int
    m0 = NE(G)
    for e in edge_list
        add!(G, e[1], e[2])
    end
    return NE(G) - m0
end

# edge deletion
"""
`delete!(G,v)` deletes vertex `v` (and any edges incident with `v`)
from the graph.

`delete!(G,v,w)` deletes the edge `(v,w)` from `G`.
"""
function SimpleGraphs.delete!(G::SimpleGraph, v, w)
    flag = false
    if has(G, v, w)
        cache_clear(G)
        flag = true
        delete!(G.E, (v, w))
        delete!(G.E, (w, v))
        if G.Nflag
            delete!(G.N[v], w)
            delete!(G.N[w], v)
        end
    end
    return flag
end

# vertex deletion
function SimpleGraphs.delete!(G::SimpleGraph, v)
    flag = false
    if has(G, v)
        cache_clear(G)
        flag = true
        Nv = G[v]
        for w in Nv
            SimpleGraphs.delete!(G, v, w)
        end
        SimpleGraphs.delete!(G.V, v)
        if G.Nflag
            SimpleGraphs.delete!(G.N, v)
        end
    end
    return flag
end

# Contract an edge in a graph. If uv is an edge of G, we add all of
# v's neighbors to u's neighborhood and then delete v. This mutates
# the graph. Usually, vertices u and v are adjacent, but we don't
# require that. If either u or v is not a vertex of this graph, or if
# u==v, we return false. Otherwise we return true to indicate success.
"""
`contract!(G,u,v)` contracts the edge `(u,v)` in the graph. The merged
vertex is named `u` (hence `contract!(G,v,u)` results in a different,
albeit isomorphic, graph).

Note: The edge `(u,v)` need not be present in the graph. If missing,
this is equivalent to first adding the edge to the graph and then
contracting it.
"""
function contract!(G::SimpleGraph, u, v)
    if !has(G, u) || !has(G, v) || u == v
        return false
    end
    cache_clear(G)
    Nv = G[v]
    for x in Nv
        add!(G, u, x)
    end
    SimpleGraphs.delete!(G, v)
    return true
end

# Given a simple graph G and a set of vertices A, form the induced
# subgraph G[A]. Note that A must be a subset of V(G).
"""
`induce(G,A)` creates the induced subgraph of `G` with vertices in the
set `A`. The graph may be either a `SimpleGraph` or a `SimpleDigraph`.
"""
function induce(G::SimpleGraph{T}, A::Set) where {T}
    # Check that A is a subset of V(G)
    for v in A
        if ~has(G, v)
            error("The set A must be a subset of V(G)")
        end
    end
    H = SimpleGraph{T}()  # place to hold the answer

    # add all the vertices in A to H
    for v in A
        add!(H, v)
    end

    # The method we choose depends on the size of A. For small A, best
    # to iterate over pairs of elements of A. For large A, best to
    # iterate over G.E
    a = length(A)

    if a < 2
        return H
    end

    # case of small A, iterate over pairs from A
    if a * (a - 1) < 2 * NE(G)
        alist = collect(A)
        for i = 1:a-1
            u = alist[i]
            for j = i+1:a
                v = alist[j]
                if has(G, u, v)
                    add!(H, u, v)
                end
            end
        end

    else  # case of large A, iterate over G.E
        for e in G.E
            if in(e[1], A) && in(e[2], A)
                add!(H, e[1], e[2])
            end
        end
    end
    return H
end

# Digraph version
function induce(G::SimpleDigraph{T}, A::Set) where {T}
    # Check that A is a subset of V(G)
    for v in A
        if ~has(G, v)
            error("The set A must be a subset of V(G)")
        end
    end
    H = SimpleDigraph{T}()  # place to hold the answer

    # add all the vertices in A to H
    for v in A
        add!(H, v)
    end

    for x in A
        for y in A
            if has(G, x, y)
                add!(H, x, y)
            end
        end
    end
    return H
end

# Create the line graph of a given graph
"""
`line_graph(G)` creates the line graph of `G`.
"""
function line_graph(G::SimpleGraph{T}) where {T}
    H = SimpleGraph{Tuple{T,T}}()

    m = NE(G)
    E = elist(G)
    for e in E
        add!(H, e)
    end

    for i = 1:m-1
        e = E[i]
        for j = i+1:m
            f = E[j]
            if e[1] == f[1] || e[2] == f[1] || e[1] == f[2] || e[2] == f[2]
                add!(H, e, f)
            end
        end
    end
    return H
end

# Create the complement of a graph.
"""
`complement(G)` creates (as a new graph) the complement of `G`.
Note that `G'` is a short cut for `complement(G)`.
"""
function complement(G::SimpleGraph{T}) where {T}
    H = SimpleGraph{T}()
    V = vlist(G)
    n = NV(G)

    for i = 1:n
        add!(H, V[i])
    end

    for i = 1:n-1
        v = V[i]
        for j = i+1:n
            w = V[j]
            if ~has(G, v, w)
                add!(H, v, w)
            end
        end
    end
    if cache_check(G, :name)
        name(H, "Complement of $(name(G))")
    end
    return H
end

# We can use G' to mean complement(G)
"""
`G'` is equivalent to `complement(G)`.
"""
adjoint(G::SimpleGraph) = complement(G)

# complement in place. Returns None.
"""
`complement!(G)` replaces `G` with its complement.
"""
function complement!(G::SimpleGraph)
    name_flag = cache_check(G, :name)
    name_hold = ""
    if name_flag
        name_hold = name(G)
    end
    cache_clear(G)
    n = NV(G)
    V = vlist(G)
    for i = 1:n-1
        v = V[i]
        for j = i+1:n
            w = V[j]
            if has(G, v, w)
                SimpleGraphs.delete!(G, v, w)
            else
                add!(G, v, w)
            end
        end
    end
    if name_flag
        name(G, "Complement of $name_hold")
    end
    nothing
end

# Create the cartesian product of two graphs
"""
`cartesian(G,H)` creates the Cartesian product of the two graphs.
This can be abbreviated as `G*H`.
"""
function cartesian(G::SimpleGraph{S}, H::SimpleGraph{T}) where {S,T}
    K = SimpleGraph{Tuple{S,T}}()
    for v in G.V
        for w in H.V
            add!(K, (v, w))
        end
    end

    for e in G.E
        (u, v) = e
        for w in H.V
            add!(K, (u, w), (v, w))
        end
    end

    for e in H.E
        (u, v) = e
        for w in G.V
            add!(K, (w, u), (w, v))
        end
    end

    return K
end

# Use G*H for Cartesian product
"""
For `SimpleGraph`s: `G*H` is equivalent to `cartesian(G,H)`.
"""
function *(G::SimpleGraph{S}, H::SimpleGraph{T}) where {S,T}
    return cartesian(G, H)
end

# The join of two graphs is formed by taking disjoint copies of the
# graphs and all possible edges between the two.

"""
`join(G,H)` is a new graph formed by taking disjoint copies of
`G` and `H` together with all possible edges between those copies.
"""
function join(G::SimpleGraph{S}, H::SimpleGraph{T}) where {S,T}
    K = disjoint_union(G, H)
    for v in G.V
        for w in H.V
            add!(K, (v, 1), (w, 2))
        end
    end
    return K
end

# Create the union of two graphs. If they have vertices or edges in
# common, that's OK.

"""
`union(G,H)` creates the union of the graphs `G` and `H`. The graphs
may (and typically do) have common vertices or edges.
"""
function union(G::SimpleGraph{S}, H::SimpleGraph{T}) where {S,T}
    if S == T
        K = SimpleGraph{S}()
    else
        K = SimpleGraph{Any}()
    end

    for v in G.V
        add!(K, v)
    end
    for v in H.V
        add!(K, v)
    end

    for e in G.E
        add!(K, e[1], e[2])
    end
    for e in H.E
        add!(K, e[1], e[2])
    end
    return K
end

# This is an unexposed helper function that takes a graph and creates
# a new graph in which the name of each vertex has an integer
# appended. For example, if the vertex type is String in the original
# graph, the new vertices are type (String, Int).
function label_append(G::SimpleGraph{S}, a::Int) where {S}
    mapper = Dict{S,Tuple{S,Int}}()
    for v in G.V
        mapper[v] = (v, a)
    end
    return relabel(G, mapper)
end

# The disjoint_union of two graphs, G and H, is a new graph consisting of
# disjoint copies of G and H.

"""
`disjoint_union(G,H)` is a new graph formed by taking disjoint copies
of `G` and `H` (and no additional edges).
"""
function disjoint_union(G::SimpleGraph{S}, H::SimpleGraph{T}) where {S,T}
    GG = label_append(G, 1)
    HH = label_append(H, 2)
    if S == T
        K = SimpleGraph{Tuple{S,Int}}()
    else
        K = SimpleGraph{Any}()
    end

    for v in GG.V
        add!(K, v)
    end
    for v in HH.V
        add!(K, v)
    end

    for e in GG.E
        add!(K, e[1], e[2])
    end
    for e in HH.E
        add!(K, e[1], e[2])
    end
    return K
end

# Repeatedly remove vertices with the given degree or less until there
# are no such vertices remaining in the graph. The default trim(G)
# simply removes all isolated vertices.

"""
`trim(G)` returns a copy of `G` with all isolated vertices removed.

`trim(G,d)` returns a copy of `G` in which we iteratively remove all
vertices of degree `d` or smaller. For example, if `G` is a tree,
`trim(G,1)` will eventually remove all vertices.
"""
function trim(G::SimpleGraph, d::Int = 0)
    H = deepcopy(G)
    while NV(H) > 0 && minimum(deg(H)) <= d
        for v in H.V
            if deg(H, v) <= d
                SimpleGraphs.delete!(H, v)
            end
        end
    end
    return H
end

# Relabel the vertics of a graph based on a dictionary mapping old
# vertex names to new

"""
`relabel(G)` returns a copy of `G` in which the vertices are renamed
`1:n`.

`relabel(G,d)` (where `d` is a `Dict`) returns a copy of `G` in which
vertex `v` is renamed `d[v]`.
"""
function relabel(G::SimpleGraph{S}, label::Dict{S,T}) where {S,T}
    H = SimpleGraph{T}()
    for v in G.V
        add!(H, label[v])
    end

    for e in G.E
        u = label[e[1]]
        v = label[e[2]]
        add!(H, u, v)
    end
    return H
end

# Relabel the vertices with the integers 1:n
function relabel(G::SimpleGraph{S}) where {S}
    verts = vlist(G)
    n = length(verts)
    label = Dict{S,Int}()
    sizehint!(label, n)

    for idx = 1:n
        label[verts[idx]] = idx
    end

    return relabel(G, label)
end

"""
`lex(G,H)` computes the lexicographic product of the graphs `G` and
`H`. The vertex set of this graph is the set of all ordered pairs
`(g,h)` (with `g` a vertex of `G` and `h` a vertex of `H` in which we
have the edge `(g,h)~(g',h')` if either `g~g'` in `G` or else `g=g`
and `h~h'`.

We can use the notation `G[H]` also to create `lex(G,H)`.
"""
function lex(G::SimpleGraph{S}, H::SimpleGraph{T}) where {S,T}
    VT = Tuple{S,T}
    K = SimpleGraph{VT}()
    # Create vertex set
    for a in G.V
        for b in H.V
            add!(K, (a, b))
        end
    end

    # Add edges of the form (a,x)~(b,y) where a~b
    for e in G.E
        a = e[1]
        b = e[2]

        for x in H.V
            for y in H.V
                add!(K, (a, x), (b, y))
            end
        end
    end

    # Add edges of the form (a,x)~(a,y) where x~y
    for a in G.V
        for e in H.E
            x = e[1]
            y = e[2]
            add!(K, (a, x), (a, y))
        end
    end

    return K
end

"""
Abbreviation for `lex(G,H)` for `SimpleGraph`s.
"""
getindex(G::SimpleGraph, H::SimpleGraph) = lex(G, H)

"""
`subdivide(G::SimpleGraph)` creates a new `SimpleGraph` by replacing
every edge with a path of length two.
"""
function subdivide(G::SimpleGraph{T}) where {T}
    TT = Tuple{T,T}
    H = SimpleGraph{TT}()

    for v in G.V
        add!(H, (v, v))
    end

    for e in G.E
        u, v = e
        add!(H, (u, u), (u, v))
        add!(H, (v, v), (u, v))
    end

    if hasxy(G)
        Gxy = getxy(G)
        Hxy = Dict{TT,Vector{Float64}}()
        for v in G.V
            Hxy[(v, v)] = Gxy[v]
        end
        for e in G.E
            u, v = e
            Hxy[e] = (Gxy[u] + Gxy[v]) / 2
        end
        embed(H, Hxy)
    end

    return H
end


text:
draw(PNG("butt-mitre.png", 140, 100),
     Axes(border=0.2),
     Paper("lightgrey"),
     Pen(0.3, cap="butt", join="mitre")) do
    move(0, 0)
    line(0.5, 0)
    line(0.5, 0.5)
    line(1, 1)
end

ext: jl


text:
# This script lets us load the compiled .so file to bintray

using HTTPClient

function load_apikey()
	fin  = open("apikey.txt")
	lines = readlines(fin)
	close(fin)
	apikey = lines[1]
	@assert(length(apikey) == 40)
	apikey
end
const APIKEY = load_apikey()::ASCIIString
const USERNAME = "tawheeler"
const DESTPATH = "generic/libsmilejl/libsmilejl_unix_x64_0_1_0/generic/"

curpath = dirname(@__FILE__() )
downloads = curpath*"/downloads"
@assert(isdir(downloads))

srcfile = downloads * "/libsmilejl.so"
println("uploading ", srcfile)
println("username: ", USERNAME)
println("destpath: ", DESTPATH)

println("curl -T $srcfile -u$(USERNAME):$(APIKEY) https://api.bintray.com/content/$(USERNAME)/$DESTPATH")
run(`curl -T $srcfile -u$(USERNAME):$(APIKEY) https://api.bintray.com/content/$(USERNAME)/$DESTPATH`)

# NOTE(tim): if this complains that you published more than 180 days ago, just delete the "version" on bintray and recreate it

println("\tTry running build.jl to make sure the file can be downloaded")
println("done")


text:
# helper functions for interface2.jl.  These functions are not interface
# functions themselves, hence this file

"""
  This function starts parallel communication, if necessary, for
  functionals.  This function also unpacks `eqn.q_vec` into `eqn.q`
  if communication is started.

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * func: the functional
   * start_comm: bool specifying if communication must be started.
                 Communication may also be started when `start_comm` is false,
                 if the parallel data of `eqn.shared_data` is different.

  **Outputs**

   * start_comm_q: bool describing if communication was started or not
"""
function startCommunicationFunctional(mesh, sbp, eqn, opts,
                                      func::AbstractFunctional, start_comm::Bool)
  # start communication
  pdata = getParallelData(func)
  println(eqn.params.f, "pdata = ", pdata)
  start_comm_q = pdata != PARALLEL_DATA_NONE && (start_comm || pdata != getParallelData(eqn.shared_data))
  println(eqn.params.f, "start_comm_q = ", start_comm_q)
  if start_comm
    array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)
    eqn.params.time.t_send += @elapsed if eqn.commsize > 1 && start_comm_q
      setParallelData(eqn.shared_data, pdata)
      startSolutionExchange(mesh, sbp, eqn, opts)
    end
  end

  return start_comm_q
end


text:
push!(LOAD_PATH, "..")

using LinearAlgebra
using JuMP, Optim
using LogisticOptTools
using Xpress, GLPK

const LOT = LogisticOptTools

SOLVER = optimizer_with_attributes(Xpress.Optimizer, "PRESOLVE" => 0,
                                   "MIPPRESOLVE" => 0)
SOLVER = () -> Xpress.Optimizer(
        PRESOLVE = 0,
        CUTSTRATEGY = 0,
        HEURSTRATEGY = 0,
        SYMMETRY = 0,
        OUTPUTLOG = 3
    )
SOLVER = GLPK.Optimizer

sparsity(s::Vector; tol=0.5) = findall(x -> x > tol, s)

function oracle(X, y, γ, s)
    pattern = sparsity(s)
    Xs = X[:, pattern]
    glm = LOT.DualLogisticRegressor(Xs, y, penalty=LOT.L2Penalty(γ))
    f, gradient!, _, _ = LOT.generate_callbacks(glm)
    algo = LBFGS()
    options = Optim.Options(iterations=250, g_tol=1e-5, show_trace=false)
    lower = LOT.lowerbound(glm.data)
    upper = LOT.upperbound(glm.data)
    x0 = 0.5 * (lower .+ upper)
    res = Optim.optimize(f, gradient!, lower, upper,
                         x0, Fminbox(algo), options)
    c = res.minimum
    α♯ = res.minimizer
    ∇c = zeros(length(s))

    @inbounds for i in eachindex(s)
        atmp = dot(X[:, i], α♯)
        ∇c[i] = - γ * atmp^2
    end
    return -c, ∇c
end

# TODO: for some reason lazy constraints are not working with Xpress yet
function benders(X, y, s0, γ, k; solver=SOLVER)
    n, p = size(X)
    c0, ∇c0 = oracle(X, y, γ, s0)

    benders_model = Model(solver)

    # Variables s and t
    @variable(benders_model, 0 <= s[i=1:p] <= 1, Bin, start=s0[i])
    @variable(benders_model, 0 <= t , start=c0)
    # Global constraints
    @constraint(benders_model, t >= c0 + dot(∇c0, s - s0))
    @constraint(benders_model, sum(s) <= k)
    # Objective
    @objective(benders_model, Min, t)

    function lazy_constraint_callback(cb_data)
        s_current = [callback_value(cb_data, s[j]) for j in eachindex(s)]
        println(s_current)
        ck, ∇ck = oracle(X, y, γ, s_current)
        new_optimality_cons = @build_constraint(ck + dot(∇ck, s - s_current) <= t)
        MOI.submit(benders_model, MOI.LazyConstraint(cb_data), new_optimality_cons)
    end

    MOI.set(benders_model, MOI.LazyConstraintCallback(), lazy_constraint_callback)
    optimize!(benders_model)
    return benders_model
end

function benders_manual(X, y, s0, γ, k; solver=SOLVER, niter=100)
    n, p = size(X)
    c0, ∇c0 = oracle(X, y, γ, s0)

    benders_model = Model(solver)

    # Variables s and t
    @variable(benders_model, 0 <= s[i=1:p] <= 1, Bin, start=s0[i])
    @variable(benders_model, 0 <= t , start=c0)
    # Global constraints
    @constraint(benders_model, t >= c0 + dot(∇c0, s - s0))
    @constraint(benders_model, sum(s) <= k)
    # Objective
    @objective(benders_model, Min, t)

    for i in 1:niter
        optimize!(benders_model)
        sk = JuMP.value.(s)
        ck, ∇ck = oracle(X, y, γ, sk)
        @constraint(benders_model, t >= ck + dot(∇ck, s - sk))
        lb = JuMP.value(t)
        ub = ck
        println("LB: ", lb,  " UB ", ub)

        if (ub - lb) / max(1.0, abs(lb)) <= 1e-6
            break
        end
    end

    return benders_model
end

γ = 1e-2
# Set number of non-null component
k = 2
DATASET = "diabetes"
#= DATASET = "covtype.libsvm.binary.bz2" =#
#= DATASET = "colon-cancer.bz2" =#
dataset = LOT.LogitData(DATASET, scale_data=true)
X, y = dataset.X, dataset.y

n, p = size(X)
# Scale penalty for dual problem
dγ = 1.0 / (4.0 * γ * n)

s0 = zeros(p)
s0[1:k] .= 1

m = benders_manual(X, y, s0, dγ, k)


text:
mutable struct DemoTrainingAlephStar
    # network weights
    w::Vector{Array{Float32,N} where N}

    # tree params
    stepc::Int64
    epsilon::Float64
    MIN_EPSILON::Float64
    fac::Float64
    gamma::Float64

    # training params
    LR::Float64
    batchsize::Int64
    train_epochs::Int64
    max_states::Int64
    train_when_min::Int64

    # for validation
    val_safety::Float64

    # experience buffer
    sensors::Vector{Array{UInt8,2}}
    qs::Vector{Vector{Float32}}
    priorities::Vector{Float64}

    # for stats
    rewards::Vector{Float64}
    ranks::Vector{Int64}
    num_of_done_leafs::Vector{Int64}
    num_of_leafs::Vector{Int64}
    num_of_car_crash::Vector{Int64}
    num_of_wall_crash::Vector{Int64}
    num_of_vel_crash::Vector{Int64}
    tree_mean_speed::Vector{Float64}
    tree_std_speed::Vector{Float64}
    val_mean_speed::Vector{Float64}
    val_std_speed::Vector{Float64}
    val_steps::Vector{Float64}
    val_rewards::Vector{Float64}
    avg_window::Int64
    weighted_nodes_threshold::Int64
    function DemoTrainingAlephStar()
        new(
            Vector{Array{Float32,N} where N}(),
            0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            Vector{Array{UInt8,2}}(),
            Vector{Vector{Float32}}(),
            Vector{Float64}(),
            Vector{Float64}(), Vector{Int64}(),
            Vector{Int64}(), Vector{Int64}(),
            Vector{Int64}(), Vector{Int64}(),
            Vector{Int64}(), Vector{Float64}(),
            Vector{Float64}(), Vector{Float64}(),
            Vector{Float64}(), Vector{Float64}(),
            Vector{Float64}(), 0, 0,
        )
    end
end

function InitializeDTAS()
    dtas = DemoTrainingAlephStar()
    dtas.w = initialize_weights(actionc_steer*actionc_accel)
    dtas.stepc = 5500
    dtas.epsilon = 0.5
    dtas.MIN_EPSILON = 0.01
    dtas.fac = 0.995
    dtas.gamma = 0.98
    dtas.LR = 0.01
    dtas.batchsize = 64
    dtas.train_epochs = 20
    dtas.max_states = 80000
    dtas.train_when_min = 50000
    dtas.val_safety = 1.0
    dtas.avg_window = 50
    dtas.weighted_nodes_threshold = 200
    dtas
end

function traindtas(dtas::DemoTrainingAlephStar, iters::Integer, fname)
    w = map(KnetArray, dtas.w)
    for i in 1:iters
        #### Accumulating tree ###################################################

        println("--------------------------------------------------------------- i="*string(i))
        flush(stdout)

        state, env = initialize_simple_road();
        tree = build_tree(w, env, state, dtas.stepc, dtas.epsilon, dtas.gamma)
        backprop_weighted_q!(tree, dtas.gamma, dtas.weighted_nodes_threshold)

        # gather some stats
        push!(dtas.num_of_car_crash, 0)
        push!(dtas.num_of_wall_crash, 0)
        push!(dtas.num_of_vel_crash, 0)
        # some crash stats
        for node in tree.nodes
            !is_leaf(node) && continue
            sim_state = tree.states[node.id]
            yolo_x = env.splx(sim_state.yolocar_dist)
            yolo_y = env.sply(sim_state.yolocar_dist)
            _, result = my_calc_reward(getcar(sim_state), env.coll, yolo_x, yolo_y, SAFETY)
            if result==1
                dtas.num_of_vel_crash[end] += 1
            elseif result==2
                dtas.num_of_wall_crash[end] += 1
            elseif result==3
                dtas.num_of_car_crash[end] += 1
            end
            @assert is_leaf(node)
            # this is possible due to floating point issues
            # result==0 && @assert !node.done
        end
        # calculate mean speed
        treev = Float64[]
        n = get_best_leaf(tree, dtas.gamma)
        times_passed = tree.states[n.id].times_passed_yolocar
        while !is_root(n)
            sim_state = tree.states[n.id]
            push!(treev, sim_state.speed)
            n = get(n.parent)
        end
        max_rank = get_tree_rank(tree, dtas.gamma)
        push!(dtas.ranks, max_rank)
        push!(dtas.tree_mean_speed, mean(treev))
        push!(dtas.tree_std_speed, std(treev))
        push!(dtas.rewards, get_accumulated_reward(tree))
        push!(dtas.num_of_done_leafs, sum([tree.dones[n.id] for n in tree.nodes if is_leaf(n)]))
        push!(dtas.num_of_leafs, sum([is_leaf(n) for n in tree.nodes]))
        # these could be slightly different due to floating point issues:
        # @assert num_of_car_crash[end] + num_of_vel_crash[end] + num_of_wall_crash[end] == num_of_done_leafs[end]

        ixs = findall([!is_leaf(n) || tree.dones[n.id] for n in tree.nodes])
        append!(dtas.qs, tree.children_qs[ixs])
        append!(dtas.sensors, tree.sensors[ixs])
        f = if length(dtas.priorities)==0
            1.0
        else
            median(dtas.priorities)
        end
        append!(dtas.priorities, f.*ones(Float32, length(ixs)))
        if length(dtas.qs) > dtas.max_states
            dtas.qs = dtas.qs[(end-dtas.max_states+1):end]
            dtas.sensors = dtas.sensors[(end-dtas.max_states+1):end]
            dtas.priorities = dtas.priorities[(end-dtas.max_states+1):end]
        end

        #### Training #######################################################

        if length(dtas.qs) >= dtas.train_when_min
            N = round(Int64, length(ixs) * dtas.train_epochs / dtas.batchsize)
            println("\n---- training for "*string(N)*" iterations"); flush(stdout)
            w,ll = trainit(N, w, dtas.sensors, dtas.qs, dtas.LR, dtas.batchsize, dtas.priorities)
            dtas.w = map(Array, w)
            if dtas.epsilon > dtas.MIN_EPSILON
                dtas.epsilon *= dtas.fac
            end
        end

        #### Validating #####################################################

        valv = Float64[]
        valr = 0.0
        spd = mean(dtas.tree_mean_speed[max(1,end-dtas.avg_window):end])
        if isnan(spd)
            spd = 0.8
        end
        state, val_env = initialize_simple_road(spd);
        for _ in 1:dtas.stepc
            push!(valv, state.speed)

            _sensors = get_sensors(val_env, state)
            _vqs = Vector{Float32}(network_predict(val_env, w, _sensors))
            action = action_ix_to_action(val_env,Int32(argmax(_vqs)))
            state, reward, done = sim!(val_env, state, action, 5, dtas.val_safety)
            valr += reward
            done && break
        end
        vpassed = state.times_passed_yolocar
        push!(dtas.val_rewards, valr)
        push!(dtas.val_mean_speed, mean(valv))
        push!(dtas.val_std_speed, std(valv))
        push!(dtas.val_steps, length(valv))

        #### Reporting #####################################################

        println()
        @show length(dtas.qs), length(ixs)
        println("      eps   = "*string(round(dtas.epsilon, digits=4))*"    gamma="*string(round(dtas.gamma, digits=4)))
        println(" tree_rnk   = "*string(max_rank)*    "    mvel    = "*string(round(mean(treev), digits=4))*"    svel = "*string(round(std(treev), digits=4)))
        println("val_steps   = "*string(length(valv))*"    mvel    = "*string(round(mean(valv),digits=4))*"    svel = "*string(round(std(valv),digits=4)))
        println("tree_passed = "*string(times_passed)*"    vpassed = "*string(vpassed))
        flush(stdout)

        if i%10 == 0
            IJulia.clear_output();
        end
        if i%100 == 0
            @save fname dtas
        end

    end
end

function plotdtas(dtas)
    p1 = plot(moving_avg(dtas.rewards, dtas.avg_window), label="max tree reward")
    plot!(p1, moving_avg(dtas.val_rewards, dtas.avg_window), lw=3, label="network rewards")

    p2 = plot(moving_avg(dtas.ranks/dtas.stepc, dtas.avg_window), label="max tree rank (perc.)")
    plot!(p2, moving_avg(dtas.val_steps/dtas.stepc, dtas.avg_window), lw=3, label="network rank")
    plot!(p2, moving_avg(dtas.num_of_done_leafs/dtas.stepc, dtas.avg_window), label="done perc.")
    plot!(p2, moving_avg(dtas.num_of_leafs/dtas.stepc, dtas.avg_window), label="leafs perc.")

    p3 = plot(moving_avg(dtas.num_of_vel_crash./dtas.num_of_done_leafs, dtas.avg_window), label="vel crash perc.", lw=3)
    plot!(p3, moving_avg(dtas.num_of_wall_crash./dtas.num_of_done_leafs, dtas.avg_window), label="wall crash perc.", lw=3)
    plot!(p3, moving_avg(dtas.num_of_car_crash./dtas.num_of_done_leafs, dtas.avg_window), label="car crash perc.", lw=3)

    p4 = plot(moving_avg(dtas.tree_mean_speed, dtas.avg_window), label="tree mean speed")
    plot!(p4, moving_avg(dtas.tree_std_speed, dtas.avg_window), label="tree std speed")
    plot!(p4, moving_avg(dtas.val_std_speed, dtas.avg_window), label="network std speed")
    plot!(p4, moving_avg(dtas.val_mean_speed, dtas.avg_window), label="network mean speed")

    plot(p1,p2,p3,p4, layout=(4,1), size=(900,1300))
end


text:
# using Pkg; for p in ("Knet","ArgParse"); haskey(Pkg.installed(),p) || Pkg.add(p); end

"""

This example learns to classify hand-written digits from the
[MNIST](http://yann.lecun.com/exdb/mnist) dataset.  There are 60000
training and 10000 test examples. Each input x consists of 784 pixels
representing a 28x28 image.  The pixel values are normalized to
[0,1]. Each output y is converted to a ten-dimensional one-hot vector
(a vector that has a single non-zero component) indicating the correct
class (0-9) for a given image.  10 is used to represent 0.

You can run the demo using `julia lenet.jl` at the command line or
`julia> LeNet.main()` at the Julia prompt.  Use `julia lenet.jl
--help` or `julia> LeNet.main("--help")` for a list of options.  The
dataset will be automatically downloaded.  By default the
[LeNet](http://yann.lecun.com/exdb/lenet) convolutional neural network
model will be trained for 10 epochs.  The accuracy for the training
and test sets will be printed at every epoch and optimized parameters
will be returned.

"""
module LeNet
using Knet,CUDA,ArgParse
include(Knet.dir("data","mnist.jl"))

function predict(w,x)
    n=length(w)-4
    for i=1:2:n
        x = pool(relu.(conv4(w[i],x;padding=0) .+ w[i+1]))
    end
    x = mat(x)
    for i=n+1:2:length(w)-2
        x = relu.(w[i]*x .+ w[i+1])
    end
    return w[end-1]*x .+ w[end]
end

loss(w,x,ygold) = nll(predict(w,x), ygold)

lossgradient = grad(loss)

function train(w, data; lr=.1, epochs=3, iters=1800)
    for epoch=1:epochs
        for (x,y) in data
            g = lossgradient(w, x, y)
            update!(w, g, lr=lr)
            if (iters -= 1) <= 0
                return w
            end
        end
    end
    return w
end

function weights(;atype=KnetArray{Float32})
    w = Array{Any}(undef,8)
    w[1] = xavier(5,5,1,20)
    w[2] = zeros(1,1,20,1)
    w[3] = xavier(5,5,20,50)
    w[4] = zeros(1,1,50,1)
    w[5] = xavier(500,800)
    w[6] = zeros(500,1)
    w[7] = xavier(10,500)
    w[8] = zeros(10,1)
    return map(a->convert(atype,a), w)
end

function main(args=ARGS)
    s = ArgParseSettings()
    s.description="lenet.jl (c) Deniz Yuret, 2016. The LeNet model on the MNIST handwritten digit recognition problem from http://yann.lecun.com/exdb/mnist."
    s.exc_handler=ArgParse.debug_handler
    @add_arg_table s begin
        ("--seed"; arg_type=Int; default=-1; help="random number seed: use a nonnegative int for repeatable results")
        ("--batchsize"; arg_type=Int; default=128; help="minibatch size")
        ("--lr"; arg_type=Float64; default=0.1; help="learning rate")
        ("--fast"; action=:store_true; help="skip loss printing for faster run")
        ("--epochs"; arg_type=Int; default=3; help="number of epochs for training")
        ("--iters"; arg_type=Int; default=typemax(Int); help="maximum number of updates for training")
        ("--gcheck"; arg_type=Int; default=0; help="check N random gradients per parameter")
        ("--atype"; default=(gpu()>=0 ? "KnetArray{Float32}" : "Array{Float32}"); help="array and float type to use")
    end
    isa(args, AbstractString) && (args=split(args))
    if in("--help", args) || in("-h", args)
        ArgParse.show_help(s; exit_when_done=false)
        return
    end
    println(s.description)
    o = parse_args(args, s; as_symbols=true)
    println("opts=",[(k,v) for (k,v) in o]...)
    o[:seed] > 0 && Knet.seed!(o[:seed])
    atype = eval(Meta.parse(o[:atype]))
    if atype <: Array; warn("CPU conv4 support is experimental and very slow."); end

    xtrn,ytrn,xtst,ytst = mnist()
    global dtrn = minibatch(xtrn, ytrn, o[:batchsize]; xtype=atype)
    global dtst = minibatch(xtst, ytst, o[:batchsize]; xtype=atype)
    w = weights(atype=atype)
    report(epoch)=println((:epoch,epoch,:trn,accuracy(w,dtrn,predict),:tst,accuracy(w,dtst,predict)))

    if o[:fast]
        @time (train(w, dtrn; lr=o[:lr], epochs=o[:epochs], iters=o[:iters]); gpu()>=0 && Knet.cudaDeviceSynchronize())
    else
        report(0)
        iters = o[:iters]
        for epoch=1:o[:epochs]
            @time train(w, dtrn; lr=o[:lr], epochs=1, iters=iters)
            report(epoch)
            if o[:gcheck] > 0
                gradcheck(loss, w, first(dtrn)...; gcheck=o[:gcheck], verbose=true)
            end
            if (iters -= length(dtrn)) <= 0; break; end
        end
    end
    return w
end

# This allows both non-interactive (shell command) and interactive calls like:
# $ julia lenet.jl --epochs 10
# julia> LeNet.main("--epochs 10")
PROGRAM_FILE == "lenet.jl" && main(ARGS)

end # module

# SAMPLE RUN 65f57ff+ Wed Sep 14 10:02:30 EEST 2016
#
# lenet.jl (c) Deniz Yuret, 2016. The LeNet model on the MNIST handwritten digit recognition problem from http://yann.lecun.com/exdb/mnist.
# opts=(:seed,-1)(:batchsize,100)(:epochs,3)(:lr,0.1)(:gcheck,0)(:fast,true)
# ..................
# 9.319163 seconds (5.84 M allocations: 277.927 MB, 7.37% gc time)


text:
# Define exit time
A = 600

# Starting amount of money (S0)
k = 500

# Define probability of winning
p = .49
q = 1 - p

# Define b vector
b = zeros(A + 1)
b[A+1] = 1

# Define M matrix
M = zeros(A+1, A+1)
M[1,1] = 1
M[A+1,A+1] = 1

for i=2:A
    global M, q, p
    M[i,i] = -1
    M[i,i-1] = q
    M[i,i+1] = p
end

# Compute exit probabilities by computing M^(-1)b
u = M\b

# Get probability of interest
using Printf
@printf("Probability of walking away with \$%g: %f", k,
 u[k+1])

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
using Hangul # Jamo Letter Consonant Vowel Syllable DecomposeError to_char to_syllable
using Test

g = Jamo('ㄱ')
@test g isa Letter{Consonant}
@test g == 'ㄱ'

o = Jamo('ㅗ')
@test o isa Letter{Vowel}
@test o == 'ㅗ'

m = Jamo('ㅁ')
@test m isa Letter{Consonant}
@test m == 'ㅁ'

@test '고' == to_char(Syllable(g, o))
@test '고' == to_char(Syllable(g, o, nothing))
@test '곰' == to_char(Syllable(g, o, m))

@test 'ㄱ' == to_char(Syllable(g, nothing, nothing))
@test 'ㅗ' == to_char(Syllable(nothing, o, nothing))
@test 'ㅁ' == to_char(Syllable(nothing, nothing, m))

@test_throws ComposeError to_char(Syllable(g, nothing, m))
@test_throws ComposeError to_char(Syllable(nothing, o, m))
@test_throws ComposeError to_char(Syllable(nothing, nothing, nothing))

@test Syllable(g, o, nothing) == to_syllable('고')
@test Syllable(g, o, m) == to_syllable('곰')
@test Syllable(g, nothing, nothing) == to_syllable('ㄱ')
@test Syllable(nothing, o, nothing) == to_syllable('ㅗ')

@test_throws DecomposeError to_syllable('ㄳ') # 종성은 DecomposeError
@test_throws DecomposeError to_syllable('a')

@test Letter{Consonant}('ㄱ') == Jamo('ㄱ') == Jamo(0x3131)
@test Letter{Vowel}('ㅗ') == Jamo('ㅗ') == Jamo(0x3157)

@test_throws Hangul.LetterError Jamo('a')


text:
"""
    shortestpath(distancemap, startpoint, sourcepoint, stepsize)

    This function traces the shortest path from `startpoint` to
    `sourcepoint` using Runge Kutta 4 in a 2D `distancemap`.

    Optionally specify the line tracing `stepsize`.

    Returns a 2D-array with the shortest path.
"""
function shortestpath(
    distancemap::AbstractArray{T},
    startpoint::AbstractArray{T},
    sourcepoint=ones(T, 2, 1),
    stepsize::T=T(0.5),
) where {T}
    gradientvolume = zeros(T, size(distancemap, 1), size(distancemap, 2), 2)
    return shortestpath!(gradientvolume, distancemap, startpoint, sourcepoint, stepsize)
end

function shortestpath!(
    gradientvolume,
    distancemap::AbstractArray{T},
    startpoint::AbstractArray{T},
    sourcepoint=ones(T, 2, 1),
    stepsize::T=0.5,
) where {T}
    distancetoend = T(Inf)
    Fx = view(gradientvolume, :, :, 1)
    Fy = view(gradientvolume, :, :, 2)
    pointmin!(distancemap, Fx, Fy)
    Fx .= -Fx
    Fy .= -Fy

    i = 0
    # Reserve a block of memory for the shortest line array
    ifree = size(distancemap, 1) * 2
    shortestline = zeros(T, ifree, ndims(distancemap))

    # Iteratively trace the shortest line
    movement = T(1)
    while true
        endpoint = rk4(startpoint, gradientvolume, stepsize)

        # Calculate the distance to the end point
        if !isempty(sourcepoint) && !isnan(sourcepoint[1])
            (distancetoend, ind) = findmin(sqrt.(sum(
                (sourcepoint .- repeat(endpoint, 1, size(sourcepoint, 2))) .^ 2
            )))
        else
            distancetoend = Inf
        end

        # Calculate the movement between current point and point 10 iterations back
        if i > 10
            movement = sqrt(sum((endpoint - shortestline[i - 10, :]) .^ 2))
        end

        # Stop if out of boundary, distance to end smaller then a pixel or
        # if we have not moved for 10 iterations
        if (endpoint[1] == 0) || (movement < stepsize)
            break
        end

        # Count the number of iterations
        i = i + 1

        # Add a new block of memory if nearly full
        if i > (ifree - 2)
            ifree = ifree * 2
            shortestline = vcat(shortestline, zeros(T, ifree, ndims(distancemap)))
        end

        # Add current point to the shortest line array
        shortestline[i, :] = endpoint

        if distancetoend < stepsize
            i = i + 1
            # Add (Last) Source point to the shortest line array
            shortestline[i, :] = sourcepoint
            break
        end

        # Current point is next Starting Point
        startpoint = endpoint
    end

    if (distancetoend > 1) && (~isempty(sourcepoint))
        @warn "The shortest path trace did not finish at the source point" distancetoend
    end

    # Remove unused memory from array
    shortestline = shortestline[1:i, :]
    return shortestline
end


text:
import Base: read, write

"""
    read(fil::String,x::MeshArray)

Read file / array into MeshArray. Methods:

```
read(fil::String,x::MeshArray) #from File
read(xx::Array,x::MeshArray) #from Array
read(xx::Array,γ::gcmgrid) #from Array
```
"""
function read(fil::String,x::MeshArray)

  (n1,n2)=x.grid.ioSize
  (nFaces,n3)=nFacesEtc(x)

  fid = open(fil)
  xx = Array{eltype(x),2}(undef,(n1*n2,n3))
  read!(fid,xx)
  xx = reshape(hton.(xx),(n1,n2,n3))
  close(fid)

  return x.grid.read(xx,x)

end

function read(xx::Array,γ::gcmgrid)
  n3=Int(round(prod(size(xx))/prod(γ.ioSize)))
  read(xx,MeshArray(γ,γ.ioPrec,n3))
end

function read(xx::Array,x::MeshArray)

  facesSize=x.grid.fSize
  (n1,n2)=x.grid.ioSize
  (nFaces,n3)=nFacesEtc(x)

  size(xx)!=(n1*n2,n3) ? xx=reshape(xx,(n1*n2,n3)) : nothing

  y=similar(x; m=x.meta)
  i0=0; i1=0;
  for iFace=1:nFaces
    i0=i1+1;
    nn=facesSize[iFace][1]; mm=facesSize[iFace][2];
    i1=i1+nn*mm;

    if n3>1 && ndims(x.f)==1
      y.f[iFace]=reshape(xx[i0:i1,:],(nn,mm,n3))
    elseif n3>1 && ndims(x.f)==2
      for i3=1:n3
        y.f[iFace,i3]=reshape(xx[i0:i1,i3],(nn,mm))
      end
    else
      y.f[iFace]=reshape(xx[i0:i1,:],(nn,mm))
    end
  end

  return y

end

"""
    write(fil::String,x::MeshArray)

Write MeshArray to binary file. Other methods:

```
write(xx::Array,x::MeshArray) #to Array
```
"""
function write(fil::String,x::MeshArray)
  y=x.grid.write(x)
  fid = open(fil,"w")
  write(fid,ntoh.(y))
  close(fid)
end

function write(x::MeshArray)

  facesSize=x.grid.fSize
  (n1,n2)=x.grid.ioSize
  (nFaces,n3)=nFacesEtc(x)

  y = Array{eltype(x),2}(undef,(n1*n2,n3))
  i0=0; i1=0;
  for iFace=1:nFaces;
    i0=i1+1;
    nn=facesSize[iFace][1];
    mm=facesSize[iFace][2];
    i1=i1+nn*mm;
    if n3>1 && ndims(x.f)==2
      for i3=1:n3
        y[i0:i1,i3]=reshape(x.f[iFace,i3],(nn*mm,1))
      end
    else
      y[i0:i1,:]=reshape(x.f[iFace],(nn*mm,n3))
    end
  end

  y=reshape(y,(n1,n2,n3));
  n3==1 ? y=dropdims(y,dims=3) : nothing

  return y

end

##

function read_tiles(fil::String,x::MeshArray)

  (n1,n2)=x.grid.ioSize
  (nFaces,n3)=nFacesEtc(x)

  fid = open(fil)
  xx = Array{eltype(x),2}(undef,(n1*n2,n3))
  read!(fid,xx)
  xx = reshape(hton.(xx),(n1,n2,n3))
  close(fid)

  return x.grid.read(xx,x)

end

function read_tiles(xx::Array,γ::gcmgrid)
  n3=Int(round(prod(size(xx))/prod(γ.ioSize)))
  read_tiles(xx,MeshArray(γ,γ.ioPrec,n3))
end

function read_tiles(xx::Array,x::MeshArray)
	tmp=similar(x)
	s=x.grid.ioSize
	(n1,n2)=x.grid.fSize[1]
	ni=Int(s[1]/n1)
	nj=Int(s[2]/n2)
	ii=[0]
	for j in 1:nj, i in 1:ni
		ii[1]+=1
		tmp[ii[1]]=xx[(i-1)*n1 .+ collect(1:n1),(j-1)*n2 .+ collect(1:n2)]
	end
	tmp
end

function write_tiles(x::MeshArray)
  tmp = Array{eltype(x),2}(undef,x.grid.ioSize)
	s=x.grid.ioSize
	(n1,n2)=x.grid.fSize[1]
	ni=Int(s[1]/n1)
	nj=Int(s[2]/n2)
	ii=[0]
	for j in 1:nj, i in 1:ni
		ii[1]+=1
		tmp[(i-1)*n1 .+ collect(1:n1),(j-1)*n2 .+ collect(1:n2)].=x[ii[1]]
	end
	tmp
end

function write_tiles(fil::String,x::MeshArray)
  y=x.grid.write(x)
  fid = open(fil,"w")
  write(fid,ntoh.(y))
  close(fid)
end

##


text:
@testset "Scalar domains" begin
    @test in_domain(sin, 10.)
    @test in_domain(cos, 10.)
    @test !in_domain(acos, 10.)
    @test !in_domain(asin, 10.)
    @test domain1(sin) == (minimum(points), maximum(points))
    @test domain1(log) == (minimum(points[points .> 0]), maximum(points))
    @test domain1(acos) == (minimum(points[points .> -1]),
                             maximum(points[points .< 1]))
    @test domain2((+)) == ((minimum(points), maximum(points)),
                            (minimum(points), maximum(points)))
    @test domain2((^)) == ((minimum(points[points .> 0]), maximum(points)),
                            (minimum(points), maximum(points)))
    @test domain2(beta) == ((minimum(points[points .> 0]), maximum(points)),
                                 (minimum(points[points .> 0]), maximum(points)))
end

@testset "Scalar" begin
    let v = 1.0, ȳ = 5.0, z̄ = 4.0, rng = MersenneTwister(123456)
        let
            @test ∇(identity)(5) === (1,)
            @test ∇(identity)(5.0) === (1.0,)
        end

        unary_check(f, x) = check_errs(f, ȳ, x, v)
        @testset "$f" for f in UNARY_SCALAR_SENSITIVITIES
            domain = domain1(f)
            domain === nothing && error("Could not determine domain for $f.")
            lb, ub = domain
            randx = () -> rand(rng) * (ub - lb) + lb

            for _ in 1:10
                @test unary_check(f, randx())
            end
        end

        @testset "$f" for f in BINARY_SCALAR_SENSITIVITIES
            if f in ONLY_DIFF_IN_SECOND_ARG_SENSITIVITIES
                # First argument is not differentiable, it is integer-valued.
                domain = domain1(y -> f(0, y))
                domain === nothing && error("Could not determine domain for $f.")
                lb, ub = domain
                randx = () -> rand(rng, 0:5)
                randy = () -> rand(rng) * (ub - lb) + lb

                for _ in 1:10
                    x = randx()
                    @test check_errs(y -> f(x, y), ȳ, randy(), v)
                end
            else  # Both arguments are differentiable
                domain = domain2(f)
                domain === nothing && error("Could not determine domain for $f.")
                (x_lb, x_ub), (y_lb, y_ub) = domain
                randx = () -> rand(rng) * (x_ub - x_lb) + x_lb
                randy = () -> rand(rng) * (y_ub - y_lb) + y_lb

                for _ in 1:10
                    @test check_errs(f, z̄, (randx(), randy()), (v, v))
                end
            end
        end

        # Test whether the exponentiation amibiguity is resolved.
        @test ∇(x -> x^2)(1) == (2.0,)
    end

    @testset "float" begin
        # Scalars
        x_ = 4
        x = Leaf(Tape(), x_)
        y = float(x)
        @test y isa Branch{Float64}
        @test unbox(y) == 4.0

        # Arrays
        X_ = [1,2,3,4]
        X = Leaf(Tape(), X_)
        Y = float(X)
        @test Y isa Branch{Vector{Float64}}
        @test unbox(Y) == Float64[1,2,3,4]

        # In expressions.
        @test ∇(x->2*float(x))(1) === (2.0,)

        # Note that actually float or no float out is the same as
        # integers always have floating point derivatives
        @test ∇(x->2x)(1) === (2.0,)

    end
end


text:
@testset "MPC: Policy for Quadruped" begin
    T = Float64

    s = get_simulation("quadruped", "flat_2D_lc", "flat")
    model = s.model
    env = s.env

    ref_traj_ = deepcopy(ContactImplicitMPC.get_trajectory(s.model, s.env,
        joinpath(module_dir(), "src/dynamics/quadruped/gaits/gait2.jld2"),
        load_type = :split_traj_alt))

    ref_traj = deepcopy(ref_traj_)

    # time
    H = ref_traj.H
    h = ref_traj.h
    N_sample = 5
    H_mpc = 10
    h_sim = h / N_sample
    H_sim = 500

    # barrier parameter
    κ_mpc = 2.0e-4

    # obj
    obj = TrackingVelocityObjective(model, env, H_mpc,
        q = [Diagonal(1e-2 * [10; 0.02; 0.25; 0.25 * ones(model.dim.q-3)]) for t = 1:H_mpc],
        v = [Diagonal(0e-2 * [10; 0.02; 0.25; 0.25 * ones(model.dim.q-3)]) for t = 1:H_mpc],
        u = [Diagonal(3e-2 * ones(model.dim.u)) for t = 1:H_mpc],
        γ = [Diagonal(1.0e-100 * ones(model.dim.c)) for t = 1:H_mpc],
        b = [Diagonal(1.0e-100 * ones(model.dim.c * friction_dim(env))) for t = 1:H_mpc])

    # linearized MPC policy
    p = ContactImplicitMPC.linearized_mpc_policy(ref_traj, s, obj,
        H_mpc = H_mpc,
        N_sample = N_sample,
        κ_mpc = κ_mpc,
        n_opts = ContactImplicitMPC.NewtonOptions(
            r_tol = 3e-4,
            max_iter = 5),
        mpc_opts = ContactImplicitMPC.LinearizedMPCOptions())

    # initial configurations
    q1_ref = copy(ref_traj.q[2])
    q0_ref = copy(copy(ref_traj.q[1]))
    q1_sim = SVector{model.dim.q}(q1_ref)
    q0_sim = SVector{model.dim.q}(copy(q1_sim - (q1_ref - q0_ref) / N_sample))
    @assert norm((q1_sim - q0_sim) / h_sim - (q1_ref - q0_ref) / h) < 1.0e-8

    # simulator
    sim = ContactImplicitMPC.simulator(s, q0_sim, q1_sim, h_sim, H_sim,
        p = p,
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
			γ_reg = 0.0,
			undercut = Inf,
            r_tol = 1.0e-8,
            κ_tol = 1.0e-8),
        sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = true))

    # simulator
    @test status = ContactImplicitMPC.simulate!(sim)
    ref_traj = deepcopy(ref_traj_)

    qerr, uerr, γerr, berr = ContactImplicitMPC.tracking_error(ref_traj, sim.traj, N_sample, idx_shift=[1])
    @test qerr < 0.0201 * 1.5 # 0.0201
    @test uerr < 0.0437 * 1.5 # 0.0437
    @test γerr < 0.374 * 1.5 # 0.374
    @test berr < 0.0789 * 1.5 # 0.0789
    qerr > 0.0201 * 1.2 && @warn "mild regression on q tracking: current tracking error = $qerr, nominal tracking error = 0.0201"
    uerr > 0.0437 * 1.2 && @warn "mild regression on u tracking: current tracking error = $uerr, nominal tracking error = 0.0437"
    γerr > 0.3740 * 1.2 && @warn "mild regression on γ tracking: current tracking error = $γerr, nominal tracking error = 0.374"
    berr > 0.0789 * 1.2 && @warn "mild regression on b tracking: current tracking error = $berr, nominal tracking error = 0.0789"
end

# @testset "MPC: Policy for Quadruped on Sinusoidal Terrain" begin
#     T = Float64

#     s_sim = get_simulation("quadruped", "sine1_2D_lc", "sinusoidal")
#     s = get_simulation("quadruped", "flat_2D_lc", "flat")
#     model = s.model
#     env = s.env

#     nq = model.dim.q

#     ref_traj_ = deepcopy(ContactImplicitMPC.get_trajectory(s.model, s.env,
#         joinpath(module_dir(), "src/dynamics/quadruped/gaits/gait2.jld2"),
#         load_type = :split_traj_alt))
#     ref_traj = deepcopy(ref_traj_)

#     # time
#     H = ref_traj.H
#     h = ref_traj.h
#     N_sample = 5
#     H_mpc = 10
#     h_sim = h / N_sample
#     H_sim = 1500

#     # barrier parameter
#     κ_mpc = 2.0e-4

#     obj = TrackingObjective(model, env, H_mpc,
#         q = [Diagonal(1e-2 * [10; 0.02; 0.25; 0.25 * ones(nq-3)]) for t = 1:H_mpc],
#         u = [Diagonal(3e-2 * ones(model.dim.u)) for t = 1:H_mpc],
#         γ = [Diagonal(1.0e-100 * ones(model.dim.c)) for t = 1:H_mpc],
#         b = [Diagonal(1.0e-100 * ones(model.dim.c * friction_dim(env))) for t = 1:H_mpc])

#     p = linearized_mpc_policy(ref_traj, s, obj,
#         H_mpc = H_mpc,
#         N_sample = N_sample,
#         κ_mpc = κ_mpc,
#         n_opts = NewtonOptions(
#             r_tol = 3e-4,
#             max_iter = 5,
#             # verbose = true,
#             ),
#         mpc_opts = LinearizedMPCOptions(
#             # live_plotting=true,
#             altitude_update = true,
#             altitude_impact_threshold = 0.05,
#             altitude_verbose = false,
#             )
#         )

#     q1_ref = copy(ref_traj.q[2])
#     q0_ref = copy(ref_traj.q[1])
#     q1_sim = SVector{model.dim.q}(q1_ref)
#     q0_sim = SVector{model.dim.q}(copy(q1_sim - (q1_ref - q0_ref) / N_sample))
#     @assert norm((q1_sim - q0_sim) / h_sim - (q1_ref - q0_ref) / h) < 1.0e-8

#     sim = simulator(s_sim, q0_sim, q1_sim, h_sim, H_sim,
#         p = p,
#         ip_opts = InteriorPointOptions(
# 			γ_reg = 0.0,
# 			undercut = Inf,
#             r_tol = 1.0e-8,
#             κ_tol = 1.0e-8),
#         sim_opts = SimulatorOptions(warmstart = true)
#         )

#     @time status = simulate!(sim)
#     ref_traj = deepcopy(ref_traj_)

#     qerr, uerr, γerr, berr = ContactImplicitMPC.tracking_error(ref_traj, sim.traj, N_sample, idx_shift=[1])
#     @test qerr < 0.0333 * 1.5 # 0.0333
#     @test uerr < 0.0437 * 1.5 # 0.0437
#     @test γerr < 0.3810 * 1.5 # 0.3810
#     @test berr < 0.0795 * 1.5 # 0.0795
#     qerr > 0.0333 * 1.2 && @warn "mild regression on q tracking: current tracking error = $qerr, nominal tracking error = 0.0333"
#     uerr > 0.0437 * 1.2 && @warn "mild regression on u tracking: current tracking error = $uerr, nominal tracking error = 0.0437"
#     γerr > 0.3810 * 1.2 && @warn "mild regression on γ tracking: current tracking error = $γerr, nominal tracking error = 0.381"
#     berr > 0.0795 * 1.2 && @warn "mild regression on b tracking: current tracking error = $berr, nominal tracking error = 0.0795"
# end

@testset "MPC Quadruped (long trajectory)" begin

	s = get_simulation("quadruped", "flat_2D_lc", "flat")
	model = s.model
	env = s.env

	ref_traj = deepcopy(ContactImplicitMPC.get_trajectory(s.model, s.env,
	    joinpath(module_dir(), "src/dynamics/quadruped/gaits/gait2.jld2"),
	    load_type = :split_traj_alt))

	# time
	H = ref_traj.H
	h = ref_traj.h
	N_sample = 5
	H_mpc = 10
	h_sim = h / N_sample
	H_sim = 9000

	# barrier parameter
	κ_mpc = 2.0e-4

	obj_mpc = TrackingObjective(model, env, H_mpc,
	    q = [Diagonal(1e-2 * [1.0; 0.02; 0.25; 0.25 * ones(model.dim.q-3)]) for t = 1:H_mpc],
	    u = [Diagonal(3e-2 * ones(model.dim.u)) for t = 1:H_mpc],
	    γ = [Diagonal(1.0e-100 * ones(model.dim.c)) for t = 1:H_mpc],
	    b = [Diagonal(1.0e-100 * ones(model.dim.c * friction_dim(env))) for t = 1:H_mpc])

	p = linearized_mpc_policy(ref_traj, s, obj_mpc,
	    H_mpc = H_mpc,
	    N_sample = N_sample,
	    κ_mpc = κ_mpc,
		mode = :configuration,
	    n_opts = NewtonOptions(
			solver = :lu_solver,
			r_tol = 3e-4,
			max_iter = 5,
			# max_time = ref_traj.h, # HARD REAL TIME
			),
	    mpc_opts = LinearizedMPCOptions(),
	    )

	q1_ref = copy(ref_traj.q[2])
	q0_ref = copy(ref_traj.q[1])
	q1_sim = SVector{model.dim.q}(q1_ref)
	q0_sim = SVector{model.dim.q}(copy(q1_sim - (q1_ref - q0_ref) / N_sample))
	@assert norm((q1_sim - q0_sim) / h_sim - (q1_ref - q0_ref) / h) < 1.0e-8

	sim = ContactImplicitMPC.simulator(s, q0_sim, q1_sim, h_sim, H_sim,
	    p = p,
	    ip_opts = ContactImplicitMPC.InteriorPointOptions(
			undercut = Inf,
			γ_reg = 0.0,
	        r_tol = 1.0e-8,
	        κ_tol = 1.0e-8,),
	    sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = true))

	status = ContactImplicitMPC.simulate!(sim, verbose = false)

	qerr, uerr, γerr, berr = tracking_error(ref_traj, sim.traj, N_sample, idx_shift = [1])
	@test qerr < 0.0202 * 1.5 # 0.0201
	@test uerr < 0.0437 * 1.5 # 0.0437
	@test γerr < 0.3780 * 1.5 # 0.3790
	@test berr < 0.0799 * 1.5 # 0.0798
	qerr > 0.0201 * 1.2 && @warn "mild regression on q tracking: current tracking error = $qerr, nominal tracking error = 0.0202"
	uerr > 0.0437 * 1.2 && @warn "mild regression on u tracking: current tracking error = $uerr, nominal tracking error = 0.0437"
	γerr > 0.3790 * 1.2 && @warn "mild regression on γ tracking: current tracking error = $γerr, nominal tracking error = 0.3780"
	berr > 0.0798 * 1.2 && @warn "mild regression on b tracking: current tracking error = $berr, nominal tracking error = 0.0799"
end

# @testset "MPC quadruped: long trajectory with structured newton solver" begin
#
# 	s = get_simulation("quadruped", "flat_2D_lc", "flat")
# 	model = s.model
# 	env = s.env
#
# 	ref_traj = deepcopy(ContactImplicitMPC.get_trajectory(s.model, s.env,
# 	    joinpath(module_dir(), "src/dynamics/quadruped/gaits/gait2.jld2"),
# 	    load_type = :split_traj_alt))
#
# 	# time
# 	H = ref_traj.H
# 	h = ref_traj.h
# 	N_sample = 5
# 	H_mpc = 10
# 	h_sim = h / N_sample
# 	H_sim = 9000
#
# 	# barrier parameter
# 	κ_mpc = 1.0e-4
#
# 	obj_mpc = quadratic_objective(model, H_mpc,
# 	    q = [Diagonal(1e-2 * [1.0; 0.02; 0.25; 0.25 * ones(model.dim.q-3)]) for t = 1:H_mpc+2],
# 	    v = [Diagonal(0.0 * ones(model.dim.q)) for t = 1:H_mpc],
# 	    u = [Diagonal(3e-2 * ones(model.dim.u)) for t = 1:H_mpc-1])
# 	# obj_mpc = TrackingObjective(model, env, H_mpc,
# 	#     q = [Diagonal(1e-2 * [1.0; 0.02; 0.25; 0.25 * ones(model.dim.q-3)]) for t = 1:H_mpc],
# 	#     u = [Diagonal(3e-2 * ones(model.dim.u)) for t = 1:H_mpc],
# 	#     γ = [Diagonal(1.0e-100 * ones(model.dim.c)) for t = 1:H_mpc],
# 	#     b = [Diagonal(1.0e-100 * ones(model.dim.c * friction_dim(env))) for t = 1:H_mpc])
#
# 	p = linearized_mpc_policy(ref_traj, s, obj_mpc,
# 	    H_mpc = H_mpc,
# 	    N_sample = N_sample,
# 	    κ_mpc = κ_mpc,
# 		mode = :configuration,
# 		ip_type = :mehrotra,
# 		newton_mode = :structure,
# 	    n_opts = NewtonOptions(
# 			solver = :lu_solver,
# 			r_tol = 3e-4,
# 			max_iter = 5,
# 			# max_time = ref_traj.h, # HARD REAL TIME
# 			),
# 	    mpc_opts = LinearizedMPCOptions(),
# 		ip_opts = InteriorPointOptions(
# 			max_iter = 100,
# 			verbose = false,
# 			r_tol = 1.0e-4,
# 			κ_tol = 1.0e-4,
# 			diff_sol = true,
# 			# κ_reg = 1e-3,
# 			# γ_reg = 1e-1,
# 			solver = :empty_solver,
# 			),
# 	    )
#
# 	q1_ref = copy(ref_traj.q[2])
# 	q0_ref = copy(ref_traj.q[1])
# 	q1_sim = SVector{model.dim.q}(q1_ref)
# 	q0_sim = SVector{model.dim.q}(copy(q1_sim - (q1_ref - q0_ref) / N_sample))
# 	@assert norm((q1_sim - q0_sim) / h_sim - (q1_ref - q0_ref) / h) < 1.0e-8
#
# 	sim = ContactImplicitMPC.simulator(s, q0_sim, q1_sim, h_sim, H_sim,
# 	    p = p,
# 	    ip_opts = ContactImplicitMPC.InteriorPointOptions(
# 	        r_tol = 1.0e-8,
# 	        κ_tol = 1.0e-6,
# 	        diff_sol = true),
# 	    sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = true))
#
# 	status = ContactImplicitMPC.simulate!(sim, verbose = true)
#
# 	qerr, uerr, γerr, berr = tracking_error(ref_traj, sim.traj, N_sample, idx_shift = [1])
# 	@test qerr < 0.0202 * 1.5 # 0.0202
# 	@test uerr < 0.0437 * 1.5 # 0.0437
# 	@test γerr < 0.3780 * 1.5 # 0.3780
# 	@test berr < 0.0799 * 1.5 # 0.0799
# 	qerr > 0.0201 * 1.2 && @warn "mild regression on q tracking: current tracking error = $qerr, nominal tracking error = 0.0202"
# 	uerr > 0.0437 * 1.2 && @warn "mild regression on u tracking: current tracking error = $uerr, nominal tracking error = 0.0437"
# 	γerr > 0.3790 * 1.2 && @warn "mild regression on γ tracking: current tracking error = $γerr, nominal tracking error = 0.3780"
# 	berr > 0.0798 * 1.2 && @warn "mild regression on b tracking: current tracking error = $berr, nominal tracking error = 0.0799"
#
# end


text:
using MakieLayout
using Makie
using Colors

begin
    scene = Scene(camera = campixel!)
    display(scene)

    g = GridLayout(scene, alignmode=Outside(30))

    rowsize!(g, 1, Auto(false))
    colsize!(g, 1, Auto(false))
    ax = g[1, 1] = LAxis(scene)

    nlines = 10
    xs = 1:0.1:10
    linearr = [lines!(ax, xs, sin.(xs .+ i/3) .* 5, color = col, linestyle = :dash)
        for (i, col) in enumerate(LCHuv.(50, 80, LinRange(0, 360, nlines + 1)[1:end-1]))]

    xs2 = 1:0.5:10
    scatarr = [scatter!(ax, xs2, sin.(xs2 .+ i/3) .* 5, color = col, strokewidth=5, marker = rand(('◀', '■', :circle, :cross, :x)))
        for (i, col) in enumerate(LCHuv.(50, 80, LinRange(0, 360, nlines + 1)[1:end-1]))]

    ll = g[1, 2] = LLegend(scene)

    sg = g[2, :] = GridLayout()

    sg[1, 1] = LText(scene, "Columns", halign=:left)
    colslider = sg[1, 2] = LSlider(scene, range=1:5, height=30)
    on(colslider.value) do v; ll.ncols = v; end;

    sg[2, 1] = LText(scene, "Patch Width", halign=:left)
    patchslider_h = sg[2, 2] = LSlider(scene, range=10:80, height=30)
    on(patchslider_h.value) do v; ll.patchsize = Base.setindex(ll.patchsize[], v, 1); end;

    sg[3, 1] = LText(scene, "Patch Height", halign=:left)
    patchslider_v = sg[3, 2] = LSlider(scene, range=10:80, height=30)
    on(patchslider_v.value) do v; ll.patchsize = Base.setindex(ll.patchsize[], v, 2); end;

    bg = sg[4, :] = GridLayout()
    addbutton = bg[1, 1] = LButton(scene, label = "Add legend entry",
        height = Auto())
    delbutton = bg[1, 2] = LButton(scene, label = "Remove legend entry",
        height = Auto())

    on(addbutton.clicks) do c
        i_randline = rand(1:length(linearr))
        ll.entries[] = [
            ll.entries[];
            LegendEntry("Line $i_randline", linearr[i_randline], scatarr[i_randline])
        ]
    end
    on(delbutton.clicks) do c
        ll.entries[] = ll.entries[][1:end-1];
    end

    nothing
end

g[1, 2] = ll
lg = gridnest!(g, 1, 2)
ll2 = lg[2, 1] = LLegend(scene)
ll2.valign = :top

rowsize!(lg, 1, Auto(true))

ll2.entries[] = LegendEntry.(["Size $i" for i in 1:5:30], [MarkerElement(color = :black, markersize = i, marker = '⚫', strokecolor = :transparent) for i in 1:5:30])

ll.entries[] = [ll.entries[];
    MakieLayout.LegendEntry(
        MakieLayout.LegendElement[
            MakieLayout.LineElement(; color = :green, linestyle = :dot)], Attributes(label = "Blibli"))]

ll.title = "Long\ntitle"
ll.entries[] = LegendEntry.(["Line $i" for i in 1:length(linearr)], linearr, scatarr)

linearr[1].color = :black
ll.bgcolor = :gray
ll.halign = :left
ll.valign = :top
ll.margin = (10, 10, 10, 10)

ll.strokecolor = :transparent
ll.patchcolor = :white
ll.titlesize = 24
ll.entries[][1].labelsize = 20
ll.entries[][1].labelcolor = :red
ll.labelcolor = :black
ll.entries[][1].label = "whoopdeedoo"
ll.labelsize = 18
ll.patchsize = (40, 40)
ll.labelhalign = :left
ll.labelvalign = :center
ll.patchlabelgap = 5
ll.rowgap = 5
ll.colgap = 10
ll.markerpoints = Point2f0.([0.2, 0.5, 0.8], [0.2, 0.5, 0.8])
ll.linepoints = [Point2f0(0, 0), Point2f0(1, 1)]
ll.entries[][1].linepoints = Point2f0.(LinRange(0, 1, 40), sin.(LinRange(0, 2pi, 40)) .* 0.5 .+ 0.5)

ll.entries[] = reverse(ll.entries[])
using Random
ll.entries[] = shuffle!(ll.entries[])
ll.entries[] = ll.entries[][1:end-1]; display(scene);
scene
ll.patchsize = 60
ll.halign = rand(setdiff([:left, :center, :right], [ll.halign[]]))
ll.valign = rand(setdiff([:bottom, :center, :top], [ll.valign[]]))
ll.margin = rand(4) .* 0

translate!(ll.scene, 0, 0, 10)

linearr[2].color = :red
scatarr[2].color = :red
scatarr[2].marker = 'm'

@macroexpand @lift($nodes[1])


text:
@testset "Misc/NumberField" begin
  @testset "issubfield" begin
    Qx, x = FlintQQ["x"]
    K, a = NumberField(x^2 + 1, "a")
    L, b = NumberField(x^4 + 1, "b")

    c, KtoL = Hecke.issubfield(K, L)
    @test c == true
    @test parent(KtoL(a)) == L

    c, KtoL = Hecke.issubfield_normal(K, L)
    @test c == true
    @test parent(KtoL(a)) == L

    OK = maximal_order(K)
    OL = maximal_order(L)
    c, KtoL = Hecke.issubfield(K, L)
    @test c == true
    @test parent(KtoL(a)) == L
  end

  @testset "isisomorphic" begin
    Qx, x = FlintQQ["x"]
    f = x^5 + 12x - 92
    K, a = NumberField(f, "a")

    g = x^5 - 172x^4 + 7024x^3 + 8656448x^2 + 55735552x + 45796197888
    K2, a2 = NumberField(g, "a2")

    c, KtoK2 = Hecke.isisomorphic(K, K2)
    @test c == true
    @test parent(KtoK2(a)) == K2

    OK = maximal_order(K)
    OK2 = maximal_order(K2)
    c, KtoK2 = Hecke.isisomorphic(K, K2)
    @test c == true
    @test parent(KtoK2(a)) == K2

    h = f - 1
    K3, a3 = NumberField(h, "a3")
    d, KtoK3 = Hecke.isisomorphic(K, K3)
    @test d == false
  end
end


text:
root_dir = "../."
using Pkg; Pkg.activate(root_dir)
import DotEnv

function main()
    cfg = DotEnv.config(path="$root_dir/.env")
    println(cfg["INPUT_DIR"])
    println(ARGS[end])
end
open("../data/opts.csv", "w") do io
    cols = join(["coin", "nslow","nfast","return"], ',')
    println(io, cols)
    for (key, value) in d
        println(io, join([key, join(value, ',')], ','))
    end
end

main()

max_stars_repo_name: lithiferous/Strathreads

max_issues_repo_name: lithiferous/Strathreads

max_forks_repo_name: lithiferous/Strathreads


text:
#=
   False coin in Turing.jl

   Inspired by Peter Winkler in
   From Celebration of Mind (Gathering For Gardner):
   "CoM Apr 2021 - Drawing from Urns - Peter Winkler"
   https://www.youtube.com/watch?v=A0GxvYPTKDk&list=PL5_D5QQxvWKVtX4Vklvhsf5KTYZBLDSwO&index=3

   Note: This is not what Peter talked about, but what I - wrongly -
   remembered.

   *  We have one fake and one proper coin:
   1. head, head
   2. head, tail

  * We take one coin randomly and flip it: it's a head
  * To identify the fake coin, is it better to flip the same
    coin or flip the other coin?

  According to this model, it's better to throw the same coin again, i.e.
    coinB_switchIsFake: 1/3
  vs
    coinB_sameIsFake: 2/3

  """
  Distributions of variable coinB_switchIsFake (num:0)
  0.00000 =>    6586  (0.658600)
  1.00000 =>    3414  (0.341400)

  Distributions of variable coinB_sameIsFake (num:0)
  1.00000 =>    6586  (0.658600)
  0.00000 =>    3414  (0.341400)
  """

  Note: The problem that Peter Winkler's talked about was a little different:
  @3:43
  """
  You have two quarters.
  Both are sitting on a table heads up.
  One of them is an ordinary quarter, but the
  other has a head on the other side as well.

  You get two flips to help you guess which is which.

  Should you flip each coin once or one coin twice.
  """

  See false_coin2.jl for this version.

=#

using Turing
include("jl_utils.jl")

@model function false_coin()
    #  The two coins
    coin1 = 1 #  fake coin: head,head
    coin2 = 2 #  proper coin: head,tail
    coins = [coin1,coin2];

    #  Sides of the coins
    head = 1
    tail = 2

    function flip_coin(cond)
        return cond ? Categorical([1.0,0.0]) :  Categorical([0.5,0.5])
    end

    #  1. Pick a coin randomly and flip it. It's head.
    coinA ~ Categorical([0.5,0.5])
    sideA ~ flip_coin(coinA == coin1)

    #  We observe that the first flip was a head.
    true ~ Dirac(sideA == head)

    # 2a: Take the other coin and flip it.
    coinB_switch ~ Dirac(coinA == coin1 ? coin2 : coin1)
    sideB_switch ~ flip_coin(coinB_switch == coin1)

    # 2b: Take the same coin and flip it.
    coinB_same ~ Dirac(coinA)
    sideB_same ~ flip_coin(coinB_same == coin1);

    # The probability of fake coin
    coinB_switchIsFake ~ Dirac(coinB_switch==coin1)
    coinB_sameIsFake ~ Dirac(coinB_same==coin1)

end

model = false_coin()
# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, IS(), 10_000)
# chns = sample(model, SMC(), 10_000)

display(chns)
show_var_dist_pct(chns,:coinB_switchIsFake)
show_var_dist_pct(chns,:coinB_sameIsFake)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Test
using LinearAlgebra: norm
using Random: seed!

using SimplexMethod

@testset "Smoke" begin
    A = Float64[
        1   2
        1   1
        3   2
    ]
    b = Float64[16, 9, 24]
    c = Float64[40, 30]
    @test simplex_method(A, b, c)[1] ≈ 330.0

    A = Float64[
        1   1   3
        2   2   5
        4   1   2
    ]
    b = Float64[30, 24, 36]
    c = Float64[3, 1, 2]
    @test simplex_method(A, b, c)[1] ≈ 28.0

    A = Float64[
        2   1
        2   3
        3   1
    ]
    b = Float64[18, 42, 24]
    c = Float64[3, 2]
    @test simplex_method(A, b, c)[1] ≈ 30.0
end

@testset "Randomized" begin
    # JuliaOpt/JuMP seems too heavyweight so roll my own test routine.
    seed!()
    N_TESTS = 10
    N_PERTURBATIONS = 5
    for _ in 1:N_TESTS
        n = rand(10:20)
        m = rand(1:(n - 1))
        A = rand(m, n)
        b = rand(m)
        c = abs.(rand(n))

        (z, x) = simplex_method(A, b, c)
        @test c' * x ≈ z
        lambda = 1.0e-3

        for _ in 1:N_PERTURBATIONS
            delta = rand(n)
            delta .*= lambda / norm(delta)
            xd = x + delta
            @test any(xd .< 0.0) || any(A * xd .> b) || c' * xd <= z
        end
    end
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
using Test
using OpenQuantumSystems
using SparseArrays, LinearAlgebra

@testset "interaction picture" begin

    D(op1::Array, op2::Array) = abs(norm(op1 - op2))
    D(x1::StateVector, x2::StateVector) = norm(x2 - x1)
    D(op1::AbstractOperator, op2::AbstractOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))
    D(op1::AbstractSuperOperator, op2::AbstractSuperOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))

    # TODO: change to macro
    mode1 = Mode(0.2, 1.0)
    mode2 = Mode(0.3, 2.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1], 3, [2.0, 200.0])
    mol2 = Molecule([mode2], 3, [3.0, 300.0])
    aggCore = AggregateCore([mol1, mol2])
    aggCore.coupling[2, 3] = 50
    aggCore.coupling[3, 2] = 50
    agg = setupAggregate(aggCore)

    Ham = agg.operators.Ham

    Ham_I = agg.operators.Ham_I
    Ham_0 = agg.operators.Ham_0
    Ham = agg.operators.Ham

    Ham_0_lambda, Ham_0_S = eigen(Ham_0.data)
    Ham_0_Sinv = inv(Ham_0_S)
    Ham_0_lambda = diagm(Ham_0_lambda)

    for t in [0.0, 1.0, 2.0]
        U_op = evolutionOperator(Ham_0, t)
        U_II_t_ref = U_op' * Ham_I * U_op

        U_II_t = getInteractionHamIPicture(Ham_0.data, Ham_I.data, t)
        @test 1e-14 > D(U_II_t_ref.data, U_II_t)
        U_II_t = getInteractionHamIPicture(Ham_0, Ham_I, t)
        @test 1e-14 > D(U_II_t_ref, U_II_t)
        U_II_t = getInteractionHamIPictureA(Ham_0, Ham_I, t)
        @test 1e-14 > D(U_II_t_ref.data, U_II_t)
        U_II_t = getInteractionHamIPictureA(Ham_I.data, Ham_0_lambda, Ham_0_S, Ham_0_Sinv, t)
        @test 1e-10 > D(U_II_t_ref.data, U_II_t)
    end

end # testset

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
module TidalFluxConfigurations

const config = Dict{Symbol,Any}()

end # module


text:
Nfile = 8

myFiles = 1:Nfile

for i = myFiles
    tic()
    fname = "data/"* string(i)

    fidRow = open(fname * "r.txt","r")
    fidCol = open(fname * "c.txt","r")
    fidVal = open(fname * "v.txt","r")

    rowStr = readline(fidRow)
    colStr = readline(fidCol)
    valStr = readline(fidVal)

    close(fidRow)
    close(fidCol)
    close(fidVal)

    A = Assoc(rowStr,colStr,1,(sum))

    pTime = toq()
    println("Sum Time: ", pTime)
    println("Edges/sec: ", string(NumStr(rowStr)/pTime))
end


text:
"""
Luby Transform Code

# References

- David J. C. MacKay. 2002. Information Theory, Inference & Learning Algorithms (pp. 589-596).
  Cambridge University Press, New York, NY, USA.
- https://franpapers.com/en/algorithmic/2018-introduction-to-fountain-codes-lt-codes-with-python/
"""
module LubyTransformCode

include("lehmer.jl")
include("soliton.jl")
include("core.jl")
include("encode.jl")
include("decode.jl")

end


text:
# V-nodes `[x_1, x_2, x_3, ..., dx_1, dx_2, ..., y_1, y_2, ...]` where `x`s are
# differential variables and `y`s are algebraic variables.
function get_vnodes(sys)
    dxvars = Operation[]
    edges = map(_->Int[], 1:length(sys.eqs))
    for (i, eq) in enumerate(sys.eqs)
        if !(eq.lhs isa Constant)
            # Make sure that the LHS is a first order derivative of a var.
            @assert eq.lhs.op isa Differential
            @assert !(eq.lhs.args[1] isa Differential) # first order

            push!(dxvars, eq.lhs)
            # For efficiency we note down the diff edges here
            push!(edges[i], length(dxvars))
        end
    end

    xvars = (first ∘ var_from_nested_derivative).(dxvars)
    algvars  = setdiff(states(sys), xvars)
    return xvars, dxvars, edges, algvars
end

function sys2bigraph(sys)
    xvars, dxvars, edges, algvars = get_vnodes(sys)
    xvar_offset = length(xvars)
    algvar_offset = 2xvar_offset
    for edge in edges
        isempty(edge) || (edge .+= xvar_offset)
    end

    for (i, eq) in enumerate(sys.eqs)
        # T or D(x):
        # We assume no derivatives appear on the RHS at this point
        vs = vars(eq.rhs)
        for v in vs
            for (j, target_v) in enumerate(xvars)
                if v == target_v
                    push!(edges[i], j)
                end
            end
            for (j, target_v) in enumerate(algvars)
                if v == target_v
                    push!(edges[i], j+algvar_offset)
                end
            end
        end
    end

    fullvars = [xvars; dxvars; algvars] # full list of variables
    vars_asso = [(1:xvar_offset) .+ xvar_offset; zeros(Int, length(fullvars) - xvar_offset)] # variable association list
    return edges, fullvars, vars_asso
end

print_bigraph(sys, vars, edges) = print_bigraph(stdout, sys, vars, edges)
function print_bigraph(io::IO, sys, vars, edges)
    println(io, "Equations:")
    foreach(x->println(io, x), [i => sys.eqs[i] for i in 1:length(sys.eqs)])
    for (i, edge) in enumerate(edges)
        println(io, "\nEq $i has:")
        print(io, '[')
        for e in edge
            print(io, "$(vars[e]), ")
        end
        print(io, ']')
    end
    return nothing
end

function match_equation!(edges, i, assign, active, vcolor=falses(length(active)), ecolor=falses(length(edges)))
    # `edge[active]` are active edges
    # i: equations
    # j: variables
    # assign: assign[j] == i means (i-j) is assigned
    #
    # color the equation
    ecolor[i] = true
    # if a V-node j exists s.t. edge (i-j) exists and assign[j] == 0
    for j in edges[i]
        if active[j] && assign[j] == 0
            assign[j] = i
            return true
        end
    end
    # for every j such that edge (i-j) exists and j is uncolored
    for j in edges[i]
        (active[j] && !vcolor[j]) || continue
        # color the variable
        vcolor[j] = true
        if match_equation!(edges, assign[j], assign, active, vcolor, ecolor)
            assign[j] = i
            return true
        end
    end
    return false
end

function matching(edges, nvars, active=trues(nvars))
    assign = zeros(Int, nvars)
    for i in 1:length(edges)
        match_equation!(edges, i, assign, active)
    end
    return assign
end

function pantelides(sys::ODESystem; kwargs...)
    edges, fullvars, vars_asso = sys2bigraph(sys)
    return pantelides!(edges, fullvars, vars_asso; kwargs...)
end

function pantelides!(edges, vars, vars_asso; maxiter = 8000)
    neqs = length(edges)
    nvars = length(vars)
    assign = zeros(Int, nvars)
    eqs_asso = fill(0, neqs)
    neqs′ = neqs
    for k in 1:neqs′
        i = k
        pathfound = false
        # In practice, `maxiter=8000` should never be reached, otherwise, the
        # index would be on the order of thousands.
        for _ in 1:maxiter
            # run matching on (dx, y) variables
            active = vars_asso .== 0
            vcolor = falses(nvars)
            ecolor = falses(neqs)
            pathfound = match_equation!(edges, i, assign, active, vcolor, ecolor)
            pathfound && break # terminating condition
            # for every colored V-node j
            for j in eachindex(vcolor); vcolor[j] || continue
                # introduce a new variable
                nvars += 1
                push!(vars_asso, 0)
                vars_asso[j] = nvars
                push!(assign, 0)
            end

            # for every colored E-node l
            for l in eachindex(ecolor); ecolor[l] || continue
                neqs += 1
                # create new E-node
                push!(edges, copy(edges[l]))
                # create edges from E-node `neqs` to all V-nodes `j` and
                # `vars_asso[j]` s.t. edge `(l-j)` exists
                for j in edges[l]
                    if !(vars_asso[j] in edges[neqs])
                        push!(edges[neqs], vars_asso[j])
                    end
                end
                push!(eqs_asso, 0)
                eqs_asso[l] = neqs
            end

            # for every colored V-node j
            for j in eachindex(vcolor); vcolor[j] || continue
                assign[vars_asso[j]] = eqs_asso[assign[j]]
            end
            i = eqs_asso[i]
        end # for _ in 1:maxiter
        pathfound || error("maxiter=$maxiter reached! File a bug report if your system has a reasonable index (<100), and you are using the default `maxiter`. Try to increase the maxiter by `pantelides(sys::ODESystem; maxiter=1_000_000)` if your system has an incredibly high index and it is truly extremely large.")
    end # for k in 1:neqs′
    return edges, assign, vars_asso, eqs_asso
end


text:
const list_section_template = mt"""
{{{description}}}

```@raw html
<div class="list-card-section">
```

{{{cards}}}

```@raw html
</div>
```
"""

const list_card_template = mt"""
```@raw html
<div class="list-card">
<table>
  <td valign="bottom"><div class="list-card-cover">
```
[![list-card-cover-image]({{{coverpath}}})](@ref {{id}})
```@raw html
  </div></td>
  <td><div class="list-card-text">
```
[{{{title}}}](@ref {{id}})
```@raw html
</div>
    <div class="list-card-description">
```
{{{description}}}
```@raw html
    </div>
  </td>
</tbody></table>
</div>
```

"""

function cardtheme(::Val{:list})
    templates = Dict(
        "card" => list_card_template,
        "section" => list_section_template
    )
    return templates, abspath(@__DIR__, "style.css")
end