text:
# AIM 3 Analysis:
# What are the relative frequencies of serious vs. non-serious adverse events (and general, all AEs)
# for different groups of vaccine distributors?

#loading librarys
using DataFrames, CSV, Missings, Pipe, StatsBase, Query,  HypothesisTests

# setup filepaths
path_source = string(@__DIR__,"\\..\\source");
path_dev = string(@__DIR__,"\\..\\dev");

# load the labelled data
df = DataFrame(CSV.read(joinpath(path_dev, "19-21VAERSCOMB.csv", DataFrame));

# rewriting the get frequency code into function
# this function takes in one column name(default), or two column names
# and calculate the overall frequency/counts for col-AE (or col1-col2-AE) pairs
function get_freq(df, col1, col2 = nothing)

    # when only one column name specified
    if isnothing(col2)

        # one observation is associated with 5 symptom columns,
        # calculate the frequency/counts for col1-AE pairs separately by symptom column
        count1 = dropmissing(rename!(combine(groupby(df, [Symbol(col1), :SYMPTOM1]),
                    nrow => :COUNT), :SYMPTOM1 => :SYMPTOM), :SYMPTOM);
        count2 = dropmissing(rename!(combine(groupby(df, [Symbol(col1), :SYMPTOM2]),
                    nrow => :COUNT), :SYMPTOM2 => :SYMPTOM), :SYMPTOM);
        count3 = dropmissing(rename!(combine(groupby(df, [Symbol(col1), :SYMPTOM3]),
                    nrow => :COUNT), :SYMPTOM3 => :SYMPTOM), :SYMPTOM);
        count4 = dropmissing(rename!(combine(groupby(df, [Symbol(col1), :SYMPTOM4]),
                    nrow => :COUNT), :SYMPTOM4 => :SYMPTOM), :SYMPTOM);
        count5 = dropmissing(rename!(combine(groupby(df, [Symbol(col1), :SYMPTOM5]),
                    nrow => :COUNT), :SYMPTOM5 => :SYMPTOM), :SYMPTOM);

        # concatenating, note that col1-AE pairs may not be unique at this step
        count_raw = reduce(vcat, [count1, count2, count3, count4, count5]);

        # recalculating frequency/counts by summing up the counts for non-unique col1-AE pairs
        count = combine(groupby(count_raw, [Symbol(col1), :SYMPTOM]), :COUNT => sum)

    # when two column names specified
    else
        count1 = dropmissing(rename!(combine(groupby(df, [Symbol(col1), Symbol(col2), :SYMPTOM1]),
                    nrow => :COUNT), :SYMPTOM1 => :SYMPTOM), :SYMPTOM);
        count2 = dropmissing(rename!(combine(groupby(df, [Symbol(col1), Symbol(col2), :SYMPTOM2]),
                    nrow => :COUNT), :SYMPTOM2 => :SYMPTOM), :SYMPTOM);
        count3 = dropmissing(rename!(combine(groupby(df, [Symbol(col1), Symbol(col2), :SYMPTOM3]),
                    nrow => :COUNT), :SYMPTOM3 => :SYMPTOM), :SYMPTOM);
        count4 = dropmissing(rename!(combine(groupby(df, [Symbol(col1), Symbol(col2), :SYMPTOM4]),
                    nrow => :COUNT), :SYMPTOM4 => :SYMPTOM), :SYMPTOM);
        count5 = dropmissing(rename!(combine(groupby(df, [Symbol(col1), Symbol(col2), :SYMPTOM5]),
                    nrow => :COUNT), :SYMPTOM5 => :SYMPTOM), :SYMPTOM);

        count_raw = reduce(vcat, [count1, count2, count3, count4, count5]);

        count = combine(groupby(count_raw, [Symbol(col1), Symbol(col2), :SYMPTOM]), :COUNT => sum)
    end
    return count
end

#get the manufacturer w/ the most frequent aderse event overall
tb1 = sort!(get_freq(df, "VAX_MANU"), :COUNT_sum, rev = true)

#get the manufacturer with most frequent non-serious adverse event overall
tb2 = sort!(filter(x -> x.SERIOUS_EVENT == 0, get_freq(df, "VAX_MANU", "SERIOUS_EVENT")), :COUNT_sum, rev = true)

#get the manufacturer with most frequent non-serious adverse event overall
tb3 = sort!(filter(x -> x.SERIOUS_EVENT == 1, get_freq(df, "VAX_MANU", "SERIOUS_EVENT")), :COUNT_sum, rev = true)

# top 10 manufacturers people got vaccines from
manu10 = sort!(combine(groupby(df, :VAX_MANU), nrow => :CT), :CT, rev=true)[1:11, 1][1:end .!= 7]

# get relative frequency of AE by manufacturer, sorted
filtered = filter(row -> row.VAX_MANU ∈ manu10, tb1);
sort(transform(groupby(filtered, :VAX_MANU), :COUNT_sum => (x -> x / sum(x)) => :freq), :freq, rev=true)

# get relative frequency of non-serious AE by manufacturer, sorted
filtered_n_serious = filter(row -> row.VAX_MANU ∈ manu10, tb2);
show(sort(transform(groupby(filtered_n_serious, :VAX_MANU), :COUNT_sum => (x -> x / sum(x)) => :freq), :freq, rev=true), allcols = true)

# get relative frequency of serious AE by manufacturer, sorted
filtered_serious = filter(row -> row.VAX_MANU ∈ manu10, tb3);
show(sort(transform(groupby(filtered_serious, :VAX_MANU), :COUNT_sum => (x -> x / sum(x)) => :freq), :freq, rev=true), allcols = true)

# hypothesis testing
# perform unequal variance t test to investigate the difference in mean # AE reports across distributors

# get counts
moderna = filter(x -> x.VAX_MANU .== "MODERNA", tb1).COUNT_sum;
pfizer = filter(x -> x.VAX_MANU .== "PFIZER\\BIONTECH", tb1).COUNT_sum;
seqirus = filter(x -> x.VAX_MANU .== "SEQIRUS, INC.", tb1).COUNT_sum;

# moderna vs. pfizer (top 2 distributors having the highest relative freq of non-serious AE)
UnequalVarianceTTest(moderna, pfizer)

# moderna vs. seqirus (top 2 distributors having the highest relative freq of serious AE)
UnequalVarianceTTest(moderna, seqirus)


text:
"""
    EconJobMarket
An interface with the [EconJobMarket](https://econjobmarket.org/pages/about) [API](https://backend.econjobmarket.org/new_api).
"""
module EconJobMarket
    using Cascadia: nodeText, parsehtml, Selector
    using CSV: CSV, Tables, normalizename, Tables.columntable
    using Dates: Date
    using HTTP: HTTP
    using JSON3: JSON3

    foreach(include,
            [ filename for filename ∈ readdir(@__DIR__)
              if occursin(r"^\d+.*\.jl$", filename) ] |>
            sort!)

    export CSV, Tables,
           Date,
           fetch_ads
end

ext: jl


text:
# lex:
#
# contains algorithms in the lex-minimizer paper
# includes iterative lex algorithm
#
# Started by xiao.shi@yale.edu on Sep 22, 2015
# Contributers:

import Laplacians.intHeap
import Laplacians.intHeapAdd!, Laplacians.intHeapPop!

LEX_DEBUG = false # debug flag for the lex module
LEX_EPS = 1e-12 # default error tolerance

function setLexDebugFlag(f::Bool)
  global LEX_DEBUG = f
end

#=
Simulate IterLex on Uniformly Weighted Graphs

numIter: number of iterations;

A: n by n adjacency matrix of the graph;

isTerm: boolean array[n] of whether each vertex is a terminal;

initVal: array[n] of initial voltage assignments;

=#
function simIterLexUnwtd{Tv, Ti}(numIter::Int64,
                                 A::SparseMatrixCSC{Tv, Ti},
                                 isTerm::Array{Bool, 1},
                                 initVal::Array{Float64, 1}, )
  n = A.n
  val = copy(initVal)
  nextVal = zeros(Float64, n)

  for t = 1:numIter
    # if the bits representation of val and nextVal are the same for
    # ever vertex, then there is no point in keeping iterating.
    progress = false
    for u = 1:n
      if (!isTerm[u])
        nbrs = A.rowval[A.colptr[u]:(A.colptr[u + 1] - 1)]
        maxNeighbor = maximum(val[nbrs])
        minNeighbor = minimum(val[nbrs])
        nextVal[u] = minNeighbor + (maxNeighbor - minNeighbor) / 2.0
        if (bits(val[u]) != bits(nextVal[u]))
          progress = true
        end
      else
        nextVal[u] = val[u]
      end
    end

    if (!progress)
      @printf("INFO: simIterLexUnwtd: terminating early after %d iterations,
              as numerical error prevents further progress.\n", t)
      return val
    end

    tmp = val
    val = nextVal
    nextVal = tmp
  end
  return val
end

function simIterLex{Tv<:Float64, Ti}(numIter::Int64,
                                     A::SparseMatrixCSC{Tv, Ti},
                                     isTerm::Array{Bool, 1},
                                     initVal::Array{Float64, 1}, )
  n = A.n
  val = copy(initVal)
  nextVal = zeros(Float64, n)

  for t = 1:numIter
    # if the bits representation of val and nextVal are the same for
    # ever vertex, then there is no point in keeping iterating.
    progress = false
    for u = 1:n
      if (!isTerm[u])
        nbrs = A.rowval[A.colptr[u]:(A.colptr[u + 1] - 1)]
        w = A.nzval[A.colptr[u]:(A.colptr[u + 1] - 1)] # weights
        deg = length(nbrs)
        maxGrad = -Inf
        maxI = 0
        maxJ = 0
        for i in 1:deg
          for j in 1:deg
            grad = (val[nbrs[i]] - val[nbrs[j]]) / (1/w[i] + 1/w[j])
            if (grad > maxGrad)
              maxGrad = grad
              maxI = i
              maxJ = j
            end
          end
        end

        wi = w[maxI]
        vi = val[nbrs[maxI]]
        wj = w[maxJ]
        vj = val[nbrs[maxJ]]
        nextVal[u] = (wi * vi + wj * vj) / (wi + wj)
        if (bits(val[u]) != bits(nextVal[u]))
          progress = true
        end
      else
        nextVal[u] = val[u]
      end
    end

    if (!progress)
      @printf("INFO: simIterLexUnwtd: terminating early after %d iterations,
              as numerical error prevents further progress.\n", t)
      return val
    end

    tmp = val
    val = nextVal
    nextVal = tmp
  end
  return val
end

#=
check the correctness of lex assignment for uniformly weighted graphs

lex: the assignment to check;

eps: absolute error tolerance;

fatal: if true, throws error and halt; if false, return false;

=#
function checkLexUnwtd{Tv, Ti}(A::SparseMatrixCSC{Tv, Ti},
                               isTerm::Array{Bool, 1},
                               initVal::Array{Float64, 1},
                               lex::Array{Float64, 1};
                               eps::Float64 = LEX_EPS,
                               fatal::Bool = true)
  n = A.n
  for i in 1:n
    if (isTerm[i])
      expected = initVal[i]
    else
      nbrs = A.rowval[A.colptr[i]:A.colptr[i+1]-1]
      maxval = maximum(lex[nbrs])
      minval = minimum(lex[nbrs])
      expected = minval + (maxval - minval) / 2.0
    end
    if (abs(lex[i] - expected) > eps)
      if (fatal)
        error(@sprintf("The average of the max and min neighbors is %f
                       at vertex %d, but got %f\n", expected, i, lex[i]))
      end
      return false
    end
  end
  return true
end

#=
check the correctness of lex assignment for uniformly weighted graphs

lex: the assignment to check;

eps: absolute error tolerance;

fatal: if true, throws error and halt; if false, return false;

=#
function checkLex{Tv, Ti}(A::SparseMatrixCSC{Tv, Ti},
                          isTerm::Array{Bool, 1},
                          initVal::Array{Float64, 1},
                          lex::Array{Float64, 1};
                          eps::Float64 = LEX_EPS,
                          fatal::Bool = true)
  n = A.n

  for i in 1:n
    correct = true
    if (isTerm[i])
      correct = (lex[i] == initVal[i])
    else
      nbrs = A.rowval[A.colptr[i]:A.colptr[i+1]-1]
      w = A.nzval[A.colptr[i]:A.colptr[i+1]-1] # weights
      deg = length(nbrs)
      grads = zeros(deg)
      for j in 1:deg
        grads[j] = (lex[nbrs[j]] - lex[i]) * w[j]
      end

      maxgrad = maximum(grads)
      mingrad = minimum(grads)

      if (maxgrad + mingrad > eps)
        correct = false
      end
    end
    if (!correct)
      if (fatal)
        error(@sprintf("The average of the max and min neighbors is %f
                       at vertex %d, but got %f\n", expected, i, lex[i]))
      end
      return false
    end
  end
  return true
end

#==========
 LEX TESTS
===========#

function simIterLexTest()
  simIterLexTestPath(4)
  simIterLexTestPath(40)
end

function simIterLexTestPath(n::Int64)
  if (n < 3)
    error("n must be at least 3.")
  end

  numIter = 1000 * n * n

  # path graph with n verticies
  Pn = pathGraph(n)
  # terminal assignments
  isTerm = zeros(Bool, n)
  isTerm[1] = true
  isTerm[n] = true

  initVal = zeros(Float64, n)
  initVal[n] = n - 1

  val = simIterLexUnwtd(numIter, Pn, isTerm, initVal)
  # correct = checkLexUnwtd(Pn, isTerm, initVal, val, fatal=false)

  correct = true
  for i in 1:n
    expected = (i - 1)
    if (abs(val[i] - expected) > LEX_EPS)
      println(val[i], expected)
      println(val)
      error("simIterLex is not correct.")
    end
  end
end

# returns the maximum gradient of an edge in an assignment;
# for ease of verifying CompInfMin
function MaxEdgeGrad(A, val)
  maxGrad = -Inf
  maxX = 0
  maxY = 0
  for x in 1:A.n
    for i in A.colptr[x]:A.colptr[x+1]-1
      y = A.rowval[i]
      grad = abs(val[x]-val[y]) / A[x,y]
      if (grad > maxGrad)
        maxGrad = grad
        maxX = x
        maxY = y
      end
    end
  end
  return maxX, maxY, maxGrad
end

#===============
 Lex Algorithms
===============#

# Shortest Paths _without_ going past terminals.
# Treat weights as reciprocals of distances.
function termFreeShortestPaths{Tv,Ti}(mat::SparseMatrixCSC{Tv, Ti},
                                      start::Ti,
                                      isTerm = zeros(Bool, mat.n), )
  n = mat.n
  visited = zeros(Bool,n)

  nh = intHeap(n)
  dists = nh.keys

  pArray = [1:n...]

  intHeapAdd!(nh, start, 0.0)
  pArray[start] = start

  while nh.nitems > 0
    v::Ti = intHeapPop!(nh)
    visited[v] = true

    dv = dists[v]
    if (isTerm[v] && v != start)
      continue
    end
    for ind in mat.colptr[v]:(mat.colptr[v+1]-1)
      nbr = mat.rowval[ind]
      if !visited[nbr]
        newdist = dv + 1.0 / mat.nzval[ind]
        if newdist < dists[nbr]
          dists[nbr] = newdist
          intHeapAdd!(nh,nbr,newdist)
          pArray[nbr] = v
        end # if
      end # if
    end # for

  end # while
  return copy(dists), pArray
end # termFreeShortestPaths

function ModDijkstra{Tv, Ti}(G::SparseMatrixCSC{Tv, Ti},
                             isTerm::Array{Bool, 1},
                             initVal::Array{Float64, 1},
                             alpha::Float64, )
  n = G.n
  finished = zeros(Bool, n)

  nh = intHeap(n)
  key = nh.keys
  val = copy(initVal)

  parent = [1:n...]
  for i = 1:n
    if (isTerm[i])
      key[i] = initVal[i]
    end
    intHeapAdd!(nh, i, key[i]);
  end

  while (nh.nitems > 0)
    u = intHeapPop!(nh)
    val[u] = key[u]
    finished[u] = true

    for ind in G.colptr[u]:(G.colptr[u+1]-1)
      nbr = G.rowval[ind]
      if (G.nzval[ind] > 0 && !finished[nbr] && !isTerm[nbr])
        newkey = val[u] + alpha / G.nzval[ind]
        if (key[nbr] > newkey)
          intHeapAdd!(nh, nbr, newkey)
          parent[nbr] = u
        end
      end
    end
  end

  return val, parent
end

function CompVLow{Tv, Ti}(G::SparseMatrixCSC{Tv, Ti},
                          isTerm::Array{Bool, 1},
                          initVal::Array{Float64, 1},
                          alpha::Float64)
  return ModDijkstra(G, isTerm, initVal, alpha)
end

function CompVHigh{Tv, Ti}(G::SparseMatrixCSC{Tv, Ti},
                           isTerm::Array{Bool, 1},
                           initVal::Array{Float64, 1},
                           alpha::Float64)
  newInitVal = -copy(initVal)
  temp, parent = ModDijkstra(G, isTerm, newInitVal, alpha)
  return -temp, parent
end

# returns a vertex set (of non terminals)
function CompHighPressGraph(G, isTerm, initVal, alpha)
  vLow, LParent = CompVLow(G, isTerm, initVal, alpha)
  vHigh, HParent = CompVHigh(G, isTerm, initVal, alpha)

  find(vHigh .> vLow * (1 + LEX_EPS))
end

# visited and val will change!
function fixPath!(A, isTerm, val, path)
  l = length(path)
  if (l <= 2)
    return
  end
  if (!isTerm[path[1]])
    error(@sprintf("Value of vertex %d is not set!", path[1]))
  end
  if (!isTerm[path[l]])
    error(@sprintf("Value of vertex %d is not set!", path[l]))
  end

  dist = Array(Float64, l-1)
  cdist = Array(Float64, l-1)
  for i in 1:l-1
    if (A[path[i], path[i+1]] == 0)
      error(@sprintf("No edge between vertices %d and %d!", path[i], path[i+1]))
    end
    dist[i] = A[path[i], path[i+1]]
  end
  cdist = cumsum!(cdist, dist)
  sum = cdist[l-1]

  v1 = val[path[1]]
  v2 = val[path[l]]
  for i in 1:l-2
    d1 = cdist[i]
    d2 = sum - cdist[i]
    val[path[i + 1]] = (d1 * v2 + d2 * v1) / (d1 + d2)
    isTerm[path[i + 1]] = true
  end
end

# check the correctness of lex assignment
function CheckLex(A, isTerm, initVal, lex, eps = 1e-14, fatal = true)
  n = A.n
  correct = true
  for i in 1:n
    if (isTerm[i])
      expected = initVal[i]
    else
      nbrs = A.rowval[A.colptr[i]:A.colptr[i+1]-1]
      maxval = maximum(lex[nbrs])
      minval = minimum(lex[nbrs])
      expected = minval + (maxval - minval) / 2.0
    end
    if (abs(lex[i]-expected) > eps)
      if (fatal)
        @printf("Expect value %f at vertex %d, but got %f\n", expected, i, lex[i])
      end
      return false
    end
  end
  return true
end

# StarSteepestPath: all terminal nodes in T connects to a single non-terminal node u.
# Assignments are v; distances from u to terminals are d.
#
# Naive version: enumerate pairs of vertices t1, t2 in T
# that maximizes (v(t1)-v(t2))/(d(t1)-d(t2)). O(|T|^2)
function NaiveStarSteepestPath(T, v, d)
  max = -Inf
  res1 = 0
  res2 = 0
  for i in 1:length(T)
    for j in (i+1):length(T)
      if (d[i] == Inf || d[j] == Inf)
        continue
      end
      grad = (v[i] - v[j]) / (d[i] + d[j])
      if (grad > max)
        max = grad
        res1 = i
        res2 = j
      end
      if (-grad > max)
        max = -grad
        res1 = j
        res2 = i
      end
    end
  end
  return T[res1], T[res2]
end

# grad error tolerance
gradTol = 1e-14

# Improved version. Expected O(|T|). (Alg 10 in Lex Paper)
function StarSteepestPath(T, v, d)
  if (LEX_DEBUG)
    println("[Star] ", T', " ", v', " ", d')
  end
  numT = length(T)
  if (numT != length(v) || numT != length(d))
    error("The dimensions of T, v, and d do not match!")
  end

  if (numT < 1)
    error("No terminal in the star graph, cannot process!")
  end
  if (numT == 1)
    return (T[1], T[1])
  end
  if (numT == 2)
    return (v[1] > v[2])? (T[1],T[2]) : (T[2], T[1])
  end

  i = rand(1:length(T))
  alpha = -Inf
  i2 = 0
  for j in 1:length(T)
    grad = abs(v[i] - v[j]) / (d[i] + d[j])
    if (grad > alpha)
      alpha = grad
      i2 = j
    end
  end

  vPlusArr = copy(v)
  vMinusArr = copy(v)
  for j in 1:length(T)
    vPlusArr[j] = v[j] + alpha * d[j]
    vMinusArr[j] = v[j] - alpha * d[j]
  end

  vLow = minimum(vPlusArr)
  vHigh = maximum(vMinusArr)

  T2 = falses(length(T))
  for j in 1:length(T)
    if (vMinusArr[j] > vLow * (1 + gradTol)
        || vPlusArr[j] < vHigh * (1 - gradTol))
      T2[j] = true
    end
  end
  if (findfirst(T2, true) == 0 || vLow >= vHigh)
    return (v[i] > v[i2])? (T[i],T[i2]) : (T[i2], T[i])
  else
    if (findfirst(T2, false) == 0)
      println(T, v, d)
      error("Infinite recursion!")
    end
    return StarSteepestPath(T[T2], v[T2], d[T2])
  end
end

# VertexSteepestPath. (Alg 9 in Lex Paper)
# does NOT work for directed graphs
function VertexSteepestPath(A, cpnts, isTerm, val, x)
  n = A.n
  d, parents = termFreeShortestPaths(A, x, isTerm)
  if (isTerm[x])
    maxGrad = -Inf
    y = x

    # at this point the graph may not be connected
    for i in 1:n
      if (isTerm[i] && cpnts[i] == cpnts[x])
        grad = abs(val[x]-val[i])/(d[i])
        if (grad > maxGrad)
          maxGrad = grad
          y = i
        end
      end
    end

    path = pathFromParents(parents, y)

    if (LEX_DEBUG)
      println(x, " ", y, " ", path)
    end

    if (val[x] > val[y])
      return reverse(path)
    else
      return path
    end
  else
    # filter out terminals in other components
    filter = copy(isTerm)
    for i in find(isTerm)
      if (d[i] == Inf)
        filter[i] = false
      end
    end
    if (LEX_DEBUG)
      println(x, "; filter = ", find(filter))
    end

    t1, t2 = StarSteepestPath(find(filter), val[filter], d[filter])
    p1 = pathFromParents(parents, t1)
    p2 = pathFromParents(parents, t2)
    if (LEX_DEBUG)
      println(t1, " ", t2, " ", p1', " ",  p2', " ", val[isTerm], " ", d[isTerm])
    end
    return vcat(p1[1:end-1], reverse(p2))
  end
end

# given a path from terminal to terminal, compute its gradient
function gradOfPath(A, isTerm, val, path)
  l = length(path)
  if (l < 2)
    return 0.0
  end
  @assert isTerm[path[1]] @sprintf("Node %d is not a terminal.", path[1])
  @assert isTerm[path[l]] @sprintf("Node %d is not a terminal.", path[l])

  dist = 0.0
  for i in 1:l-1
    # verify that it is actually a path
    @assert A[path[i], path[i + 1]] != 0.0
    @sprintf("Vertices %d and %d are not connected.", path[i], path[i + 1])

    dist += A[path[i], path[i + 1]]
  end

  dVolt = val[path[1]] - val[path[l]]
  return dVolt/dist
end # gradOfPath

function SteepestPath(A, cpnts, isTerm, val, vert=collect(1:n))

  n = A.n

  # do the following in each component, recurse at the same time

  # randomly pick an edge (x1, x2), and a vertex x3
  ind = rand(1:A.colptr[end]-1)
  x1 = findfirst(x -> x > ind, A.colptr) - 1
  x2 = A.rowval[ind]
  x3 = rand(vert)

  maxPath = VertexSteepestPath(A, cpnts, isTerm, val, x1)
  if (LEX_DEBUG)
    println(isTerm)
  end
  if (LEX_DEBUG)
    print("[SteepestPath] x1: ")
    print(x1)
    println(maxPath')
  end
  maxGrad = gradOfPath(A, isTerm, val, maxPath)

  p2 = VertexSteepestPath(A, cpnts, isTerm, val, x2)
  if (LEX_DEBUG)
    print("[SteepestPath] x2: ")
    print(x2)
    println(p2')

  end
  g2 = gradOfPath(A, isTerm, val, p2)
  if (g2 > maxGrad)
    maxPath = p2
    maxGrad = g2
  end

  p3 = VertexSteepestPath(A, cpnts, isTerm, val, x3)
  if (LEX_DEBUG)
    print("[SteepestPath] x3: ")
    print(x3)
    println(p3')
  end
  g3 = gradOfPath(A, isTerm, val, p3)

  if (g3 > maxGrad)
    maxPath = p3
    maxGrad = g3
  end

  hPresSet = CompHighPressGraph(A, isTerm, val, maxGrad)

  if (LEX_DEBUG)
    print("[SteepestPath] maxPath = ")
    println(maxPath')
    print("[SteepestPath] high pressure graph then becomes ")
    println(hPresSet')
  end

  if (length(hPresSet) == 0)
    return maxPath
  end

  # construct G' to recurse
  # TODO: make this cleaner and much faster, preferably without copying
  # by using a filter or something.
  t = copy(isTerm)
  t[hPresSet] = true
  newA = copy(A)
  for i = 1:n
    for j = 1:n
      if (!t[i] || !t[j])
        newA[i,j] = 0.0
      end
    end
  end

  newCpnts = components(newA)
  newVertices = []
  append!(newVertices, find(isTerm))
  append!(newVertices, hPresSet)
  return SteepestPath(newA, newCpnts, isTerm, val, newVertices)
end # SteepestPath

# remove (x,y) in E(G) such that both x and y are terminals;
# returns max_{x,y} grad[x,y]
function removeTerm2TermEdges!{Tv, Ti}(G::SparseMatrixCSC{Tv, Ti},
                                       isTerm::Array{Bool, 1},
                                       val::Array{Float64, 1}, )
  n = G.n
  terms = find(isTerm)
  alpha = -Inf

  for u in terms
    for v in terms
      if (G[u,v] > 0)
        grad = abs(val[u] - val[v]) / G[u, v]
        alpha = max(alpha, grad)
        # remove terminal to terminal edges
        G[u,v] = 0.0
        G[v,u] = 0.0
      end
    end
  end

  return alpha
end

# G may change (terminal to terminal edges will be removed)
# returns a new assignment
function CompInfMin{Tv, Ti}(G::SparseMatrixCSC{Tv, Ti},
                            isTerm::Array{Bool, 1},
                            val::Array{Float64, 1}, )
  alpha = removeTerm2TermEdges!(G, isTerm, val)

  cpnts = components(G)

  p = SteepestPath(G, cpnts, isTerm, val)
  grad = gradOfPath(G, isTerm, val, p)
  alpha = max(alpha, grad)

  vl, vlp = CompVLow(G, isTerm, val, alpha)
  vh, vhp = CompVHigh(G, isTerm, val, alpha)

  newVal = copy(val)
  for i in 1:G.n
    if (!isTerm[i])
      newVal[i] = vl[i] + (vh[i] - vl[i]) / 2.0
    end
  end

  return newVal
end # CompInfMin

function CompLexMin(A, isTerm, initVal)
  visited = copy(isTerm)

  B = sparse(copy(A))
  val = copy(initVal)

  # we keep going until every vertex has an assignment
  while (findfirst(visited, false) != 0)
    # remove terminal to terminal edges
    removeTerm2TermEdges!(B, visited, val)

    cpnts = components(B)
    p = SteepestPath(B, cpnts, visited, val)

    # visited will change here:
    fixPath!(B, visited, val, p)

    if (LEX_DEBUG)
      print("[CompLexMin] found steepest path: ")
      println(p')
      println(val')
    end
  end
  return val
end


text:
function exp_rand()
    # cdf of Possion+(log(2)), since P(X=0) = 1/2
    # for positive Possion, q[k] = sum(log(2)^k / k!)
    q = [
    0.6931471805599453,
    0.9333736875190459,
    0.9888777961838675,
    0.9984959252914960,
    0.9998292811061389,
    0.9999833164100727,
    0.9999985691438767,
    0.9999998906925558,
    0.9999999924734159,
    0.9999999995283275,
    0.9999999999728814,
    0.9999999999985598,
    0.9999999999999289,
    0.9999999999999968,
    0.9999999999999999,
    1.0000000000000000
    ]
    a = 0.0
    u = rand()
    while (u <= 0) | (u >= 1)
        u = rand()
    end
    while true
        u += u;
        if (u > 1.0)
            break
        end
        a += q[1]  # μ = q[1], and the times of loop is M, so a = μM
        # TODO: but why the times of loop follows a geometric (1/2) ??
    end
    u -= 1.0 # an independent uniform variate
    if u <= q[1]
        return a + u
    end
    i = 1 # one-based
    ustar = rand()
    umin = ustar
    while true
        ustar = rand()
        if (umin > ustar)
            umin = ustar
        end
        i += 1
        if u <= q[i]
            return a + umin * q[1]
        end
    end
end

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
module TempleGate

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Temple Gate (Theo, Holding Theo)" => Ahorn.EntityPlacement(
        Maple.TempleGate,
        "point",
        Dict{String, Any}(
            "sprite" => "theo",
            "type" => "HoldingTheo"
        )
    ),

    "Temple Gate (Default, Close Behind)" => Ahorn.EntityPlacement(
        Maple.TempleGate,
        "point",
        Dict{String, Any}(
            "sprite" => "default",
            "type" => "CloseBehindPlayer",
        )
    ),
    "Temple Gate (Mirror, Close Behind)" => Ahorn.EntityPlacement(
        Maple.TempleGate,
        "point",
        Dict{String, Any}(
            "sprite" => "mirror",
            "type" => "CloseBehindPlayer",
        )
    ),

    "Temple Gate (Default, Nearest Switch)" => Ahorn.EntityPlacement(
        Maple.TempleGate,
        "point",
        Dict{String, Any}(
            "sprite" => "default",
            "type" => "NearestSwitch",
        )
    ),
    "Temple Gate (Mirror, Nearest Switch)" => Ahorn.EntityPlacement(
        Maple.TempleGate,
        "point",
        Dict{String, Any}(
            "sprite" => "mirror",
            "type" => "NearestSwitch",
        )
    ),

    "Temple Gate (Default, Touch Switches)" => Ahorn.EntityPlacement(
        Maple.TempleGate,
        "point",
        Dict{String, Any}(
            "sprite" => "default",
            "type" => "TouchSwitches",
        )
    ),
    "Temple Gate (Mirror, Touch Switches)" => Ahorn.EntityPlacement(
        Maple.TempleGate,
        "point",
        Dict{String, Any}(
            "sprite" => "mirror",
            "type" => "TouchSwitches",
        )
    ),
)

textures = String["default", "mirror", "theo"]
modes = Maple.temple_gate_modes

Ahorn.editingOptions(entity::Maple.TempleGate) = Dict{String, Any}(
    "type" => modes,
    "sprite" => textures
)

function Ahorn.selection(entity::Maple.TempleGate)
    x, y = Ahorn.position(entity)
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x - 4, y, 15, height)
end

sprites = Dict{String, String}(
    "default" => "objects/door/TempleDoor00",
    "mirror" => "objects/door/TempleDoorB00",
    "theo" => "objects/door/TempleDoorC00"
)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.TempleGate, room::Maple.Room)
    sprite = get(entity.data, "sprite", "default")

    if haskey(sprites, sprite)
        Ahorn.drawImage(ctx, sprites[sprite], -4, 0)
    end
end

end


text:
using Test

using BraidChains
using BraidChains: SynchronicBraid, Hash, MetaData

import PeaceCypher: sign, hash
using PeaceCypher

using SynchronicBallot

include("basics.jl")

include("services.jl")

include("keychain.jl")


text:
function sieve_of_atkin(limit)
    primes = Int64[]

    sieve = falses(limit)

    square_root_of_limit = round(Int64, sqrt(limit))

    for x in 1:square_root_of_limit
        for y in 1:square_root_of_limit
            n = (4 * x * x) + (y * y)

            if n <= limit && (n % 12 == 1 || n % 12 == 5)
                sieve[n] = !sieve[n]
            end

            n = (3 * x * x) + (y * y)

            if n <= limit && n % 12 == 7
                sieve[n] = !sieve[n]
            end

            n = (3 * x * x) - (y * y)

            if x > y && n <= limit && n % 12 == 11
                sieve[n] = !sieve[n]
            end
        end
    end

    sieve[2] = true

    sieve[3] = true

    push!(primes, 2)

    push!(primes, 3)

    n = 5

    while n <= square_root_of_limit
        if sieve[n]
            i = n * n

            while i < limit
                sieve[i] = false

                i += n * n
            end

            push!(primes, n)
        end

        n += 2
    end

    while n < limit
        if sieve[n]
            push!(primes, n)
        end

        n += 2
    end

    primes
end

n = 10001

limit = round(Int64, (n * log(n)) + (n * (log(log(n)))))

println(sieve_of_atkin(limit)[n])

max_stars_repo_name: seaneshbaugh/rosetta-euler

max_issues_repo_name: seaneshbaugh/rosetta-euler

max_forks_repo_name: seaneshbaugh/rosetta-euler


text:
module io

export IORestoreSettings
export IOSetParam
export IOSetParamMat
export IOSetParamUMat
export IOSetParamRow
export IOSetParamCol
export IOSetParamURow
export IOSetParamUCol
export IOGetParamBool
export IOGetParamInt
export IOGetParamDouble
export IOGetParamString
export IOGetParamVectorStr
export IOGetParamVectorInt
export IOGetParamMat
export IOGetParamUMat
export IOGetParamCol
export IOGetParamRow
export IOGetParamUCol
export IOGetParamURow
export IOGetParamMatWithInfo
export IOEnableVerbose
export IODisableVerbose
export IOSetPassed

import mlpack_jll
const library = mlpack_jll.libmlpack_julia_util

# Utility function to convert 1d object to 2d.
function convert_to_2d(in::Array{T, 1})::Array{T, 2} where T
  reshape(in, length(in), 1)
end

# Utility function to convert 2d object to 1d.  Fails if the size of one
# dimension is not 1.
function convert_to_1d(in::Array{T, 2})::Array{T, 1} where T
  if size(in, 1) != 1 && size(in, 2) != 1
    throw(ArgumentError("given matrix must be 1-dimensional; but its size is " *
        "$(size(in))"))
  end

  vec(in)
end

# Utility function to convert to and return a matrix.
function to_matrix(input, T::Type)
  if isa(input, Array{T, 1})
    convert_to_2d(input)
  else
    convert(Array{T, 2}, input)
  end
end

# Utility function to convert to and return a vector.
function to_vector(input, T::Type)
  if isa(input, Array{T, 1})
    input
  else
    convert_to_1d(convert(Array{T, 2}, input))
  end
end

function IORestoreSettings(programName::String)
  ccall((:IO_RestoreSettings, library), Nothing, (Cstring,), programName);
end

function IOSetParam(paramName::String, paramValue::Int)
  ccall((:IO_SetParamInt, library), Nothing, (Cstring, Cint), paramName,
      Cint(paramValue));
end

function IOSetParam(paramName::String, paramValue::Float64)
  ccall((:IO_SetParamDouble, library), Nothing, (Cstring, Float64), paramName,
      paramValue);
end

function IOSetParam(paramName::String, paramValue::Bool)
  ccall((:IO_SetParamBool, library), Nothing, (Cstring, Bool), paramName,
      paramValue);
end

function IOSetParam(paramName::String, paramValue::String)
  ccall((:IO_SetParamString, library), Nothing, (Cstring, Cstring), paramName,
      paramValue);
end

function IOSetParamMat(paramName::String,
                        paramValue,
                        pointsAsRows::Bool)
  paramMat = to_matrix(paramValue, Float64)
  ccall((:IO_SetParamMat, library), Nothing, (Cstring, Ptr{Float64}, Csize_t,
      Csize_t, Bool), paramName, Base.pointer(paramMat), size(paramMat, 1),
      size(paramMat, 2), pointsAsRows);
end

function IOSetParamUMat(paramName::String,
                         paramValue,
                         pointsAsRows::Bool)
  paramMat = to_matrix(paramValue, Int)

  # Sanity check.
  if minimum(paramMat) <= 0
    throw(DomainError("Input $(paramName) cannot have 0 or negative values!  " *
        "Must be 1 or greater."))
  end

  m = convert(Array{Csize_t, 2}, paramMat .- 1)
  ccall((:IO_SetParamUMat, library), Nothing, (Cstring, Ptr{Csize_t}, Csize_t,
      Csize_t, Bool), paramName, Base.pointer(m), size(paramValue, 1),
      size(paramValue, 2), pointsAsRows);
end

function IOSetParam(paramName::String,
                     vector::Vector{String})
  # For this we have to set the size of the vector then each string
  # sequentially.  I am not sure if this is fully necessary but I have some
  # reservations about Julia's support for passing arrays of strings correctly
  # as a const char**.
  ccall((:IO_SetParamVectorStrLen, library), Nothing, (Cstring, Csize_t),
      paramName, size(vector, 1));
  for i in 1:size(vector, 1)
    ccall((:IO_SetParamVectorStrStr, library), Nothing, (Cstring, Cstring,
        Csize_t), paramName, vector[i], i .- 1);
  end
end

function IOSetParam(paramName::String,
                     vector::Vector{Int})
  cint_vec = convert(Vector{Cint}, vector)
  ccall((:IO_SetParamVectorInt, library), Nothing, (Cstring, Ptr{Cint},
      Csize_t), paramName, Base.pointer(cint_vec), size(cint_vec, 1));
end

function IOSetParam(paramName::String,
                     matWithInfo::Tuple{Array{Bool, 1}, Array{Float64, 2}},
                     pointsAsRows::Bool)
  ccall((:IO_SetParamMatWithInfo, library), Nothing, (Cstring, Ptr{Bool},
      Ptr{Float64}, Int, Int, Bool), paramName,
      Base.pointer(matWithInfo[1]), Base.pointer(matWithInfo[2]),
      size(matWithInfo[2], 1), size(matWithInfo[2], 2), pointsAsRows);
end

function IOSetParamRow(paramName::String,
                        paramValue)
  paramVec = to_vector(paramValue, Float64)
  ccall((:IO_SetParamRow, library), Nothing, (Cstring, Ptr{Float64}, Csize_t),
      paramName, Base.pointer(paramVec), size(paramVec, 1));
end

function IOSetParamCol(paramName::String,
                        paramValue)
  paramVec = to_vector(paramValue, Float64)
  ccall((:IO_SetParamCol, library), Nothing, (Cstring, Ptr{Float64}, Csize_t),
      paramName, Base.pointer(paramVec), size(paramVec, 1));
end

function IOSetParamURow(paramName::String,
                         paramValue)
  paramVec = to_vector(paramValue, Int)

  # Sanity check.
  if minimum(paramVec) <= 0
    throw(DomainError("Input $(paramName) cannot have 0 or negative values!  " *
        "Must be 1 or greater."))
  end
  m = convert(Array{Csize_t, 1}, paramVec .- 1)

  ccall((:IO_SetParamURow, library), Nothing, (Cstring, Ptr{Csize_t}, Csize_t),
      paramName, Base.pointer(m), size(paramValue, 1));
end

function IOSetParamUCol(paramName::String,
                         paramValue)
  paramVec = to_vector(paramValue, Int)

  # Sanity check.
  if minimum(paramVec) <= 0
    throw(DomainError("Input $(paramName) cannot have 0 or negative values!  " *
        "Must be 1 or greater."))
  end
  m = convert(Array{Csize_t, 1}, paramValue .- 1)

  ccall((:IO_SetParamUCol, library), Nothing, (Cstring, Ptr{Csize_t}, Csize_t),
      paramName, Base.pointer(m), size(paramValue, 1));
end

function IOGetParamBool(paramName::String)
  return ccall((:IO_GetParamBool, library), Bool, (Cstring,), paramName)
end

function IOGetParamInt(paramName::String)
  return Int(ccall((:IO_GetParamInt, library), Cint, (Cstring,), paramName))
end

function IOGetParamDouble(paramName::String)
  return ccall((:IO_GetParamDouble, library), Float64, (Cstring,), paramName)
end

function IOGetParamString(paramName::String)
  return ccall((:IO_GetParamString, library), Cstring, (Cstring,), paramName)
end

function IOGetParamVectorStr(paramName::String)
  local size::Csize_t
  local ptr::Ptr{String}

  # Get the size of the vector, then each element.
  size = ccall((:IO_GetParamVectorStrLen, library), Csize_t, (Cstring,),
      paramName);
  out = Array{String, 1}()
  for i = 1:size
    s = ccall((:IO_GetParamVectorStrStr, library), Cstring, (Cstring, Csize_t),
        paramName, i .- 1)
    push!(out, Base.unsafe_string(s))
  end

  return out
end

function IOGetParamVectorInt(paramName::String)
  local size::Csize_t
  local ptr::Ptr{Cint}

  # Get the size of the vector, then the pointer to it.  We will own the
  # pointer.
  size = ccall((:IO_GetParamVectorIntLen, library), Csize_t, (Cstring,),
      paramName);
  ptr = ccall((:IO_GetParamVectorIntPtr, library), Ptr{Cint}, (Cstring,),
      paramName);

  return convert(Array{Int, 1}, Base.unsafe_wrap(Array{Cint, 1}, ptr, (size),
      own=true))
end

function IOGetParamMat(paramName::String, pointsAsRows::Bool)
  # Can we return different return types?  For now let's restrict to a matrix to
  # make it easy...
  local ptr::Ptr{Float64}
  local rows::Csize_t, cols::Csize_t;
  # I suppose it would be possible to do this all in one call, but this seems
  # easy enough.
  rows = ccall((:IO_GetParamMatRows, library), Csize_t, (Cstring,), paramName);
  cols = ccall((:IO_GetParamMatCols, library), Csize_t, (Cstring,), paramName);
  ptr = ccall((:IO_GetParamMat, library), Ptr{Float64}, (Cstring,), paramName);

  if pointsAsRows
    # In this case we have to transpose, unfortunately.
    m = Base.unsafe_wrap(Array{Float64, 2}, ptr, (rows, cols), own=true)
    return m';
  else
    # Here no transpose is necessary.
    return Base.unsafe_wrap(Array{Float64, 2}, ptr, (rows, cols), own=true);
  end
end

function IOGetParamUMat(paramName::String, pointsAsRows::Bool)
  # Can we return different return types?  For now let's restrict to a matrix to
  # make it easy...
  local ptr::Ptr{Csize_t}
  local rows::Csize_t, cols::Csize_t;
  # I suppose it would be possible to do this all in one call, but this seems
  # easy enough.
  rows = ccall((:IO_GetParamUMatRows, library), Csize_t, (Cstring,), paramName);
  cols = ccall((:IO_GetParamUMatCols, library), Csize_t, (Cstring,), paramName);
  ptr = ccall((:IO_GetParamUMat, library), Ptr{Csize_t}, (Cstring,), paramName);

  if pointsAsRows
    # In this case we have to transpose, unfortunately.
    m = Base.unsafe_wrap(Array{Csize_t, 2}, ptr, (rows, cols), own=true);
    return convert(Array{Int, 2}, m' .+ 1)  # Add 1 because these are indexes.
  else
    # Here no transpose is necessary.
    m = Base.unsafe_wrap(Array{Csize_t, 2}, ptr, (rows, cols), own=true);
    return convert(Array{Int, 2}, m .+ 1)
  end
end

function IOGetParamCol(paramName::String)
  local ptr::Ptr{Float64};
  local rows::Csize_t;

  rows = ccall((:IO_GetParamColRows, library), Csize_t, (Cstring,), paramName);
  ptr = ccall((:IO_GetParamCol, library), Ptr{Float64}, (Cstring,), paramName);

  return Base.unsafe_wrap(Array{Float64, 1}, ptr, rows, own=true);
end

function IOGetParamRow(paramName::String)
  local ptr::Ptr{Float64};
  local cols::Csize_t;

  cols = ccall((:IO_GetParamRowCols, library), Csize_t, (Cstring,), paramName);
  ptr = ccall((:IO_GetParamRow, library), Ptr{Float64}, (Cstring,), paramName);

  return Base.unsafe_wrap(Array{Float64, 1}, ptr, cols, own=true);
end

function IOGetParamUCol(paramName::String)
  local ptr::Ptr{Csize_t};
  local rows::Csize_t;

  rows = ccall((:IO_GetParamUColRows, library), Csize_t, (Cstring,), paramName);
  ptr = ccall((:IO_GetParamUCol, library), Ptr{Csize_t}, (Cstring,), paramName);

  m = Base.unsafe_wrap(Array{Csize_t, 1}, ptr, rows, own=true);
  return convert(Array{Int, 1}, m .+ 1)
end

function IOGetParamURow(paramName::String)
  local ptr::Ptr{Csize_t};
  local cols::Csize_t;

  cols = ccall((:IO_GetParamURowCols, library), Csize_t, (Cstring,), paramName);
  ptr = ccall((:IO_GetParamURow, library), Ptr{Csize_t}, (Cstring,), paramName);

  m = Base.unsafe_wrap(Array{Csize_t, 1}, ptr, cols, own=true);
  return convert(Array{Int, 1}, m .+ 1)
end

function IOGetParamMatWithInfo(paramName::String, pointsAsRows::Bool)
  local ptrBool::Ptr{Bool};
  local ptrData::Ptr{Float64};
  local rows::Csize_t;
  local cols::Csize_t;

  rows = ccall((:IO_GetParamMatWithInfoRows, library), Csize_t, (Cstring,),
      paramName);
  cols = ccall((:IO_GetParamMatWithInfoCols, library), Csize_t, (Cstring,),
      paramName);
  ptrBool = ccall((:IO_GetParamMatWithInfoBoolPtr, library), Ptr{Bool},
      (Cstring,), paramName);
  ptrMem = ccall((:IO_GetParamMatWithInfoPtr, library), Ptr{Float64},
      (Cstring,), paramName);

  types = Base.unsafe_wrap(Array{Bool, 1}, ptrBool, (rows), own=true)
  if pointsAsRows
    # In this case we have to transpose, unfortunately.
    m = Base.unsafe_wrap(Array{Float64, 2}, ptr, (rows, cols), own=true)
    return (types, m');
  else
    # Here no transpose is necessary.
    return (types, Base.unsafe_wrap(Array{Float64, 2}, ptr, (rows, cols),
        own=true));
  end
end

function IOEnableVerbose()
  ccall((:IO_EnableVerbose, library), Nothing, ());
end

function IODisableVerbose()
  ccall((:IO_DisableVerbose, library), Nothing, ());
end

function IOSetPassed(paramName::String)
  ccall((:IO_SetPassed, library), Nothing, (Cstring,), paramName);
end

end # module io


text:
module Day9

export day9

function issumoftwo(window, number)
    for a in window
        if in(number-a, window)
            return true
        end
    end
    false
end

function day9(::Val{:a}, filename, prelen)
    input = Vector{Int}()
    for line in readlines(filename)
        push!(input, parse(Int, line))
    end
    for i in 1:length(input)-prelen
        window = input[i:prelen+i-1]
        new = input[prelen+i]
        if !issumoftwo(window, new)
            return new
        end
    end
end

function findwindow(input, target)
    for i in 1:length(input)-1
        current_sum = input[i]
        for j in i+1:length(input)
            current_sum += input[j]
            if current_sum == target
                return input[i:j]
            end
        end
    end
end

function day9(::Val{:b}, filename, prelen)
    input = Vector{Int}()
    for line in readlines(filename)
        push!(input, parse(Int, line))
    end
    target = day9a(filename, prelen)
    window = findwindow(input, target)
    (min, max) = reduce(window, init = (window[1],window[1])) do (min, max), n
        (n < min ? n : min, n > max ? n : max)
    end
    min + max
end

end


text:
abstract type Position end
abstract type Environment end

struct DefaultEnvironment <: Environment end
struct DefaultPosition <: Position end

struct Ancestor end

mutable struct Organism{P <: Position, E <: Environment, G}
    id::UInt64

    genotype::G
    parent_ids::Vector{UInt64}
    parent_genotypes::Vector{G}

    position::P
    environment::E

    time_birth::UInt64
    time_death::UInt64

    snapshot_id::String

    function Organism(genotype, parent_ids, parent_genotypes, time_birth, time_death) where
        {O <: Organism, A <: Ancestor, G}

        new{DefaultPosition, DefaultEnvironment, G}(-1, genotype, parent_ids, parent_genotypes, DefaultPosition(), DefaultEnvironment(), time_birth, time_death, 0)
    end
    function Organism{P, E, G}(id, genotype, parent_ids, parent_genotypes, position::P, environment::E, time_birth, time_death, snapshot_id) where
        {P <: Position, E <: Environment, G}

        new{P, E, G}(id, genotype, parent_ids, parent_genotypes, position, environment, time_birth, time_death, string(snapshot_id))
    end
end

get_id(::Ancestor) = -1
get_id(organism::Organism) = organism.id

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
module Omniscape
using Circuitscape
using DelimitedFiles
using Distributed

include("config.jl")
include("utils.jl")
include("functions.jl")
include("io.jl")
include("run_omniscape.jl")
export run_omniscape

end


text:
@testset "3_levels" begin
  @testset "EncryptionParameters" begin
    @test_nowarn EncryptionParameters(SchemeType.bfv)
  end
  enc_parms = EncryptionParameters(SchemeType.bfv)

  @testset "polynomial modulus degree" begin
    @test_nowarn set_poly_modulus_degree!(enc_parms, 8192)
  end

  @testset "coefficient modulus" begin
    @test_nowarn set_coeff_modulus!(enc_parms, coeff_modulus_create(8192, [50, 30, 30, 50, 50]))
  end

  @testset "plain modulus" begin
    @test_nowarn set_plain_modulus!(enc_parms, plain_modulus_batching(8192, 20))
  end

  @testset "SEALContext" begin
    @test_nowarn SEALContext(enc_parms)
  end
  context = SEALContext(enc_parms)

  @testset "key_context_data" begin
    @test_nowarn key_context_data(context)
  end
  context_data = key_context_data(context)

  @testset "modulus switching chain (key context data)" begin
    @test chain_index(context_data) == 4
    @test parms_id(context_data) == [0x26d0ad92b6a78b12,
                                     0x667d7d6411d19434,
                                     0x18ade70427566279,
                                     0x84e0aa06442af302]
    @test_nowarn coeff_modulus(parms(context_data))
    primes = coeff_modulus(parms(context_data))
    @test value(primes[1]) == 0x3ffffffef4001
    @test value(primes[2]) == 0x3ffe8001
    @test value(primes[3]) == 0x3fff4001
    @test value(primes[4]) == 0x3fffffffcc001
    @test value(primes[5]) == 0x3ffffffffc001
  end

  @testset "first_context_data" begin
    @test_nowarn first_context_data(context)
  end
  context_data = first_context_data(context)

  @testset "modulus switching chain (first context data)" begin
    @test chain_index(context_data) == 3
    @test parms_id(context_data) == first_parms_id(context)
    @test parms_id(context_data) == [0x211ee2c43ec16b18,
                                     0x2c176ee3b851d741,
                                     0x490eacf1dd5930b3,
                                     0x3212f104b7a60a0c]
    @test_nowarn coeff_modulus(parms(context_data))
    primes = coeff_modulus(parms(context_data))
    @test value(primes[1]) == 0x3ffffffef4001
    @test value(primes[2]) == 0x3ffe8001
    @test value(primes[3]) == 0x3fff4001
    @test value(primes[4]) == 0x3fffffffcc001
  end

  @testset "next_context_data" begin
    @test_nowarn next_context_data(context_data)
  end
  context_data = next_context_data(context_data)
  context_data = next_context_data(context_data)
  context_data = next_context_data(context_data)
  @testset "isnothing(next_context_data)" begin
    @test isnothing(next_context_data(context_data))
  end

  @testset "modulus switching chain (last context data)" begin
    @test chain_index(context_data) == 0
    @test parms_id(context_data) == last_parms_id(context)
    @test parms_id(context_data) == [0xaf7f6dac55528cf7,
                                     0x2f532a7e2362ab73,
                                     0x03aeaedd1059515e,
                                     0xa515111177a581ca]
    @test_nowarn coeff_modulus(parms(context_data))
    primes = coeff_modulus(parms(context_data))
    @test value(primes[1]) == 0x3ffffffef4001
  end

  keygen = KeyGenerator(context)
  public_key_ = PublicKey()
  create_public_key!(public_key_, keygen)
  secret_key_ = secret_key(keygen)
  relin_keys_ = RelinKeys()
  create_relin_keys!(relin_keys_, keygen)
  galois_keys_ = GaloisKeys()
  create_galois_keys!(galois_keys_, keygen)

  @testset "parms_id of generated keys" begin
    p = [0x26d0ad92b6a78b12, 0x667d7d6411d19434, 0x18ade70427566279, 0x84e0aa06442af302]
    @test parms_id(public_key_) == p
    @test parms_id(secret_key_) == p
    @test parms_id(relin_keys_) == p
    @test parms_id(galois_keys_) == p
  end

  encryptor = Encryptor(context, public_key_)
  evaluator = Evaluator(context)
  decryptor = Decryptor(context, secret_key_)

  @testset "Plaintext" begin
    @test_nowarn Plaintext("1x^3 + 2x^2 + 3x^1 + 4")
  end
  plain = Plaintext("1x^3 + 2x^2 + 3x^1 + 4")

  encrypted = Ciphertext()
  encrypt!(encrypted, plain, encryptor)

  @testset "parms_id of plain" begin
    @test parms_id(plain) == [0x0000000000000000,
                              0x0000000000000000,
                              0x0000000000000000,
                              0x0000000000000000]
  end

  @testset "parms_id of encrypted" begin
    @test parms_id(encrypted) == [0x211ee2c43ec16b18,
                                  0x2c176ee3b851d741,
                                  0x490eacf1dd5930b3,
                                  0x3212f104b7a60a0c]
  end

  @testset "modulus switching on encrypted (level 3)" begin
    context_data = first_context_data(context)
    @test chain_index(context_data) == 3
    @test parms_id(context_data) == [0x211ee2c43ec16b18,
                                     0x2c176ee3b851d741,
                                     0x490eacf1dd5930b3,
                                     0x3212f104b7a60a0c]
    @test invariant_noise_budget(encrypted, decryptor) in (131, 132, 133)
  end

  @testset "modulus switching on encrypted (level 2)" begin
    @test mod_switch_to_next_inplace!(encrypted, evaluator) == encrypted
    context_data = next_context_data(context_data)
    @test chain_index(context_data) == 2
    @test parms_id(context_data) == [0x85626ad91458073f,
                                     0xe186437698f5ff4e,
                                     0xa1e71da26dabe039,
                                     0x9b66f4ab523b9be1]
    @test invariant_noise_budget(encrypted, decryptor) in (81, 82, 83)
  end

  @testset "modulus switching on encrypted (level 1)" begin
    @test mod_switch_to_next_inplace!(encrypted, evaluator) == encrypted
    context_data = next_context_data(context_data)
    @test chain_index(context_data) == 1
    @test parms_id(context_data) == [0x73b7dc26d10a15b9,
                                     0x56ce8bdd07324dfa,
                                     0x7ff7b8ec16a6f20f,
                                     0xb80f7319f2a28ac1]
    @test invariant_noise_budget(encrypted, decryptor) in (51, 52, 53)
  end

  @testset "modulus switching on encrypted (level 3)" begin
    @test mod_switch_to_next_inplace!(encrypted, evaluator) == encrypted
    context_data = next_context_data(context_data)
    @test chain_index(context_data) == 0
    @test parms_id(context_data) == [0xaf7f6dac55528cf7,
                                     0x2f532a7e2362ab73,
                                     0x03aeaedd1059515e,
                                     0xa515111177a581ca]
    @test invariant_noise_budget(encrypted, decryptor) in (21, 22, 23)
  end

  @testset "decrypt! and check 1" begin
    @test_nowarn decrypt!(plain, encrypted, decryptor)
    @test to_string(plain) == "1x^3 + 2x^2 + 3x^1 + 4"
  end

  @testset "compute with modswitching" begin
    @test_nowarn encrypt!(encrypted, plain, encryptor)
    @test invariant_noise_budget(encrypted, decryptor) in (131, 132, 133)

    @test square_inplace!(encrypted, evaluator) == encrypted
    @test relinearize_inplace!(encrypted, relin_keys_, evaluator) == encrypted
    @test invariant_noise_budget(encrypted, decryptor) in (99, 100, 101)

    @test square_inplace!(encrypted, evaluator) == encrypted
    @test relinearize_inplace!(encrypted, relin_keys_, evaluator) == encrypted
    @test invariant_noise_budget(encrypted, decryptor) in (66, 67, 68)

    @test_nowarn mod_switch_to_next_inplace!(encrypted, evaluator)
    @test invariant_noise_budget(encrypted, decryptor) in (66, 67, 68)

    @test square_inplace!(encrypted, evaluator) == encrypted
    @test relinearize_inplace!(encrypted, relin_keys_, evaluator) == encrypted
    @test invariant_noise_budget(encrypted, decryptor) in (33, 34, 35)

    @test_nowarn mod_switch_to_next_inplace!(encrypted, evaluator)
    @test invariant_noise_budget(encrypted, decryptor) in (33, 34, 35)
  end

  @testset "decrypt! and check 2" begin
    @test_nowarn decrypt!(plain, encrypted, decryptor)
    @test to_string(plain) == "1x^24 + 10x^23 + 88x^22 + 330x^21 + EFCx^20 + 3A30x^19 + C0B8x^18 + 22BB0x^17 + 58666x^16 + C88D0x^15 + 9C377x^14 + F4C0Ex^13 + E8B38x^12 + 5EE89x^11 + F8BFFx^10 + 30304x^9 + 5B9D4x^8 + 12653x^7 + 4DFB5x^6 + 879F8x^5 + 825FBx^4 + F1FFEx^3 + 3FFFFx^2 + 60000x^1 + 10000"
  end

  @testset "context without expanded modulus chain" begin
    @test_nowarn SEALContext(enc_parms, expand_mod_chain=false)
  end
  context = SEALContext(enc_parms, expand_mod_chain=false)

  context_data = key_context_data(context)
  @testset "modulus switching chain (chain index 1)" begin
    @test chain_index(context_data) == 1
    @test parms_id(context_data) == [0x26d0ad92b6a78b12,
                                     0x667d7d6411d19434,
                                     0x18ade70427566279,
                                     0x84e0aa06442af302]
    @test_nowarn coeff_modulus(parms(context_data))
    primes = coeff_modulus(parms(context_data))
    @test value(primes[1]) == 0x3ffffffef4001
    @test value(primes[2]) == 0x3ffe8001
    @test value(primes[3]) == 0x3fff4001
    @test value(primes[4]) == 0x3fffffffcc001
  end

  context_data = next_context_data(context_data)
  @testset "modulus switching chain (chain index 0)" begin
    @test chain_index(context_data) == 0
    @test parms_id(context_data) == [0x211ee2c43ec16b18,
                                     0x2c176ee3b851d741,
                                     0x490eacf1dd5930b3,
                                     0x3212f104b7a60a0c]
    @test_nowarn coeff_modulus(parms(context_data))
    primes = coeff_modulus(parms(context_data))
    @test value(primes[1]) == 0x3ffffffef4001
    @test value(primes[2]) == 0x3ffe8001
    @test value(primes[3]) == 0x3fff4001
  end

end


text:
module MyMethods

# Dependencies =================================================================
using LinearAlgebra
using Random
using Distributions
using DataFrames

# Module export ================================================================

export myEstimator, myLS, myTSLS, mySieve
export coef, predict, inference # my methods
export get_basis # helper functions

# Module content ===============================================================

include("FUN/myEstimator.jl") # myEstimators

include("FUN/myLS.jl") # myLS

include("FUN/myTSLS.jl") # myLS

include("FUN/mySieve.jl") # mySieve

end


text:
export simplify_code, optimize_code, GreedyMethod, KaHyParBipartite, SABipartite, TreeSA, MergeGreedy, MergeVectors, uniformsize

abstract type CodeSimplifier end
Base.@kwdef struct MergeGreedy <: CodeSimplifier
    threshhold::Float64=-1e-12
end
struct MergeVectors <: CodeSimplifier end

"""
    optimize_code(eincode, size_dict, optimizer = GreedyMethod(), simplifier=nothing, permute=true)

Optimize the einsum contraction code and reduce the time/space complexity of tensor network contraction.
Returns a `NestedEinsum` instance. Input arguments are

* `eincode` is an einsum contraction code instance, one of `DynamicEinCode`, `StaticEinCode` or `NestedEinsum`.
* `size` is a dictionary of "edge label=>edge size" that contains the size information, one can use `uniformsize(eincode, 2)` to create a uniform size.
* `optimizer` is a `CodeOptimizer` instance, should be one of `GreedyMethod`, `KaHyParBipartite`, `SABipartite` or `TreeSA`. Check their docstrings for details.
* `simplifier` is one of `MergeVectors` or `MergeGreedy`.
* optimize the permutation if `permute` is true.
"""
function optimize_code(code::Union{EinCode, NestedEinsum}, size_dict::Dict, optimizer::CodeOptimizer, simplifier=nothing, permute::Bool=true)
    if simplifier === nothing
        optcode = _optimize_code(code, size_dict, optimizer)
    else
        simpl, code = simplify_code(code, size_dict, simplifier)
        optcode0 = _optimize_code(code, size_dict, optimizer)
        optcode = embed_simplifier(optcode0, simpl)
    end
    if permute
        optimize_permute(optcode, 0)
    end
end

simplify_code(code::Union{EinCode, NestedEinsum}, size_dict, ::MergeVectors) = merge_vectors(code)
simplify_code(code::Union{EinCode, NestedEinsum}, size_dict, method::MergeGreedy) = merge_greedy(code, size_dict; threshhold=method.threshhold)

function _optimize_code(code, size_dict, optimizer::KaHyParBipartite)
    recursive_bipartite_optimize(optimizer, code, size_dict)
end
function _optimize_code(code, size_dict, optimizer::GreedyMethod)
    optimize_greedy(code, size_dict; method=optimizer.method, nrepeat=optimizer.nrepeat)
end
function _optimize_code(code, size_dict, optimizer::SABipartite)
    recursive_bipartite_optimize(optimizer, code, size_dict)
end
function _optimize_code(code, size_dict, optimizer::TreeSA)
    optimize_tree(code, size_dict; sc_target=optimizer.sc_target, βs=optimizer.βs,
        ntrials=optimizer.ntrials, niters=optimizer.niters, nslices=optimizer.nslices,
        sc_weight=optimizer.sc_weight, rw_weight=optimizer.rw_weight, initializer=optimizer.initializer,
        greedy_method=optimizer.greedy_config.method, greedy_nrepeat=optimizer.greedy_config.nrepeat)
end

uniformsize(code::AbstractEinsum, size) = Dict([l=>size for l in uniquelabels(code)])


text:
##### Mixing length models

abstract type MixingLengthModel end

#####
##### Types
#####

"""
    ConstantMixingLength

TODO: add docs
"""
struct ConstantMixingLength{FT} <: MixingLengthModel
    value::FT
end

"""
    SCAMPyMixingLength

TODO: add docs
"""
struct SCAMPyMixingLength{FT} <: MixingLengthModel
    a_L::StabilityDependentParam{FT}
    b_L::StabilityDependentParam{FT}
end

"""
    IgnaciosMixingLength

TODO: add docs
"""
struct IgnaciosMixingLength{FT} <: MixingLengthModel
    a_L::StabilityDependentParam{FT}
    b_L::StabilityDependentParam{FT}
    c_K::FT
    c_ε::FT
    c_w::FT
    ω_1::FT
    ω_2::FT
    Prandtl_neutral::FT
    function IgnaciosMixingLength(
        a_L::StabilityDependentParam{FT},
        b_L::StabilityDependentParam{FT},
        c_K::FT,
        c_ε::FT,
        c_w::FT,
        ω_1::FT,
        Prandtl_neutral::FT,
    ) where {FT}
        return new{FT}(a_L, b_L, c_K, c_ε, c_w, ω_1, ω_1 + 1, Prandtl_neutral)
    end
end

struct MinimumDissipation{FT} <: MixingLengthModel
    a_L::StabilityDependentParam{FT}
    b_L::StabilityDependentParam{FT}
    c_K::FT
    c_ε::FT
    c_w::FT
    ω_1::FT
    ω_2::FT
    Prandtl_neutral::FT
    function MinimumDissipation(
        a_L::StabilityDependentParam{FT},
        b_L::StabilityDependentParam{FT},
        c_K::FT,
        c_ε::FT,
        c_w::FT,
        ω_1::FT,
        Prandtl_neutral::FT,
    ) where {FT}
        return new{FT}(a_L, b_L, c_K, c_ε, c_w, ω_1, ω_1 + 1, Prandtl_neutral)
    end
end

######
###### Helper funcs
######

"""
    compute_mixing_τ

TODO: add docs
"""
compute_mixing_τ(zi::FT, wstar::FT) where {FT} = zi / (wstar + FT(0.001))

"""
    ϕ_m

TODO: add docs
"""
ϕ_m(ξ, a_L, b_L) = (1 + a_L * ξ)^(-b_L)

######
###### Compute mixing length methods
######

"""
    compute_mixing_length!

Define mixing length

 - `aux[:l_mix, k, i]`

for all `k` and all `i`
"""
function compute_mixing_length! end

function compute_mixing_length!(
    grid::Grid{FT},
    q,
    aux,
    params,
    model::ConstantMixingLength,
) where {FT}
    gm, en, ud, sd, al = allcombinations(q)
    @inbounds for k in over_elems(grid)
        aux[:l_mix, k, gm] = model.value
    end
end

function compute_mixing_length!(
    grid::Grid{FT},
    q,
    aux,
    params,
    model::SCAMPyMixingLength{FT},
) where {FT}
    @unpack obukhov_length, zi, wstar, k_Karman = params
    gm, en, ud, sd, al = allcombinations(q)
    τ = compute_mixing_τ(zi, wstar)
    a_L = model.a_L(obukhov_length)
    b_L = model.b_L(obukhov_length)
    @inbounds for k in over_elems_real(grid)
        l1 = τ * sqrt(max(q[:tke, k, en], FT(0)))
        z = grid.zc[k]
        ξ = z / obukhov_length
        l2 = k_Karman * z * ϕ_m(ξ, a_L, b_L)
        aux[:l_mix, k, gm] = max(1 / (1 / max(l1, 1e-10) + 1 / l2), FT(1e-3))
    end
end

function compute_mixing_length!(
    grid::Grid{FT},
    q,
    aux,
    params,
    model::IgnaciosMixingLength{FT},
) where {FT}
    @unpack obukhov_length, SurfaceModel, k_Karman, param_set = params
    gm, en, ud, sd, al = allcombinations(q)
    ustar = SurfaceModel.ustar
    L = Vector(undef, 3)
    a_L = model.a_L(obukhov_length)
    b_L = model.b_L(obukhov_length)
    @inbounds for k in over_elems_real(grid)
        TKE_k = max(q[:tke, k, en], FT(0))
        θ_ρ = aux[:θ_ρ, k, gm]
        z = grid.zc[k]
        ts_dual = ActiveThermoState(param_set, q, aux, Dual(k), gm)
        θ_ρ_dual = virtual_pottemp.(ts_dual)
        ∇θ_ρ = ∇_z_flux(θ_ρ_dual, grid)
        buoyancy_freq = FT(grav(param_set)) * ∇θ_ρ / θ_ρ
        L[1] = sqrt(model.c_w * TKE_k) / buoyancy_freq
        ξ = z / obukhov_length
        κ_star = ustar / sqrt(TKE_k)
        L[2] = k_Karman * z / (model.c_K * κ_star * ϕ_m(ξ, a_L, b_L))
        S_squared =
            ∇_z_flux(q[:u, Dual(k), gm], grid)^2 +
            ∇_z_flux(q[:v, Dual(k), gm], grid)^2 +
            ∇_z_flux(q[:w, Dual(k), en], grid)^2
        R_g = aux[:∇buoyancy, k, gm] / S_squared
        if unstable(obukhov_length)
            Pr_z = model.Prandtl_neutral
        else
            discriminant1 = -4 * R_g + (1 + model.ω_2 * R_g)^2
            Pr_z =
                model.Prandtl_neutral *
                (1 + model.ω_2 * R_g - sqrt(discriminant1)) / (2 * R_g)
        end
        discriminant2 = S_squared - aux[:∇buoyancy, k, gm] / Pr_z
        L[3] =
            sqrt(model.c_ε / model.c_K) * sqrt(TKE_k) * 1 /
            sqrt(max(discriminant2, 1e-2))
        aux[:l_mix, k, gm] =
            sum([L[j] * exp(-L[j]) for j in 1:3]) / sum([exp(-L[j]) for j in 1:3])
    end
end

function compute_mixing_length!(
    grid::Grid{FT},
    q,
    aux,
    params,
    model::MinimumDissipation{FT},
) where {FT}
    @unpack obukhov_length, SurfaceModel, k_Karman, param_set = params
    gm, en, ud, sd, al = allcombinations(q)
    ustar = SurfaceModel.ustar
    L = Vector(undef, 3)
    a_L = model.a_L(obukhov_length)
    b_L = model.b_L(obukhov_length)
    gravitation = FT(grav(param_set))
    Rd::FT = FT(R_d(param_set))
    k_1 = first_interior(grid, Zmin())

    eps_vi::FT = FT(molmass_ratio(param_set))
    @inbounds for k in over_elems_real(grid)
        z = grid.zc[k]
        # precompute
        TKE_Shear =
            ∇_z_flux(q[:u, Dual(k), gm], grid)^2 +
            ∇_z_flux(q[:v, Dual(k), gm], grid)^2 +
            ∇_z_flux(q[:w, Dual(k), en], grid)^2
        Π = exner_given_pressure(param_set, aux[:p_0, k])
        lv::FT = latent_heat_vapor(param_set, aux[:t_cloudy, k]) # lv = latent_heat(t_cloudy)
        cpm = cp_m(param_set, PhasePartition(aux[:q_tot_cloudy, k]))
        TKE_k = max(q[:tke, k, en], FT(0))

        # compute L1 - static stability
        θ_ρ = aux[:θ_ρ, k, gm]
        ts_dual = ActiveThermoState(param_set, q, aux, Dual(k), en)
        θ_ρ_dual = virtual_pottemp.(ts_dual)
        ∇θ_ρ = ∇_z_flux(θ_ρ_dual, grid)
        buoyancy_freq = gravitation * ∇θ_ρ / θ_ρ
        if buoyancy_freq > 0
            L[1] = sqrt(model.c_w * TKE_k) / buoyancy_freq
        else
            L[1] = 1e-6
        end

        # compute L2 - law of the wall
        if obukhov_length < 0.0 #unstable case
            L[2] = (
                k_Karman * z / (
                    sqrt(max(q[:tke, k_1, en], FT(0)) / ustar / ustar) *
                    model.c_K
                ) * min((1 - 100 * z / obukhov_length)^0.2, 1 / k_Karman)
            )
        else # neutral or stable cases
            L[2] =
                k_Karman * z /
                (sqrt(max(q[:tke, k_1, en], FT(0)) / ustar / ustar) * model.c_K)
        end

        # I think this is an alrenative way for the same computation
        ξ = z / obukhov_length
        κ_star = ustar / sqrt(TKE_k)
        L[2] = k_Karman * z / (model.c_K * κ_star * ϕ_m(ξ, a_L, b_L))

        # compute L3 - entrainment detrainment sources
        prefactor = gravitation * (Rd * aux[:ρ_0, k] / aux[:p_0, k]) * Π
        dbdθl_dry = prefactor * (1 + (eps_vi - 1) * aux[:q_tot_dry, k])
        dbdqt_dry = prefactor * aux[:θ_dry, k] * (eps_vi - 1)
        if aux[:CF, k] > 0
            dbdθl_cloudy = (
                prefactor * (
                    1 +
                    eps_vi *
                    (1 + lh / Rv / aux[:t_cloudy, k]) *
                    aux[:q_vap_cloudy, k] - aux[:q_tot_cloudy, k]
                ) / (
                    1 +
                    lh * lh / cpm / Rv / aux[:t_cloudy, k] / aux[:t_cloudy, k] * aux[:q_vap_cloudy, k]
                )
            )
            dbdqt_cloudy =
                (lh / cpm / aux[:t_cloudy, k] * dbdθl_cloudy - prefactor) *
                aux[:θ_cloudy, k]
        else
            dbdθl_cloudy = 0
            dbdqt_cloudy = 0
        end
        # partial buoyancy gradients
        ∂b∂θl = (aux[:CF, k] * dbdθl_cloudy + (1 - aux[:CF, k]) * dbdθl_dry)
        ∂b∂qt = (aux[:CF, k] * dbdqt_cloudy + (1 - aux[:CF, k]) * dbdqt_dry)
        # chain rule
        # ∂b∂z_θl = ∇_z_flux(q[:θ_liq,k,en], grid) * ∂b∂θl
        # ∂b∂z_qt = ∇_z_flux(q[:q_tot,k,en], grid) * ∂b∂qt
        ∂b∂z_θl = ∇_z_flux(q[:θ_liq, Dual(k), en], grid) * ∂b∂θl
        ∂b∂z_qt = ∇_z_flux(q[:q_tot, Dual(k), en], grid) * ∂b∂qt
        Grad_Ri = min(
            ∂b∂z_θl / max(TKE_Shear, eps(FT)) +
            ∂b∂z_qt / max(TKE_Shear, eps(FT)),
            0.25,
        )
        if unstable(obukhov_length)
            Pr_z = model.Prandtl_neutral
        else
            Pr_z =
                model.Prandtl_neutral * (
                    2 * Grad_Ri / (
                        1 + (53 / 13) * Grad_Ri -
                        sqrt((1 + (53 / 130) * Grad_Ri)^2 - 4 * Grad_Ri)
                    )
                )
        end
        # Production/destruction terms
        a =
            model.c_ε *
            (TKE_Shear - ∂b∂z_θl / Pr_z - ∂b∂z_qt / Pr_z) *
            sqrt(TKE_k)
        # Dissipation term
        b = 0.0

        @inbounds for i in ud
            b +=
                q[:a, k, i] * q[:w, k, i] * aux[:δ_model, k, i] / q[:a, k, en] *
                (
                    (q[:w, k, i] - q[:w, k, en]) *
                    (q[:w, k, i] - q[:w, k, en]) / 2 - TKE_k
                ) -
                q[:a, k, i] *
                q[:w, k, i] *
                (q[:w, k, i] - q[:w, k, en]) *
                aux[:εt_model, k, i] *
                q[:w, k, en] / q[:a, k, en]
        end

        c_neg = model.c_ε * TKE_k * sqrt(TKE_k)
        if abs(a) > eps(FT) && 4 * a * c_neg > -b^2
            l_entdet = max(-b / 2.0 / a + sqrt(b^2 + 4 * a * c_neg) / 2 / a, 0)
        elseif abs(a) < eps(FT) && abs(b) > eps(FT)
            l_entdet = c_neg / b
        else
            l_entdet = 0.0
        end
        L[3] = l_entdet

        aux[:l_mix, k, gm] = lamb_smooth_minimum(L, 0.1, 1.5)
    end
end


text:
#####
##### `AbstractOscillator`
#####

abstract type AbstractOscillator <: AbstractSynthComponent end

function getfrequency end
function setfrequency! end

#####
##### `WaveForms`
#####

sine_wave(x) = sinpi(2x)
square_wave(x) = x < 0.5 ? 1.0 : -1.0
saw_wave(x) = 2x - 1
triangle_wave(x) = x < 0.5 ? 4x - 1 : 3 - 4x

#####
##### `Oscillator`
#####

mutable struct Oscillator{O<:Real,F,S<:AbstractFloat,T<:AbstractFloat} <: AbstractOscillator
    f::F
    frequency::S
    position::T
end

function Oscillator{O}(f, frequency::Real, position::Real = 0.0) where {O<:Real}
    frequency = float(frequency)
    position = float(position)
    return Oscillator{O,typeof(f),typeof(frequency),typeof(position)}(f, frequency, position)
end

function Oscillator(f, frequency::Real, position::Real = 0.0)
    frequency = float(frequency)
    position = float(position)
    return Oscillator{promote_op(f,typeof(position))}(f, frequency, position)
end

function next!(o::Oscillator{O})::O where {O}
    res = o.f(o.position)
    step = o.frequency / SAMPLE_RATE[]
    o.position = mod(o.position + step, true)
    return res
end

output_type(::Type{<:Oscillator{O}}) where {O} = O

getfrequency(o::Oscillator) = o.frequency
function setfrequency!(o::Oscillator{<:Any,<:Any,S}, frequency::Real) where {S}
    o.frequency = S(frequency)
end

#####
##### `FrequencyModulator`
#####

struct FrequencyModulator{
    O<:Real,
    S<:AbstractSynthComponent,
    T<:AbstractSynthComponent,
    U<:AbstractFloat,
    V<:NTuple{<:Any,<:Pair{<:AbstractOscillator,<:Any}},
} <: AbstractOscillator
    modulator::S
    output::T
    base_frequency::Ref{U}
    to_modulate::V

    function FrequencyModulator(
        modulator::S,
        output::T,
        base_frequency::Real,
        to_modulate::Pair{<:AbstractOscillator,<:Any}...,
    ) where {S<:AbstractSynthComponent,T<:AbstractSynthComponent}

        O = output_type(output)
        base_frequency = float(base_frequency)
        U = typeof(base_frequency)
        V = typeof(to_modulate)

        return new{O,S,T,U,V}(modulator, output, Ref(base_frequency), to_modulate)
    end
end

function FrequencyModulator(
    modulator::AbstractSynthComponent,
    output::AbstractOscillator,
    base_frequency::Real,
    modulation_function::Function,
)
    return FrequencyModulator(
        modulator,
        output,
        base_frequency,
        output => modulation_function,
    )
end

function FrequencyModulator(
    modulator::AbstractSynthComponent,
    output::AbstractOscillator,
    base_frequency::Real,
    weight::AbstractFloat,
)
    return FrequencyModulator(
        modulator,
        output,
        base_frequency,
        (base_frequency, modulation) -> base_frequency + weight * modulation,
    )
end

function FrequencyModulator(
    modulator::AbstractSynthComponent,
    output::AbstractOscillator,
    weight::AbstractFloat,
)
    return FrequencyModulator(
        modulator,
        output,
        getfrequency(output),
        weight,
    )
end

function next!(fm::FrequencyModulator)
    modulation = next!(fm.modulator)
    base_frequency = fm.base_frequency[]
    for (oscillator, modulation_function) in fm.to_modulate
        setfrequency!(oscillator, modulation_function(base_frequency, modulation))
    end
    return next!(fm.output)
end

output_type(::Type{<:FrequencyModulator{O}}) where {O} = O

getfrequency(fm::FrequencyModulator) = fm.base_frequency[]
setfrequency!(fm::FrequencyModulator, frequency::Real) = fm.base_frequency[] = float(frequency)


text:
using JuliaFormatter: format_text

handle("format") do data
  @destruct [text, indent || 4, margin || 92] = data

  formattedtext = format_text(text, indent = indent, margin = margin)
  Dict(:formattedtext => formattedtext)
end


text:
# TensorOperations compatiblity
#-------------------------------
function add!(α, A::SparseArray{<:Any, N}, CA::Symbol,
                β, C::SparseArray{<:Any, N}, indCinA) where {N}

    (N == length(indCinA) && TupleTools.isperm(indCinA)) ||
        throw(IndexError("Invalid permutation of length $N: $indCinA"))
    size(C) == TupleTools.getindices(size(A), indCinA) ||
        throw(DimensionMismatch("non-matching sizes while adding arrays"))

    β == one(β) || (iszero(β) ? _zero!(C) : LinearAlgebra.lmul!(β, C))
    for (IA, vA) in A.data
        IC = CartesianIndex(TupleTools.getindices(IA.I, indCinA))
        C[IC] += α* (CA == :C ? conj(vA) : vA)
    end
    C
end

function trace!(α, A::SparseArray{<:Any, NA}, CA::Symbol, β, C::SparseArray{<:Any, NC},
                indCinA, cindA1, cindA2) where {NA,NC}

    NC == length(indCinA) ||
        throw(IndexError("Invalid selection of $NC out of $NA: $indCinA"))
    NA-NC == 2*length(cindA1) == 2*length(cindA2) ||
        throw(IndexError("invalid number of trace dimension"))
    pA = (indCinA..., cindA1..., cindA2...)
    TupleTools.isperm(pA) ||
        throw(IndexError("invalid permutation of length $(ndims(A)): $pA"))

    sizeA = size(A)
    sizeC = size(C)

    TupleTools.getindices(sizeA, cindA1) == TupleTools.getindices(sizeA, cindA2) ||
        throw(DimensionMismatch("non-matching trace sizes"))
    sizeC == TupleTools.getindices(sizeA, indCinA) ||
        throw(DimensionMismatch("non-matching sizes"))

    β == one(β) || (iszero(β) ? _zero!(C) : LinearAlgebra.lmul!(β, C))
    for (IA, v) in A.data
        IAc1 = CartesianIndex(TupleTools.getindices(IA.I, cindA1))
        IAc2 = CartesianIndex(TupleTools.getindices(IA.I, cindA2))
        IAc1 == IAc2 || continue

        IC = CartesianIndex(TupleTools.getindices(IA.I, indCinA))
        C[IC] += α * (CA == :C ? conj(v) : v)
    end
    return C
end

function contract!(α, A::SparseArray, CA::Symbol, B::SparseArray, CB::Symbol,
                    β, C::SparseArray,
                    oindA::IndexTuple, cindA::IndexTuple,
                    oindB::IndexTuple, cindB::IndexTuple,
                    indCinoAB::IndexTuple, syms::Union{Nothing, NTuple{3,Symbol}} = nothing)

    pA = (oindA...,cindA...)
    (length(pA) == ndims(A) && TupleTools.isperm(pA)) ||
        throw(IndexError("invalid permutation of length $(ndims(A)): $pA"))
    pB = (oindB...,cindB...)
    (length(pB) == ndims(B) && TupleTools.isperm(pB)) ||
        throw(IndexError("invalid permutation of length $(ndims(B)): $pB"))
    (length(oindA) + length(oindB) == ndims(C)) ||
        throw(IndexError("non-matching output indices in contraction"))
    (ndims(C) == length(indCinoAB) && isperm(indCinoAB)) ||
        throw(IndexError("invalid permutation of length $(ndims(C)): $indCinoAB"))

    sizeA = size(A)
    sizeB = size(B)
    sizeC = size(C)

    csizeA = TupleTools.getindices(sizeA, cindA)
    csizeB = TupleTools.getindices(sizeB, cindB)
    osizeA = TupleTools.getindices(sizeA, oindA)
    osizeB = TupleTools.getindices(sizeB, oindB)

    csizeA == csizeB ||
        throw(DimensionMismatch("non-matching sizes in contracted dimensions"))
    TupleTools.getindices((osizeA..., osizeB...), indCinoAB) == size(C) ||
        throw(DimensionMismatch("non-matching sizes in uncontracted dimensions"))

    β == one(β) || (iszero(β) ? _zero!(C) : LinearAlgebra.lmul!(β, C))

    keysA = sort!(collect(nonzero_keys(A)),
                    by = IA->CartesianIndex(TupleTools.getindices(IA.I, cindA)))
    keysB = sort!(collect(nonzero_keys(B)),
                    by = IB->CartesianIndex(TupleTools.getindices(IB.I, cindB)))

    iA = iB = 1
    @inbounds while iA <= length(keysA) && iB <= length(keysB)
        IA = keysA[iA]
        IB = keysB[iB]
        IAc = CartesianIndex(TupleTools.getindices(IA.I, cindA))
        IBc = CartesianIndex(TupleTools.getindices(IB.I, cindB))
        if IAc == IBc
            Ic = IAc
            jA = iA
            while jA < length(keysA)
                if CartesianIndex(TupleTools.getindices(keysA[jA+1].I, cindA)) == Ic
                    jA += 1
                else
                    break
                end
            end
            jB = iB
            while jB < length(keysB)
                if CartesianIndex(TupleTools.getindices(keysB[jB+1].I, cindB)) == Ic
                    jB += 1
                else
                    break
                end
            end
            rA = iA:jA
            rB = iB:jB
            if length(rA) < length(rB)
                for kB in rB
                    IB = keysB[kB]
                    IBo = CartesianIndex(TupleTools.getindices(IB.I, oindB))
                    vB = B[IB]
                    for kA in rA
                        IA = keysA[kA]
                        IAo = CartesianIndex(TupleTools.getindices(IA.I, oindA))
                        IABo = CartesianIndex(IAo, IBo)
                        IC = CartesianIndex(TupleTools.getindices(IABo.I, indCinoAB))
                        vA = A[IA]
                        increaseindex!(C, α * (CA == :C ? conj(vA) : vA) * (CB == :C ? conj(vB) : vB), IC)
                    end
                end
            else
                for kA in rA
                    IA = keysA[kA]
                    IAo = CartesianIndex(TupleTools.getindices(IA.I, oindA))
                    vA = A[IA]
                    for kB in rB
                        IB = keysB[kB]
                        IBo = CartesianIndex(TupleTools.getindices(IB.I, oindB))
                        vB = B[IB]
                        IABo = CartesianIndex(IAo, IBo)
                        IC = CartesianIndex(TupleTools.getindices(IABo.I, indCinoAB))
                        increaseindex!(C, α * (CA == :C ? conj(vA) : vA) * (CB == :C ? conj(vB) : vB), IC)
                    end
                end
            end
            iA = jA+1
            iB = jB+1
        elseif IAc < IBc
            iA += 1
        else
            iB += 1
        end
    end
    C
end

# function contract_CSC!(α, A::SparseArray, CA::Symbol, B::SparseArray, CB::Symbol,
#                     β, C::SparseArray,
#                     oindA::IndexTuple, cindA::IndexTuple,
#                     oindB::IndexTuple, cindB::IndexTuple,
#                     indCinoAB::IndexTuple, syms::Union{Nothing, NTuple{3,Symbol}} = nothing)
#
#     pA = (oindA...,cindA...)
#     (length(pA) == ndims(A) && TupleTools.isperm(pA)) ||
#         throw(IndexError("invalid permutation of length $(ndims(A)): $pA"))
#     pB = (oindB...,cindB...)
#     (length(pB) == ndims(B) && TupleTools.isperm(pB)) ||
#         throw(IndexError("invalid permutation of length $(ndims(B)): $pB"))
#     (length(oindA) + length(oindB) == ndims(C)) ||
#         throw(IndexError("non-matching output indices in contraction"))
#     (ndims(C) == length(indCinoAB) && isperm(indCinoAB)) ||
#         throw(IndexError("invalid permutation of length $(ndims(C)): $indCinoAB"))
#
#     sizeA = size(A)
#     sizeB = size(B)
#     sizeC = size(C)
#
#     csizeA = TupleTools.getindices(sizeA, cindA)
#     csizeB = TupleTools.getindices(sizeB, cindB)
#     osizeA = TupleTools.getindices(sizeA, oindA)
#     osizeB = TupleTools.getindices(sizeB, oindB)
#
#     csizeA == csizeB ||
#         throw(DimensionMismatch("non-matching sizes in contracted dimensions"))
#     TupleTools.getindices((osizeA..., osizeB...), indCinoAB) == size(C) ||
#         throw(DimensionMismatch("non-matching sizes in uncontracted dimensions"))
#
#     β == one(β) || (iszero(β) ? _zero!(C) : LinearAlgebra.lmul!(β, C))
#     if isempty(nonzero_pairs(A)) || isempty(nonzero_pairs(B))
#         return C
#     end
#
#     # Build CSC-like representation, using internal Dict structure
#     # perform multiplication based on Gustafson's algorithm
#     Anz = nonzero_pairs(A)
#     Bnz = nonzero_pairs(B)
#     sortbyA = i->_subind(i, Anz, cindA)
#     sortbyB = i->_subind(i, Bnz, oindB)
#
#     indexA = sort!(findall(!iszero, Anz.slots); by = sortbyA)
#     indexB = sort!(findall(!iszero, Bnz.slots); by = sortbyB)
#     rowvalsA = map(i->_subind(i, Anz, oindA), indexA)
#     rowvalsB = map(i->_subind(i, Bnz, cindB), indexB)
#     colvalsA, colptrA = _uniqueranges(sortbyA, indexA)
#     colmapA = Dict(colvalsA[i]=>i for i in 1:length(colvalsA))
#     colvalsB, colptrB = _uniqueranges(sortbyB, indexB)
#     nzvalsA = Anz.vals[indexA]
#     nzvalsB = Bnz.vals[indexB]
#
#     colptrC = similar(colptrB)
#     rowvalsC = similar(rowvalsA, 0)
#     nzvalsC = Vector{eltype(C)}()
#     localrowsC = Dict{eltype(rowvalsA), Int}()
#
#     offset = 0
#     counter = 0
#     @inbounds for jB = 1:length(colvalsB)
#         IBo = colvalsB[jB]
#         colptrC[jB] = offset+1
#         rowsB = colptrB[jB]:(colptrB[jB+1]-1)
#         for kB in rowsB
#             IBc = rowvalsB[kB]
#             vB = CB == :C ? conj(nzvalsB[kB]) : nzvalsB[kB]
#             kA = get(colmapA, IBc, 0)
#             kA == 0 && continue
#             rowsA = colptrA[kA]:(colptrA[kA+1]-1)
#             for iA in rowsA
#                 IAo = rowvalsA[iA]
#                 vA = CA == :C ? conj(nzvalsA[iA]) : nzvalsA[iA]
#                 iC = get!(localrowsC, IAo, counter+1)
#                 if iC <= counter
#                     nzvalsC[offset + iC] += vA*vB
#                 else
#                     push!(nzvalsC, vA*vB)
#                     push!(rowvalsC, IAo)
#                     counter += 1
#                 end
#             end
#         end
#         offset += counter
#         counter = 0
#         empty!(localrowsC)
#     end
#     colptrC[end] = offset+1
#
#     # transfer CSC data to C
#     _sizehint!(C, 2*length(nzvalsC))
#     @inbounds for j in 1:length(colvalsB)
#         IBo = colvalsB[j]
#         for i in colptrC[j]:(colptrC[j+1]-1)
#             IAo = rowvalsC[i]
#             vC = nzvalsC[i]
#             IABo = CartesianIndex(IAo, IBo)
#             IC = CartesianIndex(TupleTools.getindices(IABo.I, indCinoAB))
#             increaseindex!(C, α * vC, IC)
#         end
#     end
#     return C
# end
#
# @inbounds _subind(i, d::Dict{<:CartesianIndex}, ind::IndexTuple) =
#     CartesianIndex(TupleTools.getindices(d.keys[i].I, ind))
#
# # for a sorted vector a, return a list of unique elements and a vector whose i'th entry is
# # the starting index of the i'th unique element in a
# function _uniqueranges(f, a::AbstractVector{T}) where T
#     isempty(a) && error("does not work on empty array")
#     i = 1
#     startindex = [1]
#     uniquevals = [f(a[1])]
#     @inbounds for j in 2:length(a)
#         v = f(a[j])
#         v == uniquevals[i] && continue
#         push!(uniquevals, v)
#         push!(startindex, j)
#         i += 1
#     end
#     push!(startindex, length(a)+1)
#     return uniquevals, startindex
# end


text:
using DynamicPolynomials
using MultivariateMoments
using Distributions
using PDMats

const MM=MultivariateMoments

include("../src/BeurlingLasso.jl")
using .BeurlingLasso

params = [[-0.8], [0.7]]
weights = [0.3, 0.7]

domain = reshape([-1, 1], (1, 2))

natoms = length(params)
ndims = length(params[1])

@polyvar x[1:ndims]
μ0 = AtomicMeasure(x, [WeightedDiracMeasure(params[i], weights[i]) for i=1:natoms])

gmm = MixtureModel(Normal[Normal(param[1], 0.5) for param in params],
                   weights)

λ = 1e-3
σ = 0.5

freq_distrib = Normal(0, 2.0)
nfreqs = 20
freqs = rand(freq_distrib, nfreqs)

function φ(x, n)
   μ = x[1]
   ω = freqs[n]
   exp(im*ω*μ - σ^2*ω^2/2)
end

y = [BeurlingLasso.expectation(μ0, x -> φ(x, n)) for n=1:nfreqs]

# function θ(t, n)
#    ω = freqs[n]
#    exp(im*ω*t)
# end

# nsamples = 1000
# samples = rand(gmm, nsamples)
# y = [sum(θ(t, n) for t in samples) / nsamples for n=1:nfreqs]

gridsize = 1000
niter = 10

prob = blasso(y, φ, λ, domain)
solver = sfwsolver(niter, gridsize)

solve!(prob, solver)


text:
using StaticArrays
using DifferentialDynamicsModels
using SimpleCarModels
using Test

for T in (Float32, Float64)
    dubins_bvp = DubinsSteering(r=T(0.1))
    reedsshepp_bvp = ReedsSheppSteering(r=T(0.1))

    dubins_bvp_v = DubinsSteering(v=T(10), r=T(0.1))
    reedsshepp_bvp_v = ReedsSheppSteering(v=T(10), r=T(0.1))

    for i in 1:1000
        q1 = rand(SE2State{T}).*SVector{3,T}(1, 1, 2π)
        q2 = rand(SE2State{T}).*SVector{3,T}(1, 1, 2π)
        dubins_cost, dubins_control = dubins_bvp(q1, q2)
        reedsshepp_cost, reedsshepp_control = reedsshepp_bvp(q1, q2)
        @test dubins_cost isa T
        @test reedsshepp_cost isa T
        @test dubins_control isa SVector{3,StepControl{2,T,VelocityCurvatureControl{T}}}
        @test reedsshepp_control isa SVector{5,StepControl{2,T,VelocityCurvatureControl{T}}}
        @test dubins_cost >= reedsshepp_cost - sqrt(eps(T))

        dubins_cost_v, dubins_control_v = dubins_bvp_v(q1, q2)
        reedsshepp_cost_v, reedsshepp_control_v = reedsshepp_bvp_v(q1, q2)
        @test dubins_cost ≈ 10*dubins_cost_v
        @test reedsshepp_cost ≈ 10*reedsshepp_cost_v

        @test propagate(dubins_bvp.dynamics, q1, dubins_control) ≈ q2
        @test propagate(reedsshepp_bvp.dynamics, q1, reedsshepp_control) ≈ q2
        @test propagate(dubins_bvp_v.dynamics, q1, dubins_control_v) ≈ q2
        @test propagate(reedsshepp_bvp_v.dynamics, q1, reedsshepp_control_v) ≈ q2
    end
end


text:
module UnicodeFun

include("symbols.jl")

include("sub_super_scripts.jl")
export to_superscript, to_subscript
export to_fraction, to_fraction_nl

include("fontstyles.jl")
export to_blackboardbold
export to_boldface
export to_italic
export to_caligraphic
export to_frakture
export to_underline
export to_overline

include("latex.jl")
export to_latex

include("roots.jl")
export to_root

end # module


text:
using RigidBodyDynamics

using RigidBodyDynamics: Bounds

const urdfdir = joinpath(@__DIR__, "..", "urdf")
const kindsldir = joinpath(@__DIR__, "..", "kindsl")
const sourcedir = joinpath(urdfdir, "source")
const simplifieddir = joinpath(urdfdir, "simplified")

const urdf2kindsl = joinpath(@__DIR__, "..", "urdf2kindsl", "urdf2kindsl.py")
const libdir = joinpath(@__DIR__, "..", "..", "libs")
const robcogen = joinpath(libdir, "robcogen", "robcogen", "exe", "robcogen.sh")
const robcogen_gen_dir = joinpath(libdir, "robcogen", "gen")

for urdf in filter(x -> endswith(x, ".urdf"), readdir(sourcedir))
    println("Processing $urdf")
    robot_name, _ = splitext(urdf)
    floatingfile = joinpath(sourcedir, "$(robot_name)_floating")
    isfile(floatingfile) || error("$floatingfile not found")
    floating = open(floatingfile) do io
        Bool(parse(Int, read(io, String)))
    end

    # Generate simplified URDF from raw URDF (remove fixed joints, change root joint type,
    # work around RBDL bug where links named "base_link" are misparsed)
    root_joint_type = floating ? QuaternionFloating{Float64}() : Fixed{Float64}()
    mechanism = parse_urdf(joinpath(sourcedir, urdf), root_joint_type=root_joint_type)
    for body in bodies(mechanism)
        if body.name == "base_link"
            i = 0
            while true
                new_name = "base_link_$i"
                if new_name ∈ map(b -> b.name, bodies(mechanism))
                    i += 1
                else
                    body.name = new_name
                    break
                end
            end
        end
    end
    simplified_urdf = joinpath(simplifieddir, "$robot_name.urdf")
    write_urdf(simplified_urdf, mechanism; robot_name=robot_name, include_root=!floating)

    # Generate kindsl files that are at least approximately the same using `urdf2kindsl`.
    # We suspect that that urdf2kindsl has a flaw that makes it produce .kindsl files that
    # aren't quite equivalent to the URDF files in some cases (notably, Atlas).
    kindsl = joinpath(kindsldir, "$robot_name.kindsl")
    run(pipeline(`$urdf2kindsl --digits=12 --pi-digits=12 $simplified_urdf`, stdout=kindsl))

    # Change `RobotBase x` to `RobotBase x floating` as appropriate
    if floating
        modified_kindsl_contents = replace(read(kindsl, String), r"(RobotBase\s*\S*)" => s"\g<0> floating")
        write(kindsl, modified_kindsl_contents)
    end

    # Generate URDFs from kindsl files
    robcogen_gen_urdf_option = 23
    robcogen_quit_option = 28
    cd(dirname(robcogen)) do # robcogen.sh needs to be run from its parent directory
        open(pipeline(`$robcogen $kindsl`, stdout=devnull), "w", stdout) do io
            println(io, robcogen_gen_urdf_option)
            println(io, robcogen_quit_option)
        end
    end
    robcogen_urdf = joinpath(robcogen_gen_dir, robot_name, "models", urdf)

    # Remove explicit floating joints
    mechanism = parse_urdf(robcogen_urdf, root_joint_type=Fixed{Float64}())
    root = root_body(mechanism)
    floating_joints = filter(isfloating, out_joints(root, mechanism))
    convert(Bool, length(floating_joints)) == floating || error()
    if !isempty(floating_joints)
        mechanism, = submechanism(mechanism, successor(floating_joints[1], mechanism))
    end

    # RBDL doesn't handle missing inertia tag for root body properly, so add one explicitly:
    root = root_body(mechanism)
    if root.inertia == nothing
        root.inertia = zero(SpatialInertia{Float64}, root_frame(mechanism))
    end

    # Pinocchio seems to misparse the 'continuous' joint type (model.nq is double what it should be)
    # so add some artificial joint limits, causing RigidBodyDynamics.jl to output 'revolute' rather
    # than 'continuous' joint types
    for joint in joints(mechanism)
        if joint_type(joint) isa Revolute
            joint.position_bounds .= Bounds(-1e5, 1e5)
        end
    end

    # Write final result to `urdfdir`
    write_urdf(joinpath(urdfdir, urdf), mechanism; robot_name=robot_name, include_root=true)
end

max_stars_repo_name: CobbledSteel/rbd-benchmarks-riscv

max_issues_repo_name: CobbledSteel/rbd-benchmarks-riscv

max_forks_repo_name: CobbledSteel/rbd-benchmarks-riscv


text:
module PMapReduceTiming

using ParallelUtilities
using Distributed

function initialize(sleeptime)
    A = Array{Int}(undef, 20, 20)
    for ind in eachindex(A)
        sleep(sleeptime)
        A[ind] = ind
    end
    return A
end

function main_mapreduce(sleeptime)
    mapreduce(x -> initialize(sleeptime), hcat, 1:20)
end

function main_pmapreduce(sleeptime)
    pmapreduce(x -> initialize(sleeptime), hcat, 1:20)
end

function compare_with_serial()
    # precompile
    main_mapreduce(0)
    main_pmapreduce(0)

    # time
    println("Tesing serial")
    A = @time main_mapreduce(5e-6)
    println("Tesing parallel")
    B = @time main_pmapreduce(5e-6)

    # check results
    println("Results match : ", A == B)
end

end


text:
module Redux

include("action.jl")
export AbstractAction, AbstractSyncAction, AbstractAsyncAction
using .ActionTypes

include("store.jl")
export AbstractStore, Store
export create_store, get_state, subscribe!, dispatch!

include("state.jl")
export AbstractState, AbstractMutableState, AbstractImmutableState

end # module


text:
using Test, Infiltrator, REPL

baz = 3
function f(x)
    x = 2
    y = [1, 2, 3]
    @infiltrate
    return x .+ y
end

function g(x)
    x *= 12
    @infiltrate x > 12
    return x
end

h(x) = [f(x) for x in x]

function i(x)
    y = [1.0 for _ in 1:x]
    @infiltrate
    return y
end

function j(x)
    xxxxx = 4
    @infiltrate
    return x*xxxxx
end

struct Bar
    xxx
    yyy
end
function k()
    zzzzz = 333
    aaaa = Bar(zzzzz, zzzzz)
    @infiltrate
    aaaa
end

module Jmod
using Random
using ..Infiltrator
jfunc() = @infiltrate
end

function multiexfiltrate()
    for i in 1:10
        @infiltrate
    end
end

function function_form_infiltration(x)
    if isdefined(Main, :Infiltrator)
        Main.infiltrate(@__MODULE__, Base.@locals, @__FILE__, @__LINE__)
    end
end

@testset "infiltration tests" begin
    if Sys.isunix() && VERSION >= v"1.1.0"
        using TerminalRegressionTests

        function run_terminal_test(func, result, commands, validation)
            TerminalRegressionTests.automated_test(joinpath(@__DIR__, validation), commands) do emuterm
            # TerminalRegressionTests.create_automated_test(joinpath(@__DIR__, validation), commands) do emuterm
                Infiltrator.end_session!()
                repl = REPL.LineEditREPL(emuterm, true)
                repl.interface = REPL.setup_interface(repl)
                repl.specialdisplay = REPL.REPLDisplay(repl)

                Infiltrator.TEST_TERMINAL_REF[] = repl.t
                Infiltrator.TEST_NOSTACK[] = true
                Infiltrator.TEST_REPL_REF[] = repl; Infiltrator.CHECK_TASK[] = false
                @test func() == result
                Infiltrator.TEST_TERMINAL_REF[] = nothing
                Infiltrator.TEST_NOSTACK[] = false
                Infiltrator.TEST_REPL_REF[] = nothing
            end
        end

        run_terminal_test(() -> f(3), [3, 4, 5],
                        ["?\n", "@trace\n", "@locals\n", "x.*y\n", "3+\n4\n", "ans\n", "baz\n", "0//0\n", "@toggle\n", "@toggle\n", "@toggle\n", "\x4"],
                        "Julia_f_$(VERSION.major).$(VERSION.minor).multiout")

        @test f(3) == [3, 4, 5] # `@toggle`d `@infiltrate` should not open a prompt

        Infiltrator.clear_disabled!()

        run_terminal_test(() -> f(3), [3, 4, 5],
                        ["@locals\n", "\x4"],
                        "Julia_f2_$(VERSION.major).$(VERSION.minor).multiout")

        run_terminal_test(() -> f(3), [3, 4, 5],
                        ["@locals x\n", "\x4"],
                        "Julia_f2_filter_$(VERSION.major).$(VERSION.minor).multiout")

        @test g(1) == 12 # conditional @infiltrate should not open a prompt

        run_terminal_test(() -> g(2), 24,
                        ["?\n", "@trace\n", "@locals\n", "x\n", "\x4"],
                        "Julia_g_$(VERSION.major).$(VERSION.minor).multiout")

        run_terminal_test(() -> h([1,2,3]), [[3,4,5], [3,4,5], [3,4,5]],
                        ["\x4", "@locals\n", "@exit\n"],
                        "Julia_h_$(VERSION.major).$(VERSION.minor).multiout")

        run_terminal_test(() -> i(1000), i(1000),
                        ["2+2\n", "@locals\n", "\x4"],
                        "Julia_i_$(VERSION.major).$(VERSION.minor).multiout")

        # scratch pad test
        run_terminal_test(() -> g(2), 24,
                        ["@locals\n", "xxxxx = 12\n", "aa, bb = ('a', 'b')\n","foo(x) = x\n", "function bar(x); 2x; end\n", "x = 2\n", "@exfiltrate xxxxx aa bb foo bar\n", "\x4"],
                        "Julia_exfil_$(VERSION.major).$(VERSION.minor).multiout")

        @test Infiltrator.store.xxxxx == 12
        @test Infiltrator.store.aa == 'a'
        @test Infiltrator.store.bb == 'b'
        @test Infiltrator.store.foo(3) == 3
        @test Infiltrator.store.bar(3) == 6
        @test !isdefined(getfield(Infiltrator.store, :store), :x)

        # proper scoping of scratch pad
        run_terminal_test(() -> j(2), 8,
                        ["@locals\n", "@exfiltrate\n", "\x4"],
                        "Julia_scoping_$(VERSION.major).$(VERSION.minor).multiout")

        @test Infiltrator.store.xxxxx == 4

        # persistent history
        run_terminal_test(() -> h([1,2,3]), [[3,4,5], [3,4,5], [3,4,5]],
                        ["y\n", "\e[A\n", "\x4", "\e[A\n", "@exit\n"],
                        "Julia_hist_$(VERSION.major).$(VERSION.minor).multiout")

        # top-level test
        run_terminal_test(() -> include(joinpath(@__DIR__, "toplevel-fixture.jl")), "success",
                        ["@exit\n"],
                        "Julia_toplevel_$(VERSION.major).$(VERSION.minor).multiout")

        # completions test
        run_terminal_test(k, Bar(333, 333),
                        ["struct Foo\n  xxx\n  yyy\nend\n", "foo = Foo(1, 2)\n", "fo\t\t\x3", "foo.xx\t\t\n", "zz\t\t\x3", "aa\t\t\x3", "aaaa.xx\t\t\n", "@exfiltrate foo nope\n", "@exit\n"],
                        "Julia_completions_$(VERSION.major).$(VERSION.minor).multiout")
        @test Infiltrator.store.foo.xxx == 1

        # imported globals
        run_terminal_test(Jmod.jfunc, nothing,
                          ["x\n", "randstring\n", "@exit\n"],
                          "Julia_imported_globals_$(VERSION.major).$(VERSION.minor).multiout")

        # safehouse should not shadow local variables
        run_terminal_test(multiexfiltrate, nothing,
                          ["i\n", "@continue\n", "i\n", "@exfiltrate\n", "@continue\n", "i\n", "safehouse.i\n", "@continue\n", "@exit\n"],
                          "Julia_multi_exfiltrate_$(VERSION.major).$(VERSION.minor).multiout")
        @test Infiltrator.store.i == 2

        # function-form infiltration
        run_terminal_test(() -> function_form_infiltration(2), nothing,
                          ["x\n", "@exit\n"],
                          "Julia_function_inf_$(VERSION.major).$(VERSION.minor).multiout")
    else
        @warn "Skipping UI tests on non unix systems"
    end
end

@testset "exfiltration tests" begin
    Infiltrator.clear_store!()
    function foo_ex(x)
        y = 3
        foo = :asd
        store = 33
        bar = function (yy) yy end
        @exfiltrate
    end
    foo_ex(55)

    @test Infiltrator.store.y == 3
    @test Infiltrator.store.foo == :asd
    @test Infiltrator.store.store == 33
    @test Infiltrator.store.bar(3) == 3
    @test_throws UndefVarError Infiltrator.store.asd

    # `@with` is basically for dynamic usage only
    @test 6 == Core.eval(@__MODULE__, :(Infiltrator.@withstore(2y)))
    @test "asd" == Core.eval(@__MODULE__, :(Infiltrator.@withstore(string(foo))))
end