text:
# ollam.jl
#
# author: Wade Shen
# swade@ll.mit.edu
# Copyright © 2009 Massachusetts Institute of Technology, Lincoln Laboratory
# version 0.1
#
# 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.
module Ollam
using Stage, DataStructures #LIBSVM
import Base: copy, start, done, next, length, dot, getindex
export LinearModel, RegressionModel, copy, score, best, train_perceptron, test_classification, test_regression,
       train_svm, train_mira, train_libsvm, lazy_map, indices,
       print_confusion_matrix, hildreth, setup_hildreth, zero_one_loss, linear_regression_loss,
       regress_perceptron, regress_mira

# ----------------------------------------------------------------------------------------------------------------
# Utilities
# ----------------------------------------------------------------------------------------------------------------
immutable Map{I}
    flt::Function
    itr::I
end
lazy_map(f::Function, itr) = Map(f, itr)

function start(m :: Map)
  s = start(m.itr)
  return s
end

function next(m :: Map, s)
  n, ns = next(m.itr, s)
  return (m.flt(n), ns)
end

done(m :: Map, s) = done(m.itr, s)
length(m :: Map) = length(m.itr)

# linear algebra helpers
indices(a :: SparseMatrixCSC) = a.rowval
indices(a :: Vector)          = 1:length(a)
getnth(a :: SparseMatrixCSC, i :: Integer) = a.nzval[i]

sqr(a::Vector) = norm(a)^2
function sqr(a::SparseMatrixCSC)
  total = 0.0
  for i in indices(a)
    total += a[i] * a[i]
  end
  return total
end
dot(a::SparseMatrixCSC, b::Vector) = dot(b, a)
dot(a::SparseMatrixCSC, b::Matrix) = dot(b, a)
function dot(a::Union{Vector, SparseMatrixCSC}, b::SparseMatrixCSC)
  total = 0.0
  for i in indices(b)
    total += a[i] * b[i]
  end
  return total
end
function dot(a::Matrix, b::SparseMatrixCSC)
  total = 0.0
  for i in indices(b)
    total += a[1, i] * b[i]
  end
  return total
end
dot(a::Matrix, b::Vector) = (a * b)[1]

function print_confusion_matrix(confmat; width = 10, logger = Log(STDERR))
  width = max(5, width)
  sfmt = "%$(width)s"
  dfmt = "%$(width)d"
  ffmt = "%$(width).$(width-3)f"
  @eval s(x) = @sprintf($sfmt, x)
  @eval d(x) = @sprintf($dfmt, x)
  @eval f(x) = @sprintf($ffmt, x)
  total, errors = 0, 0
  accs = 0.0

  str = s("")
  for t in keys(confmat)
    str *= " " * s(t)
  end
  @sep logger
  @info logger "$str " * s("N") * " " * s("cl %")

  for t in keys(confmat)
    str = s(t)
    rtotal, rerrors = 0, 0
    for h in keys(confmat)
      str *= " " * d(confmat[t][h]) # @sprintf(" %10d", confmat[t][h])
      if t != h
        rerrors += confmat[t][h]
      end
      rtotal += confmat[t][h]
    end
    errors += rerrors
    total  += rtotal
    @info logger "$str" * " " * d(rtotal) * " " * f(1.0 - rerrors/rtotal) #@sprintf(" %10d %10.7f", rtotal, 1.0 - rerrors/rtotal)
    accs += 1.0 - rerrors/rtotal
  end
  @sep logger

  @info logger "overall accuracy = $(1.0 - errors/total), average class accuracy = $(accs / length(keys(confmat)))"

end

# ----------------------------------------------------------------------------------------------------------------
# Types
# ----------------------------------------------------------------------------------------------------------------
type LinearModel{T}
  weights     :: Matrix{Float64}
  b           :: Vector{Float64}

  class_index :: Dict{T, Int32}
  index_class :: Vector{T}
end

dims(lm :: LinearModel)    = size(lm.weights, 2)
classes(lm :: LinearModel) = size(lm.weights, 1)

function LinearModel{T}(classes::Dict{T, Int32}, dims)
  index = Array(T, length(classes))
  for (k, i) in classes
    index[i] = k
  end

  return LinearModel(zeros(length(index), dims), zeros(length(index)), classes, index)
end

RegressionModel(dims) = LinearModel(Dict{AbstractString,Int32}("regressor" => 1), dims)

copy(lm :: LinearModel) = LinearModel(copy(lm.weights), copy(lm.b), copy(lm.class_index), copy(lm.index_class))
score(lm :: LinearModel, fv::Vector) = lm.weights * fv + lm.b
score(lm :: LinearModel, fv::SparseMatrixCSC) = vec(lm.weights * fv + lm.b)

function acc_update(model :: LinearModel, acc :: LinearModel)
  for x = 1:size(acc.weights, 1)
    for y = 1:size(acc.weights, 2)
      acc.weights[x, y] += model.weights[x, y]
    end
  end
end

function best{T <: AbstractFloat}(scores :: Vector{T})
  bidx = indmax(scores)
  return bidx, scores[bidx]
end

function test_classification(lm :: LinearModel, fvs, truth; record = (truth, hyp) -> nothing)
  errors = 0
  total  = 0

  for (fv, t) in zip(fvs, truth)
    scores  = score(lm, fv)
    bidx, b = best(scores)
    if lm.index_class[bidx] != t
      errors += 1
    end
    record(t, lm.index_class[bidx])
    total += 1
  end

  return errors / total
end

function test_regression(lm :: LinearModel, fvs, truth; lossfn = linear_regression_loss)
  N = 0
  total_se = 0.0
  for (fv, t) in zip(fvs, truth)
    scores  = score(lm, fv)
    bidx, b = best(scores)
    total_se += lossfn(b, t)^2
    N += 1
  end

  return total_se / N
end

# ----------------------------------------------------------------------------------------------------------------
# Loss functions
# ----------------------------------------------------------------------------------------------------------------
zero_one_loss(hyp, ref) = hyp == ref ? 0.0 : 1.0
linear_regression_loss(hyp, ref) = ref - hyp

# ----------------------------------------------------------------------------------------------------------------
# Perceptron
# ----------------------------------------------------------------------------------------------------------------
function perceptron_update(model, c, alpha, fv)
  for i in indices(fv)
    model.weights[c, i] += alpha * fv[i]
  end
end

function train_perceptron(fvs, truth, init_model; learn_rate = 1.0, average = true, iterations = 40,
                          logger = Log(STDERR), verbose = true)
  model = copy(init_model)
  acc   = LinearModel(init_model.class_index, dims(init_model))
  numfv = 0
  for fv in fvs
    numfv += 1
  end
  avg_w = iterations * numfv

  for i = 1:iterations
    fj = 1
    for (fv, t) in zip(fvs, truth)
      scores     = score(model, fv)
      bidx, b    = best(scores)
      w          = (avg_w - (numfv * (i - 1) + fj) + 1)
      if model.index_class[bidx] != t
        for c = 1:classes(model)
          sign = model.index_class[c] == t ? 1.0 : (-1.0 / (classes(model) - 1))
          perceptron_update(model, c, sign * learn_rate, fv)
          perceptron_update(acc, c, w * sign * learn_rate, fv)
        end
      end
      fj += 1
    end
    if verbose
      @info logger @sprintf("iteration %3d complete (Training error rate: %7.3f%%)", i, test_classification(model, fvs, truth) * 100.0)
    end
  end

  if average
    return acc
  else
    return model
  end
end

function regress_perceptron(fvs, truth, init_model; learn_rate = 0.01, average = true, iterations = 40, C = 0.1,
                            logger = Log(STDERR), verbose = true, lossfn = linear_regression_loss)
  model = copy(init_model)
  acc   = LinearModel(init_model.class_index, dims(init_model))
  numfv = 0
  for fv in fvs
    numfv += 1
  end
  avg_w = iterations * numfv

  for i = 1:iterations
    fj = 1
    for (fv, t) in zip(fvs, truth)
      scores  = score(model, fv)
      bidx, b = best(scores)
      w       = (avg_w - (numfv * (i - 1) + fj) + 1)
      loss    = lossfn(b, t)
      alpha   = max(min(C, loss * learn_rate), -C)
      #@debug @sprintf("loss = %10.3f, alpha = %10.7f, expected = %10.7f, ref = %10.7f fv = %s", loss, alpha, b, t, fv)
      perceptron_update(model, 1, alpha, fv)
      perceptron_update(acc, 1, w * alpha, fv)
      fj += 1
    end
    if verbose
      @info logger @sprintf("iteration %3d complete (Training RMSE: %7.3f)", i, sqrt(test_regression(model, fvs, truth)))
    end
  end

  if average
    acc.weights /= avg_w
    return acc
  else
    return model
  end
end

# ----------------------------------------------------------------------------------------------------------------
# MIRA
# ----------------------------------------------------------------------------------------------------------------
function mira_update(weights, bidx, tidx, alpha, fv :: SparseMatrixCSC)
  ind = indices(fv)
  for i in 1:length(ind)
    idx = ind[i]
    tmp = alpha * fv.nzval[i] # direct access avoids fv[idx] binary search in getindex()
    weights[bidx, idx] -= tmp
    weights[tidx, idx] += tmp
  end
end

function mira_update(weights, bidx, tidx, alpha, fv :: Array)
  for idx in indices(fv)
    tmp = alpha * fv[idx] # slow for sparse because of getindex() [see above specialization]
    weights[bidx, idx] -= tmp
    weights[tidx, idx] += tmp
  end
end

type HildrethState
  k           :: Int32
  alpha       :: Vector{Float64}
  F           :: Vector{Float64}
  kkt         :: Vector{Float64}
  C           :: Float64
  A           :: Matrix{Float64}
  is_computed :: Vector{Bool}
  EPS         :: Float64
  ZERO        :: Float64
  MAX_ITER    :: Float64
end

function setup_hildreth(;k = 5, C = 0.1, EPS = 1e-8, ZERO = 0.0000000000000001, MAX_ITER = 10000) # assumes that the number of contraints == number of distances
  alpha       = zeros(k)
  F           = zeros(k)
  kkt         = zeros(k)
  A           = zeros(k, k)
  is_computed = falses(k)
  return HildrethState(k, alpha, F, kkt, C, A, is_computed, EPS, ZERO, MAX_ITER)
end

# translated from Ryan McDonald's MST Parser
function hildreth(a, b, h)
  max_kkt = -Inf
  max_kkt_i = -1

  for i = 1:h.k
    h.A[i, i] = dot(a[i], a[i])
    h.kkt[i] = h.F[i] = b[i]
    h.is_computed[i] = false
    if h.kkt[i] > max_kkt
      max_kkt   = h.kkt[i]
      max_kkt_i = i
    end
    h.alpha[i] = 0.0
  end

  iter = 0
  while max_kkt >= h.EPS && iter < h.MAX_ITER
    diff_alpha = h.A[max_kkt_i, max_kkt_i] <= h.ZERO ? 0.0 : h.F[max_kkt_i] / h.A[max_kkt_i, max_kkt_i]
    try_alpha  = h.alpha[max_kkt_i] + diff_alpha
    add_alpha  = 0.0

    if try_alpha < 0.0
      add_alpha = - h.alpha[max_kkt_i]
    elseif try_alpha > h.C
      add_alpha = h.C - h.alpha[max_kkt_i]
    else
      add_alpha = diff_alpha
    end

    h.alpha[max_kkt_i] += add_alpha

    if !h.is_computed[max_kkt_i]
      for i = 1:h.k
	h.A[i, max_kkt_i] = dot(a[i], a[max_kkt_i])
	h.is_computed[max_kkt_i] = true
      end
    end
    for i = 1:h.k
      h.F[i]  -= add_alpha * h.A[i, max_kkt_i]
      h.kkt[i] = h.F[i]
      if h.alpha[i] > (h.C - h.ZERO)
	h.kkt[i] = -h.kkt[i]
      elseif h.alpha[i] > h.ZERO
	h.kkt[i] = abs(h.F[i])
      end
    end
    max_kkt   = -Inf
    max_kkt_i = -1
    for i = 1:h.k
      if h.kkt[i] > max_kkt
        max_kkt   = h.kkt[i]
        max_kkt_i = i
      end
    end
    iter += 1
  end

  return h.alpha
end

type ScaledVec{T}
  k :: T
  v
end
getindex(dv :: ScaledVec, i :: Integer) = dv.k * dv.v[i]
dot(v1 :: ScaledVec, v2 :: ScaledVec) = v1.k * v2.k * dot(v1.v, v2.v)

function train_mira(fvs, truth, init_model;
                    average = true, C = 0.1, k = 1, iterations = 20, lossfn = zero_one_loss,
                    logger = Log(STDERR), verbose = true)
  model = copy(init_model)
  acc   = LinearModel(init_model.class_index, dims(init_model))
  acc2  = LinearModel(init_model.class_index, dims(init_model))
  numfv = 0
  for fv in fvs
    numfv += 1
  end

  h       = setup_hildreth(k = min(k, length(model.class_index)), C = C)
  b       = Array(Float64, h.k)
  kidx    = Array(Int32, h.k)
  distvec = Array(ScaledVec, h.k) #Array(Union{SparseMatrixCSC, Vector}, h.k)
  avg_w   = iterations * numfv

  for i = 1:iterations
    fj = 1
    alpha = 0.0
    for (fv, t) in zip(fvs, truth)
      scores    = score(model, fv)
      tidx      = model.class_index[t]
      tgt_score = scores[tidx]
      w         = (avg_w - (numfv * (i - 1) + fj) + 1)

      if h.k > 1 # K-best
        sorted = sortperm(scores, rev = true)
        for n = 1:h.k
          cidx        = sorted[n]
          score       = scores[cidx]
          class       = model.index_class[cidx]
          loss        = lossfn(t, class)
          dist        = tgt_score - score

          b[n]       = loss - dist
          distvec[n] = ScaledVec(2.0, fv) # 2 * fv # slow due to allocation
          kidx[n]    = cidx
        end

        alphas = hildreth(distvec, b, h)

        # update
        for n = 1:h.k
          mira_update(model.weights, kidx[n], tidx, alphas[n], fv)
          mira_update(acc.weights, kidx[n], tidx, w * alphas[n], fv)
        end
      else # 1-best
        bidx, b_score = best(scores)
        class = model.index_class[bidx]
        loss  = lossfn(t, class)
        dist  = tgt_score - b_score
        alpha = max(min((loss - dist) / (2 * sqr(fv)), C), -C)

        #@debug logger "truth: $t -- best class $(model.index_class[bidx]) -- best score: $b_score, truth score: $tgt_score"
        #@debug logger "loss = $loss, dist = $dist [$tgt_score - $b_score], denom = $(2 * norm(fv)^2), alpha = $alpha"
        mira_update(model.weights, bidx, tidx, alpha, fv)
        mira_update(acc.weights, bidx, tidx, w * alpha, fv)
      end

      fj += 1
    end
    if verbose
      @info logger @sprintf("iteration %3d complete (Training error rate: %7.3f%%)", i, test_classification(model, fvs, truth) * 100.0)
    end
  end

  if average
    return acc
  else
    return model
  end
end

function mira_regress_update(weights, bidx, alpha, fv :: SparseMatrixCSC)
  ind = indices(fv)
  for i in 1:length(ind)
    idx = ind[i]
    tmp = alpha * fv.nzval[i] # direct access avoids fv[idx] binary search in getindex()
    weights[bidx, idx] += tmp
  end
end

function mira_regress_update(weights, bidx, alpha, fv :: Array)
  for idx in indices(fv)
    tmp = alpha * fv[idx] # slow for sparse because of getindex() [see above specialization]
    weights[bidx, idx] += tmp
  end
end

function regress_mira(fvs, truth, init_model;
                    average = true, C = 0.1, min_loss = 0.0, iterations = 20, lossfn = linear_regression_loss,
                    logger = Log(STDERR), verbose = true)
  model = copy(init_model)
  acc   = LinearModel(init_model.class_index, dims(init_model))
  numfv = 0
  for fv in fvs
    numfv += 1
  end

  avg_w   = iterations * numfv

  for i = 1:iterations
    fj = 1
    alpha = 0.0
    for (fv, t) in zip(fvs, truth)
      scores        = score(model, fv)
#      scores        =  max(min(score(model, fv),1.0),0.0)
      w             = (avg_w - (numfv * (i - 1) + fj) + 1)
      bidx, b_score = best(scores)
      loss          = lossfn(b_score, t)
      if abs(loss) < min_loss
        loss = 0.0
      end
      alpha         = max(min(loss / sqr(fv), C), -C)

      mira_regress_update(model.weights, bidx, alpha, fv)
      #if fj % 100 == 0 # abs(loss) > 7.0
      #  @debug "$fj: $(model.weights) <- $loss"
      #end
      mira_regress_update(acc.weights, bidx, w * alpha, fv)

      fj += 1
    end
    if verbose
      @info logger @sprintf("iteration %3d complete (Training RMSE: %7.3f)", i, sqrt(test_regression(model, fvs, truth)))
    end
  end

  if average
    acc.weights /= avg_w
    return acc
  else
    return model
  end
end

# ----------------------------------------------------------------------------------------------------------------
# SVM
# ----------------------------------------------------------------------------------------------------------------
#=
immutable SVMNode
  index::Int32
  value::Float64
end

immutable SVMModel
  param::LIBSVM.SVMParameter
  nr_class::Int32
  l::Int32
  SV::Ptr{Ptr{LIBSVM.SVMNode}}
  sv_coef::Ptr{Ptr{Float64}}
  rho::Ptr{Float64}
  probA::Ptr{Float64}
  probB::Ptr{Float64}
  sv_indices::Ptr{Int32}

  label::Ptr{Int32}
  nSV::Ptr{Int32}

  free_sv::Int32
end

function transfer_sv(p::Ptr{LIBSVM.SVMNode})
  ret  = (LIBSVM.SVMNode)[]
  head = unsafe_load(p)

  while head.index != -1
    push!(ret, head)
    p += 16
    head = unsafe_load(p)
  end
  return ret
end

function transfer(svm; logger = Log(STDERR))
  # unpack svm model
  ptr    = unsafe_load(convert(Ptr{SVMModel}, svm.ptr))
  nSV    = pointer_to_array(ptr.nSV, ptr.nr_class)
  xSV    = pointer_to_array(ptr.SV, ptr.l)
  SV     = (Array{LIBSVM.SVMNode, 1})[ transfer_sv(x) for x in xSV ]
  xsvc   = pointer_to_array(ptr.sv_coef, ptr.nr_class)
  svc    = (Array{Float64, 1})[ pointer_to_array(x, ptr.l) for x in xsvc ]
  labels = pointer_to_array(ptr.label, ptr.nr_class)
  rho    = pointer_to_array(ptr.rho, 1)
  @debug logger "# of SVs = $(length(SV)), labels = $labels, rho = $rho, $(svm.labels)"

  # precompute classifier weights
  start = 1
  weights = zeros(svm.nfeatures)

  for i = 1:ptr.nr_class
    for sv_offset = 0:(nSV[i]-1)
      sv = SV[start + sv_offset]
      for d = 1:length(sv)
        weights[sv[d].index] += svc[1][start + sv_offset] * sv[d].value
      end
    end
    start += nSV[i]
  end
  b = -rho[1]

  if svm.labels[1] == -1
    weights = -weights
    b       = -b
  end

  return (weights, b)
end

function train_libsvm(fvs, truth; C = 1.0, nu = 0.5, cache_size = 200.0, eps = 0.0001, shrinking = true, verbose = false, gamma = 0.5, logger = Log(STDERR))
  i = 1
  classes = Dict{Any, Int32}()

  for t in truth
    if !(t in keys(classes))
      classes[t] = i
      i += 1
    end
  end

  feats = hcat(fvs...)
  model = LinearModel(classes, size(feats, 1))

  svms  = Array(Any, length(classes))
  refs  = Array(Any, length(classes))

  for (t, ti) in classes
    @timer logger "training svm for class $t (index: $ti)" begin
      refs[ti] = @spawn begin
        svm_t = svmtrain(map(c -> c == t ? 1 : -1, truth), feats;
                         gamma = gamma, C = C, nu = nu, kernel_type = int32(0), degree = int32(1), svm_type = int32(0),
                         cache_size = cache_size, eps = eps, shrinking = shrinking, verbose = verbose)
        transfer(svm_t, logger = logger)
      end
    end
  end

  for c = 1:length(refs)
    svms[c] = fetch(refs[c])
  end

  for c = 1:length(svms)
    weights_c, b_c = svms[c]
    for i = 1:length(weights_c)
      model.weights[c, i] = weights_c[i]
    end
    model.b[c] = b_c
  end

  return model
end

function train_svm(fvs, truth; C = 0.01, batch_size = -1, norm = 2, iterations = 100, logger = Log(STDERR))
  i = 1
  classes = Dict{Any, Int32}()

  for t in truth
    if !(t in keys(classes))
      classes[t] = i
      i += 1
    end
  end

  fs = hcat(fvs...)
  feats = vcat(fs, ones(1, size(fs, 2)))

  if batch_size == -1
    batch_size = size(feats, 2)
  end

  model = LinearModel(classes, size(fs, 1))

  svms  = Array(Any, length(classes))
  refs  = Array(Any, length(classes))

  for (t, ti) in classes
    @timer logger "training svm for class $t (index: $ti)" begin
      refs[ti] = @spawn begin
        svm_t = cddual(feats, map(c -> c == t ? 1 : -1, truth);
                       C = C, norm = norm, randomized = true, maxpasses = iterations)
        (svm_t.w[1:end-1], svm_t.w[end])
      end
    end
  end

  for c = 1:length(refs)
    svms[c] = fetch(refs[c])
  end

  for c = 1:length(svms)
    weights_c, b_c = svms[c]
    for i = 1:length(weights_c)
      model.weights[c, i] = weights_c[i]
    end
    model.b[c] = b_c
  end

  return model
end
=#
end # module end


text:
using Distributed
using Base.Filesystem
using DataFrames
using CSV
using Query
using Statistics
using UnicodePlots
using ClusterManagers
using Dates
using DelimitedFiles

## load the packages by covid19abm
using Parameters, Distributions, StatsBase, StaticArrays, Random, Match, DataFrames
using covid19abm
const cv=covid19abm

#addprocs(30, exeflags="--project=.")
#@everywhere using covid19abm

addprocs(SlurmManager(500), N=17, topology=:master_worker, exeflags="--project=.")
@everywhere using covid19abm

function run(myp::ModelParameters, nsims=500, folderprefix="./", writefiles = true)
    println("starting $nsims simulations...\nwriting files: $writefiles")
    #dump(myp)
    myp.calibration && error("can not run simulation, calibration is on.")
    # will return 6 dataframes. 1 total, 4 age-specific
    cdr = pmap(1:nsims) do x
            cv.runsim(x, myp)
    end

    println("simulations finished")
    println("total size of simulation dataframes: $(Base.summarysize(cdr))")
    ## write the infectors
    writefiles && writedlm("$(folderprefix)/infectors.dat", [cdr[i].infectors for i = 1:nsims])

    ## write contact numbers
    writefiles && writedlm("$(folderprefix)/ctnumbers.dat", [cdr[i].ct_numbers for i = 1:nsims])
    ## stack the sims together
    allag = vcat([cdr[i].a  for i = 1:nsims]...)
    ag1 = vcat([cdr[i].g1 for i = 1:nsims]...)
    ag2 = vcat([cdr[i].g2 for i = 1:nsims]...)
    ag3 = vcat([cdr[i].g3 for i = 1:nsims]...)
    ag4 = vcat([cdr[i].g4 for i = 1:nsims]...)
    ag5 = vcat([cdr[i].g5 for i = 1:nsims]...)
    #mydfs = Dict("all" => all, "ag1" => ag1, "ag2" => ag2, "ag3" => ag3, "ag4" => ag4, "ag5" => ag5)
    mydfs = Dict("all" => allag)

    ## save at the simulation and time level
    ## to ignore for now: miso, iiso, mild
    c1 = Symbol.((:LAT, :ASYMP, :MILD, :MISO, :INF, :IISO, :HOS, :ICU, :DED), :_INC)
    c2 = Symbol.((:LAT, :ASYMP, :MILD, :MISO, :INF, :IISO, :HOS, :ICU, :DED), :_PREV)
    for (k, df) in mydfs
        println("saving dataframe sim level: $k")
        # simulation level, save file per health status, per age group
        for c in vcat(c1..., c2...)
            udf = unstack(df, :time, :sim, c)
            fn = string("$(folderprefix)/simlevel_", lowercase(string(c)), "_", k, ".dat")
            writefiles && CSV.write(fn, udf)
        end
        println("saving dataframe time level: $k")
        # time level, save file per age group
        yaf = compute_yearly_average(df)
        fn = string("$(folderprefix)/timelevel_", k, ".dat")
        writefiles && CSV.write(fn, yaf)

        # write the secondary (Ro) values
        secvals = unstack(df, :time, :sim, :secondarys)
        fn = string("$(folderprefix)/secondarys_", k, ".dat")
        writefiles && CSV.write(fn, secvals)
    end
    return mydfs
end

function compute_yearly_average(df)
    ya = df |> @groupby(_.time) |> @map({time=key(_), cnt=length(_),
              sus_prev=mean(_.SUS_PREV),
              lat_prev=mean(_.LAT_PREV),
              pre_prev=mean(_.PRE_PREV),
              asymp_prev=mean(_.ASYMP_PREV),
              mild_prev=mean(_.MILD_PREV),
              miso_prev=mean(_.MISO_PREV),
              inf_prev=mean(_.INF_PREV),
              iiso_prev=mean(_.IISO_PREV),
              hos_prev=mean(_.HOS_PREV),
              icu_prev=mean(_.ICU_PREV),
              rec_prev=mean(_.REC_PREV),
              ded_prev=mean(_.DED_PREV),
              sus_inc=mean(_.SUS_INC),
              lat_inc=mean(_.LAT_INC),
              pre_inc=mean(_.PRE_INC),
              asymp_inc=mean(_.ASYMP_INC),
              mild_inc=mean(_.MILD_INC),
              miso_inc=mean(_.MISO_INC),
              inf_inc=mean(_.INF_INC),
              iiso_inc=mean(_.IISO_INC),
              hos_inc=mean(_.HOS_INC),
              icu_inc=mean(_.ICU_INC),
              rec_inc=mean(_.REC_INC),
              ded_inc=mean(_.DED_INC)
              }) |> DataFrame
    return ya
end

function savestr(p::ModelParameters, custominsert="simresults", customstart="")
    ## format: /data/covid19abm/$(custominsert)/$(customstart)_$(params)
    datestr = (Dates.format(Dates.now(), dateformat"mmdd_HHMM"))
    ## setup folder name based on model parameters
    taustr = replace(string(p.τmild), "." => "")
    fstr = replace(string(p.fmild), "." => "")
    rstr = replace(string(p.β), "." => "")
    prov = replace(string(p.prov), "." => "")
    eldr = replace(string(p.eldq), "." => "")
    eldqag = replace(string(p.eldqag), "." => "")
    fpreiso = replace(string(p.fpreiso), "." => "")
    tpreiso = replace(string(p.tpreiso), "." => "")
    fsev = replace(string(p.fsevere), "." => "")
    frelasymp = replace(string(p.frelasymp), "." => "")

    if p.ctstrat == 1
        strat = "strat1sc"
    elseif p.ctstrat == 3
        strat = "strat3s" * string(p.strat3qdays)
    else
        strat = "strat0"
    end
    pct = replace(string(p.fctcapture), "." => "")
    cct = replace(string(p.fcontactst), "." => "")
    idt = replace(string(p.cidtime), "." => "")
    tback = replace(string(p.cdaysback), "." => "")
    fldrname = "/data/covid19abm/$(custominsert)/$(customstart)_$(prov)_$(strat)_pct$(pct)_cct$(cct)_idt$(idt)_tback$(tback)_fsev$(fsev)_tau$(taustr)_fmild$(fstr)_q$(eldr)_qag$(eldqag)_relasymp$(frelasymp)_tpreiso$(tpreiso)_preiso$(fpreiso)/"
    mkpath(fldrname)
end

function _calibrate(nsims, myp::ModelParameters)
    myp.calibration != true && error("calibration parameter not turned on")
    vals = zeros(Int64, nsims)
    println("calibrating with beta: $(myp.β), total sims: $nsims, province: $(myp.prov)")
    println("calibration parameters:")
    #dump(myp)
    cdr = pmap(1:nsims) do i
        h = main(myp) ## gets the entire model.
        val = sum(_get_column_incidence(h, covid19abm.LAT))
        return val
    end
    return mean(cdr), std(cdr)
end

function calibrate(beta, nsims=1000, init_inf=1, size=10000, prov=:newyork)
    myp = ModelParameters() # set up default parameters
    myp.β = beta
    myp.prov = prov
    myp.popsize = size
    myp.modeltime = 50
    myp.calibration = true
    myp.initialinf = init_inf
    m, sd = _calibrate(nsims, myp)
    println("mean R0: $(m) with std: $(sd)")
    myp.calibration = false
    return m
end

function calibrate_regression()
    ## runs through a range of beta values and fits a linear regression line
    betas = 0.0225:0.0005:0.085
    ros = similar(betas)
    for i = 1:length(betas)
        ros[i] = calibrate(betas[i])
    end
    # run regression, myb is a dataframe of ros and betas
    #linearRegressor = lm(@formula(ros ~ betas), df)  ## formula expression requires dataframe
    regressed = lm(hcat(ones(length(betas)), betas), ros)
    return regressed
    # intercept = -0.00285 + 46.1761(beta) = ro
    # (ro + 0.00285)/46.1761
end

function calibrate_robustness(beta, reps, prov=:ontario)
    #[:ontario, :alberta, :bc, :manitoba, :newbruns, :newfdland, :nwterrito, :novasco, :nunavut, :pei, :quebec, :saskat, :yukon]
    # once a beta is found based on nsims simulations,
    # see how robust it is. run calibration with same beta 100 times
    # to see the variation in R0 produced.
    nsims = [500, 1000]
    means = zeros(Float64, reps, length(nsims))
    for (i, ns) in enumerate(nsims)
        cd = map(1:reps) do x
            println("iter: $x, sims: $ns")
            mval = calibrate(beta, ns, prov)
            return mval
        end
        means[:, i] = cd
    end
    # for i in 2:nworkers()
    #     ## mf defined as: @everywhere mg() = covid19abm.p.β
    #     rpr = remotecall_fetch(mf,  i+1).prov
    #     rpr != prov && error("province didn't get set in the remote workers")
    # end
    return means
end
#provs = []
#for x in [:ontario, :alberta, :bc, :manitoba, :newbruns]
#    m = calibrate_robustness(0.101, x)
#    push!(provs, m)
#end

#@everywhere mg() = covid19abm.p.β
#@everywhere mf(b) = covid19abm.p.β = b
# function hmatrix_analysis()
#     ## function takes the output of the model
#     ## and splits it into an array of 2d matrices for plotting purposes.
#     myp = cv.ModelParameters()
#     myp.β = 0.038
#     reset_params($myp)
#     cd, ags = main(1)

#     #mats = Array{Int64, 3}(undef, 500)
#     mats = Array{Int64, 3}(undef, 100,100,500)
#     for i = 1:500
#         col = cd[:, i]
#         mats[:, :, i] = reshape(col, (100, 100))
#     end

#     colorarray = IndirectArray(mats, [colorant"green", colorant"teal", colorant"blue", colorant"azure",
#     colorant"red", colorant"lightsalmon", colorant"magenta", colorant"magenta", colorant"grey60", colorant"black"])
#     save("myanimation.gif", colorarray; fps=5)

#     idxarray = rand(1:3, 100, 100, 50);
#     colorarray = IndirectArray(idxarray, [colorant"red", colorant"purple", colorant"indigo"])

#     return mats
# end

# function makeheatmaps(bufs)
#     heatmaps = map(bufs) do buf
#         heatmap(
#             buf, scale_plot = false, show_axis = false)

#     end
#     scene = hbox(map(i-> vbox(heatmaps[i, :]), 1:size(bufs, 1))...)
#     scene, last.(heatmaps)
# end
#datarows = 500; datacols = 500
#plotrows = 4; plotcols = 4
#bufs = [fill(0.0f0, datarows, datacols) for _ in 1:plotrows, _ in 1:plotcols]

# http://juliaplots.org/MakieReferenceImages/gallery//lots_of_heatmaps/index.html
# http://juliaplots.org/MakieReferenceImages/gallery//chess_game/index.html

ext: jl


text:
"""
    evaluate_field_array(f::AbstractArray,x::AbstractArray) -> AbstractArray

Evaluates the fields in the array `f` at all the vector of points in the
array of vector of points `x` and returns the result as a lazy array.

The result is numerically equivalent to

    map(evaluate_field,a,x)
"""
function evaluate_field_array(a::AbstractArray,x::AbstractArray)
  _evaluate_field_array(a,x)
end

function _evaluate_field_array(a::AbstractArray,x::AbstractArray)
  k = Eval()
  apply(k,a,x)
end

struct Eval <: Kernel end

function kernel_cache(k::Eval,a,x)
  field_cache(a,x)
end

function apply_kernel!(cache,k::Eval,a,x)
  evaluate_field!(cache,a,x)
end

function kernel_return_type(k::Eval,a,x)
  field_return_type(a,x)
end

function _evaluate_field_array(a::AbstractArray{<:Field},x::AbstractArray)
  apply(a,x)
end

"""
    evaluate_field_arrays(f::Tuple,x::AbstractArray) -> Tuple

Equivalent to

    tuple((evaluate_field_array(fi,x) for fi in f)...)
"""
function evaluate_field_arrays(f::Tuple,x::AbstractArray)
  _evaluate_field_arrays(x,f...)
end

function _evaluate_field_arrays(x,a,b...)
  ax = evaluate_field_array(a,x)
  bx = evaluate_field_arrays(b,x)
  (ax,bx...)
end

function _evaluate_field_arrays(x,a)
  ax = evaluate_field_array(a,x)
  (ax,)
end

"""
    evaluate(a::AbstractArray{<:Field},x::AbstractArray)

Equivalent to

    evaluate_field_array(a,x)

But only for arrays `a` whose element type inherits from `Field`. If this is not the case,
use `evaluate_field_array(a,x)` instead.
"""
function evaluate(a::AbstractArray{<:Field},x::AbstractArray)
  evaluate_field_array(a,x)
end

# Optimized version for arrays of fields obtained from a kernel
# and other arrays
function evaluate_field_array(
  a::AppliedArray{T,N,F,<:Fill} where {T,N,F},x::AbstractArray)
  kernel_evaluate(a.g.value,x,a.f...)
end

"""
    kernel_evaluate(k,x,f...)

This function returns by default the array obtained with

    g = apply_to_field_array(k,f...)
    evaluate_field_array(g,x)

However, it can be rewritten for specific kernels in order to improve performance and simplify
the underlying operation tree associated in the returned object.

"""
function kernel_evaluate(k,x,f...)
  a = apply(k,f...)
  _evaluate_field_array(a,x)
end

"""
    field_array_gradient(a::AbstractArray)

Returns an array containing the gradients of the fields in the array `a`.
Numerically equivalent to

    map(field_gradient,a)
"""
function field_array_gradient(a::AbstractArray)
  _field_array_gradient(a)
end

function _field_array_gradient(a::AbstractArray)
  k = Grad()
  apply(k,a)
end

struct Grad <: Kernel end

@inline apply_kernel!(::Nothing,k::Grad,x) = field_gradient(x)

function field_array_gradient(
  a::AppliedArray{T,N,F,<:Fill} where {T,N,F})
  apply_gradient(a.g.value,a.f...)
end

"""
    apply_gradient(k,f...)

By default, it returns the array obtained as

    a = apply(k,f...)
    field_array_gradient(a)

However, it can be rewritten for specific kernels in order to improve performance and simplify
the underlying operation tree associated in the returned object.

"""
function apply_gradient(k,f...)
  a = apply(k,f...)
  _field_array_gradient(a)
end

"""
    gradient(f::AbstractArray{<:Field})

Equivalent to

    field_array_gradient(f)

but only for arrays whose element type is `<:Field`. Use function `field_array_gradient` otherwise.
"""
function gradient(f::AbstractArray{<:Field})
  field_array_gradient(f)
end

"""
    field_array_gradients(f...)

Equivalent to

    map(field_array_gradient,f)
"""
function field_array_gradients(f...)
  map(field_array_gradient,f)
end

"""
    field_array_cache(a::AbstractArray,x::AbstractArray) -> Tuple

Returns the caches needed to perform the following iteration

    ca, cfi, cx = field_array_cache(a,x)

    for i in length(a)
      fi = getindex!(ca,a,i)
      xi = getindex!(cx,x,i)
      fxi = evaluate!(cfi,fi,xi)
    end
"""
function field_array_cache(a::AbstractArray,x::AbstractArray)
  ca = array_cache(a)
  fi = testitem(a)
  xi = testitem(x)
  cfi = field_cache(fi,xi)
  cx = array_cache(x)
  (ca,cfi,cx)
end

"""
    function test_array_of_fields(
      a::AbstractArray,
      x::AbstractArray,
      v::AbstractArray,
      cmp::Function=(==);
      grad = nothing)

Function to test an array of fields `a`. The array `v` is the expected result when calling
`evaluate_field_array(a,x)`. The entries in the computed array and the expected one are compared
with the `cmp` function. The key-word argument `grad` is optional. If present, it should contain
the expected result of

    ∇a = field_array_gradient(a)
    evaluate_field_array(∇a,x)
"""
function test_array_of_fields(
  a::AbstractArray,
  x::AbstractArray,
  v::AbstractArray,
  cmp::Function=(==);
  grad = nothing)

  ax = evaluate_field_array(a,x)
  test_array(ax,v,cmp)

  ca, cfi, cx = field_array_cache(a,x)

  t = true
  for i in 1:length(a)
    fi = getindex!(ca,a,i)
    xi = getindex!(cx,x,i)
    fxi = evaluate_field!(cfi,fi,xi)
    vi = v[i]
    ti = cmp(fxi,vi)
    t = t && ti
  end
  @test t

  if grad != nothing
    g = field_array_gradient(a)
    test_array_of_fields(g,x,grad,cmp)
  end

end

"""
    apply_to_field_array(k,f::AbstractArray...)

Returns an array of fields numerically equivalent to

    map( (x...) -> apply_kernel_to_field(k,x...), f )

The evaluation and the computation of the gradient of the resulting arrays
can be optimized for particular kernels by
adding new methods to the [`kernel_evaluate`](@ref) and  [`apply_gradient`](@ref) functions
respectively.
"""
function apply_to_field_array(
  k,f::AbstractArray...)
  v = Valued(k)
  apply(v,f...)
end

struct Valued{K} <: Kernel
  k::K
  function Valued(k)
    new{typeof(k)}(k)
  end
end

@inline function apply_kernel!(cache,k::Valued,x...)
  apply_kernel_to_field(k.k,x...)
end

function kernel_evaluate(k::Valued,x,f...)
  fx = evaluate_field_arrays(f,x)
  a = apply(k.k,fx...)
end


text:
function construct( e::EssentialMatrix, 𝐅::AbstractArray,  𝐊₁::AbstractArray, 𝐊₂::AbstractArray)
    if (size(𝐊₁) != (3,3)) || (size(𝐊₂) != (3,3))
        throw(ArgumentError("Expect 3 x 3 calibration matrices."))
    end
    if (size(𝐅) != (3,3))
        throw(ArgumentError("Expect 3 x 3 fundamental matrix."))
    end
    # Equation 9.12 Chapter 9 from Hartley & Zisserman
    𝐄 = 𝐊₂'*𝐅*𝐊₁
    SMatrix{3,3,Float64,3*3}(𝐄)
end


text:
module LagrangianDofBasesTests

using Test

using Gridap.TensorValues
using Gridap.Fields
using Gridap.ReferenceFEs

x1 = Point(0,0)
x2 = Point(1,0)
x3 = Point(0,1)
x4 = Point(1,1)
x = [x1,x2,x3,x4]
np = length(x)

db = LagrangianDofBasis(Float64,x)

@test db.nodes === x
@test db.node_and_comp_to_dof == collect(1:np)
@test db.dof_to_node == collect(1:np)
@test db.dof_to_comp == fill(1,np)

v = 3.0
f = GenericField(x->v*x[1])
fx = evaluate(f,x)
test_dof_array(db,f,fx)

ndof = 8
b = fill(f,ndof)
bx = evaluate(b,x)
test_dof_array(db,b,bx)

T = VectorValue{2,Float64}
db = LagrangianDofBasis(T,x)
@test db.nodes === x
@test db.node_and_comp_to_dof == VectorValue{2,Int}[(1, 5), (2, 6), (3, 7), (4, 8)]
@test db.dof_to_node == [1, 2, 3, 4, 1, 2, 3, 4]
@test db.dof_to_comp == [1, 1, 1, 1, 2, 2, 2, 2]

v = VectorValue(1,2)
f = GenericField(x->v*x[1])
dbf = [0, 1, 0, 1, 0, 2, 0, 2]

test_dof_array(db,f,dbf)

ndof = 8
b = fill(f,ndof)
bx = evaluate(b,x)
dbb = [
  0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1; 0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1;
  0 0 0 0 0 0 0 0; 2 2 2 2 2 2 2 2; 0 0 0 0 0 0 0 0; 2 2 2 2 2 2 2 2]

test_dof_array(db,b,dbb)

end # module


text:
# test the functions in src/bootstrap.jl

exdir = joinpath(@__DIR__,"..","examples")

@testset "testing hybridBootstrapSupport" begin
bestnet = readTopology(joinpath(exdir,"fish2hyb.net"));
bootnet = readMultiTopology(joinpath(exdir,"fish3hyb_20boostrap.net"));
# issues with bootstrap networks 12, 21, 42, 96
# plot(bootnet[20], showEdgeNumber=true)
# include(string(home, "bootstrap.jl"))
resn, rese, resc, gam, edgenum = hybridBootstrapSupport(bootnet,bestnet);
#@show resn; @show rese; showall(gam); @show edgenum; resc
# plot(bestnet, showIntNodeNumber=true)

@test resn[:clade][1:2] == ["H26","H25"]
@test resn[:BS_hybrid_samesisters][1:2] == [25.0,100.0]
@test resn[:BS_hybrid] == [100.0,100,0,0,0,75,0,0,0,0,0,0,5,5,5,5,5,0,0,0]
@test resn[:BS_minor_sister] == [0.0,0,100,0,0,5,10,70,75,25,5,5,0,0,0,0,0,0,0,5]
@test resn[:BS_major_sister] == [0.0,0,0,100,100,0,70,10,0,0,5,5,0,0,0,0,0,5,5,0]
@test rese[:BS_minor][2] == 25.0  # BS of introgression for H26
@test rese[:BS_minor][4] == 100.0 # BS of introgression for H25
@test resc[:taxa]==["Xgordoni","Xmeyeri","Xcouchianus","Xvariatus","Xevelynae","Xxiphidium",
 "Xmilleri","Xandersi","Xmaculatus","Xhellerii","Xalvarezi","Xmayae","Xsignum","Xclemenciae_F2",
 "Xmonticolus","Xmontezumae","Xnezahuacoyotl","Xbirchmanni_GARC","Xmalinche_CHIC2","Xcortezi",
 "Xcontinens","Xpygmaeus","Xnigrensis","Xmultilineatus"]
@test resc[:taxa][resc[:H26]] == ["Xnezahuacoyotl"]
@test resc[:taxa][resc[:H25]] == ["Xmontezumae","Xnezahuacoyotl","Xbirchmanni_GARC","Xmalinche_CHIC2","Xcortezi","Xcontinens","Xpygmaeus","Xnigrensis","Xmultilineatus"]
@test resc[:taxa][resc[:c_minus27]] == ["Xnigrensis","Xmultilineatus"] # minor sis of H26
@test resc[:taxa][resc[:Xxiphidium]] == ["Xxiphidium"] # minor sis of H25
@test resc[:taxa][resc[:Xsignum]] == ["Xsignum"] # donor8 previously
@test resc[:taxa][resc[:c_minus24]] == ["Xcontinens","Xpygmaeus","Xnigrensis","Xmultilineatus"] # donor7
@test resc[:taxa][resc[:Xmontezumae]] == ["Xmontezumae"] # major sis of H26. Below: major sis of H25
@test resc[:taxa][resc[:c_minus12]] == ["Xhellerii","Xalvarezi","Xmayae","Xsignum","Xclemenciae_F2","Xmonticolus"]
@test gam[:,2] == [.0,.0,.192,.0,.0,.0,.0,.0,.193,.0,.184,.193,.0,.0,.0,.0,.0,.193,.0,.0]
@test gam[:,4] == [.165,.166,.165,.166,.165,.165,.166,.165,.165,.166,.164,.166,.166,.165,.165,.165,.166,.165,.166,.166]
@test edgenum ==[25,39,43,7]
end # of testset, hybridBootstrapSupport

# exdir = ""
info("testing bootsnaq from quartet CF intervals")
T=readTopology(joinpath(exdir,"startTree.txt"))
df=CSV.read(joinpath(exdir,"tableCFCI.csv"))
net1 = readTopology("(5,(((2,(1)#H7:::0.7143969494428192):1.5121337017411736,4):0.4894187322508883,3):0.519160762355313,(6,#H7:::0.2856030505571808):10.0);")
bootnet = bootsnaq(T,df,nrep=2,runs=1,seed=1234,filename="",Nfail=2,ftolAbs=1e-3,ftolRel=1e-3,
                   xtolAbs=1e-4,xtolRel=1e-3,liktolAbs=0.01)
@test size(bootnet)==(2,)
@test writeTopology(bootnet[1], round=true)=="(4,(((6,#H7:0.0::0.222):10.0,2):0.0,((5,1):0.0)#H7:0.0::0.778):0.0,3);"
@test writeTopology(bootnet[2], round=true)=="(2,3,((4,5):0.0,(1,6):0.011):0.0);"
# "(2,((5,#H9:0.0::0.298):3.927,3):1.331,(((1,6):0.019,4):0.0)#H9:0.0::0.702);"
# above: bad diamond 2, and both edges above the hybrid have estimated length of 0.0...

@testset "bootsnaq from bootstrap gene trees, multiple procs" begin
treefile = joinpath(exdir,"treefile.txt") # pretending these are bootstrap trees, for all genes
boottrees = Vector{HybridNetwork}[]
for i=1:13 push!(boottrees, readMultiTopology(treefile)) end
for i=1:13 @test size(boottrees[i])==(10,) end # 10 bootstrap trees for each of 13 "genes"
addprocs(1)
@everywhere using PhyloNetworks
bootnet = bootsnaq(T,boottrees,nrep=2,runs=2,otherNet=net1,seed=1234,prcnet=0.5,filename="",Nfail=2,ftolAbs=1e-3,ftolRel=1e-3)
rmprocs(workers())
@test size(bootnet)==(2,)
@test writeTopology(bootnet[1], round=true)=="((((2,(1)#H7:::0.597):10.0,4):0.407,(6,#H7:::0.403):0.307):0.0,3,5);"
# "((5,((2,(1)#H7:::0.629):2.374,4):0.487):0.0,(6,#H7:::0.371):1.409,3);"
# "((((2,(1)#H7:::0.678):1.774,4):0.235,3):0.899,5,(6,#H7:::0.322):10.0);"
@test writeTopology(bootnet[2], round=true)=="(5,(((2,(1)#H7:::0.691):2.881,4):1.444,3):0.0,(6,#H7:::0.309):0.198);"
# "(5,(((2,(1)#H7:::0.751):1.559,4):0.373,3):0.688,(6,#H7:::0.249):10.0);"
# "(((5,(6,#H7:::0.249):10.0):0.688,3):0.373,(2,(1)#H7:::0.751):1.559,4);"
end


text:
export SDL, sdl
export getSDLError, getWindowSize, window_paused
export to_ARGB, sdl_colors, clear

mutable struct SDL <: System
  win # window
  rnd # renderer
  bgrd # background
  wdth # width
  hght # height
  ttl # title
  function SDL()
    return new(nothing, nothing, nothing, nothing, nothing, nothing)
  end
end

const sdl = SDL()

listenFor(sdl, TICK)

function draw()
  clear()
  map(draw, scn) # TODO investigate parallelization
  SDL_RenderPresent(sdl.rnd)
end

function handleMessage(s::SDL, m::TICK)
  @debug "SDL tick"
  try
    # handle events
    event_ref = Ref{SDL_Event}()
    while Bool(SDL_PollEvent(event_ref))
      evt = event_ref[]
      sendMessage(evt)
    end

    # draw the scene
    draw()
  catch
    handleException()
  end

end

function getSDLError()
  x = SDL_GetError()
  return unsafe_string(x)
end

# Forward reference for @cfunction
function windowEventWatcher end
const window_event_watcher_cfunc = Ref(Ptr{Nothing}(0))

const window_paused = Threads.Atomic{UInt8}(0) # Whether or not the game should be running (if lost focus)

function makeWinRenderer(title::String, w::Int, h::Int)
    win = SDL_CreateWindow(title,
        Int32(SDL_WINDOWPOS_CENTERED), Int32(SDL_WINDOWPOS_CENTERED), Int32(w), Int32(h),
        UInt32(SDL_WINDOW_ALLOW_HIGHDPI|SDL_WINDOW_OPENGL|SDL_WINDOW_SHOWN));
        SDL_SetWindowMinimumSize(win, Int32(w), Int32(h))
    window_event_watcher_cfunc[] = @cfunction(windowEventWatcher, Cint, (Ptr{Nothing}, Ptr{SDL_Event}))
    SDL_AddEventWatch(window_event_watcher_cfunc[], win);

    renderer = SDL_CreateRenderer(win, Int32(-1), UInt32(SDL_RENDERER_ACCELERATED | SDL_RENDERER_PRESENTVSYNC))
    SDL_SetRenderDrawBlendMode(renderer, SDL_BLENDMODE_BLEND)
    return win,renderer
end

# This function handles all window events.
# We currently do no allow window resizes
function windowEventWatcher(data_ptr::Ptr{Cvoid}, event_ptr::Ptr{SDL_Event})::Cint
    ev = unsafe_load(event_ptr, 1)
    #ee = evt.type
    #t = UInt32(ee[4]) << 24 | UInt32(ee[3]) << 16 | UInt32(ee[2]) << 8 | ee[1]
    t = ev.type
    if (t == SDL_WindowEvent)
        event = unsafe_load( Ptr{SDL_WindowEvent}(pointer_from_objref(ev)) )
        winevent = event.event;  # confusing, but that's what the field is called.
        if (winevent == SDL_WINDOWEVENT_FOCUS_LOST || winevent == SDL_WINDOWEVENT_HIDDEN || winevent == SDL_WINDOWEVENT_MINIMIZED)
            # Stop game playing when out of focus
                window_paused[] = 1
            #end
        elseif (winevent == SDL_WINDOWEVENT_FOCUS_GAINED || winevent == SDL_WINDOWEVENT_SHOWN)
            window_paused[] = 0
        end
    end
    return 0
end

function getWindowSize()
    w,h,w_highDPI,h_highDPI = Int32[0],Int32[0],Int32[0],Int32[0]
    SDL_GetWindowSize(sdl.win, w, h)
    SDL_GL_GetDrawableSize(sdl.win, w_highDPI, h_highDPI)
    return w[],h[],w_highDPI[],h_highDPI[]
end

to_ARGB(c) = c
to_ARGB(c::ARGB) = c
to_ARGB(c::Colorant) = ARGB(c)

clear() = fill(sdl.bgrd)

sdl_colors(c::Colorant) = sdl_colors(convert(ARGB{Colors.FixedPointNumbers.Normed{UInt8,8}}, c))
sdl_colors(c::ARGB) = Int.(reinterpret.((red(c), green(c), blue(c), alpha(c))))

function Base.fill(c::Colorant)
  SDL_SetRenderDrawColor(
      sdl.rnd,
      sdl_colors(c)...,
  )
  SDL_RenderClear(sdl.rnd)
end

function Base.fill(sf::Ptr{SDL_Surface})
  texture = SDL_CreateTextureFromSurface(sdl.rnd, sf)
  SDL_RenderCopy(sdl.rnd, texture, C_NULL, C_NULL)
  SDL_DestroyTexture(texture)
end

function awake!(s::SDL)
  SDL_GL_SetAttribute(SDL_GL_MULTISAMPLEBUFFERS, 4)
  SDL_GL_SetAttribute(SDL_GL_MULTISAMPLESAMPLES, 4)
  r = SDL_Init(SDL_INIT_EVERYTHING)
  if r != 0
    error("unable to initialize SDL: $(getSDLError())")
  end

  TTF_Init()

  mix_init_flags = MIX_INIT_FLAC|MIX_INIT_MP3|MIX_INIT_OGG
  inited = Mix_Init(Int32(mix_init_flags))
  if inited & mix_init_flags != mix_init_flags
      @warn "Failed to initialise audio mixer properly. All sounds may not play correctly\n$(getSDLError())"
  end

  device = Mix_OpenAudio(Int32(22050), UInt16(MIX_DEFAULT_FORMAT), Int32(2), Int32(1024) )
  if device != 0
      @warn "No audio device available, sounds and music will not play.\n$(getSDLError())"
      Mix_CloseAudio()
  end

  # other fields set inside App constructor
  win, rnd = makeWinRenderer(sdl.ttl, sdl.wdth, sdl.hght)
  sdl.win = win
  sdl.rnd = rnd

  draw()

  return true
end

function shutdown!(s::SDL)
  SDL_DelEventWatch(window_event_watcher_cfunc[], sdl.win);
  SDL_DestroyRenderer(sdl.rnd)
  SDL_DestroyWindow(sdl.win)

  Mix_HaltMusic()
  Mix_HaltChannel(Int32(-1))
  Mix_CloseAudio()

  TTF_Quit()
  Mix_Quit()
  SDL_Quit()

  return false
end


text:
using Gadfly
using DataFrames

m_p = 1.6726218e-27 #Proton mass in kg
MeV2J = 1.6021766e-13 #1 MeV in J
eV2J = 1.6021766e-19 #1 eV in J
c = 2.998e8 #Speed of Light in a Vacuum in m/s
debug = false
Emin = 3.0e5 #Minimum energy in eV
exponent = 2.6

function j(E)
  Eg = E/1.0e9
  factor1 = 0.90/(0.85 + Eg)^exponent
  factor2 = 1.0/(1.0 + 0.01/Eg)
  return 4*π*factor1*factor2
end

function β_i(E)
  EJ = E*eV2J
  return sqrt((2.0*EJ)/m_p)/c
end

function σ_i(E)
  β =  β_i(E)
  factor1 = (1.23e-20)/(β^2)
  logfactor = log10((β^2)/(1.0-β^2))
  factor2 = 6.20 + logfactor - 0.43*β^2
  return factor1*factor2
end

function integrand(E)
  return j(E)*σ_i(3.0e5)
end

energies = Emin:5.0e7:1.0e10

flux = Array(Float64,size(energies,1))

for i in 1:size(energies,1)
  flux[i] = j(energies[i])
end

p1 = plot(x=energies/1.0e9,y=flux,Geom.line,Scale.x_log10,Scale.y_log10)

intFlux = quadgk(integrand,Emin,Inf)
intFlux = 2.0*(5.0/3.0)*intFlux[1]

println("Integrated Flux=$(intFlux)")

println("Ending script")


text:
using Base.Broadcast: BroadcastStyle, AbstractArrayStyle, DefaultArrayStyle, Broadcasted

struct StridedArrayStyle{N} <: AbstractArrayStyle{N}
end

Broadcast.BroadcastStyle(::Type{<:AbstractStridedView{<:Any,N}}) where {N} =
    StridedArrayStyle{N}()

StridedArrayStyle(::Val{N}) where {N} = StridedArrayStyle{N}()
StridedArrayStyle{M}(::Val{N}) where {M,N} = StridedArrayStyle{N}()

Broadcast.BroadcastStyle(a::StridedArrayStyle, ::DefaultArrayStyle{0}) = a
Broadcast.BroadcastStyle(::StridedArrayStyle{N}, a::DefaultArrayStyle) where {N} =
    BroadcastStyle(DefaultArrayStyle{N}(), a)
Broadcast.BroadcastStyle(::StridedArrayStyle{N}, ::Broadcast.Style{Tuple}) where {N} =
    DefaultArrayStyle{N}()

function Base.similar(bc::Broadcasted{<:StridedArrayStyle{N}}, ::Type{T}) where {N,T}
    StridedView(similar(convert(Broadcasted{DefaultArrayStyle{N}}, bc), T))
end

Base.dotview(a::AbstractStridedView{<:Any,N}, I::Vararg{SliceIndex,N}) where {N} =
    getindex(a, I...)

# Broadcasting implementation
@inline function Base.copyto!(dest::AbstractStridedView{<:Any,N},
                                bc::Broadcasted{StridedArrayStyle{N}}) where {N}
    # convert to map

    # flatten and only keep the AbstractStridedView arguments
    # promote AbstractStridedView to have same size, by giving artificial zero strides
    stridedargs = promoteshape(size(dest), capturestridedargs(bc)...)
    c = make_capture(bc)
    _mapreduce_fuse!(c, nothing, nothing, size(dest), (dest, stridedargs...))
    return dest
end

const WrappedScalarArgs = Union{AbstractArray{<:Any,0}, Ref{<:Any}}

@inline capturestridedargs(t::Broadcasted, rest...) =
    (capturestridedargs(t.args...)..., capturestridedargs(rest...)...)
@inline capturestridedargs(t::AbstractStridedView, rest...) =
    (t, capturestridedargs(rest...)...)
@inline capturestridedargs(t, rest...) = capturestridedargs(rest...)
@inline capturestridedargs() = ()

# broadcast by promoting size 1 dimensions to full size
# making the stride in that dimension zero
promoteshape(sz::Dims, a1::AbstractStridedView, As...) = (promoteshape1(sz, a1), promoteshape(sz, As...)...)
promoteshape(sz::Dims) = ()
function promoteshape1(sz::Dims{N}, a::StridedView) where {N}
    newstrides = ntuple(Val(N)) do d
        if size(a, d) == sz[d]
            stride(a, d)
        elseif size(a, d) == 1
            0
        else
            throw(DimensionMismatch("array could not be broadcasted to match destination"))
        end
    end
    return StridedView(a.parent, sz, newstrides, a.offset, a.op)
end
function promoteshape1(sz::Dims{N}, a::UnsafeStridedView) where {N}
    newstrides = ntuple(Val(N)) do d
        if size(a, d) == sz[d]
            stride(a, d)
        elseif size(a, d) == 1
            0
        else
            throw(DimensionMismatch("array could not be broadcast to match destination"))
        end
    end
    return UnsafeStridedView(a.ptr, sz, newstrides, a.offset, a.op)
end

struct CaptureArgs{F, Args<:Tuple}
    f::F
    args::Args
end
struct Arg
end

# construct CaptureArgs
@inline function make_capture(bc::Broadcasted)
    args = make_tcapture(bc.args)
    CaptureArgs(bc.f, args)
end
@inline make_tcapture(t::Tuple{}) = t
@inline make_tcapture(t::Tuple) = (make_capture(t[1]), make_tcapture(tail(t))...)
@inline make_capture(a::WrappedScalarArgs) = a[]
@inline make_capture(a::AbstractStridedView) = Arg()
@inline make_capture(a) = a

# Evaluate CaptureArgs
(c::CaptureArgs)(vals...) = consume(c, vals)[1]
@inline function consume(c::CaptureArgs{F,Args}, vals) where {F,Args}
    args, newvals = t_consume(c.args, vals)
    return c.f(args...), newvals
end
@inline consume(a::Arg, vals::Tuple) = vals[1], tail(vals)
@inline consume(a, vals) = a, vals
@inline t_consume(t::Tuple{}, vals) = t, vals
@inline function t_consume(t::Tuple, vals)
    t1, newvals1 = consume(t[1], vals)
    ttail, newvals = t_consume(tail(t), newvals1)
    return (t1, ttail...), newvals
end


text:
"""
Compute pair warping between two functions using Bayesian method

    pair_warping_expomap(f1, f2, timet; iter=20000, burnin=min(5000,iter/2),
                         alpha0=0.1, beta0=0.1, pbetas=[0.5,0.05,0.005,0.0001],
                         probs=[0.1,0.1,0.7,0.1],propvar=1.0,
                         init_coef=zeros(20),npoints=200,extrainfo=false)

    :param f1, f2: vectors describing functions
    :param timet: vector describing timing
    :param iter: number of MCMC iterations
    :param burnin: number of MCMC burnin iterations
    :param alpha0, beta0: IG parameters for prior of sigma1
    :param pbetas: mixture ratios for pCN
    :param probs: mixcture probabilities for zpCN
    :param propvar: variance of proposal distribution
    :param init_coef: initial g coefficients
    :param npoits: number of sample points to use during alignment
    :param extrainfo: T/F whether additional information is returned (accept, logl, gamma_mat, gamma_stats, xdist, ydist)

    Returns mcmc struct containing
    :return f2_warped: warped f2
    :return gamma: warping function
    :return g_coef: g_coeficients
    :return sigma1: variance samples
    :return accept: accept samples
    :return logl: log-likelihood
    :return gamma_mat: posterior warping function samples
    :return gamma_stats: posterior warping function samples 95% credible intervals
    :return xdist: phase distance of posterior warping functions
    :return ydist: amplitude distance of posterior warping functions
"""
function pair_warping_expomap(f1, f2, timet; iter=20000, burnin=min(5000,iter/2),
                              alpha0=0.1, beta0=0.1, pbetas=[0.8,0.1,0.005,0.0001],
                              probs=[0.1,0.1,0.7,0.1],propvar=1.0,init_coef=zeros(20),
                              npoints=200,extrainfo=false)

    # check inputs
    if (length(f1)!=length(f2))
        error("Length of f1 and f2 must be equal")
    end
    if (length(f1)!=length(timet))
        error("Length of f1 and f2 must be equal")
    end
    if (length(pbetas)!=length(probs))
        error("In zpcn, betas must equal length of probs")
    end
    if ((length(init_coef) % 2) !=0)
        error("Length of init_cef must be even")
    end

    SIG_GAM = 13
    Mg = length(init_coef)
    sigma1_ini = 1.0
    g_coef_ini = init_coef
    nbasis = Int(Mg/2)

    # normalize time
    timet = collect(LinRange(0,1,length(timet)))
    f1 = func(timet,f1)
    f2 = func(timet,f2)

    g_basis = basis_fourier(collect(LinRange(0,1,npoints)), nbasis, 1)

    function propose_g_coef(g_coef_curr)
        probm = [0; cumsum(probs)]
        z = rand()
        d = MvNormal(propvar./sort(repeat(1:nbasis,2)))
        ind = 0
        prop = g_coef_curr
        for i in 1:length(pbetas)
            if (z <= probm[i+1] && z > probm[i])
                g_coef_new = rand(d)
                prop = sqrt(1 - pbetas[i] ^ 2) * g_coef_curr + pbetas[i] * g_coef_new
                ind = i
            end
        end
        return(prop,ind)
    end

    # srsf transformation
    q1 = f_to_srsf(f1)
    q2 = f_to_srsf(f2)

    obs_domain = q1.x
    sim_domain = collect(LinRange(0,1,npoints))
    valid_index = Int(burnin):iter

    g_temp = f_basistofunction(g_basis.x,g_coef_ini,g_basis)
    tmp = f_exp1(g_temp)
    if (minimum(tmp.y)<0)
        error("Invalid initial value of g")
    end

    # result variables
    g_coef = Array{Float64}(undef,length(g_coef_ini),iter)
    sigma1 = Vector{Float64}(undef,iter)
    logl = Vector{Float64}(undef,iter)
    SSE = Vector{Float64}(undef,iter)
    accept = Vector{Bool}(undef,iter).*false
    accept_betas = Vector{Int64}(undef,iter)

    g_coef_curr = g_coef_ini
    sigma1_curr = sigma1_ini

    g = f_basistofunction(g_basis.x, g_coef_curr, g_basis)
    SSE_curr = f_SSEg_pw(g, q1, q2)
    logl_curr = f_logl_pw(g, sigma1_curr^2, q1, q2, SSEg=SSE_curr)

    g_coef[:,1] = g_coef_curr
    sigma1[1] = sigma1_curr
    SSE[1] = SSE_curr
    logl[1] = logl_curr
    accept[1] = true

    @showprogress 1 "Computing MCMC..." for m in 2:iter
        # update g
        g_coef_curr, SSE_curr, logl_curr, accepti = f_updateg_pw(g_coef_curr, g_basis,
                                                                sigma1_curr^2, q1, q2,
                                                                SSE_curr, propose_g_coef)

        # update sigma1
        newshape = length(q1.x) / 2 + alpha0
        newscale = 1 / 2 * SSE_curr + beta0
        d = InverseGamma(newshape, newscale)
        sigma1_curr = sqrt(rand(d))
        g1 = f_basistofunction(g_basis.x, g_coef_curr, g_basis)
        logl_curr = f_logl_pw(g1, sigma1_curr^2, q1, q2, SSEg=SSE_curr)

        g_coef[:,m] = g_coef_curr
        sigma1[m] = sigma1_curr
        SSE[m] = SSE_curr
        if (extrainfo)
            logl[m] = logl_curr
            accept[m] = accepti
        end
    end

    # calculate posterior mean of psi
    pw_sim_est_psi_matrix = Array{Float64}(undef,npoints,length(valid_index))
    cnt = 1
    for k in valid_index
        g_temp = f_basistofunction(sim_domain, g_coef[:,k], g_basis)
        psi_temp = f_exp1(g_temp)
        pw_sim_est_psi_matrix[:,cnt] = psi_temp.y
        cnt += 1
    end

    result_posterior_psi_simDomain = f_psimean(sim_domain,pw_sim_est_psi_matrix)

    # resample to same number of points
    result_posterior_psi = approx(sim_domain, result_posterior_psi_simDomain.y, q1.x)

    # transform to γ
    psi_tmp = func(obs_domain, result_posterior_psi)
    result_posterior_gamma = f_phiinv(psi_tmp)
    gam0 = result_posterior_gamma.y
    result_posterior_gamma = func(obs_domain,norm_gam(gam0))

    # warped f_2
    f2_warped = warp_f_gamma(result_posterior_gamma.x, f2.y, result_posterior_gamma.y)

    Dx = []
    Dy = []
    gamma_mat = []
    gamma_stats = []
    if (extrainfo)
        gamma_mat = zeros(length(obs_domain),size(pw_sim_est_psi_matrix,2))
        Dx = Vector{Float64}(undef, size(pw_sim_est_psi_matrix,2))
        Dy = Vector{Float64}(undef, size(pw_sim_est_psi_matrix,2))
        one_v = ones(size(pw_sim_est_psi_matrix,1))
        for i = 1:size(pw_sim_est_psi_matrix,2)
            tmp = approx(sim_domain, pw_sim_est_psi_matrix[:,i], obs_domain)
            tmp = func(obs_domain,tmp)
            gam0 = f_phiinv(tmp)
            gamma_mat[:,i] = norm_gam(gam0.y)

            v = inv_exp_map(one_v, pw_sim_est_psi_matrix[:,i])
            Dx[i] = sqrt(trapz(sim_domain, v.^2))
            q2warp = warp_q_gamma(q2.x, q2.y, gamma_mat[:,i])
            Dy[i] = sqrt(trapz(q2.x, (q1.y-q2warp).^2))
        end
        gamma_stats = statsFun(gamma_mat)
    end

    # return object
    out = mcmc_results(f2_warped, result_posterior_gamma.y, g_coef[:,valid_index], result_posterior_psi,
                       sigma1[valid_index], accept[valid_index], logl[valid_index], gamma_mat, gamma_stats, Dx, Dy)

    return out
end

function f_SSEg_pw(g::func, q1::func, q2::func)
    par_domain = g.x
    obs_domain = q1.x
    f = f_exp1(g)
    exp1g_temp = f_predictfunction(f, obs_domain)
    pt = zeros(length(obs_domain))
    for i in 2:length(pt)
        tmp = func(obs_domain[1:i],exp1g_temp.y[1:i])
        pt[i] = f_L2norm(tmp)^2
    end
    temp = f_predictfunction(q2,pt)
    vec1 = (q1.y - temp.y .* exp1g_temp.y).^2

    return sum(vec1)
end

function f_logl_pw(g::func, var1, q1::func, q2::func; SSEg=0.0)
    if (SSEg == 0)
        SSEg = f_SSEg_pw(g, q1, q2)
    end
    n = length(q1.y)
    out = n*log(1/sqrt(2*pi)) - n*log(sqrt(var1)) - SSEg/(2*var1)
    return out
end

function f_updateg_pw(g_coef_curr, g_basis, var1_curr, q1::func, q2::func,
                      SSE_curr, propose_g_coef::Function)

    g_coef_prop, ind = propose_g_coef(g_coef_curr)
    g_prop = f_basistofunction(g_basis.x, g_coef_prop, g_basis)
    tmp = f_exp1(g_prop)
    while (minimum(tmp.y)<0)
        g_coef_prop, ind = propose_g_coef(g_coef_curr)
        g_prop = f_basistofunction(g_basis.x, g_coef_prop, g_basis)
        tmp = f_exp1(g_prop)
    end

    g_curr = f_basistofunction(g_basis.x, g_coef_curr, g_basis)
    if (SSE_curr == 0)
        SSE_curr = f_SSEg_pw(g_curr, q1, q2)
    end

    SSE_prop = f_SSEg_pw(g_prop, q1, q2)

    logl_curr = f_logl_pw(g_curr, var1_curr, q1, q2, SSEg=SSE_curr)
    logl_prop = f_logl_pw(g_prop, var1_curr, q1, q2, SSEg=SSE_prop)

    ratio = min(1, exp(logl_prop-logl_curr))
    u = rand()
    if (u <= ratio)
        accept = true
        return g_coef_prop, SSE_prop, logl_prop, accept
    else
        accept = false
        return g_coef_curr, SSE_curr, logl_curr, accept
    end

end

function statsFun(vec)
    a = zeros(size(vec,1))
    b = zeros(size(vec,1))
    for i in 1:size(vec,1)
        a[i] = quantile(vec[i,:],0.025)
        b[i] = quantile(vec[i,:],0.975)
    end
    out = [a b]
    return out
end


text:
using Random

using Dates, TimeZones
using MicroLogging
using NumericalIntegration
using DataFrames, DataFramesMeta, CSV

using StatsBase
using TimeSeries
using HypothesisTests
using Mise
using Impute
using CurveFit
using StatsPlots

function run_analysis()

    @info "Start Analysis"
    flush(stdout); flush(stderr)

    seoul_codes = [
        111121,111123,111131,111141,111142,
        111151,111152,111161,111171,111181,
        111191,111201,111212,111221,111231,
        111241,111251,111261,111262,111273,
        111274,111281,111291,111291,111301,
        111301,111311]
    seoul_names = [
        "중구","종로구","용산구","광진구","성동구",
        "중랑구","동대문구","성북구","도봉구","은평구",
        "서대문구","마포구","강서구","구로구","영등포구",
        "동작구","관악구","강남구","서초구","송파구",
        "강동구","금천구","강북구","강북구","양천구",
        "양천구","노원구"]
    # construct named tuple
    seoul_stations = (; zip(Symbol.(seoul_names), seoul_codes)...)

    # NamedTuple (::Symbol -> ::DataFrame)
    df = load_data_DNN("/input/jongro_seoul.csv", seoul_stations)
    #=
    first(df, 5)
    │ Row │ stationCode │ date                      │ lat      │ long      │ SO2     │ CO      │ O3      │ NO2     │ PM10    │ PM25    │ temp    │ u       │ v       │ pres    │ prep     │ snow     │
    │     │ Int64       │ TimeZones.ZonedDateTime   │ Float64  │ Float64   │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64  │ Float64  │
    ├─────┼─────────────┼───────────────────────────┼──────────┼───────────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┤─────────┤──────────┤──────────┤
    │ 1   │ 111123      │ 2015-01-01T01:00:00+09:00 │ 37.572   │ 127.005   │ 0.004   │ 0.2     │ 0.02    │ 0.009   │ 57.0    │ 0.0     │ -7.4    │ 4.41656 │ 1.60749 │ 1011.8  │ missing  │ missing  │
    │ 2   │ 111123      │ 2015-01-01T02:00:00+09:00 │ 37.572   │ 127.005   │ 0.005   │ 0.2     │ 0.019   │ 0.008   │ 70.0    │ 3.0     │ -8.0    │ 4.22862 │ 1.53909 │ 1011.7  │ missing  │ missing  │
    │ 3   │ 111123      │ 2015-01-01T03:00:00+09:00 │ 37.572   │ 127.005   │ 0.005   │ 0.2     │ 0.02    │ 0.006   │ 92.0    │ 5.0     │ -8.4    │ 3.57083 │ 1.29968 │ 1012.1  │ missing  │ missing  │
    │ 4   │ 111123      │ 2015-01-01T04:00:00+09:00 │ 37.572   │ 127.005   │ 0.004   │ 0.2     │ 0.019   │ 0.005   │ 111.0   │ 2.0     │ -8.8    │ 4.60449 │ 1.6759  │ 1012.3  │ missing  │ missing  │
    │ 5   │ 111123      │ 2015-01-01T05:00:00+09:00 │ 37.572   │ 127.005   │ 0.005   │ 0.2     │ 0.019   │ 0.006   │ 127.0   │ 5.0     │ -9.1    │ 5.35625 │ 1.94951 │ 1011.8  │ missing  │ missing  │
    =#

    #===== start of parameter zone =====#
    total_fdate, total_tdate = get_date_range(df)
    train_fdate = ZonedDateTime(2008, 1, 1, 1, tz"Asia/Seoul")
    train_tdate = ZonedDateTime(2017, 12, 31, 23, tz"Asia/Seoul")
    test_fdate = ZonedDateTime(2018, 1, 1, 0, tz"Asia/Seoul")
    test_tdate = ZonedDateTime(2018, 12, 31, 23, tz"Asia/Seoul")

    # stations
    #=
    train_stn_names = [
        :중구,:종로구,:용산구,:광진구,:성동구,
        :중랑구,:동대문구,:성북구,:도봉구,:은평구,
        :서대문구,:마포구,:강서구,:구로구,:영등포구,
        :동작구,:관악구,:강남구,:서초구,:송파구,
        :강동구,:금천구,:강북구,:강북구,:양천구,
        :양천구,:노원구]
    =#
    #train_stn_names = [:종로구, :강서구, :송파구, :강남구]

    train_stn_names = [
        :중구,:종로구,:용산구,:광진구,:성동구,
        :중랑구,:동대문구,:성북구,:도봉구,:은평구,
        :서대문구,:마포구,:강서구,:구로구,:영등포구,
        :동작구,:관악구,:강남구,:서초구,:송파구,
        :강동구,:금천구,:강북구,:강북구,:양천구,
        :양천구,:노원구]
    train_stn_names = [:종로구]

    df = filter_raw_data(df, train_fdate, train_tdate, test_fdate, test_tdate)
    @show first(df, 10)

    features = [:PM10, :PM25]
    # If you want exclude some features, modify train_features
    # exclude :PM10, :PM25 temporarily for log transform
    train_features = [:PM10, :PM25]
    target_features = [:PM10, :PM25]

    # For GPU, change precision of Floating numbers
    eltype::DataType = Float32

    # simply collect dates, determine exact date for prediction (for 1h, 24h, and so on) later

    Base.Filesystem.mkpath("/mnt/analysis/")

    for ycol in target_features
        # DNN's windowed dataframe is already filterred with dates.
        # size(df) = (sample_size, length(features))
        @info "Analysis of $(string(ycol))"
        for name in train_stn_names
            code = seoul_stations[name]
            stn_df = filter_raw_data(df, code, train_fdate, train_tdate)
            # for imputatation
            CSV.write("/mnt/analysis/raw_data_analysis_$(string(ycol))_$((name)).csv", stn_df)
            dates = stn_df[!, :date]

            # Imputation
            # allow missing again
            DataFrames.allowmissing!(stn_df)

            zero2Missing!(stn_df, ycol)

            # Analysis : pdf
            @info "Estimate probability density function without imputation..."
            ycol_arr = collect(skipmissing(stn_df[!, ycol]))
            if ycol == :PM10
                npts = 256
            elseif ycol == :PM25
                npts = 256
            end
            ycol_pdfx, ycol_pdfy = pdf(ycol_arr, npts, method=:kernel)
            plot_anal_pdf(ycol_pdfx, ycol_pdfy, ycol,
                "PDF for $(ycol)", "/mnt/analysis/", "$(string(ycol))_$(string(name))_pdf_noimpute")
            ycol_logpdfx, ycol_logpdfy, = pdf(log10.(ycol_arr), npts, method=:kernel)
            plot_anal_pdf(ycol_logpdfx, ycol_logpdfy, ycol,
                "Log PDF for $(ycol)", "/mnt/analysis/", "$(string(ycol))_$(string(name))_logpdf_noimpute")

            # mean imputation
            #impute!(df, ycol, :mean; total_mean = total_mean)
            # simple sampling
            stn_df[!, ycol] = Impute.srs(stn_df[!, ycol])
            # knn imputation
            #nK = 10
            #ycols = [:SO2, :CO, :O3, :NO2, :PM10, :PM25, :temp, :u, :v, :pres, :humid, :prep, :snow]
            #impute!(stn_df, ycol, :knn; ycols = ycols, leafsize = nK)

            # disallow missing again
            DataFrames.disallowmissing!(stn_df)

            ta = TimeArray(dates, stn_df[!, ycol])

            @info "Estimate probability density function without imputation..."
            ycol_arr = stn_df[!, ycol]
            ycols = stn_df[!, ycol]
            logycols = log10.(stn_df[!, ycol])

            total_mean, total_std = StatsBase.mean_and_std(values(ta))

            @info "Total mean of $(string(ycol)) : ", total_mean
            @info "Total std  of $(string(ycol)) : ", total_std

            if ycol == :PM10
                npts = 200
            elseif ycol == :PM25
                npts = 200
            end
            ycol_pdfx, ycol_pdfy = pdf(ycol_arr, npts, method=:kernel)
            plot_anal_pdf(ycol_pdfx, ycol_pdfy, ycol,
                "PDF for $(ycol)", "/mnt/analysis/", "$(string(ycol))_$(string(name))_pdf_impute")
            ycol_logpdfx, ycol_logpdfy, = pdf(log10.(ycol_arr), npts, method=:kernel)
            plot_anal_pdf(ycol_logpdfx, ycol_logpdfy, ycol,
                "Log PDF for $(ycol)", "/mnt/analysis/", "$(string(ycol))_$(string(name))_logpdf_impute")
            DataFrames.insertcols!(stn_df, 3, Symbol(ycol, "_impute") => ycol_arr)
            DataFrames.insertcols!(stn_df, 3, Symbol(ycol, "_logimpute") => log10.(ycol_arr))

            # Analysis : HypothesisTests
            @info "Hypothesis Tests"
            max_lag = 30 * 24
            #@show HypothesisTests.ADFTest(values(ta), :none, max_lag)
            #@show HypothesisTests.ADFTest(values(ta), :constant, max_lag)
            #@show HypothesisTests.ADFTest(values(ta), :trend, max_lag)
            #@show HypothesisTests.ADFTest(values(ta), :squared_trend, max_lag)
            #@show HypothesisTests.LjungBoxTest(df[!, ycol], 365*24)

            # Analysis : Autocorrelation Functions]
            #=
            for _ycol in [ycol, Symbol(ycol, :_impute), Symbol(ycol, :_logimpute)]
                Base.Filesystem.mkpath("/mnt/analysis/$(string(_ycol))/")
                #max_time = Int(ceil(min(size(stn_df[!, ycol],1)-1, 10*log10(size(stn_df[!, ycol],1)))))
                max_time = 15 * 24
                @info "Autocorrelation..."
                lags = 0:max_time
                acf = StatsBase.autocor(stn_df[!, _ycol], lags)
                plot_anal_correlogram(acf, _ycol,
                    "Autocorrelation",
                    "/mnt/analysis/$(string(_ycol))/", "$(string(ycol))_$(string(name))_mt$(max_time)_acf")
                acf_df = DataFrame(time = collect(0:max_time), corr = acf)
                CSV.write("/mnt/analysis/$(string(_ycol))/$(string(ycol))_$(string(name))_mt$(max_time)_acf.csv", acf_df)
                @info "Integral Time Scale for raw data..", compute_inttscale(0:max_time, acf)

                pacf = StatsBase.pacf(stn_df[!, _ycol], lags)
                plot_anal_correlogram(pacf, _ycol,
                    "Partial Autocorrelation",
                    "/mnt/analysis/$(string(_ycol))/", "$(string(ycol))_$(string(name))_mt$(max_time)_pacf")

                # Analysis : Autocorrelation Functions
                # No need to compute pacf
                max_time = 365 * 24
                lags = 0:max_time
                acf = StatsBase.autocor(stn_df[!, _ycol], lags)
                plot_anal_correlogram(acf, ycol,
                    "Autocorrelation",
                    "/mnt/analysis/$(string(_ycol))/", "$(string(ycol))_$(string(name))_mt$(max_time)_acf")
                acf_df = DataFrame(time = collect(0:max_time), corr = acf)
                CSV.write("/mnt/analysis/$(string(_ycol))/$(string(ycol))_$(string(name))_mt$(max_time)_acf.csv", acf_df)
                @info "Integral Time Scale for raw data..", compute_inttscale(0:max_time, acf)
            end
            =#
            #=
            # Analysis : Integral Length Scale (Total)
            int_scale = integrate(1:length(acf), acf, SimpsonEven())
            @info "Integral Time Scale of $(ycol) at $(name) (total): ", int_scale

            plot_anal_autocor(acf, ycol,
                "Autocorrelation (total series)",
                "/mnt/analysis/", "$(string(ycol))_$(string(name))_total")
            =#
            # Analysis : Time difference
            #=
            @info "Time difference..."
            order = 1
            lag = 1
            diffed = TimeSeries.diff(ta, lag, padding=true, differences=order)
            val_diffed = copy(values(diffed))
            replace!(val_diffed, Inf => 0, -Inf => 0, NaN => 0)
            @info "mean of time $(order)th difference of $(ycol) at $(name): ", StatsBase.mean(val_diffed)
            @info "std  of time $(order)th difference of $(ycol) at $(name): ", StatsBase.std(val_diffed)
            plot_anal_lineplot(DateTime.(timestamp(diffed)), val_diffed, ycol,
                order, lag, "/mnt/analysis/", "$(string(ycol))_$(string(name))")

            order = 2
            lag = 1
            diffed = TimeSeries.diff(ta, lag, padding=true, differences=order)
            val_diffed = copy(values(diffed))
            replace!(val_diffed, Inf => 0, -Inf => 0, NaN => 0)
            @info "mean of time $(order)th difference of $(ycol) at $(name): ", StatsBase.mean(val_diffed)
            @info "std  of time $(order)th difference of $(ycol) at $(name): ", StatsBase.std(val_diffed)
            plot_anal_lineplot(DateTime.(timestamp(diffed)), val_diffed, ycol,
                order, lag, "/mnt/analysis/", "$(string(ycol))_$(string(name))")

            order = 1
            lag = 24
            diffed = TimeSeries.diff(ta, lag, padding=true, differences=order)
            val_diffed = copy(values(diffed))
            replace!(val_diffed, Inf => 0, -Inf => 0, NaN => 0)
            @info "mean of time $(order)th difference of $(ycol) at $(name): ", StatsBase.mean(val_diffed)
            @info "std  of time $(order)th difference of $(ycol) at $(name): ", StatsBase.std(val_diffed)
            plot_anal_lineplot(DateTime.(timestamp(diffed)), val_diffed, ycol,
                order, lag, "/mnt/analysis/", "$(string(ycol))_$(string(name))")
            =#
            # Analysis : Grouping by time
            @info "Groupby time..."
            #=
            for _ycol in [ycol, Symbol(ycol, :_impute), Symbol(ycol, :_logimpute)]
                Base.Filesystem.mkpath("/mnt/analysis/$(string(_ycol))/")
                period_df = copy(stn_df)

                # find periodicty
                pfuncs = (
                    hour = :hour,
                    week = :week,
                    month = :month)

                # Analysis : seasonality analysis, hourly, weekly, monthly
                for period_dir in [:hour, :week, :month]
                    @info "Periodic mean of $(string(_ycol)) ($(string(period_dir))ly...) "
                    pfunc = pfuncs[period_dir]
                    periodic_col = eval(pfunc).(period_df[!, :date])
                    DataFrames.insertcols!(period_df, 3, period_dir => periodic_col)

                    periodic_means, periodic_means_repeated = populate_periodic_mean(period_df, _ycol, period_dir)
                    # insert mean and fluctuation to DataFrame
                    DataFrames.insertcols!(period_df, 3, Symbol(period_dir, "_mean") => periodic_means_repeated)
                    DataFrames.insertcols!(period_df, 3, Symbol(period_dir, "_fluc") => period_df[!, ycol] .- periodic_means_repeated)

                    # violin plot mean
                    plot_anal_violin(period_df, _ycol, period_dir, periodic_means,
                        "Violin plot for $(ycol)", "/mnt/analysis/$(string(_ycol))/", "$(string(ycol))_$(string(name))")
                    # plot mean only
                    plot_anal_periodic_mean(period_df, _ycol, period_dir, periodic_means,
                        "$(uppercasefirst(string(period_dir)))ly mean for $(ycol)", "/mnt/analysis/$(string(_ycol))/", "$(string(ycol))_$(string(name))_time_mean")
                end

                # Analysis : fluctuation analysis, hourly, weekly, monthly
                # Use Split-Apply-Combine Strategy
                for period_dir in [:hour, :week, :month]
                    @info "Autocorrelation of periodic fluctuations of $(string(_ycol)) $(string(period_dir))ly..."

                    period_fluc_df = DataFrames.by(period_df, _ycol) do _df
                        (; period_dir => getproperty(_df, period_dir),
                        :date => _df.date,
                        :mean => getproperty(_df, Symbol(period_dir, "_mean")),
                        :fluc => getproperty(_df, _ycol) - getproperty(_df, Symbol(period_dir, "_mean")))
                    end

                    period_fluc_means = periodic_mean(period_fluc_df, _ycol, period_dir)
                    plot_anal_periodic_mean(period_fluc_df, _ycol, period_dir, period_fluc_means,
                        "$(uppercasefirst(string(period_dir)))ly fluctuation mean for $(ycol)", "/mnt/analysis/$(string(_ycol))/",
                        "$(string(ycol))_$(string(name))_$(string(_ycol))_period_flucmean")

                    #max_time = Int(ceil(min(size(period_fluc_df[!, _ycol],1)-1, 10*log10(size(period_fluc_df[!, _ycol],1)))))
                    max_time = 15 * 24
                    lags = 0:max_time
                    @info "$(string(period_dir))ly Fluctuation Autocorrelation..."
                    acf = StatsBase.autocor(period_fluc_df[!, :fluc], lags)
                    plot_anal_correlogram(acf, _ycol,
                        "$(uppercasefirst(string(period_dir)))ly Fluctuation Autocorrelation",
                        "/mnt/analysis/$(string(_ycol))/",
                        "$(string(ycol))_$(string(name))_mt$(max_time)_fluc_$(string(period_dir))ly_acf")

                    max_time = 365 * 24
                    lags = 0:max_time
                    acf = StatsBase.autocor(period_fluc_df[!, :fluc], lags)
                    plot_anal_correlogram(acf, _ycol,
                        "$(uppercasefirst(string(period_dir)))ly Fluctuation Autocorrelation",
                        "/mnt/analysis/$(string(_ycol))/",
                        "$(string(ycol))_$(string(name))_mt$(max_time)_fluc_$(string(period_dir))ly_acf")
                end
            end
            =#
            # decompose seasonal components
            Base.Filesystem.mkpath("/mnt/analysis/decomposition/$(string(ycol))/")
            #lee_year_sea1, lee_year_sea2, lee_year_res2, lee_day_sea1, lee_day_res1 =
            #    season_adj_lee(stn_df, Symbol(ycol, "_impute"))

            decomposed_df = decompose_seasonality(stn_df, Symbol(ycol, "_impute"))
            CSV.write("/mnt/analysis/decomposition/$(string(ycol))/" * "$(string(ycol))_$(string(name))_new_sea_decompose.csv", decomposed_df)
            plot_year = 2016
            plot_seasonality(ycols, decomposed_df, ycol, plot_year, DateTime(train_fdate, Local),
                "Seasonal Decomposition in $(plot_year)",
                "/mnt/analysis/decomposition/$(string(ycol))/", "$(string(ycol))_$(string(name))_new_sea_decompose")
            #=
            plot_seasonality(ycols, lee_year_sea1, lee_year_sea2, lee_day_sea1, lee_day_res1, ycol, plot_year,
                DateTime(stn_df[1, :date], Local),
                ["Original data", "Annual seasonality", "Daily seasonality", "Residual"],
                "Seasonal Decomposition in $(plot_year)",
                "/mnt/analysis/decomposition/$(string(ycol))/", "$(string(ycol))_$(string(name))_seasonality")
            plot_seasonality(ycols, lee_year_sea1, lee_year_sea2, lee_day_sea1, lee_day_res1, ycol, 2015, 2017,
                DateTime(stn_df[1, :date], Local),
                ["Original data", "Annual seasonality", "Daily seasonality", "Residual"],
                "Seasonal Decomposition from 2015 to 2017",
                "/mnt/analysis/decomposition/$(string(ycol))/", "$(string(ycol))_$(string(name))_seasonality_3year")
            =#
            # Autocorrelation
            max_time = 15*24
            acf_raw = StatsBase.autocor(decomposed_df[!, :raw], 0:max_time)
            acf_yres2_s = StatsBase.autocor(decomposed_df[!, :year_res_s], 0:max_time)
            acf_yres2 = StatsBase.autocor(decomposed_df[!, :year_res], 0:max_time)
            acf_dres1 = StatsBase.autocor(decomposed_df[!, :day_res], 0:max_time)
            #acf_dres1_log = StatsBase.autocor(log.(lee_day_res1), 0:15*24)

            plot_anal_correlogram(acf_raw, ycol,
                "Raw Autocorrelation",
                "/mnt/analysis/decomposition/$(string(ycol))/",
                "$(string(ycol))_$(string(name))_mt$(max_time)_acf_raw")
            plot_anal_correlogram(acf_yres2_s, ycol,
                "Annual Residual (Smoothed) Autocorrelation",
                "/mnt/analysis/decomposition/$(string(ycol))/",
                "$(string(ycol))_$(string(name))_mt$(max_time)_acf_yres2_s")
            plot_anal_correlogram(acf_yres2, ycol,
                "Annual Residual (Raw) Autocorrelation",
                "/mnt/analysis/decomposition/$(string(ycol))/",
                "$(string(ycol))_$(string(name))_mt$(max_time)_acf_yres2")
            plot_anal_correlogram(acf_dres1, ycol,
                "Daily Residual Autocorrelation",
                "/mnt/analysis/decomposition/$(string(ycol))/",
                "$(string(ycol))_$(string(name))_mt$(max_time)_acf_dres1")
            acf_df = DataFrame(time = collect(0:max_time), corr = acf_yres2)

            Base.mkpath("/mnt/analysis/$(string(ycol))")
            CSV.write("/mnt/analysis/$(string(ycol))/$(string(ycol))_$(string(name))_mt$(max_time)_acf_yres2.csv", acf_df)

            intT = compute_inttscale(0:max_time, acf_yres2)
            @info "Integral Time Scale for Annual Residual..", intT
            # not to log negative value when fitting
            positiveIntT = max(1, Int(round(intT)) - 1)
            fit = exp_fit(acf_df[1:positiveIntT, :time], acf_df[1:positiveIntT, :corr])
            @info "Coef a and b in a * exp(b * x) for Annual Residual.. to 0:$(positiveIntT)", fit[1], fit[2]

            fit = exp_fit(acf_df[1:positiveIntT, :time], acf_df[1:positiveIntT, :corr])
            fit_x = acf_df[1:positiveIntT, :time]
            fit_y = fit[1] .* exp.(fit_x .* fit[2])
            intT = compute_inttscale(fit_x, fit_y)
            @info "Integral Time Scale for Fitted Annual Residual..", intT
            @info "Coef a and b in a * exp(b * x / T) for Fitted Annual Residual.. to 0:$(positiveIntT)", fit[1], fit[2] / intT

            acf_df = DataFrame(time = collect(0:max_time), corr = acf_dres1)
            CSV.write("/mnt/analysis/$(string(ycol))/$(string(ycol))_$(string(name))_mt$(max_time)_acf_dres1.csv", acf_df)
            intT = compute_inttscale(0:max_time, acf_dres1)
            @info "Integral Time Scale for Daily Residual..", intT
            positiveIntT = max(1, Int(round(intT)) - 1)
            fit = exp_fit(acf_df[1:positiveIntT, :time], acf_df[1:positiveIntT, :corr])
            @info "Coef a and b in a * exp(b * T) for Daily Residual.. to 0:$(positiveIntT)", fit[1], fit[2]

            fit = exp_fit(acf_df[1:positiveIntT, :time], acf_df[1:positiveIntT, :corr])
            fit_x = acf_df[1:positiveIntT, :time]
            fit_y = fit[1] .* exp.(fit_x .* fit[2])
            intT = compute_inttscale(fit_x, fit_y)
            @info "Integral Time Scale for Fitted Daily Residual..", intT
            @info "Coef a and b in a * exp(b * x / T) for Fitted Daily Residual.. to 0:$(positiveIntT)", fit[1], fit[2] / intT

            max_time = 365*24
            acf_yres2_long = StatsBase.autocor(decomposed_df[!, :year_res], 0:max_time)
            acf_dres1_long = StatsBase.autocor(decomposed_df[!, :day_res], 0:max_time)
            plot_anal_correlogram(acf_yres2_long, ycol,
                "Annual Residual Autocorrelation",
                "/mnt/analysis/decomposition/$(string(ycol))/",
                "$(string(ycol))_$(string(name))_mt$(max_time)_acf_yres2_long")
            plot_anal_correlogram(acf_dres1_long, ycol,
                "Daily Residual Autocorrelation",
                "/mnt/analysis/decomposition/$(string(ycol))/",
                "$(string(ycol))_$(string(name))_mt$(max_time)_acf_dres1_long")
            acf_df = DataFrame(time = collect(0:max_time), corr = acf_yres2_long)
            CSV.write("/mnt/analysis/$(string(ycol))/$(string(ycol))_$(string(name))_mt$(max_time)_acf_yres2.csv", acf_df)
            intT = compute_inttscale(0:max_time, acf_yres2_long)
            @info "Integral Time Scale for Annual Residual..", intT
            acf_df = DataFrame(time = collect(0:max_time), corr = acf_dres1_long)
            CSV.write("/mnt/analysis/$(string(ycol))/$(string(ycol))_$(string(name))_mt$(max_time)_acf_dres1.csv", acf_df)
            intT = compute_inttscale(0:max_time, acf_dres1_long)
            @info "Integral Time Scale for Daily Residual..", compute_inttscale(0:max_time, acf_dres1_long)
        end
    end

end

run_analysis()


text:
"""
  Compare the fields of two meshes where the first has Tmsh = Float64 and
  the second one has Tmsh = Complex128
"""
function compare_meshes(mesh, mesh_c)

  # decide which fields to check
  fieldnames = [:vert_coords, :coords_bndry, :coords_interface, :dxidx,
                :jac, :nrm_face, :nrm_bndry]

  if mesh.coord_order == 1
    push!(fieldnames, :dxidx_face, :dxidx_bndry, :jac_face, :jac_bndry)
  end

  fieldnames_shared = [:coords_sharedface, :dxidx_sharedface, :jac_sharedface,
                       :nrm_sharedface]

  if mesh.coord_order == 1
    push!(fieldnames_shared, :dxidx_sharedface, :jac_sharedface)
  end

  # check that the real parts of the field are the same and the imaginary
  # parts are zero
  @testset "----- Comparing meshes -----" begin
    for i in fieldnames

      f_real = getfield(mesh, i)
      f_imag = getfield(mesh_c, i)

      if length(f_real) != 0
        @test isapprox( vecnorm(f_real - real(f_imag))/sqrt(length(f_real)), 0.0) atol=1e-14
        @test isapprox( vecnorm(imag(f_imag)), 0.0) atol=1e-14
      end
    end

    for i in fieldnames_shared

      f_real = getfield(mesh, i)
      f_imag = getfield(mesh_c, i)

      for peer=1:mesh.npeers
        f_real_peer = f_real[peer]
        f_imag_peer = f_imag[peer]

        @test isapprox( vecnorm(f_real_peer - real(f_imag_peer))/sqrt(length(f_real_peer)), 0.0) atol=1e-14
        @test isapprox( vecnorm(imag(f_imag_peer)), 0.0) atol=1e-14
      end
    end
  end

  return nothing
end


text:
using DataFrames
using DataFramesMeta
using CSV

d = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "icaps_2018","data")

solver_order = ["pomcpow" => "POMCPOW",
                "pft" => "PFT-DPW",
                "qmdp" => "QMDP",
                "pomcpdpw" => "POMCP-DPW",
                "despot" => "DESPOT",
                "d_pomcp" => "POMCP\\textsuperscript{D}",
                "d_despot" => "DESPOT\\textsuperscript{D}"]

cardinality = Dict("lasertag" => "(D, D, D)",
                   "lightdark" => "(D, D, C)",
                   "subhunt" => "(D, D, C)",
                   "vdpbarrier" => "(C, C, C)",
                   "vdptag" => "(C, C, C)")

problem_order = ["lasertag", "lightdark", "subhunt", "vdpbarrier"]

filenames = Dict("simpleld"=>"/home/zach/.julia/v0.6/ContinuousPOMDPTreeSearchExperiments/icaps_2018/data/simpleld_Monday_5_Mar_21_39.csv","subhunt"=>"/home/zach/.julia/v0.6/ContinuousPOMDPTreeSearchExperiments/icaps_2018/data/subhunt_Monday_5_Mar_22_08.csv","lasertag"=>"/home/zach/.julia/v0.6/ContinuousPOMDPTreeSearchExperiments/icaps_2018/data/lasertag_Monday_5_Mar_20_17.csv","vdpbarrier"=>"/home/zach/.julia/v0.6/ContinuousPOMDPTreeSearchExperiments/icaps_2018/data/vdpbarrier_Monday_5_Mar_23_01.csv")
filenames["lightdark"] = filenames["simpleld"]

# filenames = Dict("lasertag" => "$(Pkg.dir("ContinuousPOMDPTreeSearchExperiments"))/icaps_2018/data/lasertag_Monday_26_Feb_18_46.csv",
#                  "lightdark" => "$(Pkg.dir("ContinuousPOMDPTreeSearchExperiments"))/icaps_2018/data/simpleld_Monday_26_Feb_20_13.csv",
#                  "subhunt" => "$(Pkg.dir("ContinuousPOMDPTreeSearchExperiments"))/icaps_2018/data/subhunt_Monday_26_Feb_20_44.csv",
#                  "vdpbarrier" => "$(Pkg.dir("ContinuousPOMDPTreeSearchExperiments"))/icaps_2018/data/bdpbarrier_Monday_26_Feb_21_42.csv")

data = Dict(
    "lightdark" => Dict(
        "limits" => (-20.0, 80.0),
        "name" => "Light Dark"
    ),

    "subhunt" => Dict(
        "limits" => (0.0, 80.0),
        "name" => "Sub Hunt"
    ),

    "lasertag" => Dict(
        "limits" => (-20.0, -8.0),
        "name" => "Laser Tag"
    ),

    "vdptag" => Dict(
        "limits" => (-20.0, 40.0),
        "name" => "VDP Tag"
    ),

    "vdpbarrier" => Dict(
        "limits" => (0.0, 31.0),
        "name" => "VDP Tag"
    ),

)

for p in problem_order
    df = CSV.read(filenames[p])
    d = data[p]
    for s in unique(df[:solver])
        rs = @where(df, :solver.==s)[:reward]
        m = mean(rs)
        sem = std(rs)/sqrt(length(rs))
        if s == "ar_despot"
            s = "despot"
        elseif s == "despot_01"
            s = "despot"
        elseif p == "lasertag" && s == "pomcp"
            s = "d_pomcp"
        end
        d[s] = (m, sem)
    end
end

hbuf = IOBuffer()

for k in problem_order
    print(hbuf, "& $(data[k]["name"]) \\makebox[0pt][l]{$(cardinality[k])} & ")
end
print(hbuf, "\\\\")

tbuf = IOBuffer()
for (k, n) in solver_order
    print(tbuf, n*" ")
    for p in problem_order
        d = data[p]
        if haskey(d, k)
            m, sem = d[k]
            lo, hi = d["limits"]
            frac = (m-lo)/(hi-lo)
            @printf(tbuf, "& \\result{%.1f}{%.1f}{%d}{%.2f} ",
                    m, sem, round(Int, 100*frac), frac
                   )
        else
            print(tbuf,  "& \\noresult{} ")
        end
    end
    print(tbuf, "\\\\\n")
end

columns = "l"*"rl"^length(data)

tabletex = """
    \\begin{tabular}{$columns}
        \\toprule
            $(String(hbuf))
        \\midrule
            $(String(tbuf))
        \\bottomrule
    \\end{tabular}
"""

println(tabletex)


text:
using VersionVigilante

using Pkg
using Test

const repo_url = "https://github.com/bcbi-test/versionvigilante-integration-test-repo.git"

@testset "VersionVigilante.jl" begin
    @testset "assert.jl" begin
        @test nothing == VersionVigilante.always_assert(true, "")
        @test_throws VersionVigilante.AlwaysAssertionError VersionVigilante.always_assert(false, "")
    end
    @testset "compare_versions.jl" begin
        @testset "isprerelease" begin
            @test !VersionVigilante.isprerelease(v"1")
            @test !VersionVigilante.isprerelease(v"1-")
            @test VersionVigilante.isprerelease(v"1-dev")
            @test VersionVigilante.isprerelease(v"1-PRE")
        end
        @testset "_calculate_increment" begin
            include("test_calculate_increment.jl")
        end

    end
    @testset "main.jl" begin
        for (pass, master_branch, feature_branch, version) in [
                (false, "master-1.0.0", "feature-0.9.9", "0.9.9"), # 1.0.0 -> 0.9.9 is FAIL
                (false, "master-1.0.0", "feature-1.0.0-", "1.0.0-"), # 1.0.0 -> 1.0.0- is FAIL
                (false, "master-1.0.0", "feature-1.0.0-", "1.0.0-DEV"), # 1.0.0 -> 1.0.0-DEV is FAIL
                (false, "master-1.0.0", "feature-1.0.0", "1.0.0"), # 1.0.0 -> 1.0.0 is FAIL
                (true, "master-1.0.0", "feature-1.0.1", "1.0.1"), # 1.0.0 -> 1.0.1 is PASS
                (true, "master-1.0.0", "feature-1.1.0", "1.1.0"), # 1.0.0 -> 1.1.0 is PASS
                (true, "master-1.0.0", "feature-2.0.0", "2.0.0"), # 1.0.0 -> 2.0.0 is PASS
                (true, "master-1.0.0-dev", "feature-1.0.0-dev", "1.0.0-DEV"), # 1.0.0-DEV -> 1.0.0-DEV is PASS if `allow_unchanged_prerelease = true`
                (true, "master-1.0.0-dev", "feature-1.0.0", "1.0.0"), # 1.0.0-DEV -> 1.0.0 is PASS
                ]
            VersionVigilante.with_branch(repo_url, master_branch, feature_branch) do
                rm("Project.toml"; force = true, recursive = true)
                open("Project.toml", "w") do io
                    println(io, "version = \"$(version)\"")
                end
                if pass
                    withenv() do
                        @test nothing == VersionVigilante.main(repo_url;
                                                               master_branch = master_branch)
                    end
                else
                    withenv() do
                        @test_throws ErrorException VersionVigilante.main(repo_url;
                                                                          master_branch = master_branch)
                    end
                end
            end
        end

        for (pass, master_branch, feature_branch, version, allow_unchanged_prerelease) in [
                (false, "master-1.0.0", "feature-1.0.0", "1.0.0", true),
                (false, "master-1.0.0", "feature-1.0.0", "1.0.0", false),
                (true, "master-1.0.0-dev", "feature-1.0.0-dev", "1.0.0-DEV", true), # 1.0.0-DEV -> 1.0.0-DEV is FAIL if `allow_unchanged_prerelease = false`
                (false, "master-1.0.0-dev", "feature-1.0.0-dev", "1.0.0-DEV", false), # 1.0.0-DEV -> 1.0.0-DEV is FAIL if `allow_unchanged_prerelease = false`
                ]
            VersionVigilante.with_branch(repo_url, master_branch, feature_branch) do
                rm("Project.toml"; force = true, recursive = true)
                open("Project.toml", "w") do io
                    println(io, "version = \"$(version)\"")
                end
                if pass
                    withenv() do
                        @test nothing == VersionVigilante.main(repo_url;
                                                               master_branch = master_branch,
                                                               allow_unchanged_prerelease = allow_unchanged_prerelease)
                    end
                else
                    withenv() do
                        @test_throws ErrorException VersionVigilante.main(repo_url;
                                                                          master_branch = master_branch,
                                                                          allow_unchanged_prerelease = allow_unchanged_prerelease)
                    end
                end
            end
        end

        # (master_branch, feature_branch, version)
        _skipping_versions_integration_tests = [("master-1.0.0", "feature-$(x)", "$(x)") for x in ["1.0.2",
                                                                                "1.0.3",
                                                                                "1.2.0",
                                                                                "1.2.1",
                                                                                "1.2.2",
                                                                                "1.2.3",
                                                                                "3.0.0",
                                                                                "3.0.1",
                                                                                "3.1.0",
                                                                                "3.1.1"]]
        for (pass, master_branch, feature_branch, version, allow_skipping_versions) in vcat([(true, y..., true) for y in _skipping_versions_integration_tests],
                                                                                [(false, y..., false) for y in _skipping_versions_integration_tests])
            VersionVigilante.with_branch(repo_url, master_branch, feature_branch) do
                rm("Project.toml"; force = true, recursive = true)
                open("Project.toml", "w") do io
                    println(io, "version = \"$(version)\"")
                end
                if pass
                    withenv() do
                        @test nothing == VersionVigilante.main(repo_url;
                                                               master_branch = master_branch,
                                                               allow_skipping_versions = allow_skipping_versions)
                    end
                else
                    withenv() do
                        @test_throws ErrorException VersionVigilante.main(repo_url;
                                                                          master_branch = master_branch,
                                                                          allow_skipping_versions = allow_skipping_versions)
                    end
                end
            end
        end

        # printing output to stdout for GitHub Actions
        VersionVigilante.with_branch(repo_url,
                                     "master-1.0.0-dev",
                                     "feature-1.0.0") do
            rm("Project.toml"; force = true, recursive = true)
            open("Project.toml", "w") do io
                println(io, "version = \"1.0.0\"")
            end
            withenv("GITHUB_ACTIONS" => "true") do
                @test nothing == VersionVigilante.main(repo_url;
                                                       master_branch = "master-1.0.0-dev")
            end
            withenv("GITHUB_ACTIONS" => "false") do
                @test nothing == VersionVigilante.main(repo_url;
                                                       master_branch = "master-1.0.0-dev")
            end
        end
    end
    @testset "parse_project.jl" begin
        VersionVigilante.with_temp_dir() do
            rm("Project.toml"; force = true, recursive = true)
            open("Project.toml", "w") do io
                println(io, "version = \"1.2.3\"")
            end
            @test VersionVigilante.get_version() == v"1.2.3"
        end
    end
    @testset "utils.jl" begin
    end
end


text:
module jet_engine
using ReachabilityAnalysis
@taylorize function jet_engine!(dx, x, params, t)
    dx[1] = -x[2] - 1.5 * x[1]^2 - 0.5 * x[1]^3 - 0.5
    dx[2] = 3 * x[1] - x[2]
    return dx
end

function model(X0)
    S = @system(x' = jet_engine!(x), dim:2)
    return IVP(S, X0)
end
end # module

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary