text:
#=

Log-likelihood preparation given the data

Ignacio Quintero Mächler

t(-_-t)

Created 05 03 2020

=#

"""
    prepare_ll(X    ::Array{Array{Float64,1},1},
               p    ::Array{Float64,1},
               E0   ::Array{Float64,1},
               ns   ::Int64,
               k    ::Int64,
               h    ::Int64,
               ny   ::Int64,
               model::NTuple{3, Bool},
               abts ::Array{Float64,2},
               af!  ::Function)

Prepare **EGeoHiSSE** likelihoods using the **pruning** algorithm
given input data.
"""
function prepare_ll(X    ::Array{Array{Float64,1},1},
                    p    ::Array{Float64,1},
                    E0   ::Array{Float64,1},
                    ns   ::Int64,
                    k    ::Int64,
                    h    ::Int64,
                    ny   ::Int64,
                    model::NTuple{3, Bool},
                    abts ::Array{Float64,2},
                    af!  ::Function)

  # add extinction to likelihood vector
  for i in Base.OneTo(lastindex(X))
    append!(X[i], E0)
  end

  # make them vectors for indexing efficiency
  abts1 = abts[:,1]
  abts2 = abts[:,2]

  # make speciation events and closure
  λevent! = make_λevent(h, k, ny, false, model, af!)

  # make root likelihood conditioning
  rootll = make_rootll(h, k, ny, model, af!)

  # make ODE function
  ode_fun = make_egeohisse(Val(k), Val(h), Val(ny), Val(model), af!)

  # make integral problem
  prob = ODEProblem(ode_fun, zeros(2*h*(2^k-1)), (0.0,1.0), p)

  int = init(prob, Tsit5(),
             save_everystep  = false,
             calck           = false,
             force_dtmin     = true,
             save_start      = false,
             initialize_save = false,
             maxiters        = 100_000_000,
             verbose         = false)

  return X, int, λevent!, rootll, abts1, abts2
end

"""
    prepare_ll(cov_mod::NTuple{M,String},
               tv     ::Dict{Int64,Array{Float64,1}},
               ed     ::Array{Int64,2},
               el     ::Array{Float64,1},
               bts    ::Array{Float64,1},
               E0     ::Array{Float64,1},
               h      ::Int64;
               Eδt    ::Float64 = 0.01,
               ti     ::Float64 = 0.0) where{M}

Prepare **EGeoHiSSE** likelihoods using the **flow** algorithm given input data.
"""
function prepare_ll(p    ::Array{Float64,1},
                    bts  ::Array{Float64,1},
                    E0   ::Array{Float64,1},
                    k    ::Int64,
                    h    ::Int64,
                    ny   ::Int64,
                    ns   ::Int64,
                    ned  ::Int64,
                    model::NTuple{3, Bool},
                    Eδt  ::Float64,
                    ti   ::Float64,
                    abts ::Array{Float64,2},
                    af!  ::Function)

  # Estimate extinction at `ts` times
  tf = maximum(bts)
  ts = [ti:Eδt:tf...]

  # Make extinction integral
  egeohisse_E = make_egeohisse_E(Val(k), Val(h), Val(ny), Val(model), af!)

  # Make extinction function at times `ts`
  Et = make_Et(egeohisse_E, p, E0, ts, ti, tf)

  # estimate first Extinction probabilities at times `ts`
  Ets  = Et(p)
  nets = length(Ets) + 1

  # Make extinction approximated function
  afE! = make_af(ts,  Ets, Val(ns))

  # make likelihood integral
  ode_intf =
    make_egeohisse_M(Val(k), Val(h), Val(ny), Val(model), af!, afE!, nets)

  # push parameters as the last vector in Ets
  push!(Ets, p)

  # make Gt function
  Gt = make_Gt(ode_intf, Ets, Matrix{Float64}(I, ns, ns), bts, ti, tf)

  # sort branching times
  sort!(bts)

  # add the present `0.0`
  pushfirst!(bts, 0.0)

  nbts = lastindex(bts)

  # estimate `Gts` according to `Ets` and `p0`
  Gts = Gt(Ets)

  # connect branches with branching times
  lbts = Array{Int64,2}(undef,size(abts))

  # preallocate absolute minimum matrix of differences
  minM = Array{Float64,1}(undef, nbts)
  # find links
  for i in Base.OneTo(2ned)
    for j in Base.OneTo(nbts)
      minM[j] = abs(abts[i] - bts[j])
    end
    lbts[i] = argmin(minM)
  end

  # make λevent!
  λevent! = make_λevent(h, k, ny, true, model, af!)

  # make root likelihood conditioning
  rootll = make_rootll(h, k, ny, model, af!)

  return Gt, Et, lbts, nets, λevent!, rootll
end

"""
    prepare_ll(ode_make_lik,
               ode_ext,
               tv  ::Dict{Int64,Array{Float64,1}},
               ed  ::Array{Int64,2},
               el  ::Array{Float64,1},
               bts ::Array{Float64,1},
               p0  ::Array{Float64,1},
               E0  ::Array{Float64,1},
               k   ::Int64,
               h   ::Int64,
               ntip::Int64;
               Eδt::Float64 = 0.01,
               ti ::Float64 = 0.0)

Prepare **MuSSE** likelihoods for flow algorithm given input model
"""
function prepare_ll(ode_make_lik,
                    ode_ext,
                    tv  ::Dict{Int64,Array{Float64,1}},
                    ed  ::Array{Int64,2},
                    el  ::Array{Float64,1},
                    bts ::Array{Float64,1},
                    p0  ::Array{Float64,1},
                    E0  ::Array{Float64,1},
                    k   ::Int64,
                    h   ::Int64,
                    ntip::Int64;
                    Eδt::Float64 = 0.01,
                    ti ::Float64 = 0.0)

  # number of states
  ns = k*h

  # Estimate extinction at `ts` times
  tf = maximum(bts)
  ts = [ti:Eδt:tf...]

  # Make extinction integral
  Et = make_Et(ode_ext, p0, E0, ts, ti, tf)

  # estimate first Extinction probabilities at times `ts`
  Ets  = Et(p0)
  nets = length(Ets) + 1

  # Make extinction approximated function
  # ** this make order is crucial **
  afE! = make_af(ts,  Ets, Val(ns))

  # make likelihood integral
  ode_intf = ode_make_lik(afE!, nets)

  # push parameters as the last vector in Ets
  push!(Ets, p0)

  # make Gt function
  Gt = make_Gt(ode_intf, Ets, Matrix{Float64}(I, ns, ns), bts, ti, tf)

  # sort branching times
  sort!(bts)

  # add the present `0.0`
  pushfirst!(bts, 0.0)

  nbts = lastindex(bts)

  # estimate `Gts` according to `Ets` and `p0`
  Gts = Gt(Ets)

  ## link edges with initial and end branching times
  abts = abs_time_branches(el, ed, ntip)
  lbts = Array{Int64,2}(undef,size(abts))

  # preallocate absolute minimum matrix of differences
  minM = Array{Float64,1}(undef, nbts)
  # find links
  for i in Base.OneTo(ned*2)
    for j in Base.OneTo(nbts)
      minM[j] = abs(abts[i] - bts[j])
    end
    lbts[i] = argmin(minM)
  end

  # make internal node triads
  triads = maketriads(ed)

  # initialize X for tips for GeoHiSSE
  X = [Array{Float64,1}(undef,ns) for i in Base.OneTo(ned)]

  # initialize X for tips
  for i in Base.OneTo(ned)
    edi = ed[i,2]
    edi > ntip && continue
    X[i] = tv[edi]
  end

  return Gt, Et, X, triads, lbts, ns, ned, nets
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
function bubblesort!(arr::AbstractVector)
    for _ in 2:length(arr), j in 1:length(arr)-1
        if arr[j] > arr[j+1]
            arr[j], arr[j+1] = arr[j+1], arr[j]
        end
    end
    return arr
end

array = rand(1:100, 5)
println("Unsorted array: $array\nSorted array: ", bubblesort!(array))

ext: jl

max_stars_repo_name: russellbanks/BubbleSorts

max_issues_repo_name: russellbanks/BubbleSorts

max_forks_repo_name: russellbanks/BubbleSorts


text:
# run this script to update the pages for Documenter.jl

# to generate doc pages
using Literate

function run_literate(names, folder)
        for t in names
                fp = joinpath(folder, string(t, ".jl"))
                output_folder=joinpath(@__DIR__, "..", "docs", "src", folder)
                println(output_folder)
                Literate.markdown(fp, output_folder, documenter=true,credit=false)
#                Literate.notebook(fp, output_folder, documenter=true,credit=false)
        end
end

run_literate(["tut1", "tut2"], "")


text:
immutable S2C <: Model; rnn; fnn; params;
    function S2C(rnn::Function,fnn::Function; o...)
        r = Net(rnn; o...)
        f = Net(fnn; o...)
        p = vcat(params(r), params(f))
        new(r, f, p)
    end
end

params(m::S2C)=m.params
reset!(m::S2C;o...)=(reset!(m.rnn;o...); reset!(m.fnn;o...))

train(m::S2C, data, loss; o...)=s2c_loop(m, data, loss; trn=true, o...)
test(m::S2C, data, loss; o...)=(l=zeros(2); s2c_loop(m, data, loss; losscnt=l, o...); l[1]/l[2])

function s2c_loop(m::S2C, data, loss; gcheck=false, o...)
    reset!(m; o...)
    for item in data
        (x,ygold) = item2xy(item)
        ycell = forw(m.rnn, x...; seq=true, o...)
        ygold == nothing && continue
        ypred = forw(m.fnn, ycell; o...)
        s2c_eos(m, x, ygold, ypred, loss; gcheck=gcheck, o...)
        gcheck && break
        reset!(m; o...)
    end
end

function s2c_eos(m::S2C, x, ygold, ypred, loss; trn=false, gcheck=false, maxnorm=nothing, losscnt=nothing, gclip=0, o...)
    if losscnt != nothing
        losscnt[1] += loss(ypred, ygold)
        losscnt[2] += 1
    end
    if trn
        ygrad = back(m.fnn, ygold, loss; getdx=true, o...)
        back(m.rnn, ygrad; seq=true, o...)                      # TODO: avoid unnecessary copy of ygrad
        while m.rnn.sp > 0
            back(m.rnn; seq=true, o...)
        end
        g = (gclip > 0 || maxnorm!=nothing ? gnorm(m) : 0)
        if !gcheck
            gclip=(g > gclip > 0 ? gclip/g : 0)
            update!(m; gclip=gclip, o...)
        end
    end
    if maxnorm != nothing
        w=wnorm(m)
        w > maxnorm[1] && (maxnorm[1]=w)
        g > maxnorm[2] && (maxnorm[2]=g)
    end
end

### DEAD CODE

# Decided not to expose the Net() function.  Users only interact with
# knet functions and models.  Having an intermediate Net object is
# confusing.  To solve the problem of different initialization
# parameters, one could just define new knet functions.

# S2C(rnn::Net,fnn::Net)=S2C(rnn,fnn,vcat(params(rnn),params(fnn)))

# function train(m::S2C, data, loss; getloss=true, getnorm=true, gclip=0, o...)
#     numloss = sumloss = maxwnorm = maxgnorm = w = g = 0
#     for item in data
#         (x,ygold) = item2xy(item)
#         ycell = forw(m.rnn, x...; trn=true, seq=true, o...)     # TODO: do we need trn/seq everywhere?
#         ygold == nothing && continue                            # ygold!=nothing signals end of sequence
#         ypred = forw(m.fnn, ycell; trn=true, o...)              # TODO: avoid unnecessary copy of ycell
#         getloss && (sumloss += loss(ypred, ygold); numloss += 1)
#         ygrad = back(m.fnn, ygold, loss; getdx=true, o...)
#         back(m.rnn, ygrad; seq=true, o...)                      # TODO: avoid unnecessary copy of ygrad
#         while m.rnn.sp > 0
#             back(m.rnn; seq=true, o...)
#         end
#         (getnorm || gclip>0) && (g = gnorm(m); g > maxgnorm && (maxgnorm = g))
#         update!(m; gclip=(g > gclip > 0 ? gclip/g : 0), o...)
#         getnorm && (w = wnorm(m); w > maxwnorm && (maxwnorm = w))
#         reset!(m.rnn; o...)
#     end
#     return (sumloss/numloss, maxwnorm, maxgnorm)
# end

# function test(m::S2C, data, loss; o...)
#     sumloss = numloss = 0
#     for item in data
#         (x,ygold) = item2xy(item)
#         ycell = forw(m.rnn, x...; trn=false, seq=true, o...)
#         ygold == nothing && continue
#         ypred = forw(m.fnn, ycell; trn=false, o...)
#         sumloss += loss(ypred, ygold); numloss += 1
#         reset!(m.rnn; o...)
#     end
#     sumloss / numloss
# end


text:
# This file is part of Kpax3. License is MIT.

function processchain(x::AminoAcidData,
                      fileroot::String)
  # open and close output files to check writing permissions
  # fail early if we can't write, saving the computational time required to
  # process the chain
  # use "append" for a non destructive operation
  fp = open(string(fileroot, "_posterior_k.csv"), "a")
  close(fp)

  fp = open(string(fileroot, "_posterior_R.csv"), "a")
  close(fp)

  fp = open(string(fileroot, "_posterior_C.csv"), "a")
  close(fp)

  n = zeros(Int, 1)
  m = zeros(Int, 1)
  N = zeros(Int, 1)

  k = zeros(Int, 1)

  h = 0

  fpR = open(string(fileroot, "_row_partition.bin"), "r")

  read!(fpR, n)
  read!(fpR, m)
  read!(fpR, N)

  R = zeros(Int, n[1])

  pk = zeros(Float64, n[1])
  pR = zeros(Float64, div(n[1] * (n[1] - 1), 2))

  T = 0
  try
    while !eof(fpR)
      read!(fpR, k)
      read!(fpR, R)

      pk[k[1]] += 1

      h = 0
      for i in 1:(n[1] - 1), j in (i + 1):n[1]
        h += 1
        pR[h] += float(R[i] == R[j])
      end

      T += 1
    end
  finally
    close(fpR)
  end

  if T != N[1]
    warn(string("Expecting ", N[1], " simulations but instead found ", T, "."))
  end

  pk /= T
  pR /= T

  fp = open(string(fileroot, "_posterior_k.csv"), "w")
  write(fp, "\"k\",\"p(k | x)\"\n")
  for g in 1:n[1]
    if pk[g] > 0.0
      write(fp, string(g, ",", pk[g], "\n"))
    end
  end
  close(fp)

  fp = open(string(fileroot, "_posterior_R.csv"), "w")
  write(fp, "\"Unit i\",\"Unit j\",\"p(i ~ j | x)\"\n")
  h = 0
  for i in 1:(n[1] - 1), j in (i + 1):n[1]
    h += 1
    write(fp, string("\"", x.id[i], "\",", "\"", x.id[j], "\",", pR[h], "\n"))
  end
  close(fp)

  fpC = open(string(fileroot, "_col_partition.bin"), "r")

  read!(fpC, n)
  read!(fpC, m)
  read!(fpC, N)

  C = zeros(UInt8, m[1])

  pC = zeros(Float64, 3, m[1])

  T = 0
  try
    while !eof(fpC)
      readbytes!(fpC, C, m[1])

      for b in 1:m[1]
        pC[C[b], b] += 1
      end

      T += 1
    end
  finally
    close(fpC)
  end

  if T != N[1]
    warn(string("Expecting ", N[1], " simulations but instead found ", T, "."))
  end

  pC /= T

  # estimate the probabilities of amino acids at each site
  #
  # Use a Dirichlet(1/J, ..., 1/J) as the prior distibution for the vector of
  # probabilities, where J is the total number of unique amino acids observed
  # at the current site. If y[j] is the count of amino acid j, then
  # a[j] = y[j] + 1 / J
  # is the corresponding parameter of the posterior probability.
  #
  # Posterior expected value is p = a / sum(a) = (y + 1 / J) / (n + 1)
  p = sum(float(x.data), dims=2)

  key = 1
  idxstart = 1
  a = 0.0
  count = 0.0
  for idxend in 1:m[1]
    if x.key[idxend] != key
      idx = idxstart:(idxend - 1)
      a = 1 / length(idx)

      for i in idx
        p[i] = (p[i] + a) / (n[1] + 1)
      end

      key += 1
      idxstart = idxend
      count = 0.0
    end

    count += p[idxend]
  end

  idx = idxstart:m[1]
  a = 1 / length(idx)

  for i in idx
    p[i] = (p[i] + a) / (n[1] + 1)
  end

  fp = open(string(fileroot, "_posterior_C.csv"), "w")
  write(fp, string("\"Site\",\"AminoAcid\",\"Proportion\",\"p('noise' | x)\",",
                   "\"p('weak signal' | x)\",\"p('strong signal' | x)\"\n"))

  isuint8 = eltype(x.ref) == UInt8

  polval = if isuint8
              UInt8('.')
            else
              "."
            end

  key = 0
  b = 1
  for j in 1:length(x.ref)
    if x.ref[j] == polval
      key += 1

      while (b <= m[1]) && (x.key[b] == key)
        write(fp, string(j, ",\"",
                        uppercase(isuint8 ? string(Char(x.val[b])) : x.val[b]),
                        "\",", p[b], ",", pC[1, b], ",", pC[2, b], ",",
                        pC[3, b], "\n"))
        b += 1
      end
    else
      write(fp, string(j, ",\"",
                       uppercase(isuint8 ? string(Char(x.ref[j])) : x.ref[j]),
                       "\",", 1.0, ",", 1.0, ",", 0.0, ",", 0.0, "\n"))
    end
  end
  close(fp)

  nothing
end

function savestate!(fpR::IOStream,
                    fpC::IOStream,
                    state::AminoAcidState)
  write(fpR, state.k)
  write(fpR, state.R)

  for b in 1:size(state.C, 2)
    write(fpC, state.C[state.cl[1], b] < 0x03 ? state.C[state.cl[1], b] : 0x03)
  end

  nothing
end


text:
function download(
    emod :: ERA5Module,
    evar :: ERA5Variable,
    ereg :: ERA5Region,
    ispy :: Bool = false
)

    downloadcheckevar(evar)

    if ispy
          pythonprint(emod,evar,ereg)
    else; cdsretrieve(emod,evar,ereg)
    end

end

function downloadcheckevar(
    evar :: Union{SingleVariable,PressureVariable}
)

    @info "$(modulelog()) - The ERA5Variable ID \"$(evar)\" is a valid ERA5 Variable identifier and therefore can be used to download data from the Climate Data Store"

end

function downloadcheckevar(
    evar :: Union{SingleCustom,PressureCustom}
)

    error("$(modulelog()) - The ERA5Variable ID \"$(evar)\" is a custom ERA5 Variable identifier and therefore cannot be used to download data from the Climate Data Store, and can only be calculated/analyzed from preexisting valid ERA5 data")

end


text:
using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.Mesh.Filters
using CLIMA.MPIStateArrays
using CLIMA.DGBalanceLawDiscretizations
using Printf
using LinearAlgebra
using Logging

const ArrayTypes = (CLIMA.array_type(),)

function initialcondition!(Q, x, y, z, _)
  @inbounds Q[1] = abs(x) - 0.1
end

using Test
function run(mpicomm, dim, ArrayType, Ne, N, FT)
  brickrange = ntuple(j->range(FT(-1); length=Ne[j]+1, stop=1), dim)
  topl = StackedBrickTopology(mpicomm, brickrange,
                              periodicity=ntuple(j->true, dim))

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                         )

  spacedisc = DGBalanceLaw(grid = grid,
                           length_state_vector = 1,
                           flux! = (x...) -> (),
                           numerical_flux! = (x...) -> ())

  Q = MPIStateArray(spacedisc, initialcondition!)

  initialsumQ = weightedsum(Q)
  @test minimum(Q.realdata) < 0

  Filters.apply!(Q, 1, spacedisc.grid, TMARFilter())

  sumQ = weightedsum(Q)

  @test minimum(Q.realdata) >= 0
  @test isapprox(initialsumQ, sumQ; rtol = 10eps(FT))
end

let
  CLIMA.init()
  mpicomm = MPI.COMM_WORLD

  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
             ll == "WARN"  ? Logging.Warn  :
             ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))

  numelem = (2, 2, 2)
  polynomialorder = 4

  @testset "$(@__FILE__)" for ArrayType in ArrayTypes
    for FT in (Float64, Float32)
      for dim = 2:3
        @info (ArrayType, FT, dim)
        run(mpicomm, dim, ArrayType, numelem, polynomialorder, FT)
      end
    end
  end
end

nothing

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
using Documenter, NMRTools, Plots
ENV["GKSwstype"] = "100" # https://github.com/jheinen/GR.jl/issues/278

makedocs(;
    modules=[NMRTools],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Examples" => "examples.md",
        "Metadata" => "metadata.md",
    ],
    repo="https://github.com/chriswaudby/NMRTools.jl/blob/{commit}{path}#L{line}",
    sitename="NMRTools.jl",
    authors="Chris Waudby",
    assets=String[],
)

deploydocs(;
    repo="github.com/chriswaudby/NMRTools.jl",
)


text:
using DelimitedFiles, ProgressMeter, PyPlot, Gurobi

include("./util.jl")
include("./parsers.jl")
include("./Algorithms/TEAVAR.jl")
include("./Algorithms/FFC.jl")
include("./simulation.jl")

function network_utilization(topology,
                            demand_downscale,
                            capacity_downscale,
                            num_demands,
                            iterations,
                            cutoff,
                            availabilities;
                            k=20,
                            paths="SMORE",
                            weibull_scale=.0001,
                            plot=true,
                            dirname="./data/raw/utilization")

    env = Gurobi.Env()
    dir = nextRun(dirname)

    scenarios_all = []
    scenario_probs_all = []
    links, capacity, link_probs, nodes = readTopology(topology, downscale=1)
    for i in 1:iterations
        link_probs = weibullProbs(length(links), shape=.8, scale=.0005)
        scenarios, probs = subScenarios(link_probs, cutoff, first=true, last=true)
        push!(scenarios_all, scenarios)
        push!(scenario_probs_all, probs)
    end
    links, capacity, link_probs, nodes = readTopology(topology, downscale=1)

    min_routed = 0
    max_routed = Inf
    progress = ProgressMeter.Progress(num_demands*iterations*length(availabilities), .1, "Computing Utilization Plot for TEAVAR...", 50)

    throughputs_teavar = [[] for i in 1:length(availabilities)]
    edge_utilization_percentages = [[] for i in 1:length(availabilities)]
    for b in 1:length(availabilities)
        throughputs = []
        edge_utilizations = []
        for d in 1:num_demands
            avg_throughput = 0
            for i in 1:iterations
                scenarios = scenarios_all[i]
                probs = scenario_probs_all[i]
                demand, flows = readDemand("$(topology)/demand", length(nodes), d, scale=1, downscale=demand_downscale, matrix=true)
                T, Tf, g = getTunnels(nodes, links, capacity, flows, k, edge_disjoint=false)
                var, cvar, a = TEAVAR(env, links, capacity, flows, demand, availabilities[b], k, T, Tf, scenarios, probs)
                # beta_satisfied += sum(allowed)/sum(demand)
                bandwidth_allowed = ones(length(flows)) * (1-var)
                edge_utilization, routed = simulateUtilizationNoFailures(links, capacity, demand, flows, T, Tf, k, a, bandwidth_allowed)
                edge_utilization_percentage = edge_utilization ./ capacity
                throughput = sum(routed)
                min_routed = min(throughput, min_routed)
                max_routed = max(throughput, max_routed)

                avg_throughput += throughput
                push!(throughputs, throughput)
                push!(edge_utilizations, edge_utilization)
                edge_utilization_percentages[b] = vcat(edge_utilization_percentages[b], edge_utilization_percentage)
                ProgressMeter.next!(progress, showvalues = [(:beta, availabilities[b]), (:iterations, i), (:num_demands, d), (:throughput, throughput)])
            end
            push!(throughputs_teavar[b], avg_throughput/iterations)
        end
        writedlm("$(dir)/$(availabilities[b])_throughputs", throughputs)
        writedlm("$(dir)/$(availabilities[b])_utilization", edge_utilizations)
        writedlm("$(dir)/$(availabilities[b])_utilization_percentage", edge_utilization_percentages)
    end

    progress = ProgressMeter.Progress(num_demands, .1, "Computing Utilization Plot for FFC1...", 50)
    throughputs_ffc1 = []
    edge_utilizations_ffc1 = []
    edge_utilization_percentages_ffc1 = []
    for d in 1:num_demands
        demand, flows = readDemand("$(topology)/demand", length(nodes), d, scale=1, downscale=demand_downscale, matrix=true)
        # T, Tf, g = getTunnels(nodes, links, capacity, flows, k, edge_disjoint=false)
        T, Tf, g = getTunnels(nodes, links, capacity, flows, k, edge_disjoint=true)
        a, b = FFC(env, links, capacity, flows, demand, 1, T, Tf)
        bandwidth_allowed = ones(length(flows)) .* (b ./ demand)
        edge_utilization, routed = simulateUtilizationNoFailures(links, capacity, demand, flows, T, Tf, k, a, bandwidth_allowed)
        edge_utilization_percentage = edge_utilization ./ capacity
        throughput = sum(routed)
        min_routed = min(throughput, min_routed)
        max_routed = max(throughput, max_routed)

        push!(throughputs_ffc1, throughput)
        push!(edge_utilizations_ffc1, edge_utilization)
        # push!(edge_utilization_percentages_ffc1, edge_utilization_percentage)
        edge_utilization_percentages_ffc1 = vcat(edge_utilization_percentages_ffc1, edge_utilization_percentage)
        ProgressMeter.next!(progress, showvalues = [(:num_demands, d), (:throughput, throughput)])
    end
    writedlm("$(dir)/FFC1_throughputs", throughputs_ffc1)
    writedlm("$(dir)/FFC1_utilization", edge_utilizations_ffc1)
    writedlm("$(dir)/FFC1_utilization_percentage", edge_utilization_percentages_ffc1)

    progress = ProgressMeter.Progress(num_demands, .1, "Computing Utilization Plot for FFC2...", 50)
    throughputs_ffc2 = []
    edge_utilizations_ffc2 = []
    edge_utilization_percentages_ffc2 = []
    for d in 1:num_demands
        demand, flows = readDemand("$(topology)/demand", length(nodes), d, scale=1, downscale=demand_downscale, matrix=true)
        # T, Tf, g = getTunnels(nodes, links, capacity, flows, k, edge_disjoint=false)
        T, Tf, g = getTunnels(nodes, links, capacity, flows, k, edge_disjoint=true)
        a, b = FFC(env, links, capacity, flows, demand, 2, T, Tf)
        bandwidth_allowed = ones(length(flows)) .* (b ./ demand)
        edge_utilization, routed = simulateUtilizationNoFailures(links, capacity, demand, flows, T, Tf, k, a, bandwidth_allowed)
        edge_utilization_percentage = edge_utilization ./ capacity
        throughput = sum(routed)
        min_routed = min(throughput, min_routed)
        max_routed = max(throughput, max_routed)

        push!(throughputs_ffc2, throughput)
        push!(edge_utilizations_ffc2, edge_utilization)
        # push!(edge_utilization_percentages_ffc2, edge_utilization_percentage)
        edge_utilization_percentages_ffc2 = vcat(edge_utilization_percentages_ffc2, edge_utilization_percentage)
        ProgressMeter.next!(progress, showvalues = [(:num_demands, d), (:throughput, throughput)])
    end
    writedlm("$(dir)/FFC2_throughputs", throughputs_ffc2)
    writedlm("$(dir)/FFC2_utilization", edge_utilizations_ffc2)
    writedlm("$(dir)/FFC2_utilization_percentage", edge_utilization_percentages_ffc2)
    writedlm("$(dir)/params", [[ "topology", "demand_downscale", "capacity_downscale", "num_demands", "iterations", "cutoff", "availabilities"], [topology, demand_downscale, capacity_downscale, num_demands, iterations, cutoff, availabilities]])

    if plot
        # PyPlot.clf()
        # println(num_demands)
        # println(iterations)
        # x_vals = collect(1:num_demands)
        # for b in 1:length(availabilities)
        #     PyPlot.plot(x_vals, throughputs_teavar[b])
        # end
        # PyPlot.plot(x_vals, throughputs_ffc1)
        # PyPlot.plot(x_vals, throughputs_ffc2)
        # PyPlot.legend(vcat(map(elt -> string(elt), availabilities), ["FFC_1", "FFC_2"]), loc="upper right")
        # PyPlot.title("CVaR Graph")
        # PyPlot.xlabel("Demand Number", fontweight="bold")
        # PyPlot.ylabel("Throughput", fontweight="bold")
        # PyPlot.show()
        CDF(availabilities, edge_utilization_percentages, edge_utilization_percentages_ffc1, edge_utilization_percentages_ffc2, dir)
    end

end

# network_utilization(["TEAVAR", "FFC-1", "FFC-2"], "B4", 1, 1, 3, 1, .00001, [.99, .999, .9999])

function CDF(availabilities, edge_utilization_percentages_teavar, edge_utilization_percentages_ffc1, edge_utilization_percentages_ffc2, dir)
    PyPlot.clf()
    labels = []
    for (i,vals) in enumerate(edge_utilization_percentages_teavar)
        sorted = sort(vals)
        x_vals = sorted
        # y_vals = cumsum(sorted) ./ sum(sorted)
        y_vals = cumsum(ones(length(sorted)) ./ sum(ones(length(sorted))))
        Z = zeros(length(x_vals), 2)
        Z[:,1] = x_vals
        Z[:,2] = y_vals
        # writedlm("./gnuplot/data/guarantee/$(i)", Z)
        push!(labels, string(availabilities[i]))
        PyPlot.plot(vcat(0, x_vals, 1), vcat(0, y_vals, 1))
    end

    sorted = sort(edge_utilization_percentages_ffc1)
    x_vals = sorted
    # y_vals = cumsum(sorted) ./ sum(sorted)
    y_vals = cumsum(ones(length(sorted)) ./ sum(ones(length(sorted))))
    Z = zeros(length(x_vals), 2)
    Z[:,1] = x_vals
    Z[:,2] = y_vals
    # writedlm("./gnuplot/data/guarantee/ffc_1", Z)
    push!(labels, "FFC_1")
    PyPlot.plot(vcat(0, x_vals, 1), vcat(0, y_vals, 1))

    sorted = sort(edge_utilization_percentages_ffc2)
    x_vals = sorted
    # y_vals = cumsum(sorted) ./ sum(sorted)
    y_vals = cumsum(ones(length(sorted)) ./ sum(ones(length(sorted))))
    Z = zeros(length(x_vals), 2)
    Z[:,1] = x_vals
    Z[:,2] = y_vals
    # writedlm("./gnuplot/data/guarantee/ffc_2", Z)
    push!(labels, "FFC_2")
    PyPlot.plot(vcat(0, x_vals, 1), vcat(0, y_vals, 1))
    PyPlot.title("title")
    PyPlot.xlabel("Utilization Percentage", fontweight="bold")
    PyPlot.ylabel("CDF", fontweight="bold")
    PyPlot.legend(labels, loc="lower right")
    PyPlot.savefig("$(dir)/plot.png")
    PyPlot.show()
end


text:
#=
Author : Guy Farmer
generates a distribution of ambulances for various ambulance counts and stores to a jld file
=#
 using AmbulanceDeployment

hourly_calls = CSV.File("../test/austin-data/Full_WeekdayCalls.csv") |> DataFrame
adjacent_nbhd = CSV.File("../test/austin-data/adjacent_nbhd.csv") |> DataFrame
coverage = CSV.read("../test/austin-data/coverage_real.csv", DataFrame, header=false)
coverage = convert(Array{Bool, 2}, coverage[:, :])
incidents = CSV.File("../test/austin-data/austin_incidents.csv") |> DataFrame

regions = Int[parse(Int,string(x)) for x in names(hourly_calls[:,6:ncol(hourly_calls)])]
locations = collect(1:size(coverage,2))
adjacent = convert(Array, adjacent_nbhd[:,2:ncol(adjacent_nbhd)])[regions,regions] .> 0.5
demand = convert(Array,hourly_calls[:,6:end]);

incidents = incidents[.~ismissing.(incidents[!,:stn1_min]), :] # drop 44 calls that were "unreachable" (because all stations are reachable from each other)
incidents[!,:interarrival_seconds] = [0; incidents[!,:arrival_seconds][2:end] - incidents[!,:arrival_seconds][1:end-1]]
incidents[!,:isweekday] .= true
incidents[!,:isweekday][incidents[!,:dow] .== "Sun"] .= false
incidents[!,:isweekday][incidents[!,:dow] .== "Sat"] .= false;
incidents = incidents[incidents[!,:isweekday],:]

regions2index = Dict{Int,Int}(regions[i]=>i for i in 1:length(regions))

#= previously there was an error that didnt account for all the regions because it
 it only included regions for incidents. the regions without incidents are changed
 to 0=#
for x in incidents[!,:neighborhood]
    if(!haskey(regions2index,x))
        regions2index[x] = 0
    end
end
incidents[!,:neighborhood] = [regions2index[x] for x in incidents[!,:neighborhood]];

calls = incidents[:,[[:hour,:dow,:month,:year,:neighborhood,:interarrival_seconds];
                     Symbol[Symbol("stn$(i)_min") for i in locations]]]

DataFrames.first(calls, 6)

# We focus on emergency calls during the "peak period" (8AM - 8PM),
# with the emergency calls from the first 3 month as our training set,
# and the subsequent emergency calls from the remaining months as our test set
peak_period = (hourly_calls[!,:hour] .>= 8) .* (hourly_calls[!,:hour] .<= 20)
indices = 1:DataFrames.nrow(hourly_calls);
train_filter = (hourly_calls[!,:year] .== 2019) .* (hourly_calls[!,:month] .<= 3)
test_filter  = .~train_filter;
train_indices = indices[train_filter]
test_indices = indices[test_filter];

# Same as for the hourly calls; but this is for individual emergency calls
inc_peak_period = (calls[!,:hour] .>= 8) .* (calls[!,:hour] .<= 20)
inc_indices = 1:DataFrames.nrow(calls);

inc_train_filter = (calls[!,:year] .== 2019) .* (calls[!,:month] .<= 3)
inc_test_filter  = .~inc_train_filter

inc_train_indices = inc_indices[inc_train_filter]
inc_test_indices = inc_indices[inc_test_filter];

# we distinguish between peak and offpeak hours
train_peak = indices[peak_period .* train_filter]
train_offpeak = indices[.~peak_period .* train_filter]

test_peak = indices[peak_period .* test_filter]
test_offpeak = indices[.~peak_period .* test_filter]

p = DeploymentProblem(30, length(locations), length(regions), demand, train_indices,
                      test_indices, coverage[regions,:], Array{Bool,2}(adjacent));

test_inc_peak = inc_peak_period .* inc_test_filter
test_inc_offpeak = .~inc_peak_period .* inc_test_filter;

## make sure to include @time begin
    scenarios = Dict{Symbol, Dict{Int, Vector{Vector{Int}}}}()
    generated_deployment = Dict{Symbol, Dict{Int, Vector{Vector{Int}}}}()
    upperbounds = Dict{Symbol, Dict{Int, Vector{Float64}}}()
    lowerbounds = Dict{Symbol, Dict{Int, Vector{Float64}}}()
    upptiming = Dict{Symbol, Dict{Int, Vector{Float64}}}()
    lowtiming = Dict{Symbol, Dict{Int, Vector{Float64}}}()
    amb_deployment = Dict{Symbol, Dict{Int, Vector{Int}}}()
    for(deployment_model, name) in ((dp -> RobustDeployment(dp, α=0.0000001), :RobustCertain),
                                      (dp -> RobustDeployment(dp, α=0.5), :RobustMiddle),
                                      (dp -> RobustDeployment(dp, α=0.01), :RobustDesired),
                                      (dp -> RobustDeployment(dp, α=0.99999), :RobustUncertain))
        println("$name: ")
        amb_deployment[name] = Dict{Int, Vector{Int}}()
        scenarios[name] = Dict{Int, Vector{Vector{Int}}}()
        generated_deployment[name] = Dict{Int, Vector{Vector{Int}}}()
        upperbounds[name] = Dict{Int, Vector{Float64}}()
        lowerbounds[name] = Dict{Int, Vector{Float64}}()
        upptiming[name] = Dict{Int, Vector{Float64}}()
        lowtiming[name] = Dict{Int, Vector{Float64}}()
        for namb in 30:5:50
            println("$namb ")
            p.nambulances = namb
            model = deployment_model(p)
            set_optimizer(model.m, Gurobi.Optimizer)
            @time optimize!(model, p)
            amb_deployment[name][namb] = deployment(model)

            # for tracking purposes
            scenarios[name][namb] = model.scenarios
            generated_deployment[name][namb] = model.deployment
            upperbounds[name][namb] = model.upperbounds
            lowerbounds[name][namb] = model.lowerbounds
            upptiming[name][namb] = model.upptiming
            lowtiming[name][namb] = model.lowtiming
        end
        println
    end

    JLD.jldopen("outputs/austin_robust_test.jld", "w") do file
        write(file, "amb_deployment", amb_deployment)
        write(file, "scenarios", scenarios)
        write(file, "generated_deployment", generated_deployment)
        write(file, "upperbounds", upperbounds)
        write(file, "lowerbounds", lowerbounds)
        write(file, "upptiming", upptiming)
        write(file, "lowtiming", lowtiming)
    end

    solverstats = JLD.load("outputs/austin_robust_test.jld")
    amb_deployment = solverstats["amb_deployment"]
    json_string = JSON.json(amb_deployment)
    open("../src/outputs/austin_robust_test.json","w") do f
        write(f, json_string)
    end

max_stars_repo_name: michaelhilborn/AmbulanceDeployement

max_issues_repo_name: joshua-ong/AmbulanceDeployment

max_forks_repo_name: joshua-ong/AmbulanceDeployment


text:
"Best Effort Ranking"
module BestEffortRanking

using ..AExpressions

export get_fitness

"""p is the synthesized program
e is the list of examples to test the model on(could be a starting frame and
then a number of iterations/sequence of events and then ending frame and would
say is correct if the ending frames match)
a is the aexpr representation of the program
the closer to 0 the value the better the program"""
function get_fitness(p, e, a, distance_coefficient=3.0, size_coefficient=1.0)
  return (distance_coefficient * get_distance(p, e) +
  size_coefficient * sizeA(a))
end

"""for examples that the model gets wrong find the objects in the expected result
and take the euclidean distance of each object.
return the average of (sum of euclidean distances)/(max possible sum of euclidean distances)"""
function get_distance(p, e)
  max_distance = sqrt(15*15 + 15*15)
  overall_distance = 0
  max_possible_distance = 0

  for index = 1:length(e)
    expected, steps = e[index]
    state = p.init(steps[1]...)
    for i = 2:length(steps)
      state = p.next(state, steps[i]...)
    end
    actual_cells = p.render(state.scene)
    expected_cells = expected(p)

    overall_distance += compare_states(expected_cells, actual_cells)
    if length(actual_cells) > length(expected_cells)
      longer = length(actual_cells)
    else
      longer = length(expected_cells)
    end
    max_possible_distance += longer * max_distance
  end
  overall_distance/max_possible_distance
end

function compare_states(expected_cells, actual_cells)::Float64
  max_distance = sqrt(15*15 + 15*15)
  for cells in [expected_cells, actual_cells]
    sort(cells, by = x -> x.position.y)
    sort!(cells, by=x -> x.position.x)
    sort!(cells, by=x -> x.color)
  end
  index_expected = 1
  index_actual = 1
  distance = 0
  while index_expected <= length(expected_cells) && index_actual <= length(actual_cells)
    ex_cell = expected_cells[index_expected]
    ac_cell = actual_cells[index_actual]

    if ex_cell.color == ac_cell.color
      distance += euclidean_distance(ex_cell.position.x, ex_cell.position.y, ac_cell.position.x, ac_cell.position.y)
      index_expected += 1
      index_actual += 1
    elseif ex_cell.color > ac_cell.color
      distance += max_distance
      index_actual += 1
    else
      distance += max_distance
      index_expected += 1
    end
  end
  extra_cells = length(actual_cells) - index_actual + length(expected_cells) - index_expected + 2

  distance += abs(extra_cells) * max_distance
  distance
end

function euclidean_distance(x1, y1, x2, y2)
  sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2))
end
end


text:
module CrossValidation

    using ..ModelSelection
    using Random, Statistics, Printf
    import Base: iterate

    include("const.jl")
    include("utils.jl")
    include("strings.jl")
    include("structs/result.jl")
    include("core.jl")

    export kfoldcrossvalidation, CROSSVALIDATION_EXTRAKEY
end


text:
module DashSwitch

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict()

textures = String["default", "mirror"]
directions = Dict{String, Tuple{Type, String, Bool}}(
    "up" => (Maple.DashSwitchVertical, "ceiling", false),
    "down" => (Maple.DashSwitchVertical, "ceiling", true),
    "left" => (Maple.DashSwitchHorizontal, "leftSide", false),
    "right" => (Maple.DashSwitchHorizontal, "leftSide", true),
)

for texture in textures
    for (dir, data) in directions
        key = "Dash Switch ($(uppercasefirst(dir)), $(uppercasefirst(texture)))"
        func, datakey, val = data
        placements[key] = Ahorn.EntityPlacement(
            func,
            "rectangle",
            Dict{String, Any}(
                "sprite" => texture,
                datakey => val
            )
        )
    end
end

function Ahorn.selection(entity::Maple.DashSwitchHorizontal)
    x, y = Ahorn.position(entity)
    left = get(entity.data, "leftSide", false)

    if left
        return Ahorn.Rectangle(x, y - 1, 10, 16)

    else
        return Ahorn.Rectangle(x - 2, y, 10, 16)
    end
end

function Ahorn.selection(entity::Maple.DashSwitchVertical)
    x, y = Ahorn.position(entity)
    ceiling = get(entity.data, "ceiling", false)

    if ceiling
        return Ahorn.Rectangle(x, y, 16, 12)

    else
        return Ahorn.Rectangle(x, y - 4, 16, 12)
    end
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.DashSwitchHorizontal, room::Maple.Room)
    sprite = get(entity.data, "sprite", "default")
    left = get(entity.data, "leftSide", false)
    texture = sprite == "default" ? "objects/temple/dashButton00.png" : "objects/temple/dashButtonMirror00.png"

    if left
        Ahorn.drawSprite(ctx, texture, 20, 25, rot=pi)

    else
        Ahorn.drawSprite(ctx, texture, 8, 8)
    end
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.DashSwitchVertical, room::Maple.Room)
    sprite = get(entity.data, "sprite", "default")
    ceiling = get(entity.data, "ceiling", false)
    texture = sprite == "default" ? "objects/temple/dashButton00.png" : "objects/temple/dashButtonMirror00.png"

    if ceiling
        Ahorn.drawSprite(ctx, texture, 9, 20, rot=-pi / 2)

    else
        Ahorn.drawSprite(ctx, texture, 27, 7, rot=pi / 2)
    end
end

end


text:
using Distributed, StatsBase, Distributions, PyPlot
using TransitionalMCMC

lb  = -15
ub  = 15

prior(x) = pdf(Uniform(lb,ub), x[1,:]) .* pdf(Uniform(lb,ub), x[2,:])
priorRnd(Nsamples) = rand(Uniform(lb,ub), 2, Nsamples)

function LogLik(x)
    return log.(pdf(MvNormal([0,0],[1 -0.5; -0.5 1]), x))
end

Nsamples = 2000

samps, acc = tmcmc(LogLik, prior, priorRnd, Nsamples)

plt.scatter(samps[1,:], samps[2,:])


text:
#############################################################
# Making sure embedding is invariant under forward linear map
# of its vertices
#############################################################
import Distributions
import Statistics
import StaticArrays: SVector, MVector
import DelayEmbeddings: Dataset

"""
    invariantize(pts; noise_factor = 0.01, step = 5)

Checks whether the last point in `pts` lies inside the convex hull of the
previus points. If not, the last point is iteratively moved towards the
center of the convex hull until it reaches the interior of the convex hull (`step` percent
at a time).

This functions adds a small amount of noise to the data, equivalent to `noise_factor` times
the standard deviation for each dynamical variable.

Valid inputs are `AbstractArray{T, 2}` (the largest dimension will be treated as the index),
`Vector{Vector{T}}`, `Vector{SVector{D, T}}`, `Vector{MVector{D, T}}` (each vector is a
point) or a `Dataset` instance.
"""
function invariantize end

"""
    forwardlinearmap_invariant(pts)

Return true if the last point in `pts` lies inside the convex hull of the preceding points.

Valid inputs are `AbstractArray{T, 2}` (the largest dimension will be treated as the index),
`Vector{Vector{T}}`, `Vector{SVector{D, T}}`, `Vector{MVector{D, T}}` (each vector is a
point) or a `Dataset` instance.
"""
function forwardlinearmap_invariant end

function forwardlinearmap_invariant(pts::Vector{SVector{D, T}}) where {D, T}
    forwardlinearmap_invariant(Array(hcat(pts...,)))
end

function forwardlinearmap_invariant(pts::Vector{MVector{D, T}}) where {D, T}
    forwardlinearmap_invariant(Array(hcat(pts...,)))
end

function forwardlinearmap_invariant(pts::Dataset) where {D, T}
    forwardlinearmap_invariant(Array(hcat(pts...,)))
end

function forwardlinearmap_invariant(pts::Vector{Vector{T}}) where T
    forwardlinearmap_invariant(hcat(pts...,))
end

function forwardlinearmap_invariant(pts::AbstractArray{T, 2}) where T

    # Transpose if rows represent points
    if size(pts, 1) > size(pts, 2)
       pts = transpose(pts)
    end

    lastpoint = pts[:, end]
    dim = size(pts, 1)
    # Triangulate the embedding using all points but the last
    # Returns a vector of index vectors, one for each simplex
    t = delaunay(transpose(pts[:, 1:end-1]))
    n_simplices = length(t)
    points = pts[:, 1:end-1]

    # Loop over valid simplices and check whether the corresponding simplex
    # contains the last point.
    i = 1
    lastpoint_contained = false

    #=
    Continue checking if simplices contain the last point until some simplex
    contains it. Then the point must necessarily be inside the convex hull of
    the triangulation formed by those simplices.
    =#
    while i <= n_simplices && !lastpoint_contained
        sᵢ = transpose(points[:, t[i]])
        orientation_sᵢ = Statistics.det([ones(dim + 1, 1) sᵢ])

        # Insert the last point of the embedding in place
        # of the j-th vertex of sᵢ and check whether the
        # convex expansion coefficient β𝚥 stays non-negative.
        # If β𝚥 becomes negative, the point is not contained.
        j = 1
        β𝚥 = 1
        while j <= dim + 1 && β𝚥 >= 0
            tmp = copy(sᵢ)
            tmp[j, :] = lastpoint
            β𝚥 = Statistics.det(hcat(ones(dim + 1, 1), tmp)) * sign(orientation_sᵢ)
            j = j + 1
        end

        #=
        If the last convex expansion coefficient is positive, the last point is
        contained in the triangulation (because all the previous coefficients
        must have been nonnegative)
        =#
        if β𝚥 >= 0
            lastpoint_contained = true
        end

        i = i + 1
    end
    return lastpoint_contained
end

forwardlinearmap_invariant(E::AbstractEmbedding) = forwardlinearmap_invariant(E.points)

function invariantize(E::AbstractEmbedding{D, T};
                        noise_factor = 0.01, step = 3) where {D, T}
    pts = E.points
    dim = size(E.points, 1)

    if size(unique(pts, dims = 2), 2) < size(pts, 2)

       @warn """Embedding points not unique. Adding a little noise ($noise_factor times the maximum of the the standard deviations along each axis)"""
       # Find standard deviation along each axis
       dim = size(pts, 1)
       σ = Statistics.std(E.points, dims = 2)

       for i = 1:dim
          pts[i, :] .+= rand(Distributions.Uniform(-σ[i], σ[i])) .* noise_factor
       end
    end

    #=
    # Keep track of the embedding's center point and the original position of the
    # last point in the embedding, so we can move the last point along a line
    # from its original position towards the embedding's center, until the point
    # lies inside the convex hull of the preceding points.
    =#
    ce = sum(E.points, dims = 2)/size(E.points, 2) # embedding center
    lp = E.points[:, end] # last point of the embedding
    # What direction should we move?
    dir = ce - lp

    dir = dropdims(ce .- lp, dims = 2)

    # Points along line toward the center of the embedding.
    steps = 1:step:100
    ptsonline = [lp .+ dir .* (pct_moved/100) for pct_moved in 1:step:100]

    for i = 1:length(ptsonline)
       pt = ptsonline[i]
       P = hcat(pts[:, 1:(end - 1)], pt)
       if forwardlinearmap_invariant(P)

          embeddingdata = EmbeddingData{dim, T}(
               float.(E.embeddingdata.dataseries),  # the dataseries
               E.embeddingdata.labels, # empty labels by default
               E.embeddingdata.descriptions, # empty descriptions by default
               E.embeddingdata.in_which_pos, # the positions in which the dataseries appear
               E.embeddingdata.at_what_lags # the embedding lag for each column
               )

          return LinearlyInvariantEmbedding(
              hcat(pts[:, 1:(end-1)], pt), embeddingdata)
       end
    end
    @warn """Could not make embedding invariant. Returning unmodified $E."""
    return E
end

function invariantize(pts::AbstractArray{Float64, 2};
                        noise_factor = 0.01, step = 3)
   dim = size(pts, 1)

   if size(pts, 1) > size(pts, 2)
      error("`pts` must be an array of size dim-by-n_points, where each column is a point")
   end

   if size(unique(pts, dims = 2), 2) < size(pts, 2)

      @warn """Points not unique. Adding a little noise ($noise_factor times the maximum of the the standard deviations along each axis)"""
      # Find standard deviation along each axis
      dim = size(pts, 1)
      σ = Statistics.std(pts, dims = 2)

      for i = 1:dim
         pts[i, :] .+= rand(Distributions.Uniform(-σ[i], σ[i])) .* noise_factor
      end
   end

   #=
   # Keep track of the embedding's center point and the original position of the
   # last point in the embedding, so we can move the last point along a line
   # from its original position towards the embedding's center, until the point
   # lies inside the convex hull of the preceding points.
   =#
   ce = sum(pts, dims = 2)/size(pts, 2) # embedding center
   lp = pts[:, end] # last point of the embedding
   # What direction should we move?
   dir = ce - lp

   dir = dropdims(ce .- lp, dims = 2)

   # Points along line toward the center of the embedding.
   ptsonline = [lp .+ dir .* (pct_moved/100) for pct_moved in 1:step:100]
   for pt in ptsonline
      P = hcat(pts[:, 1:(end - 1)], pt)
      if forwardlinearmap_invariant(P)
         return P
      end
   end

   @warn """Could not make point set invariant under forward linear map. Returning unmodified points."""
   return pts
end

function invariantize(pts::Vector{VT}; step = 3) where {VT <: Union{SVector{D, T}, MVector{D, T}, Vector{T}} where {D, T}}
   pts_arr = Array(hcat(pts...,))
   inv_pts = invariantize(pts_arr, step = step)
   return [VT(inv_pts[:, i]) for i = 1:length(pts)]
end

function invariantize(pts::Vector{Vector{T}}; step = 3) where T
   pts_arr = hcat(pts...,)
   inv_pts = invariantize(pts_arr, step = step)
   return [inv_pts[:, i] for i = 1:length(pts)]
end

function invariantize(pts::Dataset; step = 3)
   inv_pts = invariantize(hcat([Vector(x) for x in pts]...,), step = step)
   Dataset([inv_pts[:, i] for i = 1:length(pts)])
end

export
forwardlinearmap_invariant,
invariantize


text:
# ---
# title: 18. 4Sum
# id: problem18
# author: zhwang
# date: 2022-03-15
# difficulty: Medium
# categories: Array, Hash Table, Two Pointers
# link: <https://leetcode.com/problems/4sum/description/>
# hidden: true
# ---
#
# Given an array `nums` of _n_ integers and an integer `target`, are there
# elements _a_ , _b_ , _c_ , and _d_ in `nums` such that _a_ \+ _b_ \+ _c_ \+
# _d_ = `target`? Find all unique quadruplets in the array which gives the sum
# of `target`.
#
# **Notice**  that the solution set must not contain duplicate quadruplets.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,0,-1,0,-2,2], target = 0
#     Output: [[-2,-1,1,2],[-2,0,0,2],[-1,0,0,1]]
#
#
# **Example 2:**
#
#
#
#     Input: nums = [], target = 0
#     Output: []
#
#
#
#
# **Constraints:**
#
#   * `0 <= nums.length <= 200`
#   * `-109 <= nums[i] <= 109`
#   * `-109 <= target <= 109`
#
#
## @lc code=start
using LeetCode

four_sum(nums::Vector{Int}, target::Int) = four_sum!(copy(nums), target)
function four_sum!(nums::Vector{Int}, target::Int)::Vector{NTuple{4,Int}}
    n, res = length(sort!(nums)), NTuple{4,Int}[]
    n <= 3 && return NTuple{4,Int}[]
    for (i, first) in enumerate(@view(nums[1:(end - 3)]))
        i > 1 && first == nums[i - 1] && continue ## avoid repetition
        sum(@view(nums[i:(i + 3)])) > target && return res ## pruning
        first + sum(@view(nums[(end - 3):end])) < target && continue ## pruning
        for j in (i + 1):(n - 2)
            j > i + 1 && nums[j] == nums[j - 1] && continue ## avoid repetition
            first + sum(@view(nums[j:(j + 2)])) > target && break ## pruning
            first + nums[j] + nums[end - 1] + nums[end] < target && continue ## pruning
            newtarget = target - first - nums[j]
            left, right = j + 1, n
            while left < right
                newsum = nums[left] + nums[right]
                if newsum < newtarget
                    left += 1
                elseif newsum > newtarget
                    right -= 1
                else
                    push!(res, (first, nums[j], nums[left], nums[right]))
                    pos = findfirst(!=(nums[left]), @view(nums[(left + 1):right]))
                    isnothing(pos) ? break : (left += pos)
                    pos = findfirst(!=(nums[right]), @view(nums[(right - 1):-1:left]))
                    isnothing(pos) ? break : (right -= pos)
                end
            end
        end
    end
    return res
end
## @lc code=end


text:
include("src/A3.jl")
using .A3

λₘₐₓ, λₘᵢₙ, λ = solve(
    Problem(
        function(i, j)
            return (-1)^(i+j) / (abs(i-j) + 1)
        end,
    ),
)

open("result", "w") do io
    println(io, "λₘₐₓ = ", λₘₐₓ)
    println(io, "λₘᵢₙ = ", λₘᵢₙ, "\n")
    println(io, "λ: ")
    for e in λ
        println(io, e)
    end
end


text:
"""
    InfiniteOpt.add_measure_variable(model::JuMP.Model,
                                     var::InfiniteOpt.PointVariable,
                                     key::Val{:TransData}
                                     )::InfiniteOpt.GeneralVariableRef

Make a `PointVariableRef` and map it to the appropriate transcription variable
and return the `GeneralVariableRef`. This is an extension of
[`add_measure_variable`](@ref InfiniteOpt.add_measure_variable(::JuMP.Model,::Any,::Any))
for `TranscriptionOpt`.
"""
function InfiniteOpt.add_measure_variable(model::JuMP.Model,
                                          var::InfiniteOpt.PointVariable,
                                          key::Val{:TransData}
                                          )::InfiniteOpt.GeneralVariableRef
    # make negative index to not conflict with the InfiniteModel
    raw_index = transcription_data(model).last_point_index -= 1
    # make the reference and map it to a transcription variable
    ivref = var.infinite_variable_ref
    pvref = InfiniteOpt.GeneralVariableRef(JuMP.owner_model(ivref), raw_index,
                                           InfiniteOpt.PointVariableIndex)
    trans_var = lookup_by_support(model, ivref, var.parameter_values)
    transcription_data(model).finvar_mappings[pvref] = trans_var
    return pvref
end

"""
    InfiniteOpt.add_measure_variable(model::JuMP.Model,
                                     var::InfiniteOpt.ReducedVariable,
                                     key::Val{:TransData}
                                     )::InfiniteOpt.GeneralVariableRef

Make a `ReducedVariableRef` and add `var` to the transcription data
and return the `GeneralVariableRef`. This is an extension of
[`add_measure_variable`](@ref InfiniteOpt.add_measure_variable(::JuMP.Model,::Any,::Any))
for `TranscriptionOpt`. Note that `internal_reduced_variable` is also extended
to be able to access the `var`.
"""
function InfiniteOpt.add_measure_variable(model::JuMP.Model,
                                          var::InfiniteOpt.ReducedVariable,
                                          key::Val{:TransData}
                                          )::InfiniteOpt.GeneralVariableRef
    # make negative index to not conflict with the InfiniteModel
    reduced_vars = transcription_data(model).reduced_vars
    raw_index = -1 * (length(reduced_vars) + 1)
    # make the reference and map it to a transcription variable
    ivref = var.infinite_variable_ref
    rvref = InfiniteOpt.GeneralVariableRef(JuMP.owner_model(ivref), raw_index,
                                           InfiniteOpt.ReducedVariableIndex)
    push!(reduced_vars, var)
    _set_reduced_variable_mapping(model, var, rvref, InfiniteOpt._index_type(ivref))
    return rvref
end

"""
    InfiniteOpt.delete_reduced_variable(model::JuMP.Model,
                                        vref::InfiniteOpt.ReducedVariableRef,
                                        key::Val{:TransData})::Nothing

This is an extension of
[`delete_reduced_variable`](@ref InfiniteOpt.delete_reduced_variable(::JuMP.Model, ::Any, ::Any))
for use in `TranscriptionOpt`. Here we do not delete reduced variables once they
have been used since there is no performance gain for this paradigm and the
memory saving is small. Note this may change in the future.
"""
function InfiniteOpt.delete_reduced_variable(model::JuMP.Model,
                                             rvref::InfiniteOpt.ReducedVariableRef,
                                             key::Val{:TransData})::Nothing
    return
end


text:
function grid_search(f, grid)
    min_value = f(grid[1])
    arg_min_value = grid[1]

    for el in grid
        if f(el) < min_value
            arg_min_value = el
        end
    end

    return arg_min_value
end


text:
"""
    mae(ŷ, y; agg=mean)

Return the loss corresponding to mean absolute error:

    agg(abs.(ŷ .- y))
"""
mae(ŷ, y; agg=mean) = agg(abs.(ŷ .- y))

"""
    mse(ŷ, y; agg=mean)

Return the loss corresponding to mean square error:

    agg((ŷ .- y).^2)
"""
mse(ŷ, y; agg=mean) = agg((ŷ .- y).^2)

"""
    msle(ŷ, y; agg=mean, ϵ=eps(ŷ))

The loss corresponding to mean squared logarithmic errors, calculated as

    agg((log.(ŷ .+ ϵ) .- log.(y .+ ϵ)).^2)

The `ϵ` term provides numerical stability.
Penalizes an under-estimation more than an over-estimatation.
"""
msle(ŷ, y; agg=mean, ϵ=epseltype(ŷ)) = agg((log.((ŷ .+ ϵ) ./ (y .+ ϵ))).^2)

"""
    huber_loss(ŷ, y; δ=1, agg=mean)

Return the mean of the [Huber loss](https://en.wikipedia.org/wiki/Huber_loss)
given the prediction `ŷ` and true values `y`.

                 | 0.5 * |ŷ - y|^2,            for |ŷ - y| <= δ
    Huber loss = |
                 |  δ * (|ŷ - y| - 0.5 * δ), otherwise
"""
function huber_loss(ŷ, y; agg=mean, δ=ofeltype(ŷ, 1))
   abs_error = abs.(ŷ .- y)
   #TODO: remove dropgrad when Zygote can handle this function with CuArrays
   temp = Zygote.dropgrad(abs_error .<  δ)
   x = ofeltype(ŷ, 0.5)
   agg(((abs_error.^2) .* temp) .* x .+ δ*(abs_error .- x*δ) .* (1 .- temp))
end

"""
    label_smoothing(y::Union{Number, AbstractArray}, α; dims::Int=1)

Returns smoothed labels, meaning the confidence on label values are relaxed.

When `y` is given as one-hot vector or batch of one-hot, its calculated as

    y .* (1 - α) .+ α / size(y, dims)

when `y` is given as a number or batch of numbers for binary classification,
its calculated as

    y .* (1 - α) .+ α / 2

in which case the labels are squeezed towards `0.5`.

α is a number in interval (0, 1) called the smoothing factor. Higher the
value of α larger the smoothing of `y`.

`dims` denotes the one-hot dimension, unless `dims=0` which denotes the application
of label smoothing to binary distributions encoded in a single number.

Usage example:

    sf = 0.1
    y = onehotbatch([1, 1, 1, 0, 0], 0:1)
    y_smoothed = label_smoothing(ya, 2sf)
    y_sim = y .* (1-2sf) .+ sf
    y_dis = copy(y_sim)
    y_dis[1,:], y_dis[2,:] = y_dis[2,:], y_dis[1,:]
    @assert crossentropy(y_sim, y) < crossentropy(y_sim, y_smoothed)
    @assert crossentropy(y_dis, y) > crossentropy(y_dis, y_smoothed)
"""
function label_smoothing(y::Union{AbstractArray,Number}, α::Number; dims::Int=1)
    if !(0 < α < 1)
        throw(ArgumentError("α must be between 0 and 1"))
    end
    if dims == 0
        y_smoothed = y .* (1 - α) .+ α*1//2
    elseif dims == 1
        y_smoothed = y .* (1 - α) .+ α* 1 // size(y, 1)
    else
        throw(ArgumentError("`dims` should be either 0 or 1"))
    end
    return y_smoothed
end

"""
    crossentropy(ŷ, y; dims=1, ϵ=eps(ŷ), agg=mean)

Return the cross entropy between the given probability distributions;
calculated as

    agg(-sum(y .* log.(ŷ .+ ϵ); dims=dims))

Cross entropy is typically used as a loss in multi-class classification,
in which case the labels `y` are given in a one-hot format.
`dims` specifies the dimension (or the dimensions) containing the class probabilities.
The prediction `ŷ` is supposed to sum to one across `dims`,
as would be the case with the output of a [`softmax`](@ref) operation.

Use [`label_smoothing`](@ref) to smooth the true labels as preprocessing before
computing the loss.

Use of [`logitcrossentropy`](@ref) is recomended over `crossentropy` for
numerical stability.

See also: [`logitcrossentropy`](@ref), [`binarycrossentropy`](@ref), [`logitbinarycrossentropy`](@ref),
[`label_smoothing`](@ref)
"""
function crossentropy(ŷ, y; dims=1, agg=mean, ϵ=epseltype(ŷ))
    agg(.-sum(xlogy.(y, ŷ .+ ϵ); dims=dims))
end

"""
    logitcrossentropy(ŷ, y; dims=1, agg=mean)

Return the crossentropy computed after a [`logsoftmax`](@ref) operation;
calculated as

    agg(.-sum(y .* logsoftmax(ŷ; dims=dims); dims=dims))

Use [`label_smoothing`](@ref) to smooth the true labels as preprocessing before
computing the loss.

`logitcrossentropy(ŷ, y)` is mathematically equivalent to
[`crossentropy(softmax(ŷ), y)`](@ref) but it is more numerically stable.

See also: [`crossentropy`](@ref), [`binarycrossentropy`](@ref), [`logitbinarycrossentropy`](@ref), [`label_smoothing`](@ref)
"""
function logitcrossentropy(ŷ, y; dims=1, agg=mean)
    agg(.-sum(y .* logsoftmax(ŷ; dims=dims); dims=dims))
end

"""
    binarycrossentropy(ŷ, y; agg=mean, ϵ=eps(ŷ))

Return the binary cross-entropy loss, computed as

    agg(@.(-y*log(ŷ + ϵ) - (1-y)*log(1-ŷ + ϵ)))

The `ϵ` term provides numerical stability.

Typically, the prediction `ŷ` is given by the output of a [`sigmoid`](@ref) activation.

Use [`label_smoothing`](@ref) to smooth the `y` value as preprocessing before
computing the loss.

Use of `logitbinarycrossentropy` is recomended over `binarycrossentropy` for numerical stability.

See also: [`crossentropy`](@ref), [`logitcrossentropy`](@ref), [`logitbinarycrossentropy`](@ref),
[`label_smoothing`](@ref)
"""
function binarycrossentropy(ŷ, y; agg=mean, ϵ=epseltype(ŷ))
    agg(@.(-xlogy(y, ŷ+ϵ) - xlogy(1-y, 1-ŷ+ϵ)))
end
# Re-definition to fix interaction with CuArrays.
# CUDA.@cufunc binarycrossentropy(ŷ, y; ϵ=eps(ŷ)) = -y*log(ŷ + ϵ) - (1 - y)*log(1 - ŷ + ϵ)

"""
    logitbinarycrossentropy(ŷ, y; agg=mean)

Mathematically equivalent to
[`binarycrossentropy(σ(ŷ), y)`](@ref) but is more numerically stable.

Use [`label_smoothing`](@ref) to smooth the `y` value as preprocessing before
computing the loss.

See also: [`crossentropy`](@ref), [`logitcrossentropy`](@ref), [`binarycrossentropy`](@ref), [`label_smoothing`](@ref)
"""
function logitbinarycrossentropy(ŷ, y; agg=mean)
    agg(@.((1-y)*ŷ - logσ(ŷ)))
end
# Re-definition to fix interaction with CuArrays.
# CUDA.@cufunc logitbinarycrossentropy(ŷ, y) = (1 - y)*ŷ - logσ(ŷ)

"""
    kldivergence(ŷ, y; agg=mean)

Return the
[Kullback-Leibler divergence](https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence)
between the given probability distributions.

KL divergence is a measure of how much one probability distribution is different
from the other.
It is always non-negative and zero only when both the distributions are equal
everywhere.
"""
function kldivergence(ŷ, y; dims=1, agg=mean, ϵ=epseltype(ŷ))
  entropy = agg(sum(xlogx.(y), dims=dims))
  cross_entropy = crossentropy(ŷ, y; dims=dims, agg=agg, ϵ=ϵ)
  return entropy + cross_entropy
end

"""
    poisson_loss(ŷ, y)

# Return how much the predicted distribution `ŷ` diverges from the expected Poisson
# distribution `y`; calculated as `sum(ŷ .- y .* log.(ŷ)) / size(y, 2)`.
REDO
[More information.](https://peltarion.com/knowledge-center/documentation/modeling-view/build-an-ai-model/loss-functions/poisson_loss).
"""
poisson_loss(ŷ, y; agg=mean) = agg(ŷ .- xlogy.(y, ŷ))

"""
    hinge_loss(ŷ, y; agg=mean)

Return the [hinge_loss loss](https://en.wikipedia.org/wiki/Hinge_loss) given the
prediction `ŷ` and true labels `y` (containing 1 or -1); calculated as
`sum(max.(0, 1 .- ŷ .* y)) / size(y, 2)`.

See also: [`squared_hinge_loss`](@ref)
"""
hinge_loss(ŷ, y; agg=mean) = agg(max.(0, 1 .-  ŷ .* y))

"""
    squared_hinge_loss(ŷ, y)

Return the squared hinge_loss loss given the prediction `ŷ` and true labels `y`
(containing 1 or -1); calculated as `sum((max.(0, 1 .- ŷ .* y)).^2) / size(y, 2)`.

See also: [`hinge_loss`](@ref)
"""
squared_hinge_loss(ŷ, y; agg=mean) = agg((max.(0, 1 .- ŷ .* y)).^2)

"""
    dice_coeff_loss(ŷ, y; smooth=1)

Return a loss based on the dice coefficient.
Used in the [V-Net](https://arxiv.org/abs/1606.04797) image segmentation
architecture.
Similar to the F1_score. Calculated as:

    1 - 2*sum(|ŷ .* y| + smooth) / (sum(ŷ.^2) + sum(y.^2) + smooth)
"""
dice_coeff_loss(ŷ, y; smooth=ofeltype(ŷ, 1.0)) = 1 - (2*sum(y .* ŷ) + smooth) / (sum(y.^2) + sum(ŷ.^2) + smooth) #TODO agg

"""
    tversky_loss(ŷ, y; β=0.7)

Return the [Tversky loss](https://arxiv.org/abs/1706.05721).
Used with imbalanced data to give more weight to false negatives.
Larger β weigh recall more than precision (by placing more emphasis on false negatives)
Calculated as:
    1 - sum(|y .* ŷ| + 1) / (sum(y .* ŷ + β*(1 .- y) .* ŷ + (1 - β)*y .* (1 .- ŷ)) + 1)
"""
function tversky_loss(ŷ, y; β=ofeltype(ŷ, 0.7))
    #TODO add agg
    num = sum(y .* ŷ) + 1
    den = sum(y .* ŷ + β*(1 .- y) .* ŷ + (1 - β)*y .* (1 .- ŷ)) + 1
    1 - num / den
end


text:
module WallTimeProgress

using ArgCheck: @argcheck

export WallTimeTracker

"""
    WallTimeTracker(period; output = stdout, count = 0)

Create an object that can be used to track progress through an *ex ante* unknown
number of records.

# Arguments

- `period`: incrementing the counter will display progress information after
  every `period`

- `output`: a stream for displaying progress information

- `count`: the initial count

# Usage

```julia
w = WallTimeTracker(10)
for _ in 1:100
   WallTimeProgress.increment!(w) # will print updates every 10 steps
end
```

See [`increment!`](@ref), also [`reset!`].

`Base.count` can be used to query the current count.
"""
mutable struct WallTimeTracker{Tcount <: Integer,
                               Toutput <: IO}
    period::Tcount
    output::Toutput
    count::Tcount
    next_report::Tcount
    last_count::Tcount
    last_time_ns::UInt64
    function WallTimeTracker(period::Tcount; output::Toutput = stdout,
                             count::Tcount = 0) where {Tcount, Toutput}
        @argcheck period > 0
        new{Tcount, Toutput}(period, output, 0,
                             next_period(count, period), 0, time_ns())
    end
end

"""
    next_period(count, period)

Return the smallest integer that is

1. strictly greater than `count`,

2. an integer multiple of `period`.

NOTE: we could extend `Base.ceil`, but that would be type piracy.

NOTE: for internal use, not exported.
"""
next_period(count::Integer, period::Integer) = (div(count, period) + 1) * period

function reset!(wtt::WallTimeTracker)
    wtt.last_count = wtt.count
    wtt.last_time_ns = time_ns()
    nothing
end

"""
    _with_underscores(x)

Return a string representing the integer `x`, with `_` as the thousands
separator.
"""
function _with_underscores(x::Integer)
    result = IOBuffer()
    ds = digits(x)
    n = length(ds)
    for (i, d) in enumerate(ds)
        write(result, '0' + d)
        i % 3 == 0 && i < n && write(result, '_')
    end
    String(reverse!(take!(result)))
end

function Base.show(io::IO, wtt::WallTimeTracker)
    if wtt.count > 0
        print(io, "processed $(_with_underscores(wtt.count)) records")
    else
        print(io, "(no records yet)")
    end
    Δ = (wtt.count - wtt.last_count)
    if Δ > 0
        avg = (time_ns() - wtt.last_time_ns) / (1e9 * Δ)
        print(io, ", $avg s/record")
    end
    println(io)
end

Base.count(wtt::WallTimeTracker) = wtt.count

function increment!(wtt::WallTimeTracker, step = 1)
    @argcheck step > 0
    wtt.count += step
    if wtt.count ≥ wtt.next_report
        show(wtt.output, wtt)
        reset!(wtt)
        wtt.next_report = next_period(wtt.count, wtt.period)
    end
    nothing
end

end # module


text:
mutable struct Future
    pool::BaseExecutor
    data::Union{Nothing, Some, Exception}
    continuation::Union{Nothing, Function}
    function Future(f::Future)
        f
    end

    function Future(pool::BaseExecutor, data::Any = nothing)
        new(pool, data, nothing)
    end
end

function isfulfilled(f::Future)
    f.data !== nothing
end

function _then!(f::Future, continuation::Function)
    @assert f.continuation === nothing

    f.continuation = continuation
    if isfulfilled(f)
        f.continuation(isa(f.data, Exception) ? f.data : something(f.data))
    end
end

finally!(f::Future, continuation::Function) = _then!(f, continuation)

function then!(f::Future, continuation::Function)
    f2 = Future(f.pool)
    _then!(f,
      x -> begin
        val = if isa(x, Exception)
          x
        else
          continuation(x)
        end

        fulfill!(f2, val)
      end
    )
    f2
end

function then!(f::Future, continuation::RemoteFunction)
    f2 = Future(f.pool)
    _then!(f,
      x -> begin
        if isa(x, Exception)
          fulfill!(f2, x)
        else
          rf = submit!(f.pool, continuation, x)
          _then!(rf, x -> fulfill!(f2, x))
        end
      end)
    f2
end

function fulfill!(f::Future, value::Any)
    @assert !isfulfilled(f)

    f.data = !isa(value, Exception) ? Some(value) : value
    if f.continuation !== nothing
        f.continuation(value)
    end
    f
end

function get!(f::Future)
    while !isfulfilled(f)
        run_until!(f.pool)
    end

    if isa(f.data, Exception)
      throw(f.data)
    end

    something(f.data)
end

function whenall!(futs::Vector{Future}, ::Type{T} = Any) where T
    n = length(futs)
    @assert n > 0

    results = Vector{T}(undef, n)
    r = Future(futs[1].pool)
    function post(i::Int, val::Any)
        if n > 0
          if !isa(val, Exception)
            results[i] = val
            n = n - 1

            if n == 0
              fulfill!(r, results)
            end
          else
            fulfill!(r, val)
            n = 0
          end
        end
    end

    for (i,f) in enumerate(futs)
        _then!(f, x -> post(i, x))
    end
    r
end


text:
mode = get(ENV, "MAKIE_COMPILE", "")

function is_ci()
    get(ENV, "TRAVIS", "") == "true" ||
    get(ENV, "APPVEYOR", "") == "true" ||
    get(ENV, "CI", "") == "true"
end

if is_ci()
    Pkg.clone("https://github.com/SimonDanisch/AbstractNumbers.jl.git")
    Pkg.checkout("GLAbstraction")
    Pkg.checkout("GeometryTypes")
    Pkg.checkout("GLVisualize")
    Pkg.checkout("MeshIO")
    Pkg.add("VisualRegressionTests")
    Pkg.checkout("GLWindow")
end

if isempty(mode)
    # we're not compiling, so we do a reference image test run
    include("visual_regression.jl")
else
    # when snoop compiling, we simply just execute all sample code directly
    include("all_samples.jl")
end


text:
function eachcons(xs::AbstractVector, n)
    (view(xs, i:i+n-oneunit(n)) for i ∈ firstindex(xs):lastindex(xs)-n+oneunit(n))
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
## Soft Evaluation
## ===============
Cassette.@context SoftExCtx
Cassette.metadatatype(::Type{<:SoftExCtx}, ::Type{Bool}) = Omega.SoftBool

soft(::typeof(>)) = softgt
soft(::typeof(>=)) = softgt
soft(::typeof(<)) = softlt
soft(::typeof(<=)) = softlt
soft(::typeof(==)) = softeq

function soften(f, ctx, args...)
  println(f, args)
  Cassette.tag(f(args...), ctx, soft(f)(args...))
end

function softboolop(f, ctx, args...)
  args_ = Cassette.untag.(args, ctx)
  tags = Cassette.metadata.(args, ctx)
  Cassette.tag(f(args_...), ctx, f(tags...))
end

Softable{T} = Union{Number, Array{T}} where T <: Number

Cassette.execute(ctx::SoftExCtx, ::typeof(Base.:>), x::Softable, y::Softable) = soften(>, ctx, x, y)
Cassette.execute(ctx::SoftExCtx, ::typeof(Base.:>=), x::Softable, y::Softable) = soften(>=, ctx, x, y)
Cassette.execute(ctx::SoftExCtx, ::typeof(Base.:<), x::Softable, y::Softable) = soften(<, ctx, x, y)
Cassette.execute(ctx::SoftExCtx, ::typeof(Base.:<=), x::Softable, y::Softable) = soften(<=, ctx, x, y)
# Cassette.execute(ctx::SoftExCtx, ::typeof(Base.:(==)), x::Softable, y::Softable) = soften(==, ctx, x, y)

Cassette.execute(ctx::SoftExCtx{T}, ::typeof(Base.:!), x::Cassette.Tagged{T}) where {T} = softboolop(!, ctx, x)
Cassette.execute(ctx::SoftExCtx{T}, ::typeof(Base.:&), x::Cassette.Tagged{T}, y::Cassette.Tagged{T}) where {T} = softboolop(&, ctx, x, y)
Cassette.execute(ctx::SoftExCtx{T}, ::typeof(Base.:|), x::Cassette.Tagged{T}, y::Cassette.Tagged{T}) where {T} = softboolop(|, ctx, x, y)
Cassette.execute(ctx::SoftExCtx{T}, ::typeof(Base.:⊻), x::Cassette.Tagged{T}, y::Cassette.Tagged{T}) where {T} = softboolop(⊻, ctx, x, y)

function softapply(f, args...)
  ctx = Cassette.withtagfor(SoftExCtx(), f)
  res = Cassette.overdub(ctx, f, args...)
  Cassette.metadata(res, ctx)
end


text:
module RenyiEntropyOfMagic

using ITensors
using ThreadsX
using Distributed
using Serialization
using LinearAlgebra
using ProgressLogging

export Consts, transfer_matrix, boundary_vector, expectation,
    sample_paulistring, renyi_entropy, exact_renyi_entropy,
    enumerate_paulistring

include("dmrg.jl")
include("entropy.jl")
include("cli.jl")

end


text:
"""
    SingleBridgeOptimizer{BT<:AbstractBridge, MT<:MOI.ModelLike, OT<:MOI.ModelLike} <: AbstractBridgeOptimizer

The `SingleBridgeOptimizer` bridges any constraint supported by the bridge `BT`.
This is in contrast with the [`LazyBridgeOptimizer`](@ref) which only bridges the constraints that are unsupported by the internal model, even if they are supported by one of its bridges.
"""
struct SingleBridgeOptimizer{BT<:AbstractBridge, MT<:MOI.ModelLike, OT<:MOI.ModelLike} <: AbstractBridgeOptimizer
    model::OT
    bridged::MT
    bridges::Dict{CI, BT}
end
function SingleBridgeOptimizer{BT, MT}(model::OT) where {BT, MT, OT <: MOI.ModelLike}
    SingleBridgeOptimizer{BT, MT, OT}(model, MT(), Dict{CI, BT}())
end

is_bridged(b::SingleBridgeOptimizer, ::Type{<:MOI.AbstractFunction}, ::Type{<:MOI.AbstractSet}) = false
bridge_type(b::SingleBridgeOptimizer{BT}, F::Type{<:MOI.AbstractFunction}, S::Type{<:MOI.AbstractSet}) where BT = BT

# :((Zeros, SecondOrderCone)) -> (:(MOI.Zeros), :(MOI.SecondOrderCone))

"""
macro bridge(modelname, bridge, scalarsets, typedscalarsets, vectorsets, typedvectorsets, scalarfunctions, typedscalarfunctions, vectorfunctions, typedvectorfunctions)

Creates a type `modelname` implementing the MOI model interface and bridging the `scalarsets` scalar sets `typedscalarsets` typed scalar sets, `vectorsets` vector sets, `typedvectorsets` typed vector sets, `scalarfunctions` scalar functions, `typedscalarfunctions` typed scalar functions, `vectorfunctions` vector functions and `typedvectorfunctions` typed vector functions.
To give no set/function, write `()`, to give one set `S`, write `(S,)`.

### Examples

The optimizer layer bridging the constraints `ScalarAffineFunction`-in-`Interval` is created as follows:
```julia
@bridge SplitInterval MOIB.SplitIntervalBridge () (Interval,) () () () (ScalarAffineFunction,) () ()
```
Given an optimizer `optimizer` implementing `ScalarAffineFunction`-in-`GreaterThan` and `ScalarAffineFunction`-in-`LessThan`, the optimizer
```
bridgedmodel = SplitInterval(model)
```
will additionally support `ScalarAffineFunction`-in-`Interval`.
"""
macro bridge(modelname, bridge, ss, sst, vs, vst, sf, sft, vf, vft)
    bridged_model_name = Symbol(string(modelname) * "Instance")
    bridged_funs = :(Union{$((sf.args)...), $((sft.args)...), $((vf.args)...), $((vft.args)...)})
    bridged_sets = :(Union{$((ss.args)...), $((sst.args)...), $((vs.args)...), $((vst.args)...)})

    esc(quote
        $MOIU.@model $bridged_model_name $ss $sst $vs $vst $sf $sft $vf $vft
        const $modelname{T, OT<:$MOI.ModelLike} = $MOIB.SingleBridgeOptimizer{$bridge{T}, $bridged_model_name{T}, OT}
        $MOIB.is_bridged(::$modelname, ::Type{<:$bridged_funs}, ::Type{<:$bridged_sets}) = true
        $MOIB.supports_bridging_constraint(::$modelname, ::Type{<:$bridged_funs}, ::Type{<:$bridged_sets}) = true
    end)
end


text:
#=
The tests here merely check that the interface functions work.
Testing of the implementation is done in selinv.jl.
=#

using Test
using Random
using SparseArrays
using LinearAlgebra
using SparseFactorizations

@testset "selinv_api" begin
    @testset "selinv" begin
        @testset for (conj,adj) in ((conj,adjoint),(identity,transpose))
            Random.seed!(1)
            A = sprand(ComplexF64, 10,10,0.1)
            A = I + A + adj(A)
            F = simplicial_ldlt(A; conj=conj)
            B = selinv(F)
            Bref = inv(Matrix(A)) .* ldlt_structure(A)
            @test B ≈ Bref
        end
    end
end

max_stars_repo_name: ettersi/SparseFactorizations

max_issues_repo_name: ettersi/SparseFactorizations

max_forks_repo_name: ettersi/SparseFactorizations


text:
export AbstractExplorer

using Flux

"""
    (p::AbstractExplorer)(x)
    (p::AbstractExplorer)(x, mask)

Define how to select an action based on action values.
"""
abstract type AbstractExplorer end

function (p::AbstractExplorer)(x) end
function (p::AbstractExplorer)(x, mask) end

"""
    get_prob(p::AbstractExplorer, x) -> AbstractDistribution

Get the action distribution given action values.
"""
function RLBase.get_prob(p::AbstractExplorer, x) end

"""
    get_prob(p::AbstractExplorer, x, mask)

Similart to `get_prob(p::AbstractExplorer, x)`, but here only the `mask`ed elements are considered.
"""
function RLBase.get_prob(p::AbstractExplorer, x, mask) end

# see discussion https://github.com/hill-a/stable-baselines/issues/819
Flux.testmode!(p::AbstractExplorer, mode = true) =
    @warn "trainmode/testmode will not take effect on explorer, you may consider switching to GreedyExplorer in testmode" maxlog =
        1


text:
function assemble_hermite_mass_matrix!(
    sysmatrix,
    basis,
    quad,
    mesh,
    dofspernode,
)
    ncells = number_of_cells(mesh)
    detjac = determinant_jacobian(mesh)
    cellmatrixvals = vec(mass_matrix(basis, quad, 1, detjac))

    for cellid = 1:ncells
        nodeids = nodal_connectivity(mesh, cellid)
        assemble_cell_matrix!(sysmatrix, nodeids, dofspernode, cellmatrixvals)
    end
end

function assemble_hermite_linear_form!(
    sysrhs,
    rhsfunc,
    basis,
    quad,
    mesh,
    dofspernode,
)
    numqp = length(quad)
    p, w = points(quad), weights(quad)
    valmatrix = mapslices(basis, p, dims = 1)
    ncells = number_of_cells(mesh)
    detjac = determinant_jacobian(mesh)

    for cellid = 1:ncells
        cellmap = cell_map(mesh, cellid)
        funcvals = detjac * [rhsfunc(cellmap(p[:, i])) * w[i] for i = 1:numqp]
        vals = valmatrix * funcvals
        nodeids = nodal_connectivity(mesh, cellid)
        assemble_cell_rhs!(sysrhs, nodeids, dofspernode, vals)
    end
end


text:
using DitherPunk
using Images
using TestImages

# # Test image gallery
# This gallery uses images from [*TestImages.jl*](https://testimages.juliaimages.org).
function preprocess(img)
    img = Gray.(img)
    return imresize(img; ratio=1 / 2)
end

file_names = [
    "cameraman", "lake_gray", "house", "fabio_gray_512", "mandril_gray", "peppers_gray"
]
imgs = [preprocess(testimage(file)) for file in file_names]
mosaic(imgs; ncol=3)

# Our test function `test_on_images` just runs a dithering algorithm on all six images
# in linear color space (`to_liner=true`).
function test_on_images(alg; to_linear=false)
    dithers = [dither(img, alg; to_linear) for img in imgs]
    return mosaic(dithers; ncol=3)
end

# ## Threshold dithering
# #### `ConstantThreshold`
test_on_images(ConstantThreshold())

# #### `WhiteNoiseThreshold`
test_on_images(WhiteNoiseThreshold())

# ## Ordered dithering
# #### Bayer matrices
# **Level 1**
test_on_images(Bayer())
# **Level 2**
test_on_images(Bayer(; level=2))
# **Level 3**
test_on_images(Bayer(; level=3))
# **Level 4**
test_on_images(Bayer(; level=4))

# #### `ClusteredDots`
test_on_images(ClusteredDots())

# #### `CentralWhitePoint`
test_on_images(CentralWhitePoint())

# #### `BalancedCenteredPoint`
test_on_images(BalancedCenteredPoint())

# #### `Rhombus`
test_on_images(Rhombus())

# #### ImageMagick methods
test_on_images(IM_checks())
#
test_on_images(IM_h4x4a())
#
test_on_images(IM_h6x6a())
#
test_on_images(IM_h8x8a())
#
test_on_images(IM_h4x4o())
#
test_on_images(IM_h6x6o())
#
test_on_images(IM_h8x8o())
#
test_on_images(IM_c5x5())
#
test_on_images(IM_c6x6())
#
test_on_images(IM_c7x7())

# ## Error diffusion
# #### `SimpleErrorDiffusion`
test_on_images(SimpleErrorDiffusion())

# #### `FloydSteinberg`
test_on_images(FloydSteinberg())

# #### `JarvisJudice`
test_on_images(JarvisJudice())

# #### `Stucki`
test_on_images(Stucki())

# #### `Burkes`
test_on_images(Burkes())

# #### `Sierra`
test_on_images(Sierra())

# #### `TwoRowSierra`
test_on_images(TwoRowSierra())

# #### `SierraLite`
test_on_images(SierraLite())

# #### `Fan93()`
test_on_images(Fan93())

# #### `ShiauFan()`
test_on_images(ShiauFan())

# #### `ShiauFan2()`
test_on_images(ShiauFan2())

# #### `Atkinson()`
test_on_images(Atkinson())


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module GVTree

export plotTree

using auxfuncs

import MCTSdpw: DPW, State, StateNode

plotTree(dpw::DPW,filename::AbstractString) = plotTree(dpw.s,filename)

function plotTree(d::Dict{State,StateNode},filename::AbstractString)
  println("Size of sdict: ", length(d))
  outputFile = open(filename,"w")
  write(outputFile,"digraph gvtree {\n")
  for (s,sn) in d
    actions = collect(keys(d[s].a))
    for a in actions
      sp = collect(keys(d[s].a[a].s))[1] #assumes deterministic transition
      @printf(outputFile,"\n%s",hash(s))
      write(outputFile," -> ")
      @printf(outputFile,"%s;",hash(sp))
      @printf(outputFile,"\n%s",hash(sp)) #there are overlaps, but that's ok
      write(outputFile,""" [label="",shape=circle, style=solid, color=black]""")
    end
    @printf(outputFile,"\n%s",hash(s)) #there are overlaps, but that's ok
    write(outputFile,""" [label="",shape=circle, style=solid, color=black]""")
  end
  write(outputFile,"}")
  close(outputFile)
end

end

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
# this file tests various aspects of the
using MagmaThermoKinematics
using Plots
using LinearAlgebra
using SpecialFunctions
using Test

const CreatePlots = false      # easy way to deactivate plotting throughout

function test_SolidFraction()
  # test solid fraction

  # Model parameters
  H                     =   75e3;                              # Height

  # Define  grid
  Nz                    =   33;                                 # resolution of coarse grid
  dz                    =   H/(Nz-1);                           # grid size [m]
  z                     =   0:dz:((Nz-1)*dz);                   # 1D coordinate arrays
  coords                =   collect(Iterators.product(z))       # generate coordinates from 1D coordinate vectors
  Z                     =   (x->x[1]).(coords);                 # transfer coords to 3D arrays
  Grid, Spacing         =   (z), (dz);

  # Define T
  T                     =   Z./1e3.*20;
  T[ (Z.>50) .& (Z.< 60)]     .=   900;

  Phi_o                 = 0*T;
  dPhi_dt               = 0*T;
  dt                    = 1.0;
  Phi                   = 0*T;

  Phi, dPhi_dt = SolidFraction(T, Phi_o, dt); # call routine

  if CreatePlots
      p1          =   plot(T[:,1], Phi[:,1],  ylabel="Solid fraction ", xlabel="Temperature [C]",  dpi=150)

      plot(p1);
      png("MeltingRelationship")
  end

  out = norm(Phi[:],2)

  return out;
end

# ===================================================================================================

@testset "MeltingRelationShips" begin
  @test test_SolidFraction() ≈  4.132199893765409  atol=1e-8;
end;


text:
flip(obj::Vector{Float64}) = reverse(obj)
flip(obj::Vector{Vector{Float64}}) = map(flip, obj)
flip(obj::Vector{Vector{Vector{Float64}}}) = map(flip, obj)

flip(obj::GeoJSON.Point) = GeoJSON.Point(flip(obj.coordinates))
flip(obj::GeoJSON.MultiPoint) = GeoJSON.MultiPoint(flip(obj.coordinates))
flip(obj::GeoJSON.LineString) = GeoJSON.LineString(flip(obj.coordinates))
flip(obj::GeoJSON.MultiLineString) = GeoJSON.MultiLineString(map(flip,obj.coordinates))
flip(obj::GeoJSON.Polygon) = GeoJSON.Polygon(map(flip,obj.coordinates))
flip(obj::GeoJSON.MultiPolygon) = GeoJSON.MultiPolygon(map(flip,obj.coordinates))
flip(obj::GeoJSON.GeometryCollection) = GeoJSON.GeometryCollection(map(flip, obj.geometries)) # be careful of type here
flip(obj::GeoJSON.Feature) = GeoJSON.Feature(flip(obj.geometry))
flip(obj::GeoJSON.FeatureCollection) = GeoJSON.FeatureCollection(map(flip, obj.features))

explode(obj::Vector{Float64}) = GeoJSON.Feature[GeoJSON.Feature(GeoJSON.Point(obj))]
explode(obj::Vector{Vector{Float64}}) = vcat(map(explode, obj)...)
explode(obj::Vector{Vector{Vector{Float64}}}) = vcat(map(explode, obj)...)

explode(obj::GeoJSON.Point) = GeoJSON.FeatureCollection(explode(obj.coordinates))
explode(obj::GeoJSON.MultiPoint) = GeoJSON.FeatureCollection(explode(obj.coordinates))
explode(obj::GeoJSON.LineString) = GeoJSON.FeatureCollection(explode(obj.coordinates))
explode(obj::GeoJSON.MultiLineString) = GeoJSON.FeatureCollection(explode(obj.coordinates))
explode(obj::GeoJSON.Polygon) = GeoJSON.FeatureCollection(explode(obj.coordinates))
explode(obj::GeoJSON.MultiPolygon) = GeoJSON.FeatureCollection(vcat(map(explode, obj.coordinates)...))
explode(obj::GeoJSON.Feature) = explode(obj.geometry)

function explode(obj::GeoJSON.GeometryCollection)
    fc = GeoJSON.Feature[]
    for geom in obj.geometries
        fc = vcat(fc, explode(geom.coordinates)) # inefficient?
    end
    GeoJSON.FeatureCollection(fc)
end

function explode(obj::GeoJSON.FeatureCollection)
    fc = GeoJSON.Feature[]
    for feature in obj.features
        fc = vcat(fc, explode(feature.geometry))
    end
    GeoJSON.FeatureCollection(fc)
end

combine(obj::Vector{GeoJSON.Point}) = GeoJSON.MultiPoint(map(GeoJSON.coordinates, obj))
combine(obj::Vector{GeoJSON.LineString}) = GeoJSON.MultiLineString(map(GeoJSON.coordinates, obj))
combine(obj::Vector{GeoJSON.Polygon}) = GeoJSON.MultiPolygon(map(GeoJSON.coordinates, obj))
combine(obj::GeoJSON.GeometryCollection) = combine(obj.geometries)
combine(obj::GeoJSON.Feature) = combine(obj.geometry::GeoJSON.GeometryCollection)
combine(obj::GeoJSON.FeatureCollection) = combine(typeof(obj.features[1].geometry)[feature.geometry for feature in obj.features])

function isclockwise(ring::Vector{Vector{Float64}}) # not sure if this works
    sum = 0.0
    for i=2:length(ring)
        prev = ring[i-1]; cur = ring[i]
        sum += (cur[1] - prev[1]) * (cur[2] + prev[2])
    end
    sum > 0
end

function intersection(line1startx::Float64, line1starty::Float64, line1endx::Float64, line1endy::Float64,
                      line2startx::Float64, line2starty::Float64, line2endx::Float64, line2endy::Float64)
    online1 = false
    online2 = false

    denom = ((line2endy - line2starty)*(line1endx - line1startx)) - ((line2endx - line2startx)*(line1endy - line1starty))

    a = line1starty - line2starty
    b = line1startx - line2startx
    num1 = ((line2endx - line2startx) * a) - ((line2endy - line2starty) * b)
    num2 = ((line1endx - line1startx) * a) - ((line1endy - line1starty) * b)
    a = num1 / denom
    b = num2 / denom

    # if we cast these lines infinitely in both directions, they intersect here:
    x = line1startx + (a * (line1endx - line1startx))
    y = line1starty + (a * (line1endy - line1starty))

    # they intersect if:
    online1 = a > 0 && a < 1 # if line1 is a segment and line2 is infinite
    online2 = b > 0 && b < 1 # if line2 is a segment and line1 is infinite

    # if line1 and line2 are segments, they intersect if both of the above are true
    if online1 && online2
        return Float64[x, y]
    else
        return Float64[]
    end
end

function kinks(obj::GeoJSON.Polygon)
    results = GeoJSON.FeatureCollection()
    for ring1 in obj.coordinates
        for ring2 in obj.coordinates
            for i=1:length(ring1)-1
                for k=1:length(ring2)-1
                    intersect = intersection(ring1[i][1], ring1[i][2], ring1[i+1][1], ring1[i+1][2],
                                             ring2[k][1], ring2[k][2], ring2[k+1][1], ring2[k+1][2])
                    if length(intersect) == 2
                        push!(results.features, GeoJSON.Feature(GeoJSON.Point(intersect)))
                    end
                end
            end
        end
    end
    results
end
kinks(obj::GeoJSON.Feature) = kinks(obj.geometry)


text:
@testset "Multi-threaded runs test, with init pop" begin

@info "Commencing test of multi-threaded parallel simulation runs, with initial population. This test may take a few minutes, depending on the user's hardware and Julia settings."
tStart = now()
rm( "parallel/tmps", recursive=true, force=true )
rm( "parallel/simDB2.sqlite", force=true )

mrs = MRS( mpSim )
nRuns = 50
setMRSruns!( mrs, nRuns )
setMRSdatabaseName!( mrs, "parallel/simDB2.sqlite" )
setMRSmaxThreads!( mrs, 10 )
uploadSnapshot( mrs, "snapshot/test foto 1", "Sheet1", ["B", "C", "D"],
    ("A", "G", "F", "E", "H") )
@test verifySimulation!( mrs )

run( mrs )
tElapsed = ( now() - tStart ).value / 1000
@info string( "Running ", nRuns, " simulations in parallel on ",
    min( 10, Threads.nthreads() ), " threads took ", tElapsed,
    " seconds. Thank you for your patience." )

db = SQLite.DB( "parallel/simDB2.sqlite" )
@test length( SQLite.tables( db )[:name] ) == 3 * nRuns + 5
@test all( 1:nRuns ) do ii
    queryCmd = string( "SELECT * FROM `", mrs.mpSim.transDBname, ii, "`",
        "\n    WHERE transition IS 'Init'")
    initpop = DataFrame(DBInterface.execute( db, queryCmd ))
    size( initpop, 1 ) == 50
end  # all( 1:nRuns ) do ii

println()

end  # @testset "Multi-threaded runs test, with init pop"


text:
## based on functions in the tests directory of the MUMPS package
using Base.Test
using SuperLU

ddx(n) = spdiagm((-ones(n), ones(n)),[0,1],n,n+1)

function getDivGrad(n1,n2,n3)
    Div = hcat(kron(speye(n3),kron(speye(n2),ddx(n1))),
               kron(speye(n3),kron(ddx(n2),speye(n1))),
               kron(ddx(n3),kron(speye(n2),speye(n1))))
    Div*Div';
end

dgcube(n) = getDivGrad(n,n,n)

for n in (8,16,24)
    A = dgcube(n)
    b = randn(size(A,1))
    xumf = A\b
    xslu = NCMat(A)\b
    @test_approx_eq (A*xumf) b
    @test_approx_eq (A*xslu) b
    gc()
    umftimes = [@elapsed A\b for i in 1:10]
    gc()
    slutimes = [@elapsed NCMat(A)\b for i in 1:10]
    println("Times for solving divgrad system of size $(size(A)), nnz = $(nfilled(A))")
    showcompact(umftimes)
    println()
    showcompact(slutimes)
    println()
    println("UMFPACK: mean = $(mean(umftimes)), median = $(median(umftimes))")
    println("SuperLU: mean = $(mean(slutimes)), median = $(median(slutimes))")
end


text:
#=
    The problems are taken from http://www.hakank.org/julia/constraints/

    Thanks a lot to Håkan Kjellerstrand
=#

function cumulative(model, start, duration, resource, limit)
    tasks = [i for i in 1:length(start) if resource[i] > 0 && duration[i] > 0]
    num_tasks = length(tasks)
    times_min_a = round.(Int,[JuMP.lower_bound(start[i]) for i in tasks])
    times_min = minimum(times_min_a)
    times_max_a = round.(Int,[JuMP.upper_bound(start[i])+duration[i] for i in tasks])
    times_max = maximum(times_max_a)
    for t in times_min:times_max
        b  = @variable(model, [1:num_tasks], Bin)
        for i in tasks
            # is this task active during this time t?
            @constraint(model, b[i] := { start[i] <= t && t < start[i]+duration[i]})
        end
        # Check that there's no conflicts in time t
        @constraint(model,sum(b[i]*resource[i] for i in tasks) <= limit)
    end
end

#
# no_overlap(model, begins,durations)
#
# Ensure that there is no overlap between the tasks.
#
function no_overlap(model, begins,durations)
    n = length(begins)
    for i in 1:n, j in i+1:n
        b = @variable(model,[1:2], Bin)
        @constraint(model,b[1] := {begins[i] + durations[i] <= begins[j]})
        @constraint(model,b[2] := {begins[j] + durations[j] <= begins[i]})
        @constraint(model, sum(b) >= 1)
    end
end

# From: http://www.hakank.org/julia/constraints/furniture_moving.jl
function furniture_moving()

    model = Model(optimizer_with_attributes(CS.Optimizer,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            "branch_split"=>:InHalf, # <-

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    # Furniture moving problem
    n = 4
    # [piano, chair, bed, table]
    durations = [30,10,15,15]
    resources = [3,1,3,2] # people needed per task
    @variable(model, 0 <= start_times[1:n] <= 60, Int)
    @variable(model, 0 <= end_times[1:n]   <= 60, Int)
    @variable(model, 1 <= limit <= 3, Int)
    @variable(model, 0 <= max_time <= 60, Int)
    @constraint(model, end_times .<= max_time)

    for i in 1:n
        @constraint(model,end_times[i] == start_times[i] + durations[i])
    end
    cumulative(model, start_times, durations, resources, limit)

    # @objective(model,Min,max_time)

    optimize!(model)

    status = JuMP.termination_status(model)
    @assert status == MOI.OPTIMAL
    @assert JuMP.value(limit) ≈ 3
end

# from http://www.hakank.org/julia/constraints/organize_day.jl
function organize_day(problem,all_solutions=true)
    # cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    # glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    # ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                        ))
    tasks = problem[:tasks]
    durations = problem[:durations]
    precedences = problem[:precedences]
    start_time = problem[:start_time]
    end_time = problem[:end_time]

    n = length(tasks)
    @variable(model, start_time <= begins[1:n] <= end_time, Int)
    @variable(model, start_time <= ends[1:n] <= end_time, Int)

    for i in 1:n
        @constraint(model,ends[i] == begins[i] + durations[i])
    end

    no_overlap(model,begins,durations)

    # precedences
    for (a,b) in eachrow(precedences)
        @constraint(model, ends[a] <= begins[b])
    end

    @constraint(model,begins[1] >= 11)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    @assert status == MOI.OPTIMAL
    if all_solutions
        @assert MOI.get(model, MOI.ResultCount()) == 5
    end
end