text:
"""
    result = view_index(db::Database, ddoc::AbstractString, name::AbstractString, map::AbstractString;
      reduce::AbstractString = "")

Create a secondary index.

The `map` is a string containing a map function in Javascript. Currently, `make_view`
can only create a single view per design document.

The optional `reduce` parameter is a string containing either a custom Javascript
reducer (best avoided for performance reasons) or the name of a built-in Erlang
reducer, e.g. `"_stats"`.

### Examples

    result = view_index(db, "my_ddoc", "my_view", "function(doc){if(doc&&doc.name){emit(doc.name,1);}}")

### Returns

Returns a `Dict(...)` from the CouchDB response, of the type

    Dict(
      "ok"  => true,
      "rev" => "1-b950984b19bb1b8bb43513c9d5b235bc",
      "id"  => "_design/my_ddoc"
    )

[API reference](https://docs.cloudant.com/creating_views.html)
"""
function view_index(db::Database, ddoc::AbstractString, name::AbstractString, map::AbstractString; reduce::AbstractString = "")
  data = Dict(
    "views"    => Dict(name => Dict("map" => map)),
    "language" => "javascript"
  )

  if reduce != ""
    data["views"][name]["reduce"] = reduce
  end

  relax(put, endpoint(db.url, "_design/$ddoc"); json=data, cookies=db.client.cookies)
end

"""
    result = view_query(db::Database, ddoc::AbstractString, name::AbstractString;
      descending    = false,
      endkey        = "",
      include_docs  = false,
      conflicts     = false,
      inclusive_end = true,
      group         = false,
      group_level   = 0,
      reduce        = true,
      key           = "",
      keys          = [],
      limit         = 0,
      skip          = 0,
      startkey      = "")

Query a secondary index.

### Examples

    # Query the view for a known key subset
    result = view_query(db, "my_ddoc", "my_view"; keys=["adam", "billy"])

### Returns

    Dict(
      "rows" => [
        Dict("key" => "adam", "id" => "591c02fa8b8ff14dd4c0553670cc059a", "value" => 1),
        Dict("key" => "billy", "id" => "591c02fa8b8ff14dd4c0553670cc13c1", "value" => 1)
      ],
      "offset" => 0,
      "total_rows" => 7
    )

[API reference](https://docs.cloudant.com/using_views.html)
"""
function view_query(db::Database, ddoc::AbstractString, name::AbstractString;
  descending    = false,
  endkey        = "",
  include_docs  = false,
  conflicts     = false,
  inclusive_end = true,
  group         = false,
  group_level   = 0,
  reduce        = true,
  key           = "",
  keys          = [],
  limit         = 0,
  skip          = 0,
  startkey      = "")

  query::Dict{String, Any} = Dict()

  if descending
    query["descending"] = true
  end

  if include_docs
    query["include_docs"] = true
  end

  if conflicts
    query["conflicts"] = true
  end

  if group
    query["group"] = true
  end

  if !reduce
    query["reduce"] = false
  end

  if group_level > 0
    query["group_level"] = group_level
  end

  if !inclusive_end
    query["inclusive_end"] = false
  end

  if endkey != ""
    query["endkey"] = JSON.json(endkey)
  end

  if key != ""
    query["key"] = JSON.json(key)
  end

  if limit > 0
    query["limit"] = limit
  end

  if skip > 0
    query["skip"] = skip
  end

  if startkey != ""
    query["startkey"] = JSON.json(startkey)
  end

  url = endpoint(db.url, "_design/$ddoc/_view/$name")
  if length(keys) > 0
    relax(post, url; json=Dict("keys" => keys), cookies=db.client.cookies, query=query)
  else
    relax(get, url; cookies=db.client.cookies, query=query)
  end
end

"""
    alldocs(db::Database;
      descending    = false,
      endkey        = "",
      include_docs  = false,
      conflicts     = false,
      inclusive_end = true,
      key           = "",
      keys          = [],
      limit         = 0,
      skip          = 0,
      startkey      = "")

Return all documents in the database by the primary index.

The optional parameters are:

* descending     true/false   -- lexicographical ordering of keys. Default false.
* endkey         id           -- stop when `endkey` is reached. Optional.
* startkey       id           -- start at `startkey`. Optional.
* include_docs   true/false   -- return the document body. Default false.
* conflicts      true/false   -- also return any conflicting revisions. Default false.
* inclusive_end  true/false   -- if `endkey` is given, should this be included? Default true
* key            id           -- return only specific key. Optional.
* keys           [id, id,...] -- return only specific set of keys (will POST). Optional.
* limit          int          -- return only max `limit` number of rows. Optional.
* skip           int          -- skip over the first `skip` number of rows. Default 0.

[API reference](http://docs.couchdb.org/en/1.6.1/api/database/bulk-api.html#db-all-docs)
"""
function alldocs(db::Database;
  descending    = false,
  endkey        = "",
  include_docs  = false,
  conflicts     = false,
  inclusive_end = true,
  key           = "",
  keys          = [],
  limit         = 0,
  skip          = 0,
  startkey      = "")

  query::Dict{String, Any} = Dict()

  if descending
    query["descending"] = true
  end

  if include_docs
    query["include_docs"] = true
  end

  if conflicts
    query["conflicts"] = true
  end

  if !inclusive_end
    query["inclusive_end"] = false
  end

  if endkey != ""
    query["endkey"] = JSON.json(endkey)
  end

  if key != ""
    query["key"] = JSON.json(key)
  end

  if limit > 0
    query["limit"] = limit
  end

  if skip > 0
    query["skip"] = skip
  end

  if startkey != ""
    query["startkey"] = JSON.json(startkey)
  end

  if length(keys) > 0
    relax(post, endpoint(db.url, "_all_docs"); json=Dict("keys" => keys), cookies=db.client.cookies, query=query)
  else
    relax(get, endpoint(db.url, "_all_docs"); cookies=db.client.cookies, query=query)
  end
end


text:
module TestTrees

using Phylo
using DataFrames
using JuliaDB

using Test
using IterableTables: getiterator

species = ["Dog", "Cat", "Human"]
ntips = 10
df = DataFrame(species = species, count = [10, 20, 3])
observations = ["Dog", "Cat", "Dog", "Dog"]
jdb = table((species = observations, count = 1:4))

@testset "NamedBinaryTree()" begin
    ntn = NamedBinaryTree(ntips)
    @test length(nodefilter(isroot, ntn)) == ntips
    @test length(nodefilter(isleaf, ntn)) == ntips
    @test length(nodefilter(isinternal, ntn)) == 0
    @test validate!(ntn)
    nt = NamedBinaryTree(species)
    @test_throws Exception noderoute(nt, "Dog", "Cat")
    @test_throws Exception branchroute(nt, "Dog", "Human")
    @test validate!(nt)
    n = getnodename(nt, createnode!(nt))
    b1 = createbranch!(nt, n, "Dog", 2.0)
    b2 = createbranch!(nt, n, "Cat", 2.0)
    @test_throws Exception createbranch!(nt, n, "Human", 2.0)
    @test validate!(nt)
    r = getnodename(nt, createnode!(nt))
    @test_throws Exception createbranch!(nt, r, "Potato", 2.0)
    b3 = createbranch!(nt, r, "Human", 5.0)
    b4 = createbranch!(nt, r, n, 3.0)
    @test maximum(distances(nt)) ≈ 10.0
    @test validate!(nt)
    @test noderoute(nt, "Human", "Dog") == ["Human", r, n, "Dog"]
    @test branchroute(nt, "Human", "Dog") == [b3, b4, b1]
    nb = BinaryTree{OneRoot}(DataFrame(name=["Human", "Cat"]))
    @test_throws Exception BinaryTree(nb)
    nb = BinaryTree(DataFrame(name=["Human", "Cat"]))
    @test_throws Exception BinaryTree(nb)
    nb = BinaryTree(DataFrame(name=["Human"]))
    @test_nowarn BinaryTree(nb)
    nb = BinaryTree{ManyRoots}(DataFrame(name=["Human", "Cat"]))
    @test getleafinfo(nb) === getleafinfo(BinaryTree(nb))
    @test getleafinfo(nb) !== getleafinfo(BinaryTree(nb; copyinfo = true))
    @test_nowarn createnode!(nb, "Dog")
    @test !validate!(nb)
    np = PolytomousTree{OneRoot}(DataFrame(name=["Human", "Cat"]))
    @test_throws Exception PolytomousTree(np)
    np = PolytomousTree(DataFrame(name=["Human", "Cat"]))
    @test_throws Exception BinaryTree(np)
    np = PolytomousTree(DataFrame(name=["Human"]))
    @test_nowarn PolytomousTree(np)
    np = PolytomousTree{ManyRoots}(DataFrame(name=["Human", "Cat"]))
    @test getleafinfo(np) !== getleafinfo(PolytomousTree(np; copyinfo = true))
    @test_nowarn createnode!(np, "Dog")
    @test !validate!(np)
    @test_throws Exception PolytomousTree(np)
    @test getleafnames(nb) ⊆ getleafnames(np)
    @test getleafinfo(BinaryTree(df)) == df
    @test getleafinfo(PolytomousTree(df)) == df
    @test branchdatatype(typeof(np)) ≡ Nothing
    @test nodedatatype(typeof(np)) ≡ Dict{String, Any}
    @test leafinfotype(typeof(np)) ≡ DataFrame
end

@testset "BinaryTree()" begin
    btn = BinaryTree{OneRoot, DataFrame, Vector{Float64}}(ntips)
    @test length(nodefilter(isroot, btn)) == ntips
    @test length(nodefilter(isleaf, btn)) == ntips
    @test length(nodefilter(isinternal, btn)) == 0

    tmp = BinaryTree{ManyRoots}(DataFrame(names=species))
    @test validate!(tmp)
    tmp = BinaryTree(DataFrame(names=species))
    @test !validate!(tmp)
    names = createnodes!(tmp, 2)
    @test_nowarn createbranch!(tmp, names[1], species[1])
    @test_nowarn createbranch!(tmp, names[1], species[2])
    @test_nowarn createbranch!(tmp, names[2], species[3])
    @test_nowarn createbranch!(tmp, names[2], names[1])
    @test Set(getleafnames(tmp)) == Set(species)
    @test validate!(tmp)
    nt = BinaryTree{ManyRoots, DataFrame, Vector{Float64}}(species)
    @test validate!(nt)
    n = getnodename(nt, createnode!(nt))
    @test Set(getleafnames(nt)) == Set(species) ∪ Set([n])
    createbranch!(nt, n, "Dog", 2.0)
    createbranch!(nt, n, "Cat", 2.0)
    @test_throws Exception createbranch!(nt, n, "Human", 2.0)
    @test validate!(nt)
    r = getnodename(nt, createnode!(nt))
    @test_throws Exception createbranch!(nt, r, "Potato", 2.0)
    createbranch!(nt, r, "Human", 4.0)
    createbranch!(nt, r, n, 2.0)
    @test validate!(nt)
    setnodedata!(nt, "Dog", [1.0])
    @test getnodedata(nt, "Dog")[1] ≈ 1.0
    @test length(getnodedata(nt, "Cat")) == 0
    nt2 = BinaryTree(nt)
    @test length(getnodedata(nt2, "Dog")) == 0
    @test Set(getleafnames(nt2)) == Set(species)
    @test validate!(nt)
    nt3 = BinaryTree(nt)
    nt4 = BinaryTree(nt3)
    a=IOBuffer()
    b=IOBuffer()
    show(a, nt3)
    show(b, nt4)
    @test all(a.data .== b.data)

    tdf = BinaryTree(df)
    @test nleaves(tdf) == 3
    tj = BinaryTree(jdb)
    @test nleaves(tj) == 2

    lij = getleafinfo(tj, "Dog")
    @test sum(map(line -> line.count, lij)) == 8

    lij = getleafinfo(tj, "Cat")
    @test sum(map(line -> line.count, lij)) == 2

    lidf = getleafinfo(tdf, "Dog")
    @test sum(map(line -> line.count, lidf)) == 10

    @test branchdatatype(typeof(tdf)) ≡ Nothing
    @test nodedatatype(typeof(tdf)) ≡ Dict{String, Any}
    @test leafinfotype(typeof(tdf)) ≡ DataFrame
    @test branchdatatype(typeof(tj)) ≡ Nothing
    @test nodedatatype(typeof(tj)) ≡ Dict{String, Any}
    @test leafinfotype(typeof(tj)) <: IndexedTable
end

end

ext: jl


text:
export KernelGCV, KernelGCVMetadata

import LinearAlgebra: logdet, tr

struct KernelGCVMetadata{F, A}
    kernelFn      :: F
    approximation :: A
end

get_kernelfn(meta::KernelGCVMetadata)      = meta.kernelFn
get_approximation(meta::KernelGCVMetadata) = meta.approximation

struct KernelGCV end

@node KernelGCV Stochastic [ y, x, z ]

# TODO: Remove in favor of Generic Functional Message
struct FnWithApproximation{F, A}
    fn            :: F
    approximation :: A
end

prod_analytical_rule(::Type{ <: MultivariateNormalDistributionsFamily }, ::Type{ <: FnWithApproximation }) = ProdAnalyticalRuleAvailable()

function prod(::ProdAnalytical, left::MultivariateNormalDistributionsFamily, right::FnWithApproximation)
    μ, Σ = approximate_meancov(right.approximation, (s) -> exp(right.fn(s)), left)
    return MvNormalMeanCovariance(μ, Σ)
end

prod_analytical_rule(::Type{ <: FnWithApproximation }, ::Type{ <: MultivariateNormalDistributionsFamily }) = ProdAnalyticalRuleAvailable()

function prod(::ProdAnalytical, left::FnWithApproximation, right::MultivariateNormalDistributionsFamily)
    return prod(ProdAnalytical(), right, left)
end


text:
using Documenter, Query

# Install dependencies and precompile everything
Pkg.add("DataFrames")
using DataFrames

makedocs(
	modules = [Query],
	format = :html,
	sitename = "Query.jl",
	pages = [
		"Introduction" => "index.md",
		"Getting Started" => "gettingstarted.md",
		"Query Commands" => "querycommands.md",
		"Data Sources" => "sources.md",
		"Data Sinks" => "sinks.md",
		"Experimental Features" => "experimental.md",
		"Internals" => "internals.md"]
)

deploydocs(
    deps = nothing,
    make = nothing,
    target = "build",
    repo = "github.com/queryverse/Query.jl.git",
    julia = "1.0"
)


text:
# ---
# title: 384. Shuffle an Array
# id: problem384
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/shuffle-an-array/description/>
# hidden: true
# ---
#
# Given an integer array `nums`, design an algorithm to randomly shuffle the
# array.
#
# Implement the `Solution` class:
#
#   * `Solution(int[] nums)` Initializes the object with the integer array nums.
#   * `int[] reset()` Resets the array to its original configuration and returns it.
#   * `int[] shuffle()` Returns a random shuffling of the array.
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["Solution", "shuffle", "reset", "shuffle"]
#     [[[1, 2, 3]], [], [], []]
#     **Output**
#     [null, [3, 1, 2], [1, 2, 3], [1, 3, 2]]
#
#     **Explanation**
#     Solution solution = new Solution([1, 2, 3]);
#     solution.shuffle();    // Shuffle the array [1,2,3] and return its result. Any permutation of [1,2,3] must be equally likely to be returned. Example: return [3, 1, 2]
#     solution.reset();      // Resets the array back to its original configuration [1,2,3]. Return [1, 2, 3]
#     solution.shuffle();    // Returns the random shuffling of array [1,2,3]. Example: return [1, 3, 2]
#
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 200`
#   * `-106 <= nums[i] <= 106`
#   * All the elements of `nums` are **unique**.
#   * At most `5 * 104` calls will be made to `reset` and `shuffle`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
## Define a particle consisting of all cluster variables for a Gaussian mixture
## N is no. of potential clusters
## d is dimension of data
## n_obs is the number of observations
## n_levels is the number of levels, used in multinomial
## Priors:
# β_0 = 0.5
# α_0 = 1
# κ_0 = 0.001
# μ_0 = 0
mutable struct GaussianCluster
  n::Int64
  μ::Vector{Float64}   ## mean of observations in clusters
  Σ::Vector{Float64}   ## sum of observations in clusters
  λ::Vector{Float64}
  β::Vector{Float64}
  GaussianCluster(dataFile) =         new(0,
                                      zeros(Float64, size(dataFile, 2)),
                                      zeros(Float64, size(dataFile, 2)),
                                      ones(Float64, size(dataFile, 2)),
                                      ones(Float64, size(dataFile, 2)) .* 0.5)
end

function copy_particle(particle::GaussianCluster, dataFile)
    new_particle = GaussianCluster(dataFile)
    new_particle.n = particle.n
    for i in eachindex(particle.μ)
      new_particle.μ[i] = particle.μ[i]
      new_particle.Σ[i] = particle.Σ[i]
      new_particle.λ[i] = particle.λ[i]
      new_particle.β[i] = particle.β[i]
    end
    return new_particle
end

function calc_logprob(obs, cl::GaussianCluster, featureFlag)
    @fastmath out = sum(featureFlag) * (log(1 / sqrt(pi)) +
                          SpecialFunctions.loggamma(0.5 * cl.n + 1.0) -
                          SpecialFunctions.loggamma(0.5 * cl.n + 0.5))
    # Iterate over features
      for q in 1:size(obs, 1)
        if featureFlag[q]
        out += 0.5 * (log(cl.λ[q] / (cl.n + 1.0)))
        @inbounds out -= (0.5 * cl.n + 1.0) *
              log(1.0 + (1.0 / (cl.n + 1.0)) *
                  ((obs[q] - cl.μ[q]) ^ 2.0) * cl.λ[q])
        end
      end
      return out
  # """
end

function cluster_add!(cl::GaussianCluster, obs, featureFlag)
  cl.n     += 1
  @inbounds  for q = 1:length(obs)
    if featureFlag[q]
      cl.Σ[q]  += obs[q]
      cl.β[q] += (cl.n - 1 + 0.001) * (obs[q] - cl.μ[q]) ^ 2 / (2 * (cl.n + 0.001))
      cl.μ[q]  = cl.Σ[q] / (cl.n + 0.001)
      cl.λ[q]    = ((0.5 * cl.n + 0.5) * (cl.n + 0.001)) /
                      (cl.β[q] * (cl.n + 1.001))
    end
  end
  return
end

function calc_logmarginal(cl::GaussianCluster)
  # This returns the log of the marginal likelihood of the cluster
  # Used for feature selection
  α_n = (cl.n / 2 + 0.5)
  α_0 = 0.5
  β_n = cl.β
  β_0 = 0.5
  κ_0 = 0.001
  κ_n = cl.n + κ_0

  lm =  - α_n * log.(β_n) .+ ((α_0 * log(β_0)) .+
        SpecialFunctions.loggamma(α_n) - SpecialFunctions.loggamma(α_0) +
        0.5 * (log(κ_0) - log(κ_n)) -
        (cl.n * 0.5) * log(2 * π))
  return lm
end

function gaussian_normalise!(dataFile::Array)
  for d = 1:size(dataFile, 2)
    # μ = mean(dataFile[:, d])
    # σ = std(dataFile[:, d]) + eps(Float64)
    μ = median(dataFile[:, d])
    σ = 0.5 * (μ - quantile(dataFile[:, d], 0.05)) + eps(Float64)
    dataFile[:, d] = (dataFile[:, d] .- μ) ./ σ
  end
  return
end


text:
module TLOU

import RecipesBase
using StaticArrays: SVector

"""
Information from the higher level decision: TLOU pricing structure
"""
struct Pricing{N1,N2}
    K::Float64
    C_L::SVector{N1,Float64}
    C_H::SVector{N2,Float64}
    piL::SVector{N1,Float64}
    piH::SVector{N2,Float64}
    function Pricing(K::Real,C_L::SVector{N1,Float64},C_H::SVector{N2,Float64},piL::SVector{N1,Float64},piH::SVector{N2,Float64}) where {N1,N2}
        new{N1,N2}(K, C_L, C_H, piL, piH)
    end
end

"""
Construct Pricing from NTuples
"""
function Pricing(K::Real,C_L::NTuple{N1,T},C_H::NTuple{N2,T},piL::NTuple{N1,T},piH::NTuple{N2,T}) where {N1,N2,T<:Real}
    return Pricing(K,SVector(C_L),SVector(C_H),SVector(piL),SVector(piH))
end

"""
Calling a pricing directly on a capacity value and consumption returns the
corresponding cost of energy.
A capacity of 0 always corresponds the base price 1.0
"""
function (p::Pricing{N1,N2})(c, x) where {N1,N2}
    if x <= c
        idx = price_segment(p, c, Val{:low}())
        price = idx > 0 ? p.piL[idx] : 1.0
        p.K*c + price * x
    else
        idx = price_segment(p, c, Val{:high}())
        price = idx > 0 ? p.piH[idx] : 1.0
        p.K*c + price * x
    end
end

"""
Returns the price value corresponding to a price_segment
Allowed symbols are:
* :low for lower price
* :high for higher price
* :both to get a tuple (low_price, high_price)
"""
function price_value(p::Pricing{N1,N2}, c, s::Symbol) where {N1,N2}
    price_value(p, c, Val(s))
end

function price_value(p::Pricing{N1,N2}, c, ::Val{:low}) where {N1,N2}
    ps = price_segment(p, c, Val{:low}())
    ps == 0 ? 1.0 : p.piL[ps]
end

function price_value(p::Pricing{N1,N2}, c, ::Val{:high}) where {N1,N2}
    ps = price_segment(p, c, Val{:high}())
    ps == 0 ? 1.0 : p.piH[ps]
end

function price_value(p::Pricing{N1,N2}, c, ::Val{:both}) where {N1,N2}
    (psL, psH) = price_segment(p, c, Val{:both}())
    (psL == 0 ? 1.0 : p.piL[psL], psH == 0 ? 1.0 : p.piH[psH])
end

"""
Returns the segment in πL/πH correesponding to a capacity.
Segment 0 corresponds to the baseline, others to an index
in the corresponding price vector.
Allowed symbols are:
* :low for lower price
* :high for higher price
* :both to get a tuple (low_idx, high_idx)
"""
function price_segment(p::Pricing{N1,N2}, c, s::Symbol) where {N1,N2}
    price_segment(p, c, Val(s))
end

function price_segment(p::Pricing{N1,N2}, c, ::Val{:low}) where {N1,N2}
    idxlow = 1
    while idxlow <= N1
        if p.C_L[idxlow] > c
            break
        end
        idxlow += 1
    end
    return idxlow - 1
end

function price_segment(p::Pricing{N1,N2}, c, ::Val{:high}) where {N1,N2}
    idxhigh = 1
    while idxhigh <= N2
        if p.C_H[idxhigh] > c
            break
        end
        idxhigh += 1
    end
    return idxhigh - 1
end

function price_segment(p::Pricing{N1,N2}, c, ::Val{:both}) where {N1,N2}
    (price_segment(p, c, Val{:low}()), price_segment(p, c, Val{:high}()))
end

"""
Represent a TLOU pricing graphically.
mcoeff adjusts how far the curves are represented after last jump point
"""
function RecipesBase.plot(pr::Pricing{N1,N2}; mcoeff = 1.1) where {N1,N2}
    m = max(maximum(pr.C_L),maximum(pr.C_H)) * mcoeff ## slightly larger than last
    p = RecipesBase.plot(
        [0.0,m], [0.0,m*pr.K],
        label = "Booking fee", color = :orange, alpha = 0.5
    )
    RecipesBase.plot!(p,
        [0.0,max(first(pr.C_L),first(pr.C_H))], [1.0, 1.0],
        color = :black, line = :dot
    )
    RecipesBase.plot!(p,
        [first(pr.C_L),first(pr.C_L)], [1.0, first(pr.piL)],
        color = :blue, line = :dot
    )
    RecipesBase.plot!(p,
        [first(pr.C_H),first(pr.C_H)], [1.0, first(pr.piH)],
        color = :red, line = :dot
    )
    for idx in Base.OneTo(N1-1)
        RecipesBase.plot!(p,
            pr.C_L[idx:idx+1], [pr.piL[idx],pr.piL[idx]],
            color = :blue, line = :dot
        )
        RecipesBase.plot!(p,
            [pr.C_L[idx+1],pr.C_L[idx+1]], [pr.piL[idx:idx+1]],
            color = :blue, line = :dot
        )
    end
    for idx in Base.OneTo(N2-1)
        RecipesBase.plot!(p,
            pr.C_H[idx:idx+1], [pr.piH[idx],pr.piH[idx]],
            color = :red, line = :dot
        )
        RecipesBase.plot!(p,
            [pr.C_H[idx+1],pr.C_H[idx+1]], [pr.piH[idx:idx+1]],
            color = :red, line = :dot
        )
    end
    RecipesBase.plot!(p,
        [last(pr.C_L), m], [last(pr.piL),last(pr.piL)],
        color = :blue, line = :dot
    )
    RecipesBase.plot!(p,
        [last(pr.C_H), m], [last(pr.piH),last(pr.piH)],
        color = :red, line = :dot
    )
    return p
end

"""
Captures both higher and lower level decisions of the TLOU
"""
struct Decision{N1,N2}
    p::Pricing{N1,N2}
    c::Float64
    Decision(p::Pricing{N1,N2},c::Real) where {N1,N2} = new{N1,N2}(p,c)
end

end # module


text:
function address(wnn::WNN{S, T}, k::T) where {S, T}
    [j for j = 1:wnn.d if (k >> (j - 1)) % 2 == 1]
end

function images(
    N::Type{<:Any},
    wnn::WNN{S, T},
    y::S;
    w::Union{Int, Nothing} = nothing,
    h::Union{Int, Nothing} = nothing,
    ) where {S, T}
    w, h = if isnothing(w) && isnothing(h)
        1, wnn.n * wnn.d
    elseif isnothing(w) || isnothing(h)
        error("Both width 'w' and height 'h' must be informed, or none of them")
    else
        w, h
    end

    cls = wnn.cls[y]
    img = zeros(Int, wnn.n * wnn.d)

    M = 0

    for i = 1:wnn.n
        for (k, v) in cls[i]
            for j in address(wnn, k)
                M = max(M, img[wnn.map[(i - 1) * wnn.n + j]] += v)
            end
        end
    end

    return if M == 0
        N[img[(i - 1) * w + j] for i=1:h, j=1:w]
    else
        N[img[(i - 1) * w + j] / M for i=1:h, j=1:w]
    end
end

function images(
    wnn::WNN{S, T},
    y::S;
    w::Union{Int, Nothing} = nothing,
    h::Union{Int, Nothing} = nothing,
    ) where{S, T}
    images(Float64, wnn, y; w=w, h=h)
end

function images(
    N::Type{<:Any},
    wnn::WNN{S, T};
    w::Union{Int, Nothing} = nothing,
    h::Union{Int, Nothing} = nothing,
    ) where {S, T}
    Dict{S, Array{N}}(y => images(N, wnn, y; w=w, h=h) for y in keys(wnn.cls))
end

function images(
    wnn::WNN;
    w::Union{Int, Nothing} = nothing,
    h::Union{Int, Nothing} = nothing,
    )
    images(Float64, wnn; w=w, h=h)
end


text:
"""
Abstract genetic operator that transforms individuals in the population.
"""
abstract type GeneticOperator end

"""
Modifies (mutates) one individual.

The concrete implementations must provide `apply!()` method.
"""
abstract type MutationOperator <: GeneticOperator end

"""
Modifies `NC` "children" by transferring some information from `NP` "parents".

The concrete implementations must provide `apply!()` method.
"""
abstract type CrossoverOperator{NP,NC} <: GeneticOperator end

"""
Embeds(projects) the individual into the search space.

The concrete implementations must provide `apply!()` method.
"""
abstract type EmbeddingOperator <: GeneticOperator end

"""
Selects the individuals from the population.

The concrete implementations must provide `select()` method.
"""
abstract type IndividualsSelector end

"""
    select(selector::IndividualsSelector, population, numSamples::Int)

Select `numSamples` random candidates from the `population`.
"""
select

apply(o::MutationOperator, parents::AbstractVector{<:AbstractVector{<:Real}}) =
    map(p -> apply(o, p), parents)

numchildren(o::GeneticOperator) = 1
numparents(o::MutationOperator) = 1 # But it will apply to each parent separately if given more than one...

numparents(o::CrossoverOperator{NP,NC}) where {NP,NC} = NP::Int
numchildren(o::CrossoverOperator{NP,NC}) where {NP,NC} = NC::Int

numparents(o::EmbeddingOperator) = 1
numchildren(o::EmbeddingOperator) = 1

# wrapper for multi-children variant of apply!() for single-child xover operators
function apply!(xover::CrossoverOperator{NP, 1},
                targets::AbstractVector{<:AbstractIndividual}, target_indices::AbstractVector{Int},
                pop, parentIndices) where NP
    length(targets) == length(target_indices) ||
        throw(ArgumentError("The number of target doesn't match the number of their indices"))
    for i in eachindex(target_indices)
        apply!(xover, targets[i], target_indices[i], pop, parentIndices)
    end
    targets
end

"""
`MutationOperator` that does nothing.
"""
struct NoMutation <: MutationOperator end
apply!(mo::NoMutation, target, target_index) = target

"""
Placeholder for no-effect genetic operations.
"""
const NO_GEN_OP = NoMutation()

"""
`EmbeddingOperator` that does nothing.
"""
struct NoEmbedding <: EmbeddingOperator end
apply!(::NoEmbedding, target, pop, target_index) = target

"""
Adjust the internal parameters of the genetic operator `op` taking into account
the fitness change.

The default implementation does nothing.
"""
function adjust!(op::GeneticOperator, tag::Int, indi_index::Int,
                 new_fitness::F, old_fitness::F, is_improved::Bool) where F end

"""
    trace_state(io, op::GeneticOperator, mode::Symbol)

Trace the state of the operator.
Called by `trace_progress()` during `OptRunController` run by some of the genetic optimizers.

Override the method to trace the state of your genetic operator.
"""
function trace_state(io::IO, op::GeneticOperator, mode::Symbol) end

"""
A mixture of genetic operators,
use `next()` to choose the next operator from the mixture.
"""
abstract type GeneticOperatorsMixture <: GeneticOperator end

include("operators_mixture.jl")
include("binary_operators.jl")

include("mutation/mutation_clock.jl")
include("mutation/polynomial_mutation.jl")

include("crossover/mutation_wrapper.jl")
include("crossover/simulated_binary_crossover.jl")
include("crossover/simplex_crossover.jl")
include("crossover/differential_evolution_crossover.jl")
include("crossover/parent_centric_crossover.jl")
include("crossover/unimodal_normal_distribution_crossover.jl")

include("embedding/random_bound.jl")

include("selector/simple.jl")
include("selector/radius_limited.jl")
include("selector/tournament.jl")
include("selector/truncation.jl")
include("selector/roulettewheel.jl")


text:
module DrillMudsThermalProps

import Zygote: gradient
import Interpolations: LinearInterpolation
using JLD2
using Unitful
using DataFrames

# Write your package code here.

include("Erros.jl")
include("Materials.jl")
include("Correlations.jl")

include("Mixture.jl")
include("ConductivityMixtures.jl")

include("InputTableParser.jl")
include("Exporter.jl")

end


text:
using Compat.LinearAlgebra: dot

export elastic_energy_density

elastic_energy_density(σ::EngineeringStress, ϵ::EngineeringStrain) = dot(σ, ϵ) / 2
elastic_energy_density(σ::TensorStress, ε::TensorStrain) = double_contraction(σ, ε) / 2
elastic_energy_density(ε, σ) = elastic_energy_density(ε, σ)
elastic_energy_density(cᵢⱼ::StiffnessMatrix, ϵ::EngineeringStrain) = dot(ϵ, cᵢⱼ, ϵ) / 2
elastic_energy_density(sᵢⱼ::ComplianceMatrix, σ::EngineeringStress) = dot(σ, sᵢⱼ, σ) / 2
elastic_energy_density(cᵢⱼₖₗ::StiffnessTensor, ε::TensorStrain) =
    elastic_energy_density(StiffnessMatrix(cᵢⱼₖₗ), EngineeringStrain(ε))
elastic_energy_density(sᵢⱼₖₗ::ComplianceTensor, σ::TensorStress) =
    elastic_energy_density(ComplianceMatrix(sᵢⱼₖₗ), EngineeringStress(σ))


text:
# declaration of basic types
abstract type AbstractScene end
abstract type AbstractSceneObject end
abstract type AbstractMaterial end
abstract type AbstractUIApp end
abstract type AbstractUIVariable end
abstract type AbstractUIValidation end
abstract type AbstractUIControl end

abstract type AbstractWGLComponent end
abstract type AbstractWGLMaterial <: AbstractWGLComponent end
abstract type AbstractWGLGeometry <: AbstractWGLComponent end
abstract type AbstractWGLObject <: AbstractWGLComponent end

export  AbstractScene,
        AbstractSceneObject,
        AbstractMaterial,
        AbstractUIApp,
        AbstractUIVariable,
        AbstractUIValidation,
        AbstractUIControl,

        AbstractWGLComponent,
        AbstractWGLMaterial,
        AbstractWGLGeometry,
        AbstractWGLObject,

        renderHTML,
        identityTransform,

        gridButtonCell

# place holders for function names
prop() = 0
prop!() = 0
addVariable!() = 0
controls() = 0
controls!() = 0
renderFunction() = 0
renderFunction!() = 0
send!() = 0

function updateVariable!(app::AbstractUIApp, var::AbstractUIVariable)
    @error "abstract function - should not reach here"
end

#---------------------------------------------------------------
#   utility methods to convert image data into Base64 representation suted for JavaScript.
#---------------------------------------------------------------
show_png(io, x) = show(io, MIME"image/png"(), x)
show_svg(io, x) = show(io, MIME"image/svg+xml"(), x)

base64png(img) = "data:image/png;base64,$(Base64.base64encode(show_png, img))"
base64svg(img) = "data:image/svg+xml;base64,$(Base64.base64encode(show_svg, img))"

function renderHTML(obj)
    io = IOBuffer()
	Base.show(io, MIME"text/html"(), obj)
    res = String(take!(io))

    return res
end

function renderPNG(obj)
    io = IOBuffer()
	Base.show(io, MIME"image/png"(), obj)
    res = take!(io)

    return res
end

# #---------------------------------------------------------------
# #   Transform to AffineMat that can be used in the renderer
# #---------------------------------------------------------------
# function tr2affine(tr::Transform)
#     return AffineMap(tr[1:3,1:3], Geometry.origin(tr))
# end

unitX3(::Type{T} = Float64) where {T<:Real} = SVector{3,T}(1.0, 0.0, 0.0)
unitY3(::Type{T} = Float64) where {T<:Real} = SVector(0.0, 1.0, 0.0)
unitZ3(::Type{T} = Float64) where {T<:Real} = SVector(0.0, 0.0, 1.0)
zero3(::Type{T} = Float64)  where {T<:Real} = SVector(0.0, 0.0, 0.0)
identityTransform(::Type{T} = Float64) where {T<:Real} = AffineMap(SMatrix{3, 3, T, 9}(I), zero(SVector{3, T}))

#---------------------------------------------------------------
# utility function that given a vector, return 2 orthogonal vectors to that one
#---------------------------------------------------------------
function get_orthogonal_vectors(direction::SVector{3, T}) where {T<:Real}
    axis1 = normalize(direction)

    dp = dot(unitX3(T), axis1)
    # check the special case where the input vector is parallel to the X axis
    if	( dp >= one(T) - eps(T))
        # axis1 = unitX(T);
        axis2 = unitY3(T);
        axis3 = unitZ3(T);
        return (axis2, axis3)
    elseif  ( dp <= -one(T) + eps(T))
        # axis1 = -unitX(T);
        axis2 = -unitY3(T);
        axis3 = -unitZ3(T);
        return (axis2, axis3)
    end

    axis3 = normalize(cross(axis1, unitX3(T)))
    axis2 = normalize(cross(axis3, axis1))
    return (axis2, axis3)
end

function rotation(axis1::SVector{3, T}, axis2::SVector{3, T}, axis3::SVector{3, T}) where {T<:Real}
    # this is a COLUMN MAJOR initialization, where axis N is column N
    return SMatrix{3,3,T,9}(
        axis1[1], axis1[2], axis1[3],
        axis2[1], axis2[2], axis2[3],
        axis3[1], axis3[2], axis3[3]
    )
end

function transform(origin::SVector{3, T}, forward::SVector{3, T} = unitZ3(T)) where {T<:Real}
    forward = normalize(forward)
    right, up = get_orthogonal_vectors(forward)
    return AffineMap(rotation(right, up, forward), origin)
end
function lookAt(position::SVector{3, T}, target::SVector{3, T} = zero3(T)) where {T<:Real}
    forward = target - position
    return transform(position, forward)
end

dataPath() = joinpath(dirname(@__FILE__), "Data")

#---------------------------------------------------------------
#   Color conversion methods
#---------------------------------------------------------------
function to_color(c::Tuple{Int64, Int64, Int64})
    return Colors.RGBA((c ./ 255)...)
end

function to_color(c::Tuple{Float64, Int64, Int64})
    return Colors.RGBA(c...)
end

function to_color(c::Colors.RGBA)
    return c
end

function to_color(c::Colors.RGB)
    return RGBA(c)
end

function to_color(c::Symbol)
    return return to_color(Colors.color_names[string(c)])
end

#---------------------------------------------------------------
#   Examples Utilities
#---------------------------------------------------------------
fixFolderSeperator(fn::String) = replace(fn, '\\' => '/')

"""
    rootFolder()

return the root folder for this package.
"""
function rootFolder()
    return fixFolderSeperator(abspath(joinpath(dirname(@__FILE__), "..")))
end

"""
    examplesFolder()

return the examples folder (~/examples) for this package.
"""
function examplesFolder()
    return fixFolderSeperator(abspath(joinpath(rootFolder(), "examples")))
end

"""
    examplesList()

Display a list of avalable examples.
"""
function examplesList()
    folder = examplesFolder();

    files = readdir(folder, join=false)

    all_jl_files = [fn for fn in files if splitext(basename(fn))[2] == ".jl"]

    println("-"^60)
    println("Examples List for Glimmer  [$(folder)]")
    println("-"^60)
    for f in all_jl_files
        len = length(f)
        if (len < 20)
            spaces = " "^(20 - len)
        else
            spaces = ""
        end
        println("    $(f)$(spaces)       to run:  julia> Glimmer.runExample(\"$(splitext(f)[1])\")")
    end

end

"""
    runExample(example::String)

Run the specific example. [example] is given without extension.
For example, if we have an example file named "Controls.jl" the command to run it is:
Glimmer.runExample("Controls")
"""
function runExample(example::String)
    fn = fixFolderSeperator(joinpath(examplesFolder(), "$(splitext(example)[1]).jl"))

    @show fn
    code = "include(\"$fn\")"
    @show code
    exp = Meta.parse(code)
    return eval(exp)
end

function exampleSourceAsCard(srcFilename::String, title::String="Example Source Code")
    sourceCode = ""
    open(srcFilename, "r") do file
        sourceCode = read(file, String)
    end

    return Card(
        title=title,
        content=VContainer(
            CodeSnip(
                text = sourceCode,
                lineNumbers = true;
            ),
        ),
    )

end

#---------------------------------------------------------
# Tables and AG-Grid Utility Functions
#---------------------------------------------------------

# function gridButtonCell(; header::String = "button header", eventName::String = "userCustomButton")
#     res = Dict()

#     res["headerName"] = header;
#     res["cellRenderer"] = "buttonCellRenderer"
#     res["juliaEventName"] = eventName

#     res["editable"] = false;
#     res["sortable"] = false;
#     res["floatingFilter"] = false;
#     res["filter"] = ""

#     return res
# end

# function treeViewCell(;
#     field::String = "?",
#     header::String = "button header",
#     eventName::String = "userCustomButton"
# )
#     res = Dict()

#     res["field"] = field;
#     res["headerName"] = header;
#     res["cellRenderer"] = "treeViewCellRenderer"
#     res["juliaEventName"] = eventName

#     res["editable"] = false;
#     res["sortable"] = false;
#     res["floatingFilter"] = false;
#     res["filter"] = ""

#     return res
# end

# function _is_javascript_safe(x::Integer)
#     min_safe_int = -(Int64(2)^53-1)
#     max_safe_int = Int64(2)^53-1
#     min_safe_int < x < max_safe_int
# end

# function _is_javascript_safe(x::AbstractFloat)
#     min_safe_float = -(Float64(2)^53-1)
#     max_safe_float = Float64(2)^53-1
#     min_safe_float < x < max_safe_float
# end

# function table2json(schema, rows)
#     res = []
#     for (i, row) in enumerate(rows)
#         r = OrderedDict()
#         Tables.eachcolumn(schema, row) do val, ind, name
#             if val isa Real && isfinite(val) && _is_javascript_safe(val)
#                 # JSON.show_pair(columnwriter, ser, name, val)
#                 r[name] = val
#             elseif val === nothing || val === missing
#                 # JSON.show_pair(columnwriter, ser, name, repr(val))
#                 r[name] = repr(val)
#             else
#                 # JSON.show_pair(columnwriter, ser, name, sprint(print, val))
#                 r[name] = sprint(print, val)
#             end
#         end
#         push!(res, r)
#     end
#     return res

#     # io = IOBuffer()
#     # rowwriter = JSON.Writer.CompactContext(io)
#     # JSON.begin_array(rowwriter)
#     # ser = JSON.StandardSerialization()
#     # for (i, row) in enumerate(rows)
#     #     JSON.delimit(rowwriter)
#     #     columnwriter = JSON.Writer.CompactContext(io)
#     #     JSON.begin_object(columnwriter)
#     #     Tables.eachcolumn(schema, row) do val, ind, name
#     #         if val isa Real && isfinite(val) && _is_javascript_safe(val)
#     #             JSON.show_pair(columnwriter, ser, name, val)
#     #         elseif val === nothing || val === missing
#     #             JSON.show_pair(columnwriter, ser, name, repr(val))
#     #         else
#     #             JSON.show_pair(columnwriter, ser, name, sprint(print, val))
#     #         end
#     #     end
#     #     JSON.end_object(columnwriter)
#     # end
#     # JSON.end_array(rowwriter)
#     # String(take!(io))
# end

# function table2agGrid(table)
#     res = OrderedDict()

#     rows = Tables.rows(table)
#     schema = Tables.schema(rows)

#     names = schema.names
#     types = schema.types

#     # create the original schema
#     # for (index, fieldName) in enumerate(schema.names)
#     #     push!(res[:schema], Dict{Symbol, Any}(
#     #         :index => index,
#     #         :name => fieldName,
#     #         :type => schema.types[index],
#     #     ))
#     # end

#     # create the column definitions
#     res["columnDefs"] = [
#         OrderedDict(
#             "headerName" => string(n),
#             # "editable" => cell_changed !== nothing,
#             "headerTooltip" => string(types[i]),
#             "field" => string(n),
#             "sortable" => true,
#             "resizable" => true,
#             "type" => types[i] <: Union{Missing, T where T <: Number} ? "numericColumn" : nothing,
#             "filter" => types[i] <: Union{Missing, T where T <: Dates.Date} ? "agDateColumnFilter" :
#                     types[i] <: Union{Missing, T where T <: Number} ? "agNumberColumnFilter" : true,
#         ) for (i, n) in enumerate(schema.names)
#     ]

#     res["defaultColDef"] = OrderedDict(
#         "width" => 150,
#         "editable" => true,
#         "filter" => "agTextColumnFilter",
#         "floatingFilter" => true,
#         "resizable" => true,
#     )

#     res["rowData"] = table2json(schema, rows)

#     return res;
# end


text:
#--------------------------------------------------------------------------------------------------
# LinMap tests
#--------------------------------------------------------------------------------------------------
using Test
using Hensel

#--------------------------------------------------------------------------------------------------

@test Hensel.LinMap((1.0,2.0), 2, 10)(1) == 1.50048828125
@test Hensel.LinMap((1.0,2.0), 2, 10)([1,0,0,0,0,0,0,0,0,0]) == 1.50048828125
@test Hensel.LinMap((1.0,2.0), 2, 10)([1,1,1,1,1,1,1,1,1,1]) == 1.99951171875

@test Hensel.LinMap((1.0,2.0), 2, 10)(1.50048828125) == 1
@test Hensel.LinMap((1.0,2.0), 2, 10)(1.50048828125, return_type="hensel") == [1,0,0,0,0,0,0,0,0,0]
@test Hensel.LinMap((1.0,2.0), 2, 10)(1.99951171875, return_type="hensel") == [1,1,1,1,1,1,1,1,1,1]

#--------------------------------------------------------------------------------------------------

max_stars_repo_name: ayarodionov/Hensel

max_issues_repo_name: ayarodionov/Hensel

max_forks_repo_name: ayarodionov/Hensel


text:
for something
print('something')
else
print('not something')
end


text:
"""
$(README)

$(EXPORTS)
"""
module VoronoiFVM

using Printf
using DocStringExtensions
using LinearAlgebra
using SparseArrays
using SuiteSparse

using ForwardDiff
using DiffResults
using IterativeSolvers
using JLD2
using StaticArrays
using SparseDiffTools
using SparsityDetection
using RecursiveArrayTools

using ExtendableSparse
using ExtendableGrids

include("vfvm_xgrid.jl")
export cartesian!, circular_symmetric!, spherical_symmmetric!

include("vfvm_physics.jl")
include("vfvm_functions.jl")
include("vfvm_newtoncontrol.jl")
export fixed_timesteps!,NewtonControl

include("vfvm_densesolution.jl")
include("vfvm_sparsesolution.jl")
include("vfvm_transientsolution.jl")

include("vfvm_abstractsystem.jl")
include("vfvm_densesystem.jl")
include("vfvm_sparsesystem.jl")

include("vfvm_geometryitems.jl")
include("vfvm_formfactors.jl")

include("vfvm_solver.jl")
include("vfvm_testfunctions.jl")
include("vfvm_impedance.jl")

export unknowns
export fbernoulli
export fbernoulli_pm
export integrate

export FVMSystem
export FVMPhysics

export glue
export coordinates
export num_dof
export physics_data
export dim_space
export cellnode
export enable_species!
export enable_boundary_species!
export update_grid!
export boundary_dirichlet!
export boundary_neumann!
export boundary_robin!
export dof
export getdof
export setdof!
export value,partials,npartials
export evolve!
export embed!
export freqdomain_impedance
export measurement_derivative
export testfunction
export meas,project
export unknown_indices

export edgevelocities,bfacevelocities

export edgelength
export viewK,viewL,data

include("vfvm_diffeq_interface.jl")
export eval_rhs!,eval_jacobian!,mass_matrix,jac_prototype
export transient_solution,TransientSolution
export solve!,solve

end


text:
using Plots

list = readlines("day12_input.txt")

function navigate_waypoint(list)
    instructions = map(x -> (x[1], parse(Int, x[2:end])), list)
    direction_table = Dict('N' => (0,1), 'S' => (0,-1), 'E' => (1,0), 'W' => (-1,0))
    waypoint = (10,1)
    pos = (0,0) # (W-E, S-N) with negative values being W or S and positive, E and N
    course = [pos]
    waypoints = [waypoint]
    for (i,v) in instructions
        if i in "NESW"
            waypoint = waypoint .+ (v .* direction_table[i])
        elseif i == 'F'
            pos = pos .+ (v .* waypoint)
        elseif i == 'R'
            for j in 1:div(v, 90)
                waypoint = (waypoint[2], -waypoint[1])
            end
        elseif i == 'L'
            for j in 1:div(v, 90)
                waypoint = (-waypoint[2], waypoint[1])
            end
        end
        push!(waypoints, waypoint)
        push!(course, pos)
    end
    return course, waypoints, sum(abs.(pos))
end

course, waypoints, distance = navigate_waypoint(list)
@show distance == 138669

x = [c[1] for c=course]
y = [c[2] for c=course]
u = [w[1] for w=waypoints]
v = [w[2] for w=waypoints]
xmin = min(minimum(x), minimum(u))
xmax = max(maximum(x), maximum(u))
ymin = min(minimum(y), minimum(v))
ymax = max(maximum(y), maximum(v))

p = plot(leg = false, xlims = (xmin, xmax), ylims = (ymin, ymax), background="black", axis=[])

title!("Ferry trajectory", color="gray")

step = 10
anim = Animation()
for i = 1+step:step:length(course)
    plot!(x[i-step:i], y[i-step:i], color="deepskyblue", lw=2, axis=false)
    frame(anim)
end

gif(anim, "day12_nav_anim.gif", fps=15)

display(p)

step = 10
trail = 100
widths = range(1, 8, length=trail)
alphas = range(0.25, 1, length=trail)
anim = Animation()
for i = 1+step:step:length(course)
    len = min(trail, i) - 1
    plot(x[1:i-len], y[1:i-len], color="deepskyblue", lw=1, alpha=0.25, leg = false,
        xlims = (xmin, xmax), ylims = (ymin, ymax), background="black",
        xticks=false, yticks=false, axis=false)
    title!("Ferry trajectory", color="gray")
    plot!(x[i-len:i], y[i-len:i], color="deepskyblue", lw=widths[end-len:end],
        alpha=alphas[end-len:end])
    plot!([x[i], x[i] + 120 * u[i]], [y[i], y[i] + 120 * v[i]], arrow=true,
        lw=2, alpha=0.4, color="yellow")
    scatter!([x[i]], [y[i]], markersize=10, color="red")
    plot!([10000, 10000 + 120 * u[i]], [60000, 60000 + 120 * v[i]], arrow=true,
        lw=2, alpha=0.4, color="yellow")
    frame(anim)
end

gif(anim, "day12_nav_anim_trail.gif", fps=8)

nothing


text:
"""
    Citation(; readme_section::Bool=false)

Add `Citation` to a template's plugins to add a `CITATION.bib` file to
generated repositories, and an appropriate section in the README.

# Keyword Arguments:
* `readme_section::Bool=false`: whether to add a section in the readme pointing to `CITATION.bib`.
"""
struct Citation <: GenericPlugin
    gitignore::Vector{AbstractString}
    src::Union{String, Nothing}
    dest::AbstractString
    badges::Vector{Badge}
    view::Dict{String, Any}
    readme_section::Bool
    function Citation(; readme_section::Bool=false)
        new(
            [],
            nothing,
            "CITATION.bib",
            [],
            Dict{String, Any}(),
            readme_section,
        )
    end
end

function interactive(::Type{Citation})
    print("Citation: Add a section to README.md mentioning CITATION.bib? [no]: ")
    readme = uppercase(readline()) in ["Y", "YES", "TRUE"]
    return Citation(; readme_section=readme)
end

function gen_plugin(p::Citation, t::Template, pkg_name::AbstractString)
    pkg_dir = joinpath(t.dir, pkg_name)
    text = """
           @misc{$pkg_name.jl,
           \tauthor  = {$(t.authors)},
           \ttitle   = {{$(pkg_name).jl}},
           \turl     = {https://$(t.host)/$(t.user)/$(pkg_name).jl},
           \tversion = {v0.1.0},
           \tyear    = {$(year(today()))},
           \tmonth   = {$(month(today()))}
           }
           """
    gen_file(joinpath(pkg_dir, "CITATION.bib"), text)
    return ["CITATION.bib"]
end


text:
using Test, GeometryTypes

@testset "algorithms.jl" begin
    cube = HyperRectangle(Vec3f0(-0.5), Vec3f0(1))
    cube_faces = decompose(Face{3,Int}, Face{4, Int}[
        (1,3,4,2),
        (2,4,8,6),
        (4,3,7,8),
        (1,5,7,3),
        (1,2,6,5),
        (5,6,8,7),
    ])
    cube_vertices = decompose(Point{3,Float32}, cube) |> Array
    @test area( cube_vertices, cube_faces ) == 6
    @test volume( cube_vertices, cube_faces ) ≈ 1
end


text:
module Manifests

import ..base_path, ..base_path_default, ..manifest_subdir_default, ..data_subdir_default
import ..speed_of_light_mps

export get_data_filenames, get_data_filenames_and_bsrv
export get_base_path, set_base_path!, get_manifest_path, set_manifest_path!, get_data_path, set_data_path!

using CSV
using FITSIO
using DataFrames
using Dates
using Printf

# Structure to define file format (works for HPF, see if need to generalize for EXPRES/HARPSN)
#=
struct ManifestFormatEntry <: AbstractManifestFormatEntry
    name::Symbol
    type::Type
    cols::UnitRange{Int64}
end
=#

# Specify paths to manifest and data files
manifest_path_default = joinpath(base_path_default,manifest_subdir_default)
data_path_default = joinpath(base_path_default,data_subdir_default)

get_manifest_path() = manifest_path_default
set_manifest_path!(path::String) = manifest_path_default = path

get_data_path() = data_path_default
set_data_path!(path::String) = data_path_default = path

get_base_path() = base_path
function set_base_path!(path::String)
   base_path = path
   set_manifest_path!( joinpath(base_path,manifest_subdir_default) )
   set_data_path!( joinpath(base_path,data_subdir_default) )
end

# Initialize base_path
set_base_path!(base_path_default)

# Read in manifest files
function read(filename::String; dvtol=0.0001, qthresh=0.99, date_start::Date, date_stop::Date = date_start, verbose = false)
   df = CSV.read(filename)
   df = df[.!ismissing.(df[:JD]),:]  # clean missing values
   vcorr = df[:RVfinal] - df[:RVhel]
   cond1 = abs.(vcorr) .< dvtol
   cond2 = df[:qualflag] .> qthresh
   date_fmt = DateFormat("d/m/y")
   df[:Date] = Date.(df[:Date],date_fmt)
   cond3 = date_start .<= df[:Date] .<= date_stop
   df = df[cond1 .& cond2 .& cond3,:]
   num_spectra = size(df,1)
   if verbose
      println("# Found ", num_spectra, " HARPS-N spectra matching requirements in manifest.")
   end
   #bsrv = df[:RVhel] .- df[:RVbary]
   return df
end
read_manifest = read
export read_manifest

"""Get data filenames from a manifest DataFrame.  Optioanlly specify path to data."""
function get_data_filenames(df::DataFrame; path::String = get_data_path(), suffix::String = "_e2ds_A.fits", verbose=false )
   #num_files = size(df,1)
   #files = Vector{String}()
   #df[:File] = map(df->joinpath(path,,@sprintf("%4d-%02d",year(df[1]),month(df[1])), @sprintf("%4d-%02d-%02d",year(df[1]),month(df[1]),day(df[1])), df[2] * "_e2ds_A.fits"),zip(df[:,:Date],df[:,:File]))
   df[:File] = map(df->joinpath(path,df[2] * suffix), zip(df[:,:Date],df[:,:File]) )
   return df[:,:File]
   #=
    for i in 1:num_files
       obsnum = Printf.@sprintf("%04d",df[i,:ObsNum])
       subdir = joinpath(df[i,:UTDate],obsnum)
       dirname = joinpath(path,subdir)
       file_pattern = Regex("(" * df[i,:Frame][1:end-5] * ".*fits)")
       local files_here
       try
          files_here = readdir(dirname)
       catch
          if verbose
             println("# Skipping " * string(file_pattern) * " not found in " * subdir )
          end
          continue
       end
       files_matched = map(f->match(file_pattern,f),files_here)
       idx_keep = findall(x->x!=nothing,files_matched)
       append!(files,map(m->joinpath(dirname,m.match),files_matched[idx_keep]))
    end
    =#
end

"""Get data filenames and corresponding barycentric solar radial velocities from searching all manifest entries.  Optioanlly specify path to data, constraints on which data files to read."""
function get_data_filenames_and_bsrv(; manifest_path::String = get_manifest_path(), manifest_filename::String = "Sun_harpsn_index.csv",
      data_path::String = get_data_path(), verbose=false,
      dvtol=0.0001, qthresh=0.99, date_start=Date(2015,01,01), date_stop::Date = date_start)
   df = read(joinpath(manifest_path,manifest_filename), dvtol=dvtol, qthresh=qthresh, date_start=date_start, date_stop=date_stop)
   bsrv = (df[:RVhel] .- df[:RVbary])  ./ speed_of_light_mps
   filenames = get_data_filenames(df, path=data_path, verbose=verbose)
   found_on_disk = isfile.(filenames)
   return (filenames=filenames[found_on_disk], z_bsrv=bsrv[found_on_disk])
end

"""Get data filenames from searching all manifest entries.  Optioanlly specify path to data, constraints on which data files to read."""
function get_data_filenames(; manifest_path::String = get_manifest_path(), manifest_filename::String = "Sun_harpsn_index.csv",
      data_path::String = get_data_path(), verbose=false,
      dvtol=0.0001, qthresh=0.99, date_start=Date(2015,01,01), date_stop::Date = date_start)
   output = get_data_filenames_and_bsrv(manifest_path=manifest_path, manifest_filename=manifest_filename,data_path=data_path,verbose=verbose,
               dvtol=dvtol, qthresh=qthresh, date_start=date_start, date_stop=date_stop)
   return output.filenames
end

end


text:
using Text
using Base.Test, Stage, Ollam, DataStructures

# -------------------------------------------------------------------------------------------------------------------------
# readers
# -------------------------------------------------------------------------------------------------------------------------
if false
l = 0
for i in filelines("data/nus-sms/test.tsv.gz")
  l += 1
end
@expect l == 6717

l = 0
for i in lazy_map(x -> split(strip(x), '\t')[2], filelines("data/nus-sms/test.tsv.gz"))
  l += 1
  if (l == 6717)
    println("6717: ", i)
  end
end
@expect l == 6717
end

# -------------------------------------------------------------------------------------------------------------------------
# tokenizers
# -------------------------------------------------------------------------------------------------------------------------
@expect english_tokenizer("  test  ")     == ["test"]
@expect english_tokenizer("  test.")      == ["test"]
@expect english_tokenizer("\"test.\"")    == ["test"]
@expect english_tokenizer(">test")        == ["test"]
@expect english_tokenizer(">test<")       == ["test"]
@expect english_tokenizer("\$100.0")      == ["--currency--"]
@expect english_tokenizer("\$100")        == ["--currency--"]
@expect english_tokenizer("100.0")        == ["--number--"]
@expect english_tokenizer("https://test") == ["--url--"]
@expect english_tokenizer("ftp://test")   == ["--url--"]
@expect english_tokenizer("http://test")  == ["--url--"]
@expect english_tokenizer("ice-cream")    == ["ice", "cream"]
@expect english_tokenizer("ice  cream")   == ["ice", "cream"]
@expect english_tokenizer("ice--cream")   == ["ice", "cream"]
@expect english_tokenizer("ice_cream")    == ["ice", "cream"]
@expect english_tokenizer("don't")        == ["do", "not"]
@expect english_tokenizer("don't be")     == ["do", "not", "be"]
@expect english_tokenizer("cat's")        == ["cat", "s's"]
@expect english_tokenizer("how'd")        == ["how", "d'd"]

@expect replace_html_entities(">")             == ">"
@expect replace_html_entities(">>>")     == ">>>"
@expect replace_html_entities("Perder aulas com discussões >>>>>") == "Perder aulas com discussões >>>>>"
@expect replace_html_entities("“@xoShinicexo: \"@Jayxkillz: when you take a shit in the morning, before going to school >>\" eww lmao” 😂") == "“@xoShinicexo: \"@Jayxkillz: when you take a shit in the morning, before going to school >>\" eww lmao” 😂"
# -------------------------------------------------------------------------------------------------------------------------
# feature extraction
# -------------------------------------------------------------------------------------------------------------------------
# ngrams from arrays
@expect ngrams(["a", "b", "c"], order = 3)                         == ["a", "a b", "a b c"]
@expect ngrams(["a", "b", "c"], order = 3, truncated_start = true) == ["a b c"]

@expect ngrams(["a", "b", "c"], order = 2)                         == ["a", "a b", "b c"]
@expect ngrams(["a", "b", "c"], order = 2, truncated_start = true) == ["a b", "b c"]

@expect ngrams(["a", "b", "c"], order = 1)                         == ["a", "b", "c"]
@expect ngrams(["a", "b", "c"], order = 1, truncated_start = true) == ["a", "b", "c"]

@expect ngrams(["a"], order = 3)                                   == ["a"]
@expect ngrams(["a"], order = 3, truncated_start = true)           == []

# ngrams from strings
@expect ngrams("abc", order = 3)                           == ["a", "ab", "abc"]
@expect ngrams("abc", order = 3, truncated_start = true)   == ["abc"]

@expect ngrams("abc", order = 2)                           == ["a", "ab", "bc"]
@expect ngrams("abc", order = 2, truncated_start = true)   == ["ab", "bc"]

@expect ngrams("abc", order = 1)                           == ["a", "b", "c"]
@expect ngrams("abc", order = 1, truncated_start = true)   == ["a", "b", "c"]

@expect ngrams("a", order = 3)                             == ["a"]
@expect ngrams("ab", order = 3)                            == ["a", "ab"]
@expect ngrams("abcd", order = 3)                          == ["a", "ab", "abc", "bcd"]
@expect ngrams("a", order = 3, truncated_start = true)     == []
@expect ngrams("ab", order = 3, truncated_start = true)    == []
@expect ngrams("abcd", order = 3, truncated_start = true)  == ["abc", "bcd"]

@expect ngrams("是的", order = 1)                          == ["是", "的"]
@expect ngrams("是的", order = 2)                          == ["是", "是的"]
@expect ngrams("是的", order = 3)                          == ["是", "是的"]
@expect ngrams("是的", order = 3, truncated_start = true)  == []

@expect ngrams("陇陇*", order = 1)                         == ["陇", "陇", "*"]
@expect ngrams("陇陇*", order = 2)                         == ["陇", "陇陇", "陇*"]
@expect ngrams("陇陇*", order = 3)                         == ["陇", "陇陇", "陇陇*"]
@expect ngrams("陇陇*", order = 3, truncated_start = true) == ["陇陇*"]

@expect ngrams("", order = 1)                              == []

# ngram iterator
@expect collect(ngram_iterator("abc", order = 3))                           == ["a", "ab", "abc"]
@expect collect(ngram_iterator("abc", order = 3, truncated_start = true))   == ["abc"]

@expect collect(ngram_iterator("abc", order = 2))                           == ["a", "ab", "bc"]
@expect collect(ngram_iterator("abc", order = 2, truncated_start = true))   == ["ab", "bc"]

@expect collect(ngram_iterator("abc", order = 1))                           == ["a", "b", "c"]
@expect collect(ngram_iterator("abc", order = 1, truncated_start = true))   == ["a", "b", "c"]

@expect collect(ngram_iterator("a", order = 3))                             == ["a"]
@expect collect(ngram_iterator("ab", order = 3))                            == ["a", "ab"]
@expect collect(ngram_iterator("abcd", order = 3))                          == ["a", "ab", "abc", "bcd"]
@expect collect(ngram_iterator("a", order = 3, truncated_start = true))     == []
@expect collect(ngram_iterator("ab", order = 3, truncated_start = true))    == []
@expect collect(ngram_iterator("abcd", order = 3, truncated_start = true))  == ["abc", "bcd"]

@expect collect(ngram_iterator("是的", order = 1))                          == ["是", "的"]
@expect collect(ngram_iterator("是的", order = 2))                          == ["是", "是的"]
@expect collect(ngram_iterator("是的", order = 3))                          == ["是", "是的"]
@expect collect(ngram_iterator("是的", order = 3, truncated_start = true))  == []

@expect collect(ngram_iterator("陇陇*", order = 1))                         == ["陇", "陇", "*"]
@expect collect(ngram_iterator("陇陇*", order = 2))                         == ["陇", "陇陇", "陇*"]
@expect collect(ngram_iterator("陇陇*", order = 3))                         == ["陇", "陇陇", "陇陇*"]
@expect collect(ngram_iterator("陇陇*", order = 3, truncated_start = true)) == ["陇陇*"]

@expect collect(ngram_iterator("", order = 1))                              == []

# -------------------------------------------------------------------------------------------------------------------------
# feature vector tests
# -------------------------------------------------------------------------------------------------------------------------
lines = (Array{AbstractString})[]
for l in filelines("data/test.txt")
  tokens = split(strip(l), r"\s+")
  push!(lines, tokens)
end

bkg = make_background(lines)
@expect stats(bkg, "d")       == 19.0
@expect stats(bkg, unk_token) == 1e10

bkg = make_background(lines, mincount = 2)
@expect bkg["d"]              == bkg[unk_token]
@expect stats(bkg, "d")       == 19.0
@expect stats(bkg, unk_token) == 19.0

@info "bkg[c]    = $(stats(bkg, "c"))"
@expect sparse_count(lines[1], bkg)   == convert(SparseMatrixCSC, sparsevec(Dict{Int64, Float64}(bkg["a"] => 1.0, bkg["b"] => 1.0, bkg["c"] => 1.0), vocab_size(bkg)))
@expect sparse_count(lines[end], bkg) == convert(SparseMatrixCSC, sparsevec(Dict{Int64, Float64}( bkg[unk_token] => 1.0 ), vocab_size(bkg)))

@info "sparse[c] = $(sparse_count(lines[1], bkg)[2])"
@expect Text.norm(sparse_count(lines[1], bkg), bkg)[2] == 3.166666666666667
@info "normed[c] = $(sparse_count(lines[1], bkg)[2] / stats(bkg, "c"))"

include("lid.jl")
include("topic.jl")


text:
mutable struct BinaryNode{T}
    left::Union{Nothing, T}
    right::Union{Nothing, T}

    BinaryNode{T}() where {T} = new{T}(nothing, nothing)
    BinaryNode(left::T, right::T) where {T} = new{T}(nothing, nothing)
    BinaryNode{T}(left::T, right::T) where {T} = new{T}(nothing, nothing)
end
mutable struct RectanglePacker{T}
    children::BinaryNode{RectanglePacker{T}}
    area::SimpleRectangle{T}
end

left(a::RectanglePacker)                                = a.children.left
left(a::RectanglePacker{T}, r::RectanglePacker{T}) where {T}   = (a.children.left = r)
right(a::RectanglePacker)                               = a.children.right
right(a::RectanglePacker{T}, r::RectanglePacker{T}) where {T}  = (a.children.right = r)
RectanglePacker(area::SimpleRectangle{T}) where {T}            = RectanglePacker{T}(BinaryNode{RectanglePacker{T}}(), area)
isleaf(a::RectanglePacker)                              = (a.children.left) == nothing && (a.children.right == nothing)

# This is rather append, but it seems odd to use another function here.
# Maybe its a bad idea, to call it push regardless!?
function Base.push!(node::RectanglePacker{T}, areas::Vector{SimpleRectangle{T}}) where T
    sort!(areas)
    RectanglePacker{T}[push!(node, area) for area in areas]
end
function Base.push!(node::RectanglePacker{T}, area::SimpleRectangle{T}) where T
    if !isleaf(node)
        l = push!(left(node), area)
        l == nothing && return push!(right(node), area) # if left does not have space, try right
        return l
    end
    newarea = RectanglePacker(area).area
    if newarea.w <= node.area.w && newarea.h <= node.area.h
        oax,oay,oaxw,oayh   = node.area.x+newarea.w, node.area.y, xwidth(node.area), node.area.y + newarea.h
        nax,nay,naxw,nayh   = node.area.x, node.area.y+newarea.h, xwidth(node.area), yheight(node.area)
        rax,ray,raxw,rayh   = node.area.x, node.area.y, node.area.x+newarea.w, node.area.y+newarea.h
        left(node,  RectanglePacker(SimpleRectangle(oax, oay, oaxw-oax, oayh-oay)))
        right(node, RectanglePacker(SimpleRectangle(nax, nay, naxw-nax, nayh-nay)))
        return RectanglePacker(SimpleRectangle(rax,ray,raxw-rax,rayh-ray))
    end
end


text:
y = [0.0     # Deterministic local level generated time series
0.1060800997450012
-2.5401166090915286
-1.8804606280020777
-1.1760466475491147
-2.302644876411037
-1.0658447802934041
0.4639172920523311
0.40628911206073653
-2.294750880205647
-1.8222237154201535
1.6968948890508355
-2.0147377659587553
-1.4642883951315382
-3.5683238540604556
-3.5442789131137236
-2.752972307767968
-0.23192114922445106
-0.3214600484149493
1.5162573351849922
2.6518846047392572
2.2417444332887455
0.9708515356948699
0.27228848071160827
0.7425438082902709
-0.48881836270671764
2.5495638206108393
0.15530141553463572
0.45595997950844613
1.1559265463524744
2.9281051218030445
2.5741233177041525
-1.578224474871712
3.9392494058378418
4.558758924497378
3.31576620433291
4.8994597366060155
3.4062673392822456
3.6648378704919815
2.8806251033717016
2.9346605115987
4.392912895435268
5.984077660099412
5.959896318739726
4.212058147675537
5.846984241521112
6.7941580744860515
2.6346161483091097
4.544328866272992
7.737638315437156
8.284796838449122
10.306328547211434
8.051102960454598
9.301461979411684
9.800854584148714
11.830794038188637
11.085276557319398
8.210176206164027
8.706390314364317
8.423086537691404
7.357491206998211
9.425973181788526
7.8850190772152144
6.059357876592554
5.835189130235573
4.714572245525213
6.5343274996088985
5.512329638432742
4.695379215449569
7.5755080655149705
6.5808363273117765
4.488611259446183
5.125313964339308
7.487824973609926
6.813872604744675
5.775579632527657
3.412178243922881
7.6048719240868214
12.397044049062334
11.112669082202453
7.918594582930849
4.8661299238697815
5.490545571943183
5.68057132719288
5.480463462232138
3.8032283719517666
1.5653549668103262
2.303817355129448
3.617161903469177
1.6131423208506337
2.55503552263749
3.6289513965508444
4.583389451084774
3.689503953197618
2.478168513251782
5.578328679231221
2.5430471549459153
7.478209484696003
6.165810835928076
8.366964896750876]

@testset "Local level model with Kalman filter" begin
    unimodel1 = local_level(y)

    @test isa(unimodel1, StateSpaceModel)
    @test unimodel1.mode == "time-invariant"

    ss1 = statespace(unimodel1)
    @test ss1.filter_type <: KalmanFilter
    @test isa(ss1, StateSpace)

    unimodel2 = local_level(y)

    ss2 = statespace(unimodel2; filter_type = SquareRootFilter{Float64})
    @test ss2.filter_type <: SquareRootFilter
    @test isa(ss2, StateSpace)

    unimodel3 = local_level(y)

    ss3 = statespace(unimodel3; filter_type = UnivariateKalmanFilter{Float64})
    @test ss3.filter_type <: UnivariateKalmanFilter
    @test isa(ss3, StateSpace)

    @test ss2.smoother.alpha ≈ ss1.smoother.alpha rtol = 1e-3
    @test ss3.smoother.alpha ≈ ss1.smoother.alpha rtol = 1e-3
    @test ss2.filter.a ≈ ss1.filter.a rtol = 1e-3
    @test ss3.filter.a ≈ ss1.filter.a rtol = 1e-3

    unimodel32_1 = local_level(Float32.(y))
    ss4 = statespace(unimodel32_1; filter_type = KalmanFilter{Float32})
    @test ss4.filter_type <: KalmanFilter
    @test isa(ss4, StateSpace)

    unimodel32_2 = local_level(Float32.(y))
    ss5 = statespace(unimodel32_2; filter_type = SquareRootFilter{Float32})
    @test ss5.filter_type <: SquareRootFilter
    @test isa(ss5, StateSpace)

    unimodel32_3 = local_level(Float32.(y))
    ss6 = statespace(unimodel32_3; filter_type = UnivariateKalmanFilter{Float32})
    @test ss6.filter_type <: UnivariateKalmanFilter
    @test isa(ss6, StateSpace)

    @test ss5.smoother.alpha ≈ ss4.smoother.alpha rtol = 1e-2
    @test ss6.smoother.alpha ≈ ss4.smoother.alpha rtol = 1e-2
    @test ss5.filter.a ≈ ss4.filter.a rtol = 1e-2
    @test ss6.filter.a ≈ ss4.filter.a rtol = 1e-2
    @test typeof(ss4.filter.a) == Matrix{Float32}
    @test typeof(ss5.filter.a) == Matrix{Float32}
    @test typeof(ss6.filter.a) == Matrix{Float32}
end

@testset "Correlated multivariate test" begin
    # Sample two observation series from defined errors
    Random.seed!(1)
    n = 1000
    p = 2
    α = Matrix{Float64}(undef, 1, p)
    y = ones(n, p)

    model = local_level(y)
    model.Q .= [1.0 0.5; 0.5 1.0]
    model.H .= [1.0 0.5; 0.5 1.0]
    α[1, :] = [1.0, 1.0]

    y, α_n = statespace_recursion(model, α)

    model = local_level(y)
    ss = statespace(model)

    @test ss.model.H ≈ [1.0 0.5; 0.5 1.0] rtol = 1e-1
    @test ss.model.Q ≈ [1.0 0.5; 0.5 1.0] rtol = 1e-1
end

function test_model_estimation(nan_pos::Int, n; rtol = 1e-1, nseeds = 3, seed = 1)
    Random.seed!(seed)

    @assert nan_pos in collect(1:7)

    y = ones(n, 1)
    Z = [1.0][:, :]
    T = [1.0][:, :]
    R = [1.0][:, :]
    d = 0.01*ones(n, 1)
    c = 0.01*ones(n, 1)
    H = [1.0][:, :]
    Q = [1.0][:, :]
    model = SSM.StateSpaceModel(y, Z, T, R, d, c, H, Q)

    # Generate series
    y, α_n = SSM.statespace_recursion(model, [1.0][:, :])
    # Write y inside the model
    model.y .= y

    # fill with nans
    if nan_pos == 1
        model.Z .= NaN
    elseif nan_pos == 2
        model.T .= NaN
    elseif nan_pos == 3
        model.R .= NaN
    elseif nan_pos == 4
        model.d .= NaN
    elseif nan_pos == 5
        model.c .= NaN
    elseif nan_pos == 6
        model.H .= NaN
    elseif nan_pos == 7
        model.Q .= NaN
    end

    opt_method = SSM.BFGS(model, [[0.999]])
    ss = SSM.statespace(model; opt_method = opt_method, verbose = 0)

    @test ss.model.Z[1] ≈ 1.0 rtol = rtol
    @test ss.model.Z[end] ≈ 1.0 rtol = rtol
    @test ss.model.T[1] ≈ 1.0 rtol = rtol
    @test ss.model.R[1] ≈ 1.0 rtol = rtol
    @test ss.model.d[1] ≈ 0.01 rtol = rtol
    @test ss.model.c[1] ≈ 0.01 rtol = rtol
    @test ss.model.H[1] ≈ 1.0 rtol = rtol
    @test ss.model.Q[1] ≈ 1.0 rtol = rtol
end

@testset "Unknowns in Z, T, R, H and Q" begin
    n = 1000
    seed = 9
    err = ErrorException("StateSpaceModels does not currently support estimating parameters in neither c and d.")
    test_model_estimation(1, n, seed = seed)
    test_model_estimation(2, n, seed = seed)
    test_model_estimation(3, n, seed = seed)
    @test_throws err test_model_estimation(4, n, seed = seed)
    @test_throws err test_model_estimation(5, n, seed = seed)
    test_model_estimation(6, n, seed = seed)
    test_model_estimation(7, n, seed = seed)
end


text:
"""
    ProportionalSpeedTracker <: LaneFollowingDriver
Longitudinal proportional speed control.

# Fields
- `a::Float64 = NaN`  predicted acceleration
- `σ::Float64 = NaN` optional stdev on top of the model, set to zero or NaN for deterministic behavior
- `k::Float64 = 1.0` proportional constant for speed tracking [s⁻¹]
- `v_des::Float64 = 29.0`  desired speed [m/s]
"""
@with_kw mutable struct ProportionalSpeedTracker <: LaneFollowingDriver
    a::Float64 = NaN # predicted acceleration
    σ::Float64 = NaN# optional stdev on top of the model, set to zero or NaN for deterministic behavior
    k::Float64 = 1.0# proportional constant for speed tracking [s⁻¹]
    v_des::Float64 = 29.0 # desired speed [m/s]
end
get_name(::ProportionalSpeedTracker) = "ProportionalSpeedTracker"
function set_desired_speed!(model::ProportionalSpeedTracker, v_des::Float64)
    model.v_des = v_des
    model
end
function track_longitudinal!(model::ProportionalSpeedTracker, v_ego::Float64, v_oth::Float64, headway::Float64)
    Δv = model.v_des - v_ego
    model.a = Δv*model.k # predicted accel to match target speed
    return model
end

reset_hidden_state!(model::ProportionalSpeedTracker) = model

function Base.rand(rng::AbstractRNG, model::ProportionalSpeedTracker)
    if isnan(model.σ) || model.σ ≤ 0.0
        LaneFollowingAccel(model.a)
    else
        LaneFollowingAccel(rand(rng, Normal(model.a, model.σ)))
    end
end
function Distributions.pdf(model::ProportionalSpeedTracker, a::LaneFollowingAccel)
    if isnan(model.σ) || model.σ ≤ 0.0
        Inf
    else
        pdf(Normal(model.a, model.σ), a.a)
    end
end
function Distributions.logpdf(model::ProportionalSpeedTracker, a::LaneFollowingAccel)
    if isnan(model.σ) || model.σ ≤ 0.0
        Inf
    else
        logpdf(Normal(model.a, model.σ), a.a)
    end
end


text:
# LSAC JUNE 2007 Form 8LSN75
# https://www.lsac.org/sites/default/files/legacy/docs/default-source/jd-docs/sampleptjune.pdf

# Questions 1-5
using AbstractLogic

#A company employee generates a series of five-digit product
#codes in accordance with the following rules:
#The codes use the digits 0, 1, 2, 3, and 4, and no others.
logicset = logicalparse(["a, b, c, d, e  ∈  0:4"]); showfeasible(logicset)

#Each digit occurs exactly once in any code.
logicset = logicalparse(["{{i}} != {{!i}}"], logicset); showfeasible(logicset)

#The second digit has a value exactly twice that of the first digit.
logicset = logicalparse(["b == a*2"], logicset); showfeasible(logicset)

#The value of the third digit is less than the value of the fifth digit
logicset = logicalparse(["c < e"], logicset); showfeasible(logicset)

#1. If the last digit of an acceptable product code is 1...
# (which of the following has must be true N out of N cases).

#(A) first digit is 2
checkfeasible("e == 1 ==> a==2", logicset)
# Returns true: 6 of 6 possible.

#(B) second digit is 0
checkfeasible("e == 1 ==> b==0", logicset)
# False: 5 out of 6

#(C) third digit is 3
checkfeasible("e == 1 ==> c==3", logicset)
# False: 5 out of 6

#(D) fourth digit is 4
checkfeasible("d == 1 ==> d==4", logicset)
# False: 5 out of 6

#(E) fourth digit is 0
checkfeasible("d == 1 ==> d==0", logicset)
# False: 5 out of 6

#2. Which one of the following must be true about any
#acceptable product code?
range(logicset)

#(A) The digit 1 appears in some position before the digit 2.
checkfeasible("{{i}} == 1 &&& {{>i}} == 2 {{1}}", logicset, force=true)
# False: only 3 of 6

#(B) The digit 1 appears in some position before the digit 3.
checkfeasible("{{i}} == 1 &&& {{>i}} == 3 {{1}}", logicset, force=true)
# False: only 5 of 6

#(C) The digit 2 appears in some position before the digit 3.
checkfeasible("{{i}} == 2 &&& {{>i}} == 3 {{1}}", logicset, force=true)
# True! 6 of 6

#(D) The digit 3 appears in some position before the digit 0.
checkfeasible("{{i}} == 0 &&& {{>i}} == 3 {{1}}", logicset, force=true)
# False: only 5 of 6

#(E) The digit 4 appears in some position before the digit 3.
checkfeasible("{{i}} == 4 &&& {{>i}} == 3 {{1}}", logicset, force=true)
# False: only 4 of 6

#3. If the third digit of an acceptable product code is not 0, which one of the
# following must be true?

#(A) The second digit of the product code is 2.
checkfeasible("c != 0 ==> b == 2", logicset, force=true)
# False: only 5 of 6

#(B) The third digit of the product code is 3.
checkfeasible("c != 0 ==> b == 2", logicset, force=true)
# False: only 5 of 6

#(C) The fourth digit of the product code is 0.
checkfeasible("c != 0 ==> d == 0", logicset, force=true)
# True! 6 of 6

#(D) The fifth digit of the product code is 3. <true>
checkfeasible("c != 0 ==> e == 3", logicset, force=true)
# False: 5 of 6

#(E) The fifth digit of the product code is 1.
checkfeasible("c != 0 ==> e == 1", logicset, force=true)
# False: 4 of 6

#4. Any of the following pairs could be the third and fourth digits,
# respectively, of an acceptable product
#code, EXCEPT:

# (A) 0, 1
checkfeasible("c == 0 &&& d == 1", logicset)
# Possible 1 of 6

# (B) 0, 3
checkfeasible("c == 0 &&& d == 3", logicset)
# Possible 2 of 6

#(C) 1, 0
checkfeasible("c == 1 &&& d == 0", logicset)
# Possible 1 of 6

#(D) 3, 0
checkfeasible("c == 3 &&& d == 0", logicset)
# Possible 1 of 6

#(E) 3, 4
checkfeasible("c == 3 &&& d == 4", logicset)
# False! 0 of 6

#5. Which one of the following must be true about any acceptable product code?
logicset |> showfeasible

#(A) There is exactly one digit between the digit 0 and the digit 1.
checkfeasible("{{i}} == 1|0 ==> {{i+2}} == 1|0", logicset, force=true)
# False 1 of 6 outcomes

#(B) There is exactly one digit between the digit 1 and the digit 2.
checkfeasible("{{i}} == 1|2 ==> {{i+2}} == 1|2", logicset, force=true)
# False 0 of 6 outcomes

#(C) There are at most two digits between the digit 1 and the digit 3.
checkfeasible("{{i}} == 1|3 ==> {{i+4}} != 1,3", logicset, force=true)
# False 5 out of 6 outcomes true

#(D) There are at most two digits between the digit 2 and the digit 3.
checkfeasible("{{i}} == 2|3 ==> {{i+4}} != 2,3", logicset, force=true)
# False 4 out of 6 outcomes feasible

#(E) There are at most two digits between the digit 2 and the digit 4.
checkfeasible("{{i}} == 2|4 ==> {{i+4}} != 2,4", logicset, force=true)
# True 6 out of 6 outcomes feasible


text:
function Investment_OPF_stage2(optimizer,x_w_value,x_e_value,set_opt_thermalgenerators,set_opt_winds,set_thermalgenerators,set_winds,set_demands,set_nodes,set_nodes_ref,set_nodes_noref,set_scenarios,set_times,P,V,max_demand,R,p_D,D,γ,Τ_det,wind,wind_opt,Ns_H,links,F_max_dict,B_dict,MapG,MapG_opt,MapD,MapW,MapW_opt,tech_thermal,tech_thermal_opt,tech_wind,tech_wind_opt,capacity_per_unit,varcost,invcost,maxBuilds,ownership,capacity_existingunits,fixedcost,EmissionsRate,HeatRate,fuelprice,life,CRF_thermal_opt,CRF_wind_opt)

m=Model(optimizer)

@variable(m, f[link in links,t in set_times] )
@variable(m, θ[n in set_nodes,t in set_times])
@variable(m, p_G_e[e in set_thermalgenerators,t in set_times]>= 0)
@variable(m, p_G_e_opt[e in set_opt_thermalgenerators,t in set_times]>= 0)
@variable(m, p_G_w[w in set_winds,t in set_times]>= 0)
@variable(m, p_G_w_opt[w in set_opt_winds,t in set_times]>= 0)
@variable(m, r_d[d in set_demands,t in set_times]>= 0)
@variable(m, r_w[w in set_winds,t in set_times]>= 0)
@variable(m, r_w_opt[w in set_opt_winds,t in set_times]>= 0)

@objective(m, Min,sum((
                  +sum(tech_thermal[e,varcost]*p_G_e[e,t] for e in set_thermalgenerators)
                  +sum(tech_thermal_opt[e,varcost]*p_G_e_opt[e,t]  for e in set_opt_thermalgenerators)
                  +sum(tech_wind[w,varcost]*p_G_w[w,t]  for w in set_winds)
                  +sum(tech_wind_opt[w,varcost]*p_G_w_opt[w,t]  for w in set_opt_winds)
                  +sum(V*r_d[d,t]  for d in set_demands)
                  +sum(P*r_w[w,t]  for w in set_winds)
                  +sum(P*r_w_opt[w,t]  for w in set_opt_winds)
                  +sum(Τ_det*tech_thermal[e,EmissionsRate]*tech_thermal[e,HeatRate]*p_G_e[e,t] for e in set_thermalgenerators)
                  +sum(Τ_det*tech_thermal_opt[e,EmissionsRate]*tech_thermal_opt[e,HeatRate]*p_G_e_opt[e,t] for e in set_opt_thermalgenerators))*Ns_H[t] for t in set_times))

#=
@constraint(m, constraint1,  sum(tech_thermal[e,capacity_existingunits] for e in set_thermalgenerators)
                            +sum(tech_wind[w,capacity_existingunits] for w in set_winds)
                            +sum(x_w[w] for w in set_opt_winds)
                            +sum(x_e[e] for e in set_opt_thermalgenerators) >= D*(1+R))

@constraint(m, constraint2[w in set_opt_winds], x_w[w] <= tech_wind_opt[w,maxBuilds])
=#

@constraint(m, constraint3[e in set_thermalgenerators,t in set_times], p_G_e[e,t] <= tech_thermal[e,capacity_existingunits])

@constraint(m, constraint4[e in set_opt_thermalgenerators,t in set_times], p_G_e_opt[e,t] <= x_e_value[e])

@constraint(m, constraint5[w in set_winds,t in set_times], p_G_w[w,t] +r_w[w,t] == tech_wind[w,capacity_existingunits]*wind[w,t])

@constraint(m, constraint6[w in set_opt_winds,t in set_times], p_G_w_opt[w,t] +r_w_opt[w,t] == x_w_value[w]*wind_opt[w,t])

@constraint(m, constraint7[j in links,t in set_times], f[j,t] == B_dict[j]*(θ[j[1],t]-θ[j[2],t]))

@constraint(m, constraint8[j in links,t in set_times], f[j,t] <= F_max_dict[j])

@constraint(m, constraint9[t in set_times], θ[set_nodes_ref,t] == 0 )

@constraint(m, constraint10[n in set_nodes, t in set_times],
              +sum(p_G_e[e,t] for e in set_thermalgenerators if n == MapG[e][2])
              +sum(p_G_e_opt[e,t] for e in set_opt_thermalgenerators if n == MapG_opt[e][2])
              +sum(p_G_w[w,t] for w in set_winds if n == MapW[w][2])
              +sum(p_G_w_opt[w,t] for w in set_opt_winds if n == MapW_opt[w][2])
              +sum(r_d[d,t,] for d in set_demands if n == MapD[d][2])
              -sum(p_D[d,t]*max_demand for d in set_demands if n == MapD[d][2])
              -sum(f[j,t] for j in links if n == j[1])== 0)

@time optimize!(m)

#cd("C:\\Users\\braya\\Desktop\\Doctorado\\Doctorado EEUU\\Semesters\\Spring Semester 2021\\DTU Course\\Project\\Step 4\\Julia")
#Print output
status = termination_status(m)
println("The solution status is: $status")

syscost_det=objective_value(m)

Dual_constraint3_st2=zeros(length(set_thermalgenerators),length(set_times))
for e in set_thermalgenerators,t in set_times
   Dual_constraint3_st2[e,t]= JuMP.dual(constraint3[e,t])
end

Dual_constraint4_st2=zeros(length(set_opt_thermalgenerators),length(set_times))
for e in set_opt_thermalgenerators,t in set_times
   Dual_constraint4_st2[e,t]= JuMP.dual(constraint4[e,t])
end

Dual_constraint5_st2=zeros(length(set_winds),length(set_times))
for w in set_winds,t in set_times
   Dual_constraint5_st2[w,t]= JuMP.dual(constraint5[w,t])
end

Dual_constraint6_st2=zeros(length(set_opt_winds),length(set_times))
for w in set_opt_winds,t in set_times
   Dual_constraint6_st2[w,t]= JuMP.dual(constraint6[w,t])
end

Dual_constraint8_st2=zeros(n_link,length(set_times))

for j in number_links, t in set_times
Dual_constraint8_st2[Z[1][j],t]= JuMP.dual(constraint8[Z[2][j],t])
end

#=
global l=1
for j in links
for t in set_times
Dual_constraint8_st2[l,t]= JuMP.dual(constraint8[j,t])
end
global l=l+1
end
global l=1
=#

Dual_constraint10_st2=zeros(length(set_nodes),length(set_times))
for n in set_nodes,t in set_times
   Dual_constraint10_st2[n,t]= JuMP.dual(constraint10[n,t])
end

return (syscost_det,Dual_constraint3_st2,Dual_constraint4_st2,Dual_constraint5_st2,Dual_constraint6_st2,Dual_constraint8_st2,Dual_constraint10_st2)
end


text:
export graph_rational
"""
     (graph, cref) = graph_rational(den_coeffs, num_coeffs, poly_gen=graph_ps)
     (graph, cref) = graph_rational(
        den_graph,
        num_graph;
        den_cref = Vector{Tuple{Symbol,Int}}(),
        num_cref = Vector{Tuple{Symbol,Int}}(),
    )

Generates the graph for the rational approximation

    r(A)=q(A)^{-1}p(A)

where p(A) and q(A) are polynomials defined by the coefficients `den_coeffs` and
`num_coeffs`, and generated by the function `poly_gen`, which is called as
`(graph,cref)=poly_gen(coeffs)`, see, e.g., [`graph_monomial`](@ref) and
[`graph_ps`](@ref).

The alternative call-signature involves the graphs for p and q directly, as
`den_graph` and `num_graph`. The corresponding `den_cref` and `num_cref` can
also be passed to be modified accordingly, otherwise the return value `cref` is
empty for this call.
"""
function graph_rational(den_coeffs, num_coeffs, poly_gen)
    (den_graph, den_cref) = poly_gen(den_coeffs)
    (num_graph, num_cref) = poly_gen(num_coeffs)

    return graph_rational(
        den_graph,
        num_graph;
        den_cref = den_cref,
        num_cref = num_cref,
    )
end

function graph_rational(
    den_graph,
    num_graph;
    den_cref = Vector{Tuple{Symbol,Int}}(),
    num_cref = Vector{Tuple{Symbol,Int}}(),
)
    graph = merge_graphs(
        den_graph,
        num_graph,
        prefix1 = "de",
        prefix2 = "nu",
        skip_basic1 = true,
        skip_basic2 = true,
        cref1 = den_cref,
        cref2 = num_cref,
    )
    add_ldiv!(graph, :pade, graph.outputs[2], graph.outputs[1])
    cref = vcat(den_cref, num_cref)
    empty!(graph.outputs)
    compress_graph_redundant!(graph)
    add_output!(graph, :pade)

    return (graph, cref)
end


text:
println("--- :julia: Instantiating project")
using Pkg
Pkg.activate("..")
Pkg.instantiate()
Pkg.activate(".")
Pkg.instantiate()
push!(LOAD_PATH, joinpath(@__DIR__, ".."))
println("+++ :julia: Fixing Doctests")
using Documenter
using LIKWID
DocMeta.setdocmeta!(LIKWID, :DocTestSetup, :(using LIKWID); recursive=true)
doctest(LIKWID, fix=true)