text:
using LinearAlgebra
using Combinatorics, DataStructures

function row_echelon_linsolve(A::Matrix{T}, b::Vector{T}) where {T}
    j = 0
    x = map(1:size(A, 1)) do i
        while j < size(A, 2)
            j += 1
            if isone(A[i, j]) && all(k -> k == i || iszero(A[k, j]), 1:size(A, 1))
                return b[j]
            end
        end
        error("Not in row_echelon_form, cannot find for `$i`th entry.")
    end
    @assert transpose(A) * x == b
    return x
end

function ordered_block_diag(As, d)
    U = block_diag(As, d)
    U === nothing && return nothing
    iU = U'
    @assert iU ≈ inv(U)
    Bs = [iU * A * U for A in As]
    @assert all(Bs) do B
        isblockdim(B, d)
    end
    refs = [B[1:d, 1:d] for B in Bs]
    offset = d
    for offset in d:d:(size(U, 1) - d)
        I = offset .+ (1:d)
        Cs = [B[I, I] for B in Bs]
        σ_ok = nothing
        for σ in permutations(1:d)
            if all(zip(refs, Cs)) do refC
                ref, C = refC
                isapprox(ref[σ, σ], C, rtol=1e-8)
            end
                σ_ok = σ
                break
            end
        end
        if σ_ok === nothing
            error("No permutation can make $refs and $Cs match")
        end
        U[:, I] = U[:, I[σ_ok]]
        offset += d
    end
    ordered_block_check(U, As, d)
    return U
end

function ordered_block_check(U, As, d)
    iU = U'
    @assert iU ≈ inv(U)
    @assert all(As) do A
        is_ordered_blockdim(iU * A * U, d)
    end
end

function block_diag(As, d)
    for A in As
        #T = eigen(A).vectors
        T = schur(A).vectors
        iT = T'
        @assert iT ≈ inv(T)
        n = LinearAlgebra.checksquare(A)
        union_find = IntDisjointSets(n)
        Bs = [iT * A * T for A in As]
        for B in Bs
            merge_sparsity!(union_find, B)
        end
        blocks = Dict{Int,Vector{Int}}()
        for i in 1:n
            r = find_root!(union_find, i)
            blocks[r] = push!(get(blocks, r, Int[]), i)
        end
        if length(blocks) > 1
            U = similar(T)
            offset = 0
            for v in values(blocks)
                @assert iszero(length(v) % d)
                if length(v) == d
                    V = T[:, v]
                else
                    Cs = [B[v, v] for B in Bs]
                    V = block_diag(Cs, d)
                    V === nothing && break
                    V *= transpose(T[:, v])
                end
                U[:, offset .+ eachindex(v)] = V
                offset += length(v)
            end
            if offset == n
                return U
            end
        end
    end
end

function merge_sparsity!(union_find::IntDisjointSets, A, tol=1e-8)
    for I in CartesianIndices(A)
        i, j = I.I
        if abs(A[I]) > tol
            union!(union_find, i, j)
        end
    end
end

function is_ordered_blockdim(A, d, tol=1e-8)
    B = A[1:d, 1:d]
    return isblockdim(A, d, tol) &&
        all(d:d:(size(A, 1) - d)) do offset
            I = offset .+ (1:d)
            isapprox(B, A[I, I], rtol=tol)
        end
end
function isblockdim(A, d, tol=1e-8)
    for I in CartesianIndices(A)
        i, j = I.I
        if abs(i - j) >= d && abs(A[I]) > tol
            return false
        end
    end
    return true
end


text:
function makepangramchecker(alphabet)
    alphabet = Set(uppercase.(alphabet))
    function ispangram(s)
        lengthcheck = length(s) ≥ length(alphabet)
        return lengthcheck && all(c in uppercase(s) for c in alphabet)
    end
    return ispangram
end

const tests = ["Pack my box with five dozen liquor jugs.",
                "The quick brown fox jumps over a lazy dog.",
                "The quick brown fox jumps\u2323over the lazy dog.",
                "The five boxing wizards jump quickly.",
                "This sentence contains A-Z but not the whole alphabet."]

is_english_pangram = makepangramchecker('a':'z')

for s in tests
    println("The sentence \"", s, "\" is ", is_english_pangram(s) ? "" : "not ", "a pangram.")
end

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using ProximalOperators
using Test

@testset "IndPolyhedral" begin

# set dimensions

m, n = 25, 10

# pick random (nonempty) polyhedron

xmin = -ones(n)
xmax = +ones(n)
x0 = min.(xmax, max.(xmin, 10 .* rand(n) .- 5.0))
A = randn(m, n)
u = A*x0 .+ 0.1
l = A*x0 .- 0.1

# pick random point

x = 10 .* randn(n)
p = similar(x)

@testset "valid" for constr in [
    () -> IndPolyhedral(l, A),
    () -> IndPolyhedral(l, A, xmin, xmax),
    () -> IndPolyhedral(A, u),
    () -> IndPolyhedral(A, u, xmin, xmax),
    () -> IndPolyhedral(l, A, u),
    () -> IndPolyhedral(l, A, u, xmin, xmax),
]
    f = constr()
    @test ProximalOperators.is_convex(f) == true
    @test ProximalOperators.is_set(f) == true
    fx = call_test(f, x)
    p, fp = prox_test(f, x)
end

@testset "invalid" for constr in [
    () -> IndPolyhedral(l, A, xmax, xmin),
    () -> IndPolyhedral(A, u, xmax, xmin),
    () -> IndPolyhedral(l, A, u, xmax, xmin),
]
    @test_throws ErrorException constr()
end

end


text:
# Represents a stop time, which is only referenced in the transit network within trips
struct StopTime
    # NB possible optimization: stop is not even needed, as it's in the pattern
    # but that might not actually help, because by iterating over stop times, we are keeping everything memory locality
    stop::Int64
    # NB possible optimization - stop sequence is only needed before the graph is built
    stop_sequence::Int32
    arrival_time::Int32
    departure_time::Int32
end


text:
mutable struct NeuralNetworkClassifier{B,F,O,L} <: MLJModelInterface.Probabilistic
    builder::B
    finaliser::F
    optimiser::O    # mutable struct from Flux/src/optimise/optimisers.jl
    loss::L         # can be called as in `loss(yhat, y)`
    epochs::Int          # number of epochs
    batch_size::Int # size of a batch
    lambda::Float64 # regularization strength
    alpha::Float64  # regularizaton mix (0 for all l2, 1 for all l1)
    optimiser_changes_trigger_retraining::Bool
    acceleration::AbstractResource       # Train on gpu
end

function NeuralNetworkClassifier(; builder::B   = Short()
                                 , finaliser::F = Flux.softmax
                                 , optimiser::O = Flux.Optimise.ADAM()
                                 , loss::L      = Flux.crossentropy
                                 , epochs       = 10
                                 , batch_size   = 1
                                 , lambda       = 0
                                 , alpha        = 0
                                 , optimiser_changes_trigger_retraining = false
                                 , acceleration = CPU1()
                                 ) where {B,F,O,L}

    model = NeuralNetworkClassifier{B,F,O,L}(builder
                                             , finaliser
                                             , optimiser
                                             , loss
                                             , epochs
                                             , batch_size
                                             , lambda
                                             , alpha
                                             , optimiser_changes_trigger_retraining
                                             , acceleration
                                             )

   message = clean!(model)
   isempty(message) || @warn message

    return model
end

function MLJModelInterface.fit(model::NeuralNetworkClassifier,
                               verbosity::Int,
                               X,
                               y)

    data = collate(model, X, y)

    levels = MLJModelInterface.classes(y[1])
    n_output = length(levels)
    n_input = Tables.schema(X).names |> length

    chain = Flux.Chain(build(model.builder, n_input, n_output),
                       model.finaliser)

    optimiser = deepcopy(model.optimiser)

    chain, history = fit!(chain,
                          optimiser,
                          model.loss,
                          model.epochs,
                          model.lambda,
                          model.alpha,
                          verbosity,
                          model.acceleration,
                          data[1],
                          data[2])

    cache = (deepcopy(model), data, history, n_input, n_output, optimiser)
    fitresult = (chain, levels)
    report = (training_losses=history, )

    return fitresult, cache, report
end

function MLJModelInterface.predict(model::NeuralNetworkClassifier,
                                   fitresult,
                                   Xnew_)
    chain, levels = fitresult
    Xnew = MLJModelInterface.matrix(Xnew_)
    probs = vcat([chain(tomat(Xnew[i, :]))' for i in 1:size(Xnew, 1)]...)
    return MLJModelInterface.UnivariateFinite(levels, probs)
end

function MLJModelInterface.update(model::NeuralNetworkClassifier,
                                  verbosity::Int,
                                  old_fitresult,
                                  old_cache,
                                  X,
                                  y)

    old_model, data, old_history, n_input, n_output, optimiser = old_cache
    old_chain, levels = old_fitresult

    optimiser_flag = model.optimiser_changes_trigger_retraining &&
        model.optimiser != old_model.optimiser

    keep_chain = !optimiser_flag && model.epochs >= old_model.epochs &&
        MLJModelInterface.is_same_except(model, old_model, :optimiser, :epochs)

    if keep_chain
        chain = old_chain
        epochs = model.epochs - old_model.epochs
    else
        chain = Flux.Chain(build(model.builder, n_input, n_output),
                           model.finaliser)
        data = collate(model, X, y)
        epochs = model.epochs
    end

    # we only get to keep the optimiser "state" carried over from
    # previous training if we're doing a warm restart and the user has not
    # changed the optimiser hyper-parameter:
    if !keep_chain ||
        !MLJModelInterface._equal_to_depth_one(model.optimiser,
                                              old_model.optimiser)
        optimiser = deepcopy(model.optimiser)
    end

    chain, history = fit!(chain,
                          optimiser,
                          model.loss,
                          epochs,
                          model.lambda,
                          model.alpha,
                          verbosity,
                          model.acceleration,
                          data[1],
                          data[2])
    if keep_chain
        # note: history[1] = old_history[end]
        history = vcat(old_history[1:end-1], history)
    end

    fitresult = (chain, levels)
    cache = (deepcopy(model), data, history, n_input, n_output, optimiser)
    report = (training_losses=history, )

    return fitresult, cache, report

end

MLJModelInterface.fitted_params(::NeuralNetworkClassifier, fitresult) =
    (chain=fitresult[1],)

MLJModelInterface.metadata_model(NeuralNetworkClassifier,
                                 input=Table(Continuous),
                                 target=AbstractVector{<:Finite},
                                 path="MLJFlux.NeuralNetworkClassifier",
                                 descr="A neural network model for making "*
                                 "probabilistic predictions of a "*
                                 "`Multiclass` or `OrderedFactor` target, "*
                                 "given a table of `Continuous` features. ")


text:
__precompile__(true)
module AxisArrayPlots

using AxisArrays
using RecipesBase

@recipe function plot{T}(arr::AxisArray{T, 1})
    xlabel --> axisnames(arr)[1]
    xrotation --> -45
    xs = first(axisvalues(arr))
    xs, arr.data
end

@recipe function plot{T}(arr::AxisArray{T, 2})
    xlbl, ylbl = axisnames(arr)
    xlabel --> xlbl
    ylabel --> ylbl
    xrotation --> -45
    axisvalues(arr)..., arr.data'
end

export Diff
struct Diff
    aa1::AxisArray{T,1} where {T}
    aa2::AxisArray{T,1} where {T}
    normalize::Bool
    label1::String
    label2::String
    label_diff::String
end
function Diff(aa1,aa2;normalize=false,label1="y1",label2="y2",label_diff="diff")
    Diff(aa1,aa2,normalize,label1,label2,label_diff)
end

@recipe function plot(diff::Diff)
    title --> "Difference Plot"

    ft1 = diff.aa1
    ft2 = diff.aa2
    if diff.normalize
        # ft1 = broadcast!(similar(ft1), *, ft1, (1/maximum(ft1)) )
        # ft2 = broadcast!(similar(ft2), *, ft2, (1/maximum(ft2)) )
        ft1 = Base.normalize(ft1, Inf)
        ft2 = Base.normalize(ft2, Inf)
    end
    difference = ft1 - ft2
    @series begin
        label := diff.label1
        ft1
    end
    @series begin
        label := diff.label2
        ft2
    end
    @series begin
        label := diff.label_diff
        difference
    end
end

end


text:
using Test
using ARMA
using HDF5

function compare(m1::ARMAModel, m2::ARMAModel)
    @test m1.p == m2.p
    @test m1.q == m2.q
    @test m1.covarIV ≈ m2.covarIV atol=3e-4 * (2^m1.p)
    @test model_covariance(m1, 50) ≈ model_covariance(m2, 50) atol=3e-4 * (2^m1.p)
end

@testset "hdf5 save/load" begin
    ntests = 5
    for i = 1:ntests
        nbases = rand(2:6)
        bases = rand(nbases)
        ampls = 10*rand(nbases)
        covarIV = [3*sum(abs.(ampls))]
        model = ARMAModel(bases, ampls, covarIV)

        fname1 = tempname()*".hdf5"
        ARMA.hdf5save(fname1, model)
        model1 = ARMA.hdf5load(fname1)
        compare(model, model1)

        fname2 = tempname()*".hdf5"
        h5open(fname2, "w") do f
            g1 = create_group(f, "top")
            g2 = create_group(g1, "mid")
            g3 = create_group(g2, "low")
            ARMA.hdf5save(g3, model)
        end
        h5open(fname2, "r") do f
            model2 = ARMA.hdf5load(f["top/mid/low"])
            compare(model, model2)
        end
    end

    # These specific bases/amplitudes gave me problems before
    bases = [0.4501949250518989 - 0.8928250258588964im, 0, .6922098389057583, 0.9713167531173473, 0.9834786370247353]
    bases[2] = conj(bases[1])
    ampls = [5.7138922321940395 + 0.9683112174170238im, 0, 5.883653745922931, 30.940176242772605, 12.743102625756544]
    ampls[2] = conj(ampls[1])
    covarIV = [195.135]
    model = ARMAModel(bases, ampls, covarIV)

    fname1 = tempname()*".hdf5"
    ARMA.hdf5save(fname1, model)
    model1 = ARMA.hdf5load(fname1)
    compare(model, model1)

    fname2 = tempname()*".hdf5"
    h5open(fname2, "w") do f
        g1 = create_group(f, "top")
        g2 = create_group(g1, "mid")
        g3 = create_group(g2, "low")
        ARMA.hdf5save(g3, model)
    end
    h5open(fname2, "r") do f
        model2 = ARMA.hdf5load(f["top/mid/low"])
        compare(model, model2)
    end
end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((a+b*sin(x))^-1*(A+B*cos(x)), x) == :(B*b^-1*log(a+b*sin(x))+2*A*(a^2+-1*b^2)^-1//2*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*x))))
@test integrate((1+sin(x))^-1*(A+B*cos(x)), x) == :(B*log(1+sin(x))+-1*A*(1+sin(x))^-1*cos(x))
@test integrate((1+-1*sin(x))^-1*(A+B*cos(x)), x) == :(-1*B*log(1+-1*sin(x))+A*(1+-1*sin(x))^-1*cos(x))
@test integrate((a+b*sin(x))^-1*(b+c+cos(x)), x) == :(b^-1*log(a+b*sin(x))+(a^2+-1*b^2)^-1//2*(2b+2c)*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*x))))
@test integrate((a+-1*b*sin(x))^-1*(b+c+cos(x)), x) == :(-1*b^-1*log(a+-1*b*sin(x))+-1*(a^2+-1*b^2)^-1//2*(2b+2c)*arctan((a^2+-1*b^2)^-1//2*(b+-1*a*tan((1/2)*x))))
@test integrate((a+b*sin(x))^-1*(A+B*tan(x)), x) == :(-1*B*(-2b+2a)^-1*log(1+sin(x))+-1*B*(2a+2b)^-1*log(1+-1*sin(x))+2*A*(a^2+-1*b^2)^-1//2*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*x)))+B*a*(a^2+-1*b^2)^-1*log(a+b*sin(x)))
@test integrate((a+b*sin(x))^-1*(A+B*cot(x)), x) == :(B*a^-1*log(sin(x))+-1*B*a^-1*log(a+b*sin(x))+2*A*(a^2+-1*b^2)^-1//2*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*x))))
@test integrate((a+b*sin(x))^-1*(A+B*sec(x)), x) == :(B*(-2b+2a)^-1*log(1+sin(x))+-1*B*(2a+2b)^-1*log(1+-1*sin(x))+2*A*(a^2+-1*b^2)^-1//2*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*x)))+-1*B*b*(a^2+-1*b^2)^-1*log(a+b*sin(x)))
@test integrate((a+b*sin(x))^-1*(A+B*csc(x)), x) == :(-1*B*a^-1*arctanh(cos(x))+a^-1*(a^2+-1*b^2)^-1//2*(-2*B*b+2*A*a)*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*x))))

max_stars_repo_path: test/4 Trig functions/4.1 Sine/4.1.8 (a+b sin)^m (c+d trig)^n.jl


text:
module ImageClassification

using ...NaiveGAflux
using ..AutoFlux: fit
using NaiveGAflux: GlobalPool
using NaiveGAflux: shapetrace, squashshapes, fshape, ndimsout, check_apply
using NaiveGAflux: StatefulGenerationIter
using NaiveNASlib.Advanced, NaiveNASlib.Extend
import Flux
using Flux: Dense, Conv, ConvTranspose, DepthwiseConv, CrossCor, LayerNorm, BatchNorm, InstanceNorm, GroupNorm,
            MaxPool, MeanPool, Dropout, AlphaDropout, GlobalMaxPool, GlobalMeanPool, cpu, gpu
using Flux: Descent, Momentum, Nesterov, ADAM, NADAM, ADAGrad, WeightDecay
import Functors
using Functors: fmap
using Random
import Logging
using Statistics
import IterTools: ncycle

# To store program state for pause/resume
using Serialization

export ImageClassifier, fit, TrainSplitAccuracy, TrainAccuracyVsSize, AccuracyVsSize, TrainIterConfig, BatchedIterConfig, ShuffleIterConfig, GlobalOptimizerMutation, EliteAndSusSelection, EliteAndTournamentSelection

modelname(c::AbstractCandidate) = NaiveGAflux.model(modelname, c)
modelname(g::CompGraph) = split(name(g.inputs[]),'.')[1]

include("strategy.jl")
include("archspace.jl")

"""
    ImageClassifier
    ImageClassifier(popinit, popsize, seed)
    ImageClassifier(;popsize=50, seed=1, newpop=false, mdir=defaultdir("ImageClassifier"); insize, outsize)

Type to make `AutoFlux.fit` train an image classifier using initial population size `popsize` using random seed `seed`.

Load models from `mdir` if directory contains models. If persistence is used (e.g. by providing `cb=persist` to `fit`) candidates will be stored in this directory.

If `newpop` is `true` the process will start with a new population and existing state in the specified directory will be overwritten.
"""
struct ImageClassifier{F}
    popinit::F
    popsize::Int
    seed::Int
end
function ImageClassifier(;popsize=50, seed=1, newpop=false, mdir=defaultdir("ImageClassifier"), insize, outsize)
    popinit = () -> generate_persistent(popsize, newpop, mdir, insize, outsize)
    return ImageClassifier(popinit, popsize, seed)
end

"""
    fit(c::ImageClassifier, x, y; cb, fitnesstrategy, evolutionstrategy, stopcriterion)

Return a population of image classifiers fitted to the given data.

# Arguments
- `c::ImageClassifier`: Type of models to train. See [`ImageClassifier`](@ref).

- `x`: Input data. Must be a 4D array.

- `y`: Output data. Can either be an 1D array in which case it is assumed that `y` is the raw labes (e.g. `["cat", "dog", "cat", ...]`) or a 2D array in which case it is assumed that `y` is one-hot encoded.

- `cb=identity`: Callback function. After training and evaluating each generation but before evolution `cb(population)` will be called where `population` is the array of candidates. Useful for persistence and plotting.

- `fitnesstrategy::AbstractFitnessStrategy=`[`TrainSplitAccuracy()`](@ref): Strategy for fitness from data.

- `evolutionstrategy::AbstractEvolutionStrategy=`[`EliteAndTournamentSelection(popsize=c.popsize)`](@ref): Strategy for evolution.

- `stopcriterion`: Takes the current population and returns true if fitting shall stop. Candidate fitness is available by calling `fitness(c)` where `c` is a member of the population.

"""
function AutoFlux.fit(c::ImageClassifier, x::AbstractArray, y::AbstractArray;
                            cb=identity,
                            fitnesstrategy::AbstractFitnessStrategy=TrainSplitAccuracy(),
                            evolutionstrategy::AbstractEvolutionStrategy=EliteAndTournamentSelection(popsize=c.popsize),
                            stopcriterion = pop -> generation(pop) > 100)
    ndims(x) == 4 || error("Must use 4D data, got $(ndims(x))D data")

    inshape = size(x)[1:2]
    return fit(c, fitnessfun(fitnesstrategy, x, y), evostrategy(evolutionstrategy, inshape), stopcriterion; cb)
end

"""
    fit(c::ImageClassifier, fitnesstrategy::AbstractFitness, evostrategy::AbstractEvolution, stopcriterion; cb)

Return a population of image classifiers fitted to the given data.

Lower level version of `fit` to use when `fit(c::ImageClassifier, x, y)` doesn't cut it.

# Arguments
- `c::ImageClassifier`: Type of models to train. See [`ImageClassifier`](@ref).

- `fitnessstrategy`: An `AbstractFitness` used to compute the fitness metric for a candidate.

- `evostrategy::AbstractEvolution`: Evolution strategy to use. Population `p` will be evolved through `p = evolve(evostrategy, p)`.

- `stopcriterion`: Takes the current population and returns true if fitting shall stop. Candidate fitness is available by calling `fitness(c)` where `c` is a member of the population.

- `cb=identity`: Callback function. After training and evaluating each generation but before evolution `cb(population)` will be called where `population` is the array of candidates. Useful for persistence and plotting.
"""
function AutoFlux.fit(c::ImageClassifier, fitnesstrategy::AbstractFitness, evostrategy::AbstractEvolution, stopcriterion; cb = identity)
    Random.seed!(NaiveGAflux.rng_default, c.seed)
    @info "Start training with baseseed: $(c.seed)"

    population = c.popinit()

    # If experiment was resumed we should start by evolving as population is persisted right before evolution
    population = generation(population) > 1 ? evolve(evostrategy, population) : population

    logfitness = LogFitness(;fitnesstrategy)

    return evolutionloop(population, evostrategy, logfitness, stopcriterion, cb)
end

function evolutionloop(population, evostrategy, fitnesstrategy, stop, cb)
    while true
        @info "Begin generation $(generation(population))"

        fittedpopulation = fitness(population, fitnesstrategy)
        cb(fittedpopulation)
        stop(fittedpopulation) && return fittedpopulation

        population = evolve(evostrategy, fittedpopulation)
    end
end

function generate_persistent(nr, newpop, mdir, insize, outsize, cwrap=identity, archspace = initial_archspace(insize[1:2], outsize))
    if newpop
        rm(mdir, force=true, recursive=true)
    end

    iv(i) = conv2dinputvertex(join(["model", i, ".input"]), insize[3])
    return Population(PersistentArray(mdir, nr, i -> create_model(join(["model", i]), archspace, iv(i), cwrap)))
end
function create_model(name, as, in, cwrap)
    optselect = optmutation(1.0)
    opt = optselect(Descent(rand() * 0.099 + 0.01))
    # Always cache even if not strictly needed for all fitnessfunctions because consequence is so bad if one forgets
    # it when needed. Users who know what they are doing can unwrap if caching is not wanted.
    cwrap(CandidateOptModel(opt, CompGraph(in, as(name, in))))
end

end


text:
@testset "smarts.smarts" begin

@testset "chain" begin
    nullmol = parse(SMARTS, "")
    @test nodecount(nullmol) == 0
    @test edgecount(nullmol) == 0

    aliphatic = parse(SMARTS, "C")
    @test nodeattr(aliphatic, 1).query == (
        :and => (:atomsymbol => :C, :isaromatic => false)
    )
    @test edgecount(aliphatic) == 0

    carbonyl = parse(SMARTS, "[CX3]=[OX1]")
    @test nodeattr(carbonyl, 1).query == (:and => (
        :and => (:atomsymbol => :C, :isaromatic => false),
        :connectivity => 3
    ))
    @test nodeattr(carbonyl, 2).query == (:and => (
        :and => (:atomsymbol => :O, :isaromatic => false),
        :connectivity => 1
    ))

    ether = parse(SMARTS, "[#6][OD2][#6]")
    @test nodeattr(ether, 1).query == (:atomsymbol => :C)
    @test nodeattr(ether, 2).query == (:and => (
        :and => (:atomsymbol => :O, :isaromatic => false),
        :nodedegree => 2
    ))

    notH = parse(SMARTS, "[!#1]")
    @test nodeattr(notH, 1).query == (:not => (:atomsymbol => :H))
end

end # smarts.SMARTS


text:
module Tools

using ProgressMeter, Base.Threads

"""
`tfor(f, rg; verbose=true, msg="tfor", costs = ones(Int, length(rg)))`

Multi-threaded for loop. At each iteration the function f(n) is executed,
where n loops over `rg`.
"""
function tfor(f, rg; verbose=true, msg="tfor", costs = ones(Int, length(rg)),
                     maxnthreads=nthreads())
   p = Progress(sum(costs))
   p_ctr = 0
   t0 = time_ns()
   nthr = max(1, min(nthreads(), maxnthreads))
   if nthr == 1
      verbose && println("$msg in serial")
      dt = verbose ? 1.0 : Inf
      for (n, c) in zip(rg, costs)
         f(n)
         if verbose
            p_ctr += c
            ProgressMeter.update!(p, p_ctr)
         end
      end
   else
      if verbose
         @info("$msg with $(nthreads()) threads")
         ProgressMeter.update!(p, 0)  # not sure this is useful/needed
         p_lock = SpinLock()
      end
      # sort the tasks by cost: do the expensive ones first
      Isort = sortperm(costs, rev=true)
      # remember what the last job was
      last_job = 0
      last_job_lock = SpinLock()
      # start a simple loop over nthreads() just to split this
      # into parallel "tasks"
      rgidx = Vector{Int}(undef, nthreads())
      @threads for i = 1:nthreads()
         while true
            # acquire a new job index
            tid = threadid()
            lock(last_job_lock)
            last_job += 1
            if last_job > length(Isort)
               unlock(last_job_lock)
               break # break the while loop and then wait for the
                     # other threads to finish
            end
            rgidx[tid] = Isort[last_job]
            unlock(last_job_lock)
            # do the job
            f(rg[rgidx[tid]])
            # submit progress
            if verbose
               lock(p_lock)
               p_ctr += costs[rgidx[tid]]
               ProgressMeter.update!(p, p_ctr)
               unlock(p_lock)
            end
         end
      end
   end
   t1 = time_ns()
   verbose && @info("Elapsed: $(round((t1-t0)*1e-9, digits=1))s")
   return nothing
end

function analyse_include_exclude(set, include, exclude)
   if include != nothing && exclude != nothing
      error("only one of `include`, `exclude` may be different from `nothing`")
   end
   if include != nothing
      if !issubset(include, set)
         error("`include` can only contain elements of `set`")
      end
      # do nothing - just keep `include` as is to return
   elseif exclude != nothing
      if !issubset(exclude, set)
         error("`exclude` can only contain config types that are in `set`")
      end
      include = setdiff(set, exclude)
   else
      # both are nothing => keep all config_types
      include = set
   end
   return include
end

macro def(name, definition)
    return quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end
end

end


text:
using ComoniconTestUtils
using ComoniconZSHCompletion: emit_zshcompletion
using Test
using Random

Random.seed!(42)

@testset "test completion" for _ in 1:5
    cmd = rand_command()
    script = emit_zshcompletion(cmd)
    @test occursin("#compdef _$(cmd.root.name) $(cmd.root.name)", script)
    @test occursin("function _$(cmd.root.name)() {", script)
end


text:
# meos.jl
################################################################################
# Note: The following is based on the nomenclature in the paper
# "Revised Release on the IAPWS Formulation 1995 for the Thermodynamic
# Properties of Ordinary Water Substance for General and Scientific Use"
# by IAPWS, the International Association for the Properties of Water and Steam.
#
# It is more or less as follows:
#
# ϕ : Dimensionless Helmholtz energy as a fucntion of density and temperature
# δ : Reduced density (= ρ/ρc)
# τ : Inverse reduced temperature (= Tc/T)
# n, d, t, α, β, γ, ϵ, A, B, C, D : these are parameters used in fitting the
# equation of state to data. As far as we're concerned, they are constants with
# no real physical significance.
#
# Note that although the code was written by consulting the IAPWS paper, it is
# applicable to any multiparameter equation of state with the same form (notably
# Span and Wagner's equation of state for carbon dioxide, which is also provided
# by this package; see R. Span and W. Wagner, J. Phys. Chem. Ref. Data Vol 25,
# No. 6, 1996 Titled "A New Equation of State for Carbon Dioxide Covering the
# Fluid Region from the Triple-Point Temperature to 1100 K at Pressures up to
# 800 MPa")
################################################################################

include("substances.jl")

__precompile__()
module meos
#import ForwardDiff
import Calculus
import pr
#using Optim
using Roots.fzero
struct multiparameos{T<:AbstractFloat}
    Tc::T
    ρc::T
    pc::T
    ω::T
    R::T

    #parameters for ideal gas portion
    n₀::Vector{T}
    γ₀::Vector{T}

    # parameters for residual term 1
    n₁::Vector{T}
    d₁::Vector{T}
    t₁::Vector{T}

    # parameters for residual term 2
    n₂::Vector{T}
    d₂::Vector{T}
    t₂::Vector{T}
    c₂::Vector{T}

    # parameters for residual term 3
    n₃::Vector{T}
    d₃::Vector{T}
    t₃::Vector{T}
    α₃::Vector{T}
    β₃::Vector{T}
    γ₃::Vector{T}
    ϵ₃::Vector{T}

    # parameters for residual term 4
    n₄::Vector{T}
    a₄::Vector{T}
    b₄::Vector{T}
    β₄::Vector{T}
    A₄::Vector{T}
    B₄::Vector{T}
    C₄::Vector{T}
    D₄::Vector{T}

    # Functions needed to get density guess values
    psat_guess::Function
    PR_::pr.pengRobinson

end

function Σ(expr)
    return sum(eval(@. expr))
end

function ϕ₀(eos::multiparameos, δ, τ)
    first = log(δ) + eos.n₀[1] + eos.n₀[2]*τ + eos.n₀[3]*log(τ)
    end_ = sum(eos.n₀[4:8].*log.(1-exp.(-eos.γ₀.*τ)))
    return first + end_
end

function res_1(eos::multiparameos, δ, τ)
    return sum(eos.n₁ .* (δ .^ eos.d₁) .* (τ .^ eos.t₁))
end

function res_2(eos::multiparameos, δ, τ)
    return sum(eos.n₂ .* (δ .^ eos.d₂) .* (τ .^ eos.t₂) .* exp.(-δ .^ eos.c₂))
end

function res_3(eos::multiparameos, δ, τ)
    return sum(eos.n₃ .* (δ .^ eos.d₃) .* (τ .^ eos.t₃) .* exp.(-eos.α₃ .* (δ - eos.ϵ₃).^2 - eos.β₃ .* (τ - eos.γ₃).^2))
end

function θ(eos::multiparameos, δ, τ)
    return sum((1-τ) + eos.A₄.*((δ - 1).^2).^(1 ./(2 .* eos.β₄)))
end

function Δ(eos::multiparameos, δ, τ)
    return sum(θ(eos, δ, τ).^2 + eos.B₄ .* ((δ - 1).^2) .* eos.a₄)
end

function Ψ(eos::multiparameos, δ, τ)
    return sum(exp(-eos.C₄ .* (δ - 1)^2 - eos.D₄ .* (τ - 1)^2))
end

function res_4(eos::multiparameos, δ, τ)
    return sum(eos.n₄ .* Δ(eos,δ,τ).^eos.b₄ .* δ * Ψ(eos,δ,τ))
end

function ϕr(eos::multiparameos, δ, τ)
    return res_1(eos, δ, τ)+res_2(eos, δ, τ)+res_3(eos, δ, τ)+res_4(eos, δ, τ)
end

function ϕ₀_δ(eos::multiparameos, δ, τ)
    f(δ_) = ϕ₀(eos, δ_, τ)
    return Calculus.derivative(f, δ)
end
function ϕ₀_δδ(eos::multiparameos, δ, τ)
    f(δ_) = ϕ₀_δ(eos, δ_, τ)
    return Calculus.derivative(f, δ)
end
function ϕ₀_τ(eos::multiparameos, δ, τ)
    f(τ_) = ϕ₀(eos, δ, τ_)
    return Calculus.derivative(f, τ)
end
function ϕ₀_ττ(eos::multiparameos, δ, τ)
    f(τ_) = ϕ₀_τ(eos, δ, τ_)
    return Calculus.derivative(f, τ)
end
function ϕ₀_δτ(eos::multiparameos, δ, τ)
    return 0.0
end

function ϕr_δ(eos::multiparameos, δ, τ)
    f(δ_) = ϕr(eos, δ_, τ)
    return Calculus.derivative(f, δ)
end
function ϕr_δδ(eos::multiparameos, δ, τ)
    f(δ_) = ϕr_δ(eos, δ_, τ)
    return Calculus.derivative(f, δ)
end
function ϕr_τ(eos::multiparameos, δ, τ)
    f(τ_) = ϕr(eos, δ, τ_)
    return Calculus.derivative(f, τ)
end
function ϕr_ττ(eos::multiparameos, δ, τ)
    f(τ_) = ϕr_τ(eos, δ, τ_)
    return Calculus.derivative(f, τ)
end
function ϕr_δτ(eos::multiparameos, δ, τ)
    f(τ_) = ϕr_δ(eos, δ, τ_)
    return Calculus.derivative(f, τ)
end

function p(eos::multiparameos, ρ, T)
    δ = ρ/eos.ρc
    τ = eos.Tc/T
    p_nondim = 1 + δ*ϕr_δ(eos, δ, τ)
    return p_nondim * ρ*eos.R*T
end

"""
This function is arguably the most critical piece of the whole equation of
state. It allows all other functions to be expressed in terms of pressure and
temperature rather than density and temperature.
"""
function ρ_from_p(eos::multiparameos, p_, T)
  pr_ = p_/eos.pc
  Tr = T/eos.Tc
  if abs(pr_-1) < 0.05 || abs(Tr-1) < 0.05
    # Interpolated guess value near critical point
    ρg = eos.ρc + (2*(p_-eos.pc)/eos.pc - (T-eos.Tc)/eos.Tc)
  elseif T > eos.Tc
    # Guess value for supercritical fluid
    ρg = eos.ρc
  elseif p_ > psat(eos, T)
    # Guess value for liquid region
    ρg = 1/pr.vl(eos.PR_,p_, T)
  else
    # Guess value for vapor region
    ρg = 1/pr.vv(eos.PR_,p_, T)
  end
  return ρ_(eos, p_, T, ρg, 2)
end

function ρ_(eos::multiparameos, p_, T, ρg, order_=0)
    function solveit(ρ_guess)
        return abs(p_ - p(eos, ρ_guess, T))
    end
    ans = 0
    try
      a = fzero(solveit, ρg, order=order_)
      ans = a
    catch ex
      lastval = float(ex.reason[18:end])
      #@printf "\nLast value: %f" lastval
      pcalc = p(eos, lastval,T)
      #@printf "\np calculated at %f with last value of ρ" pcalc
      #@printf "\nError: %f" abs(p_ - pcalc)
      ans = lastval
    end
    #print(ans)
    return ans
end

function psat(eos::multiparameos, T)
  if T > eos.Tc
    throw(DomainError("Can't have a vapor pressure above the critical point!"))
  end

  guess = eos.psat_guess(T)
  function solveit(pg)
    ρl_g = 1/pr.vl(eos.PR_,pg, T)
    ρv_g = 1/pr.vv(eos.PR_,pg, T)
    Gl = G_(eos, ρl_g, T)
    Gv = G_(eos, ρv_g, T)
    return Gl-Gv
  end

  return fzero(solveit, guess)
end

function u(eos::multiparameos, p_, T)
    ρ = ρ_from_p(eos, p_, T)
    δ = ρ/eos.ρc
    τ = eos.Tc/T
    u_nondim = τ * (ϕ₀_τ(eos, δ, τ) + ϕr_τ(eos, δ, τ))
    return u_nondim * T*eos.R
end

function s(eos::multiparameos, p_, T)
    ρ = ρ_from_p(eos, p_, T)
    δ = ρ/eos.ρc
    τ = eos.Tc/T
    s_nondim = τ * (ϕ₀_τ(eos, δ, τ) + ϕr_τ(eos, δ, τ)) - (ϕ₀(eos, δ, τ) + ϕr(eos, δ, τ))
    return s_nondim * eos.R
end

function H(eos::multiparameos, p_, T)
    ρ = ρ_from_p(eos, p_, T)
    δ = ρ/eos.ρc
    τ = eos.Tc/T
    H_nondim = 1 + τ * (ϕ₀_τ(eos, δ, τ) + ϕr_τ(eos, δ, τ)) + δ*ϕr_δ(eos, δ, τ)
    return H_nondim * T*eos.R
end

function G(eos::multiparameos, p_, T)
    ρ = ρ_from_p(eos, p_, T)
    δ = ρ/eos.ρc
    τ = eos.Tc/T
    G_nondim = 1 + (ϕ₀(eos, δ, τ) + ϕr(eos, δ, τ)) + δ*ϕr_δ(eos, δ, τ)
    return G_nondim * T*eos.R
end

function G_(eos::multiparameos, ρ, T)
    δ = ρ/eos.ρc
    τ = eos.Tc/T
    G_nondim = 1 + (ϕ₀(eos, δ, τ) + ϕr(eos, δ, τ)) + δ*ϕr_δ(eos, δ, τ)
    return G_nondim * T*eos.R
end

function Cp(eos::multiparameos, p_, T)
    ρ = ρ_from_p(eos, p_, T)
    δ = ρ/eos.ρc
    τ = eos.Tc/T
    part1 = -τ^2 * (ϕ₀_ττ(eos, δ, τ) + ϕr_ττ(eos, δ, τ))
    p2_top = (1 + δ*ϕr_δ(eos, δ, τ) - δ*τ*ϕr_δτ(eos, δ, τ))^2
    p2_bot = 1 + 2δ*ϕr_δ(eos, δ, τ) + δ^2*ϕr_δδ(eos, δ, τ)
    Cp_nondim = part1 + p2_top/p2_bot
    return Cp_nondim * eos.R
end

end # ends module meos

module IAPWS
using meos
using pr.pengRobinson
function gen_iapws()
    h2o_n₀ = [-8.3204464837497, 6.6832105275932, 3.00632, 0.012436, 0.97315, 1.27950,
              0.96956, 0.24873]
    h2o_γ₀ = [1.28728967, 3.53734222, 7.74073708, 9.24437796, 27.5075105]

    h2o_n₁ = [0.12533547935523e-1,
              0.78957634722828e1,
              -0.87803203303561e1,
              0.31802509345418,
              -0.26145533859358,
              -0.78199751687981e-2,
              0.88089493102134e-2]
    h2o_d₁ = [1, 1, 1, 2, 2, 3, 4]
    h2o_t₁ = [-0.5, 0.875, 1, 0.5, 0.75, 0.375, 1]

    h2o_n₂ = [-0.66856572307965,
              0.20433810950965,
              -0.66212605039687e-4,
              -0.19232721156002,
              -0.25709043003438,
              0.16074868486251,
              -0.40092828925807e-1,
              0.39343422603254e-6,
              -0.75941377088144e-5,
              0.56250979351888e-3,
              -0.15608652257135e-4,
              0.11537996422951e-8,
              0.36582165144204e-6,
              -0.13251180074668e-11,
              -0.62639586912454e-9,
              -0.10793600908932,
              0.17611491008752e-1,
              0.22132295167546,
              -0.40247669763528,
              0.58083399985759,
              0.49969146990806e-2,
              -0.31358700712549e-1,
              -0.74315929710341,
              0.47807329915480,
              0.20527940895948e-1,
              -0.13636435110343,
              0.14180634400617e-1,
              0.83326504880713e-2,
              -0.29052336009585e-1,
              0.38615085574206e-1,
              -0.20393486513704e-1,
              -0.16554050063734e-2,
              0.19955571979541e-2,
              0.15870308324157e-3,
              -0.16388568342530e-4,
              0.43613615723811e-1,
              0.34994005463765e-1,
              -0.76788197844621e-1,
              0.22446277332006e-1,
              -0.62689710414685e-4,
              -0.55711118565645e-9,
              -0.19905718354408,
              0.31777497330738,
              -0.11841182425981]
    h2o_c₂ = [1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 2,2,2,2,2,2,2,2, 2,2,2,2,2,2,2,2,2,
              2,2,2, 3,3,3,3, 4,6,6,6,6]
    h2o_d₂ = [1,1,1,2,2,3,4,4,5,7,9,10,11,13,15,1,2,2,2,3,4,4,4,5,6,6,7,9,9,9,9,9,
              10,10,12,3,4,4,5,14,3,6,6,6]
    h2o_t₂ = [4,6,12,1,5,4,2,13,9,3,4,11,4,13,1,7,1,9,10,10,3,7,10,10,6,10,10,1,2,3,
              4,8,6,9,8,16,22,23,23,10,50,44,46,50]

    h2o_d₃ = [3, 3, 3]
    h2o_t₃ = [0,1,4]
    h2o_n₃ = [-0.31306260323435e2, 0.31546140237781e2, -0.25213154341695e4]
    h2o_α₃ = [20,20,20]
    h2o_β₃ = [150,150,150]
    h2o_γ₃ = [1.21,1.21,1.25]
    h2o_ϵ₃ = [1,1,1]

    h2o_a₄ = [3.5,3.5]
    h2o_b₄ = [0.85,0.95]
    h2o_β₄ = [0.3,0.3]
    h2o_n₄ = [-0.14874640856724, 0.31806110878444]
    h2o_A₄ = [0.32,0.32]
    h2o_B₄ = [0.2,0.2]
    h2o_C₄ = [28,32]
    h2o_D₄ = [700,800]

    Tc = 647.096
    pc = 22.064e6
    ρc = 322
    ω = 0.344
    R = 0.46151805*1000
    function psat_dippr(T)
      return exp(73.649 - 7258.2/T - 7.3037*log(T) + 4.1653e-6T^2)
    end
    pr_thing = pengRobinson(pc, Tc, ω)

    return meos.multiparameos(Tc,ρc,pc,R,ω,h2o_n₀,h2o_γ₀,h2o_n₁,h2o_d₁,h2o_t₁,
      h2o_n₂,h2o_d₂,h2o_t₂,h2o_c₂,h2o_n₃,h2o_d₃,h2o_t₃,h2o_α₃,h2o_β₃,h2o_γ₃,
      h2o_ϵ₃,h2o_n₄,h2o_a₄,h2o_b₄,h2o_β₄,h2o_A₄,h2o_B₄,h2o_C₄,h2o_D₄, psat_dippr, pr_thing)
end

end # ends module IAPWS

module co2
import meos
using pr.pengRobinson
using subst.substance
function gen_co2()
    co2_n₀= [ 8.37304456, -3.70454304,  2.5       ,  1.99427042,  0.62105248,
        0.41195293,  1.04028922,  0.08327678]

    co2_γ₀ = [3.15163,   6.1119 ,   6.77708,
            11.32384,  27.08792]

    co2_n₁ = [ 0.38856823,  2.93854759, -5.58671885, -0.767532  ,  0.31729004,
            0.54803316,  0.12279411]
    co2_d₁ = [ 1.,  1.,  1.,  1.,  2.,  2.,  3.]
    co2_t₁ = [ 0.  ,  0.75,  1.  ,  2.  ,  0.75,  2.  ,  0.75]
    co2_n₂ = [  2.16589615e+00,   1.58417351e+00,  -2.31327054e-01,
             5.81169164e-02,  -5.53691372e-01,   4.89466159e-01,
            -2.42757398e-02,   6.24947905e-02,  -1.21758602e-01,
            -3.70556853e-01,  -1.67758797e-02,  -1.19607366e-01,
            -4.56193625e-02,   3.56127893e-02,  -7.44277271e-03,
            -1.73957049e-03,  -2.18101213e-02,   2.43321666e-02,
            -3.47401334e-02,   1.43387158e-01,  -1.34919691e-01,
            -2.31512251e-02,   1.23631255e-02,   2.10583220e-03,
            -3.39585190e-04,   5.59936518e-03,  -3.03351181e-04]

    co2_d₂ = [  1.,   2.,   4.,   5.,   5.,   5.,   6.,   6.,   6.,   1.,   1.,
             4.,   4.,   4.,   7.,   8.,   2.,   3.,   3.,   5.,   5.,   6.,
             7.,   8.,  10.,   4.,   8.]

    co2_c₂ = [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  2.,  2.,  2.,  2.,
            2.,  2.,  2.,  3.,  3.,  3.,  4.,  4.,  4.,  4.,  4.,  4.,  5.,  6.]

    co2_t₂ = [  1.5,   1.5,   2.5,   0. ,   1.5,   2. ,   0. ,   1. ,   2. ,
             3. ,   6. ,   3. ,   6. ,   8. ,   6. ,   0. ,   7. ,  12. ,
            16. ,  22. ,  24. ,  16. ,  24. ,   8. ,   2. ,  28. ,  14. ]

    co2_n₃ = [  -213.65488688,  26641.56914927, -24027.21220456,  -283.41603424,
              212.472844  ]
    co2_n₄ = [-0.66642277,  0.72608632,  0.05506867]
    co2_d₃ = [2.,2.,2.,3.,3.]
    co2_t₃ = [1.,0.,1.,3.,3.]
    co2_α₃ = [25.,25.,25.,15.,20.]
    co2_β₃ = [325.,300.,300.,275.,275.]
    co2_γ₃ = [1.16,1.19,1.19,1.25,1.22]
    co2_ϵ₃ = [1.,1.,1.,1.,1.]

    co2_a₄ = [3.5,3.5,3.0]
    co2_b₄ = [0.875,0.925,0.875]
    co2_β₄ = [0.3,0.3,0.3]
    co2_A₄ = [0.7,0.7,0.7]
    co2_B₄ = [0.3,0.3,1.0]
    co2_C₄ = [10.0,10.0,12.5]
    co2_D₄ = [275.,275.,275.]

    Tc = 304.1282
    ρc = 467.6
    pc = 7.38e6
    R = 188.9241
    ω = 0.228
    function psat_dippr(T)
      return exp(47.0169 - 2839/T - 3.86388*log(T) + 2.81115e-16T^6)
    end
    pr_thing = pengRobinson(pc, Tc, ω)
    return meos.multiparameos(Tc,ρc,pc,ω,R,co2_n₀,co2_γ₀,co2_n₁,co2_d₁,co2_t₁,
      co2_n₂,co2_d₂,co2_t₂,co2_c₂,co2_n₃,co2_d₃,co2_t₃,co2_α₃,co2_β₃,co2_γ₃,
      co2_ϵ₃,co2_n₄,co2_a₄,co2_b₄,co2_β₄,co2_A₄,co2_B₄,co2_C₄,co2_D₄,psat_dippr,pr_thing)
end

function co2_viscosity(eos::meos.multiparameos, p, T)

  i1 = [0,1,2,3,4]
  a = [0.235156,
       -0.491266,
       5.211155e-2,
       5.347906e-2,
       -1.537102e-2]
  i2 = [1,2,6,8,8]
  j = [0,0,3,0,1]
  d = [0.407119e-2,
       0.7198037e-4,
       0.211697e-16,
       0.2971072e-22,
       -0.1627888e-22]
  ρ = meos.ρ_from_p(eos, p, T)
  Tr = T / 251.196
  rhor = ρ/eos.ρc
  top = 1.00697 * T ^ 0.5
  C_star = 0
  for i in i1
    C_star += a[i+1] * log(Tr) ^ i
  end
  bottom = exp(C_star)
  eta0 = top / bottom
  etar = 0
  for (num, i) in enumerate(i2)
    etar += d[num] * (rhor ^ i) / (Tr ^ j[num])
  end
  eta = eta0 + etar
  return eta * 1e-6
end

function co2_thermalConductivity(eos::meos.multiparameos, p, T)
  ik1 = [1,2,3,4,5,6,7,8,9,10]
  g = [0.0,0.0,1.5,0.0,1.0,1.5,1.5,1.5,3.5,5.5]
  h = [1.0,5.0,1.0,1.0,2.0,0.0,5.0,9.0,0.0,0.0]
  n = [7.69857587,
       0.159885811,
       1.56918621,
       -6.73400790,
       16.3890156,
       3.69415242,
       22.3205514,
       66.1420950,
       -0.171779133,
       0.00433043347,
       0.775547504]
  a = [3.0,
       6.70697,
       0.94604,
       0.30,
       0.30,
       0.39751,
       0.33791,
       0.77963,
       0.79857,
       0.9,
       0.02,
       0.2]
  ρ = meos.ρ_from_p(eos, p, T)
  rhor = ρ / eos.ρc
  Tr = T / eos.Tc
  thermal_cond = 0
  for i in ik1
    if i != 0
      index = i
      term = n[index] * Tr ^ g[index] * rhor ^ h[index]
      if i >= 4
        term *= exp(-5 * rhor ^ 2)
      end
      thermal_cond += term
    end
  end
  alpha = 1 - a[10] * acosh(1 + a[11] * ((1 - Tr) ^ 2) ^ a[12])
  crit_correction = rhor * exp((-rhor ^ a[1]) / a[1] - (a[2] * (Tr - 1)) ^ 2 - (a[3] * (rhor - 1)) ^ 2) /
                      ((((1 - 1 / Tr) + a[4] * ((rhor - 1) ^ 2) ^ (1 / (2 * a[5]))) ^ 2) ^ a[6]
                       + ((a[7] * (rhor - alpha)) ^ 2) ^ a[8]) ^ a[9]
  nc = n[9]
  thermal_cond += nc * crit_correction
  kc = 4.81314e-3
  return thermal_cond * kc  # Should return units of W / m *K
end

function init_co2()
  general_eos = gen_co2()
  function co2_density(p, T)
    return meos.ρ_from_p(general_eos, p, T)
  end
  function co2_enthalpy(p, T)
    return meos.H(general_eos, p, T)
  end
  function co2_Cp(p, T)
    return meos.Cp(general_eos, p, T)
  end
  function co2_μ(p, T)
    return co2_viscosity(general_eos, p, T)
  end
  function co2_k(p, T)
    return co2_thermalConductivity(general_eos, p, T)
  end
  return substance("carbon dioxide", [co2_density, co2_enthalpy, co2_Cp,
                                      co2_μ, co2_k])
end

end #ends module co2


text:
"""

Functions that wrap the coordinates

They modify the coordinates of the input vector.

```
function wrap!(x::AbstractVector{T}, sides::T, center::T) where T <: AbstractVector
```

Wrap to a given center of coordinates

"""
function wrap!(x::AbstractVector{T},
               sides::T,
               center::T) where T <: AbstractVector
  for i in eachindex(x)
    x[i] = wrapone(x[i],sides,center)
  end
  nothing
end

@inline function wrapone(x::T, sides::T, center::T) where T <: AbstractVector
  s = @. (x-center)%sides
  s = @. wrapx(s,sides) + center
  return s
end

@inline function wrapx(x,s)
  if x > s/2
    x = x - s
  elseif x < -s/2
    x = x + s
  end
  x
end

"""

```
wrap!(x::AbstractVector{T}, sides::T) where T <: AbstractVector
```

Wrap to origin.

"""
function wrap!(x::AbstractVector{T}, sides::T) where T <: AbstractVector
  for i in eachindex(x)
    x[i] = wrapone(x[i],sides)
  end
  nothing
end

@inline function wrapone(x::T, sides::T) where T <: AbstractVector
  s = @. x%sides
  s = @. wrapx(s,sides)
  return s
end


text:
using Test
using TimeWarp

@testset "Sequences" begin

    # 10-dimensional time series
    X = randn(10,3)
    Sx = Sequence(X)

    @test size(Sx) == (3,)
    @test length(Sx) == 3

    for i = 1:length(Sx)
        @test Sx[i] == X[:,i]
    end

    # (10x9)-dimensional time series
    Y = randn(10,9,3)
    Sy = Sequence(Y)

    @test size(Sy) == (3,)
    @test length(Sy) == 3

    for i = 1:length(Sy)
        @test Sy[i] == Y[:,:,i]
    end

end

@testset "Basic Dynamic Time Warping" begin
    a=[1,1,1,2,4,6,5,5,5,4,4,3,1,1,1]
    b=[1,1,2,4,6,6,6,5,4,4,4,3,3,3,1]
    cost, match1, match2 = dtw(a,b)
    @test cost==0
    @test match1==[1,2,3,4,5,6,6,6,7,8,9,10,10,11,12,12,12,13,14,15]
    @test match2==[1,1,2,3,4,5,6,7,8,8,8, 9,10,11,12,13,14,15,15,15]
    @test evaluate(DTWDistance(),a,b) == cost

    a[end] += 2
    cost, match1, match2 = dtw(a,b)
    @test cost==4
    @test evaluate(DTWDistance(),a,b) == cost

    a=collect(1:10)
    b=a .+ 1
    cost, match1, match2 = dtw(a,b)
    @test cost==2
    @test evaluate(DTWDistance(),a,b) == cost

    a=zeros(Int,6)
    b=1 .+ a
    cost, match1, match2 = dtw(a,b)
    @test cost==length(a)
    @test evaluate(DTWDistance(),a,b) == cost

    # Verify that a tie prefers diagonal moves
    a=[1,1,1]
    b=[1,1,1]
    cost, pa, pb = dtw(a,b)
    @test cost==0
    @test pa==[1,2,3]
    @test pb==[1,2,3]
    @test evaluate(DTWDistance(),a,b) == cost

    # Verify that trackback ends properly if it reaches an edge before reaching [1,1]
    # Also check that trackback prefers diagonal moves
    a=[0,1,1,1]
    b=[0,0,1,1]
    cost, pa, pb = dtw(a,b)
    @test cost==0
    @test pa==[1,1,2,3,4]
    @test pb==[1,2,3,3,4]
    @test evaluate(DTWDistance(),a,b) == cost

    # test the distance api with different distances
    a,b = randn(10),randn(10)
    cost, = dtw(a,b,Euclidean())
    @test evaluate(DTWDistance(Euclidean()),a,b) == cost
    cost, = dtw(a,b,Cityblock())
    @test evaluate(DTWDistance(Cityblock()),a,b) == cost
    cost, = dtw(a,b,Chebyshev())
    @test evaluate(DTWDistance(Chebyshev()),a,b) == cost
end

@testset "DTW with windows" begin
    # Verify that a tie prefers diagonal moves
    a=[1,1,1]
    b=[1,1,1]
    cost, pa, pb = dtw(a,b,[1,1,1],[3,3,3])
    @test cost==0
    @test pa==[1,2,3]
    @test pb==[1,2,3]

    # Verify that trackback ends properly if it reaches an edge before reaching [1,1]
    # Also check that trackback prefers diagonal moves
    a=[0,1,1,1]
    b=[0,0,1,1]
    cost, pa, pb = dtw(a,b,[1,1,1,1],[4,4,4,4])
    @test cost==0
    @test pa==[1,1,2,3,4]
    @test pb==[1,2,3,3,4]

    # First do the windowed test w/o windows
    a=[0,1,2,3,4,4,4,4]
    b=[0,0,1,2,2,2,3,4]
    best_pa = [1,1,2,3,3,3,4,5,6,7,8]
    best_pb = [1,2,3,4,5,6,7,8,8,8,8]
    cost, pa, pb = dtw(a,b)
    @test cost == 0
    @test pa == best_pa
    @test pb == best_pb

    # Wide window, not touching optimal path
    rmin = [1,1,1,2,3,4,5,6]
    rmax = [4,6,7,8,8,8,8,8]
    cost, pa, pb = dtw(a,b,rmin,rmax)
    @test cost == 0
    @test pa == best_pa
    @test pb == best_pb

    # Bottom of window is optimal path
    rmin = [1,3,4,7,8,8,8,8]
    rmax = [4,6,7,8,8,8,8,8]
    cost, pa, pb = dtw(a,b,rmin,rmax)
    @test cost == 0
    @test pa == best_pa
    @test pb == best_pb

    # Top of window is optimal path
    rmin = [1,1,1,2,3,4,5,6]
    rmax = [2,3,6,7,8,8,8,8]
    cost, pa, pb = dtw(a,b,rmin,rmax)
    @test cost == 0
    @test pa == best_pa
    @test pb == best_pb

    # Top and bottom of window are optimal path
    rmin = [1,3,4,7,8,8,8,8]
    rmax = [2,3,6,7,8,8,8,8]
    cost, pa, pb = dtw(a,b,rmin,rmax)
    @test cost == 0
    @test pa == best_pa
    @test pb == best_pb

    # Now top of window cuts into optimal path
    rmin = [1,1,1,2,3,4,5,6]
    rmax = [4,4,5,6,7,8,8,8]
    cost, pa, pb = dtw(a,b,rmin,rmax)
    @test cost == 2
    @test pa == [1,1,2,3,3,4,5,6,7,8]
    @test pb == [1,2,3,4,5,6,7,8,8,8]

    # Compare windowed and regular on non-square geometry:
    seq1=collect(1:16)
    seq2=seq1[1:2:end]
    cost,pa,pb = dtw(seq1,seq2)
    n1,n2 = length(seq1), length(seq2)
    cost2,qa,qb = dtw(seq1, seq2, fill(1,n1), fill(n2,n1))
    @test cost==cost2
    @test pa==qa
    @test pb==qb

    seq1=[1.2051,1.7890,2.1314,2.2423,2.2101,2.0234,1.6898,1.3422,1.1530,1.2197,1.5435,2.07147,2.7359,3.469,4.2033,4.8704,5.4056,5.75243,5.8693,5.7362,5.35974,4.7761,4.0482,3.2548,2.4764]
    seq2=[1.4927,2.1859,2.1182,1.5186, 1.1858,1.8034,3.0968,4.5317,5.5763,5.8038,5.0724,3.6576,2.4823]

    cost,pa,pb = dtw(seq1,seq2)
    n1,n2 = length(seq1), length(seq2)
    cost2,qa,qb = dtw(seq1, seq2, fill(1,n1), fill(n2,n1))
    @test cost==cost2
    @test pa==qa
    @test pb==qb
end

@testset "FastDTW compression" begin
    compress2 = TimeWarp.compress2
    s=collect(0:2:98)
    s1 = compress2(s)
    s2 = compress2(s1)
    @test s1==float(collect(1:4:97))
    @test s2==float(vcat(collect(3:8:91),[97]))

    s = [1]
    s1 = compress2(s)
    @test s1==[1.0]
end

@testset "Window Computations" begin
    computewindow = TimeWarp.computewindow

    # Simplest path (along the diagonal)
    p=collect(1:8)
    rmin,rmax = computewindow(p,p,1)
    @test rmin==[1,1,1,2,3,4,5,6]
    @test rmax==[3,4,5,6,7,8,8,8]

    rmin,rmax = computewindow(p,p,2)
    @test rmin==[1,1,1,1,1,2,3,4]
    @test rmax==[5,6,7,8,8,8,8,8]

    # A warpy path
    pa=[1,1,2,3,4,5,6,7,8,8,8]
    pb=[1,2,3,3,3,4,4,5,6,7,8]
    rmin,rmax = computewindow(pa, pb, 1)
    @test pa[end]==length(rmin)
    @test pa[end]==length(rmax)
    @test rmin==[1,1,2,2,2,3,3,4]
    @test rmax==[4,4,4,5,5,6,8,8]

    rmin,rmax = computewindow(pa, pb, 2)
    @test pa[end]==length(rmin)
    @test pa[end]==length(rmax)
    @test rmin==[1,1,1,1,1,1,2,2]
    @test rmax==[5,5,6,6,7,8,8,8]

    rmin,rmax = computewindow(pa, pb, 20)
    @test pa[end]==length(rmin)
    @test pa[end]==length(rmax)
    @test rmin==fill(1,8)
    @test rmax==fill(8,8)

    # Extreme path: follows left then upper edge
    pa=[1,1,1,1,1,1,1,1,2,3,4,5,6,7,8]
    pb=[1,2,3,4,5,6,7,8,8,8,8,8,8,8,8]
    rmin,rmax = computewindow(pa, pb, 1)
    @test pa[end]==length(rmin)
    @test pa[end]==length(rmax)
    @test rmin==[1,1,7,7,7,7,7,7]
    @test rmax==[8,8,8,8,8,8,8,8]

    rmin,rmax = computewindow(pa, pb, 2)
    @test pa[end]==length(rmin)
    @test pa[end]==length(rmax)
    @test rmin==[1,1,1,6,6,6,6,6]
    @test rmax==[8,8,8,8,8,8,8,8]

    # More columns than rows
    pa=[1,2,3,4,5,6,7,8]
    pb=[1,2,3,4,4,4,4,4]
    rmin,rmax = computewindow(pa, pb, 1)
    @test rmin==[1,1,1,2,3,3,3,3]
    @test rmax==[3,4,4,4,4,4,4,4]

    rmin,rmax = computewindow(pa, pb, 2)
    @test rmin==[1,1,1,1,1,2,2,2]
    @test rmax==fill(4,8)

    rmin,rmax = computewindow(pa, pb, 3)
    @test rmin==fill(1,8)
    @test rmax==fill(4,8)

    rmin,rmax = computewindow(pa, pb, 4)
    @test rmin==fill(1,8)
    @test rmax==fill(4,8)

    rmin,rmax = computewindow(pa, pb, 47)
    @test rmin==fill(1,8)
    @test rmax==fill(4,8)

    # More rows than columns
    pa,pb = pb,pa
    rmin,rmax = computewindow(pa, pb, 1)
    @test rmin==[1,1,1,2]
    @test rmax==[3,4,8,8]

    rmin,rmax = computewindow(pa, pb, 2)
    @test rmin==fill(1,4)
    @test rmax==[5,8,8,8]

    rmin,rmax = computewindow(pa, pb, 3)
    @test rmin==fill(1,4)
    @test rmax==fill(8,4)

    rmin,rmax = computewindow(pa, pb, 4)
    @test rmin==fill(1,4)
    @test rmax==fill(8,4)

    rmin,rmax = computewindow(pa, pb, 47)
    @test rmin==fill(1,4)
    @test rmax==fill(8,4)
end

@testset "DTW and FastDTW agreement" begin
    t=collect(1:1600)
    pktimes=[100,300,1000,1300]
    x = 1*exp.(-0.5*((t .- pktimes[1])/100).^2);
    x+= 2*exp.(-0.5*((t .- pktimes[2])/150).^2);
    x+= 3*exp.(-0.5*((t .- pktimes[3])/250).^2);
    x+= 4*exp.(-0.5*((t .- pktimes[4])/250).^2);
    y = x[1:2:end];
    cost,px,py = dtw(x,y)
    cost1,qx,qy = fastdtw(x, y, 15)
    @test isapprox(cost,cost1)
    @test px == qx
    @test py == qy
end

@testset "DBA" begin
    x = Sequence([1., 2., 2., 3., 3., 4.])
    y = Sequence([1., 3., 4.])
    z = Sequence([1., 2., 2., 4.])
    avg, _ = dba([x, y, z], init_center=z)
    @test avg == [1.0, 1.75, 2.75, 4.0]
end


text:
# Copyright 2014 Mitchell Kember. Subject to the MIT License.
# Google Code Jam Qualification Round Africa 2010
# Problem C: T9 Spelling

# Converts a character into a string representing a keypad instruction.
function charseq(c)
	if c == ' '; return "0"; end
	q, r = divrem(int(c)-int('a'), 3)
	# Account for keys with 4 letters.
	if q > 5
		r -= 1
		if r == -1
			q -= 1
			r = (q == 5) ? 3 : 2
		elseif q == 8
			q = 7; r = 3
		end
	end
	# ABC begins at the second key.
	string(q+2)^(r+1)
end

# Inserts spaces to represent pauses between keypad instructions whenever the
# last character of one is identical to the first character of the next. Returns
# the modified instruction list.
function addpauses!(list)
	i = 2
	while i <= length(list)
		if list[i-1][end] == list[i][1]
			insert!(list, i, " ")
			i += 1
		end
		i += 1
	end
	list
end

# Convert a string of alphabetic characters and spaces into a sequence that can
# be entered using an old-fashioned cellphone keypad.
keypad(str) = join(addpauses!(map(charseq, collect(str))))

function main()
	n = int(readline())
	for i = 1:n
		str = keypad(chomp(readline()))
		println("Case #$i: ", str)
	end
end

main()


text:
using Element: generate_lebedev

# generate lebedev angular grids
generate_lebedev()

using BinDeps
using Compat

@BinDeps.setup

libxc = library_dependency("libxc")

@osx_only begin
	using Homebrew
	provides(Homebrew.HB, "libxc", libxc, os = :Darwin)
end

provides(AptGet, "libxc", libxc)
provides(Yum, "libxc-devel", libxc)

julia_usrdir = normpath(JULIA_HOME*"/../") # This is a stopgap, we need a better built-in solution to get the included libraries
libdirs = AbstractString["$(julia_usrdir)/lib"]
includedirs = AbstractString["$(julia_usrdir)/include"]

provides(
	Sources,
	URI("http://www.tddft.org/programs/octopus/down.php?file=libxc/libxc-2.2.2.tar.gz"),
	libxc
)

provides(
	BuildProcess,
	Autotools(lib_dirs = libdirs, include_dirs = includedirs),
	libxc
)

@compat @BinDeps.install Dict(:libxc => :libxc)

max_stars_repo_name: johncsnyder/Element


text:
using RoughSurface, Documenter
using Test

@testset "RoughSurface.jl" begin
    # Write your tests here.
    file_path = ["test.txt", "test.csv", "test.jld2"]
    surf = collect(reshape(1:3*3,3,3))
    @testset for i = 1:length(file_path)
        write_surf(file_path[i], surf)
        surf_read = read_surf(file_path[i])
        @test surf_read == surf
        rm(file_path[i])
    end # for
    # 文档
    doctest(RoughSurface)
end


text:
"""
This struct contains the reserve down market parameters
for price forcasts using CEM and endogeneous Economic Dispatch.
The reserve down market assumes inelastic demand.
"""

struct ReserveDownMarket{Z, T} #

    demand::AxisArrays.AxisArray{Float64, 2}         # MW
    price_cap::AxisArrays.AxisArray{Float64, 1}       # $/MWh

    function ReserveDownMarket{}(demand::AxisArrays.AxisArray{Float64, 2}, price_cap::AxisArrays.AxisArray{Float64, 1})
        @assert all(demand .>= 0)
        @assert all(price_cap .>= 0)
        Z = size(demand, 1)
        T = size(demand, 2)
        new{Z, T}(demand, price_cap)
    end

end


text:
abstract type AbstractPushInterpreter end

# Default is to have no stacks, override for types for which your vm has a stack.
stacktypes(i::AbstractPushInterpreter) = DataType[]
@inline hasstack(i::AbstractPushInterpreter, T::DataType) = in(T, stacktypes(i))

function execute!(interp::AbstractPushInterpreter, program::AbstractArray)
    for instr in program
        execute!(instr, interp)
    end
end

# The default stacks follow a certain naming structure/format for speed and simplicity.
# For this to work you must name the fields for your stacks accordingly!!
@inline stack(i::AbstractPushInterpreter, ::Type{Int}) = i.stack_Int
@inline stack(i::AbstractPushInterpreter, ::Type{Float64}) = i.stack_Float64
@inline stack(i::AbstractPushInterpreter, ::Type{String}) = i.stack_String
@inline stack(i::AbstractPushInterpreter, ::Type{Bool}) = i.stack_Bool
@inline execstack(i::AbstractPushInterpreter) = i.exec_stack

@inline stackhas(i::AbstractPushInterpreter, T::DataType, n::Int) = length(stack(i, T)) >= n
@inline pop!(i::AbstractPushInterpreter, I::DataType) = pop!(stack(i, I))
@inline pop2!(i::AbstractPushInterpreter, I::DataType) = pop2!(stack(i, I))
@inline push!(i::AbstractPushInterpreter, val::O) where {O} = push!(stack(i, O), val)

# Default execution behavior is to just push the value. This takes care of all literals.
@inline execute!(instr, i::AbstractPushInterpreter) = push!(i, instr)
# But if the instruction is an array we run it as a program
execute!(instr::Vector, i::AbstractPushInterpreter) = execute!(i, instr)

abstract type PushInstruction end

# This is the default type for things that can be part of a Push program.
const PushLiteralOrInstruction = Union{Int, Float64, Bool, String, PushInstruction}
const PushProgram = Union{Int, Float64, Bool, String, PushInstruction, Vector}

# The standard/default PushVM has stacks for all the common things and executes via
# the Exec stack.
struct PushInterpreter{T} <: AbstractPushInterpreter
    exec_stack::Stack{T}
    stack_Int::Stack{Int}
    stack_Float64::Stack{Float64}
    stack_Bool::Stack{Bool}
    stack_String::Stack{String}
    PushInterpreter{T}() where {T} = begin
        new(Stack{T}(),
            Stack{Int}(), Stack{Float64}(),
            Stack{Bool}(), Stack{String}())
    end
end
PushInterpreter() = PushInterpreter{PushProgram}()
stacktypes(i::PushInterpreter) = DataType[Int, Float64, String, Bool]

function execute!(interp::PushInterpreter, program::AbstractArray)
    # Push in reverse order
    for instr in Iterators.reverse(program)
        push!(interp.exec_stack, instr)
    end
    # Then run from exec_stack
    run_from_execstack!(interp)
end

function run_from_execstack!(interp::PushInterpreter)
    while length(execstack(interp)) > 0
        i = pop!(execstack(interp))
        execute!(i, interp)
    end
end


text:
"""
    ArrayIndex{N}

Subtypes of `ArrayIndex` represent series of transformations for a provided index to some
buffer which is typically accomplished with square brackets (e.g., `buffer[index[inds...]]`).
The only behavior that is required of a subtype of `ArrayIndex` is the ability to transform
individual index elements (i.e. not collections). This does not guarantee bounds checking or
the ability to iterate (although additional functionallity may be provided for specific
types).
"""
abstract type ArrayIndex{N} end

const MatrixIndex = ArrayIndex{2}

const VectorIndex = ArrayIndex{1}

Base.ndims(::Type{<:ArrayIndex{N}}) where {N} = N

struct BidiagonalIndex <: MatrixIndex
    count::Int
    isup::Bool
end

struct TridiagonalIndex <: MatrixIndex
    count::Int# count==nsize+nsize-1+nsize-1
    nsize::Int
    isrow::Bool
end

Base.firstindex(i::Union{BidiagonalIndex,TridiagonalIndex}) = 1
Base.lastindex(i::Union{BidiagonalIndex,TridiagonalIndex}) = i.count
Base.length(i::Union{BidiagonalIndex,TridiagonalIndex}) = lastindex(i)

"""
    StrideIndex(x)

Subtype of `ArrayIndex` that transforms and index using stride layout information
derived from `x`.
"""
struct StrideIndex{N,R,C,S,O} <: ArrayIndex{N}
    strides::S
    offsets::O

    function StrideIndex{N,R,C}(s::S, o::O) where {N,R,C,S,O}
        return new{N,R::NTuple{N,Int},C,S,O}(s, o)
    end
    function StrideIndex{N,R,C}(a::A) where {N,R,C,A}
        return StrideIndex{N,R,C}(strides(a), offsets(a))
    end
    function StrideIndex(a::A) where {A}
        return StrideIndex{ndims(A),known(stride_rank(A)), known(contiguous_axis(A))}(a)
    end
end

## getindex
@propagate_inbounds Base.getindex(x::ArrayIndex, i::CanonicalInt, ii::CanonicalInt...) = x[NDIndex(i, ii...)]
@propagate_inbounds function Base.getindex(ind::BidiagonalIndex, i::Int)
    @boundscheck 1 <= i <= ind.count || throw(BoundsError(ind, i))
    if ind.isup
        ii = i + 1
    else
        ii = i + 1 + 1
    end
    convert(Int, floor(ii / 2))
end

@propagate_inbounds function Base.getindex(ind::TridiagonalIndex, i::Int)
    @boundscheck 1 <= i <= ind.count || throw(BoundsError(ind, i))
    offsetu = ind.isrow ? 0 : 1
    offsetl = ind.isrow ? 1 : 0
    if 1 <= i <= ind.nsize
        return i
    elseif ind.nsize < i <= ind.nsize + ind.nsize - 1
        return i - ind.nsize + offsetu
    else
        return i - (ind.nsize + ind.nsize - 1) + offsetl
    end
end

@inline function Base.getindex(x::StrideIndex{N}, i::AbstractCartesianIndex) where {N}
    return _strides2int(offsets(x), strides(x), Tuple(i)) + static(1)
end
@generated function _strides2int(o::O, s::S, i::I) where {O,S,I}
    N = known_length(S)
    out = :()
    for i in 1:N
        tmp = :(((getfield(i, $i) - getfield(o, $i)) * getfield(s, $i)))
        out = ifelse(i === 1, tmp, :($out + $tmp))
    end
    return Expr(:block, Expr(:meta, :inline), out)
end


text:
# Type System

const uid = Ref(UInt64(0))

abstract type Shape end

struct AnyT <: Shape end

struct Domain <: Shape
  id::UInt64
end

Domain() = Domain(uid[] += 1)

struct ArrayT <: Shape
  xs::Vector{Any}
end

ArrayT(xs::Domain...) = ArrayT([xs...])

struct TupleT <: Shape
  xs::Vector{Any}
end

TupleT(xs...) = TupleT([xs...])

struct Arrow <: Shape
  xs::Vector{Any}
end

Arrow(xs...) = Arrow([xs...])

tostring(d::Domain, ctx) = "d$(d.id)"
tostring(d::ArrayT, ctx) = string("[", join(map(x -> tostring(x, ctx), d.xs), ","), "]")
tostring(d::TupleT, ctx) = string("(", join(map(x -> tostring(x, ctx), d.xs), ","), ")")
tostring(d::Arrow, ctx) = join(map(x -> tostring(x, ctx), d.xs), " → ")
tostring(d::AnyT, ctx) = "*"

tostring(x) = tostring(x, (Dict(), ['a':'z'...]))
Base.show(io::IO, x::Shape) = print(io, tostring(x))

shape(::Type{<:AbstractArray{T,N}}) where {T,N} = ArrayT([Domain() for i = 1:N])
shape(T::Type{<:Tuple}) = TupleT([shape(x) for x in T.parameters])

# Type Inference

using DataFlow.Interpreter
using DataFlow: Call

function lower(m::Associative, x::Domain)
  while haskey(m, x)
    x = m[x]
  end
  return x
end

lower(m::Associative, x::Shape) = typeof(x)(lower.(m, x.xs)...)

function lower(m::Associative, v::IVertex)
  prewalkλ(v) do v
    isconstant(v) && v.value.value isa Shape ? constant(lower(m, v.value.value)) : v
  end
end

function unify(ctx::Context, a::Domain, b::Domain)
  a, b = lower.(ctx[:typemap], (a, b))
  a ≠ b && (ctx[:typemap][a] = b)
  return
end

withtype(v, T) = vertex(TypeAssert(), v, constant(T))
withtype(v, ::Void) = v
vtype(v::IVertex) = v.value isa TypeAssert ? v[2].value.value : nothing

function iarray(f, ctx::Context, ::Call, T::Type{<:AbstractArray})
  ds = [Domain() for _ = 1:ndims(T)]
  ctor = vertex(Call(), constant(T), constant.(ds)...)
  withtype(ctor, ArrayT(ds...))
end

function iarray(f, ctx::Context, ::Call, T::Type{<:Ref}, x)
  ctor = vertex(Call(), constant(T), constant(x))
  withtype(ctor, ArrayT())
end

iarray(f, ctx::Context, v, args...) = f(ctx, v, args...)

function iloop(f, ctx::Context, ::Loop, l, args...)
  is = [Domain() for i = 1:l.value.args]
  vars = vtype.(l.inputs)
  inputs = [withtype(inputnode(i), sh) for (i, sh) in enumerate([vars..., is...])]
  body = interpret(ctx, l.value.body, inputs...)
  withtype(
    vertex(Loop(), vertex(Lambda(l.value.args, body), l.inputs...), args..., constant.(is)...),
    vtype(args[1]))
end

iloop(f, ctx::Context, v, args...) = f(ctx, v, args...)

function iindex(f, ctx::Context, ::Call, ::typeof(getindex), xs, is...)
  foreach(d -> unify(ctx, d...), zip(vtype(xs).xs, vtype.(is)))
  vertex(Call(), constant(getindex), xs, is...)
end

function iindex(f, ctx::Context, ::Call, ::typeof(setindex!), xs, v, is...)
  foreach(d -> unify(ctx, d...), zip(vtype(xs).xs, vtype.(is)))
  vertex(Call(), constant(setindex!), xs, v, is...)
end

iindex(f, ctx::Context, v, args...) = f(ctx, v, args...)

function iinline(f, ctx::Context, ::Call, v, args...)
  isconstant(v) && v.value.value isa Func || return f(ctx, Call(), v, args...)
  interpret(ctx, v.value.value.graph, args...)
end

iinline(f, ctx::Context, v, args...) = f(ctx, v, args...)

interp(ctx::Context, f, args...) = vertex(f, constant.(args)...)

function infer(v::IVertex, ts::Shape...)
  inputs = [vertex(TypeAssert(), inputnode(i), constant(ts[i])) for i = 1:length(ts)]
  typemap = Dict()
  ctx = Context(mux(iline, iinline, iconst, iargs, iarray, iloop, iindex, interp),
                typemap = typemap)
  v = interpret(ctx, v, inputs...)
  v = lower(typemap, v)
  Arrow(lower.(typemap, ts)..., vtype(v)), v
end

infer(v::IVertex, ts::Type...) = infer(v, shape.(ts)...)

infer(f::Func, ts...) = infer(f.graph, ts...)

lower(x, ts...) = infer(x, ts...)[2]

# Post-processing

istype(v) = v.value isa DataFlow.TypeAssert
notype(v) = istype(v) ? notype(v[1]) : v

striptypes(v::IVertex) = prewalkλ(notype, v)

domainv(sh::Domain, d, v) = sh == d ? v : nothing

domainv(sh::ArrayT, d, v) =
  get(filter(x -> x ≠ nothing,
             domainv.(sh.xs, d,
                      vertex.(domain, v, constant.(1:length(sh.xs))))),
      1, nothing)

domainv(sh::Arrow, d) =
  filter(x -> x ≠ nothing,
         domainv.(sh.xs, d, inputnode.(1:length(sh.xs))))[1]

function insert_domains(v::IVertex, t::Shape)
  prewalkλ(v) do v
    isconstant(v) && v.value.value isa Shape || return v
    domainv(t, v.value.value)
  end
end

max_stars_repo_name: JuliaTagBot/Tokamak

max_issues_repo_name: JuliaTagBot/Tokamak

max_forks_repo_name: FluxML/Tokamak


text:
using Test
using Transit
using JSON

import DataStructures.OrderedDict

# Convert the value given to transit and then read
# back in the resulting JSON. Not a round trip, since
# we are just reading back in the raw JSON.
function square_trip(x, verbose=false)
    s = sprint() do io
        e = Encoder(io, verbose)
        encode(e, x, false)
    end
    JSON.parse(s)
end

@test square_trip(1) == 1
@test square_trip(1.0) == 1.0
@test square_trip("hello") == "hello"
@test square_trip("~hello") == "~~hello"
@test square_trip(true) == true
@test square_trip(false) == false
@test square_trip(2//3) == Any["~#ratio", Any[2, 3]]
@test square_trip([1,2,3]) == [1,2,3]
@test square_trip((1,2,"hello")) == Any["~#list", Any[1,2,"hello"]]

@test square_trip([:hello, :hello, :hello]) == ["~:hello", "^0", "^0"]
@test square_trip([:aaaa, :bbbb, :aaaa, :bbbb]) == ["~:aaaa", "~:bbbb", "^0", "^1"]

symbols = [Transit.TSymbol("hello"), Transit.TSymbol("hello"), Transit.TSymbol("hello")]

@test square_trip(symbols) == ["~\$hello", "^0", "^0"]
@test square_trip([:aaaa, :bbbb, :aaaa, :bbbb]) == ["~:aaaa", "~:bbbb", "^0", "^1"]

@test square_trip([:aaaa, :bbbb, :aaaa, :bbbb], true) == ["~:aaaa", "~:bbbb", "~:aaaa", "~:bbbb"]

@test square_trip(Dict([("aaa", 11), ("bbb", 22)]), true) == Dict([("aaa", 11), ("bbb", 22)])


text:
@rosimport std_msgs.msg: ColorRGBA
@rosimport geometry_msgs.msg: Point
@rosimport visualization_msgs.msg: Marker
@rosimport safe_traffic_weaving.msg: VehicleTrajectory, path, XYThV
@rosimport auto_messages.msg: from_autobox, to_autobox
rostypegen(@__MODULE__)
import .std_msgs.msg: ColorRGBA
import .geometry_msgs.msg: Point
import .visualization_msgs.msg: Marker
import .safe_traffic_weaving.msg: VehicleTrajectory, path, XYThV
import .auto_messages.msg: from_autobox, to_autobox

TrajectoryTube(p::path) = TrajectoryTube{Float64}(
    invcumtrapz(p.Ux_des_mps, p.s_m), p.s_m, p.Ux_des_mps, p.Ax_des_mps2,
    p.posE_m, p.posN_m, p.Psi_rad, p.k_1pm, p.grade_rad, 0*p.grade_rad, p.edge_L_m, p.edge_R_m
)
TrajectoryTube(p::VehicleTrajectory) = TrajectoryTube{Float64}(
    p.t, p.s, p.V, p.A, p.E, p.N, p.heading, p.curvature, p.grade, p.bank, p.edge_L, p.edge_R
)

### DEBUG
const show_loginfo = fill(false)

### HJI visualization
include("rviz.jl")

### /des_path and /des_traj
const latest_trajectory = fill(straight_trajectory(30., 5.))
const tracking_mode = fill(:path)
function nominal_trajectory_callback(msg::path, mpc=X1DMPC)                 # /des_traj
    latest_trajectory[] = TrajectoryTube(msg)
    tracking_mode[] = :path
    mpc.time_offset = NaN
    mpc.solved = false
end
function nominal_trajectory_callback(msg::VehicleTrajectory, mpc=X1CMPC)    # /des_path
    latest_trajectory[] = TrajectoryTube(msg)
    tracking_mode[] = :traj
    mpc.time_offset = convert(Float64, msg.header.stamp)
    mpc.solved = false
end

### /to_autobox
const to_autobox_msg = to_autobox()

### /from_autobox
const use_HJI_policy = fill(false)
function from_autobox_callback(msg::from_autobox, to_autobox_pub, HJI_values_pub, HJI_contour_pub, path_mpc=X1DMPC, traj_mpc=X1CMPC)
    mpc = (tracking_mode[] == :path ? path_mpc : traj_mpc)
    mpc.current_state = BicycleState(msg.E_m, msg.N_m, msg.psi_rad, msg.ux_mps, msg.uy_mps, msg.r_radps)
    # mpc.current_control = BicycleControl(msg.delta_rad, msg.fxf_N, msg.fxr_N)
    mpc.current_control = BicycleControl(to_autobox_msg.delta_cmd_rad, to_autobox_msg.fxf_cmd_N, to_autobox_msg.fxr_cmd_N)
    mpc.trajectory = latest_trajectory[]

    ###### TODO: very messy
    relative_state = HJIRelativeState(mpc.current_state, traj_mpc.other_car_state)
    V, ∇V = mpc.HJI_cache[relative_state]
    show_loginfo[] && RobotOS.loginfo("Pigeon MPC: HJI value function = $V")
    # show_loginfo[] && RobotOS.loginfo("Pigeon MPC: HJI relative state = $(relative_state)")
    try
        update_HJI_values_marker!(HJI_values_marker, relative_state)
        update_HJI_contour_marker!(HJI_contour_marker, relative_state)
        publish(HJI_values_pub, HJI_values_marker)
        publish(HJI_contour_pub, HJI_contour_marker)
    catch err
        RobotOS.logwarn("HJI Visualization Error: $err\n$(stacktrace(catch_backtrace()))")
    end
    ######

    if msg.pre_flag == 0
        show_loginfo[] && RobotOS.loginfo("Pigeon MPC: /from_autobox pre_flag == 0, MPC inactive")
        return
    end
    if isnan(mpc.time_offset)
        show_loginfo[] && RobotOS.loginfo("Pigeon MPC: time_offset not set, running in path tracking mode")
        _, _, t = path_coordinates(mpc.trajectory, mpc.current_state)
    else
        t = convert(Float64, msg.header.stamp) - mpc.time_offset
        if t < 0 || t > mpc.trajectory.t[end]
            show_loginfo[] && RobotOS.loginfo("Pigeon MPC: current time $t outside trajectory interval [0, $(mpc.trajectory.t[end])]")
            return
        end
    end
    if mpc.current_state.Ux < 1
        show_loginfo[] && RobotOS.loginfo("Pigeon MPC: current speed < 1, pausing MPC while X1 is stopped")
        return
    end
    MPC_steps_missed = msg.header.seq - (mpc.heartbeat + 1)
    if MPC_steps_missed != 0
        RobotOS.logwarn("Pigeon MPC: $(MPC_steps_missed) from_autobox messages lost")
        mpc.heartbeat = msg.header.seq - 1
    end

    t_elapsed = @elapsed begin
        try
            compute_time_steps!(mpc, t)
            compute_linearization_nodes!(mpc)
            update_QP!(mpc)
            solve!(mpc)
        catch err
            RobotOS.logwarn("Pigeon MPC Error: $err\n$(stacktrace(catch_backtrace()))")
        end
    end

    if t_elapsed > 0.01
        RobotOS.logwarn("Pigeon MPC: OSQP took $(1000*t_elapsed) ms at heartbeat $(mpc.heartbeat)")
    else
        show_loginfo[] && RobotOS.loginfo("Pigeon MPC: OSQP took $(1000*t_elapsed) ms at heartbeat $(mpc.heartbeat)")
    end
    # show_loginfo[] && RobotOS.loginfo("Pigeon MPC: $(mpc.model.optimizer.results.info)")
    # show_loginfo[] && RobotOS.loginfo("deltas: $(value.(mpc.model, mpc.variables.δ))")
    mpc.heartbeat += 1

    s, e, _ = path_coordinates(mpc.trajectory, mpc.current_state)
    if tracking_mode[] == :traj && use_HJI_policy[] && V <= mpc.HJI_ϵ
        u_next = BicycleControl(mpc.dynamics.longitudinal_params, optimal_control(mpc.dynamics, relative_state, ∇V))
        RobotOS.logwarn("Pigeon MPC: HJI stepping in to save the day (with a hammer)")
        RobotOS.logwarn("Pigeon MPC: HJI value function = $V")
    else
        if tracking_mode[] == :traj && V <= mpc.HJI_ϵ
            RobotOS.logwarn("Pigeon MPC: HJI stepping in to save the day (with a feather)")
            RobotOS.logwarn("Pigeon MPC: HJI value function = $V")
        end
        u_next = get_next_control(mpc)
    end
    to_autobox_msg.header.stamp  = RobotOS.now()
    to_autobox_msg.post_flag     = 1    # TODO: check for OSQP failure
    to_autobox_msg.heartbeat     = mpc.heartbeat
    to_autobox_msg.s_m           = s
    to_autobox_msg.e_m           = e
    to_autobox_msg.delta_cmd_rad = u_next.δ
    to_autobox_msg.fxf_cmd_N     = u_next.Fxf
    to_autobox_msg.fxr_cmd_N     = u_next.Fxr
    if isnan(to_autobox_msg.delta_cmd_rad) || isnan(to_autobox_msg.fxf_cmd_N) || isnan(to_autobox_msg.fxr_cmd_N)
        show_loginfo[] && RobotOS.loginfo("Pigeon MPC: OSQP returned NaNs " *
                                          "($(to_autobox_msg.delta_cmd_rad), $(to_autobox_msg.fxf_cmd_N), $(to_autobox_msg.fxr_cmd_N))" *
                                          "; using previous control "*
                                          "($(mpc.current_control.δ), $(mpc.current_control.Fxf), $(mpc.current_control.Fxr))")
        to_autobox_msg.delta_cmd_rad = mpc.current_control.δ
        to_autobox_msg.fxf_cmd_N     = mpc.current_control.Fxf
        to_autobox_msg.fxr_cmd_N     = mpc.current_control.Fxr
        publish(to_autobox_pub, to_autobox_msg)
        to_autobox_msg.delta_cmd_rad = 0    # ensures that if we get 2 NaNs in a row, we return to 0 control
        to_autobox_msg.fxf_cmd_N     = 0
        to_autobox_msg.fxr_cmd_N     = 0
        Parametron.initialize!(mpc.model)   # reinitialize model (otherwise NaNs might stick)
        mpc.solved = false    # any warm starting should ignore the previous solution
    else
        publish(to_autobox_pub, to_autobox_msg)
    end
end

function other_car_callback(msg::XYThV, mpc=X1CMPC)
    mpc.other_car_state = SimpleCarState(msg.x, msg.y, msg.th - pi/2, msg.v)
end

### ROS node init
function start_ROS_node()
    init_node("pigeon", anonymous=false)
    other_car = RobotOS.get_param("human", "/xbox_car")
    to_autobox_pub = Publisher{to_autobox}("/to_autobox", queue_size=10)
    HJI_values_pub = Publisher{Marker}("/HJI_values", queue_size=1)
    HJI_contour_pub = Publisher{Marker}("/HJI_contour", queue_size=1)
    Subscriber{path}("/des_path", nominal_trajectory_callback, queue_size=1)
    Subscriber{VehicleTrajectory}("/des_traj", nominal_trajectory_callback, queue_size=1)
    Subscriber{from_autobox}("/from_autobox", from_autobox_callback, (to_autobox_pub, HJI_values_pub, HJI_contour_pub), queue_size=1)
    Subscriber{XYThV}("$(other_car)/xythv", other_car_callback, queue_size=1)    # TODO: abstract with a republisher
    @spawn spin()
end

max_stars_repo_name: Michelle-NYX/CollisionAvoidance

max_issues_repo_name: Michelle-NYX/CollisionAvoidance

max_forks_repo_name: Michelle-NYX/CollisionAvoidance


text:
const RowLike = Union{NamedTuple,Tables.Row,Tables.AbstractRow}

if isdefined(Base, :ComposedFunction) # Julia >= 1.6.0-DEV.85
    using Base: ComposedFunction
else
    const ComposedFunction = let h = identity ∘ convert
        @assert h.f === identity
        @assert h.g === convert
        getfield(parentmodule(typeof(h)), nameof(typeof(h)))
    end
    @assert identity ∘ convert isa ComposedFunction
end

const GenericInitializer = Union{typeof(Transducers.Init),Transducers.InitOf}

const InitializerFor{OP} = Union{GenericInitializer,InitialValues.GenericInitialValue{OP}}

# Just like `Function` but for defining some common methods.
abstract type _Function <: Function end

# Avoid `Function` fallbacks:
@nospecialize
Base.show(io::IO, ::MIME"text/plain", f::_Function) = show(io, f)
Base.print(io::IO, f::_Function) = show(io, f)
@specialize


text:
# jhu.jl
#######################

# Load packages.
using Plots, CSV, Dates, DataFrames

# Download data from the MHLW web site.
download("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv","./csv/jhu_cases.csv");
download("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv","./csv/jhu_deaths.csv");
download("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv","./csv/jhu_recovered.csv");

Acsv=DataFrame(CSV.File("./csv/jhu_cases.csv", header=false, delim=','));
Bcsv=DataFrame(CSV.File("./csv/jhu_deaths.csv", header=false, delim=','));
Ccsv=DataFrame(CSV.File("./csv/jhu_recovered.csv", header=false, delim=','));
(pw,qw)=size(Acsv);

D=Int64(qw-4);
d0=Date(2020,1,22);
d1=d0+Day(D-1);
l0=string(d0);
l1=string(d0+Day(Int(floor((D-1)/4))));
l2=string(d0+Day(Int(floor((D-1)/2))));
l3=string(d0+Day(Int(floor(3*(D-1)/4))));
l4=string(d0+Day(D-1));
DD=Int64(D-526);
dd0=Date(2021,7,1);
dd1=d0+Day(DD-1);
ll0=string(dd0);
ll1=string(dd0+Day(Int(floor((DD-1)/2))));
ll2=string(dd0+Day(DD-1));

# Argentina
PARG=45.670451;
AARG=parse.(Float64,Array(Acsv[8,5:qw]))/PARG;
NARG=zeros(D);
for j=1:7
    NARG[j]=AARG[j]/7
end
for j=8:D
    NARG[j]=(AARG[j]-AARG[j-7])/7
end
BARG=parse.(Float64,Array(Bcsv[8,5:qw]))/PARG;
NDARG=zeros(D);
for j=1:7
    NDARG[j]=BARG[j]/7
end
for j=8:D
    NDARG[j]=(BARG[j]-BARG[j-7])/7
end

# Bangladesh
PBGD=166.567052;
ABGD=parse.(Float64,Array(Acsv[22,5:qw]))/PBGD;
NBGD=zeros(D);
for j=1:7
    NBGD[j]=ABGD[j]/7
end
for j=8:D
    NBGD[j]=(ABGD[j]-ABGD[j-7])/7
end
BBGD=parse.(Float64,Array(Bcsv[22,5:qw]))/PBGD;
NDBGD=zeros(D);
for j=1:7
    NDBGD[j]=BBGD[j]/7
end
for j=8:D
    NDBGD[j]=(BBGD[j]-BBGD[j-7])/7
end

# Brazil
PBRA=214.289417;
ABRA=parse.(Float64,Array(Acsv[32,5:qw]))/PBRA;
NBRA=zeros(D);
for j=1:7
    NBRA[j]=ABRA[j]/7
end
for j=8:D
    NBRA[j]=(ABRA[j]-ABRA[j-7])/7
end
BBRA=parse.(Float64,Array(Bcsv[32,5:qw]))/PBRA;
NDBRA=zeros(D);
for j=1:7
    NDBRA[j]=BBRA[j]/7
end
for j=8:D
    NDBRA[j]=(BBRA[j]-BBRA[j-7])/7
end

# Colombia
PCOL=51.503463;
ACOL=parse.(Float64,Array(Acsv[94,5:qw]))/PCOL;
NCOL=zeros(D);
for j=1:7
    NCOL[j]=ACOL[j]/7
end
for j=8:D
    NCOL[j]=(ACOL[j]-ACOL[j-7])/7
end
BCOL=parse.(Float64,Array(Bcsv[94,5:qw]))/PCOL;
NDCOL=zeros(D);
for j=1:7
    NDCOL[j]=BCOL[j]/7
end
for j=8:D
    NDCOL[j]=(BCOL[j]-BCOL[j-7])/7
end

# France Italy Spain
PFRA=65.439014;
PITA=60.359899;
PESP=46.775535;
PFIE=PFRA+PITA+PESP;
CFRA=parse.(Float64,Array(Acsv[132,5:qw]));
CITA=parse.(Float64,Array(Acsv[155,5:qw]));
CESP=parse.(Float64,Array(Acsv[239,5:qw]));
AFIE=(CFRA+CITA+CESP)/PFIE;
NFIE=zeros(D);
for j=1:7
    NFIE[j]=AFIE[j]/7
end
for j=8:D
    NFIE[j]=max(0,(AFIE[j]-AFIE[j-7])/7)
end
DFRA=parse.(Float64,Array(Bcsv[132,5:qw]))
DITA=parse.(Float64,Array(Bcsv[155,5:qw]));
DESP=parse.(Float64,Array(Bcsv[239,5:qw]));
BFIE=(DFRA+DITA+DESP)/PFIE;
NDFIE=zeros(D);
for j=1:7
    NDFIE[j]=BFIE[j]/7
end
for j=8:D
    NDFIE[j]=(BFIE[j]-BFIE[j-7])/7
end

# India
PIND=1395.531433;
AIND=parse.(Float64,Array(Acsv[149,5:qw]))/PIND;
NIND=zeros(D);
for j=1:7
    NIND[j]=AIND[j]/7
end
for j=8:D
    NIND[j]=(AIND[j]-AIND[j-7])/7
end
BIND=parse.(Float64,Array(Bcsv[149,5:qw]))/PIND;
NDIND=zeros(D);
for j=1:7
    NDIND[j]=BIND[j]/7
end
for j=8:D
    NDIND[j]=(BIND[j]-BIND[j-7])/7
end

# Indonesia
PIDN=276.833206;
AIDN=parse.(Float64,Array(Acsv[150,5:qw]))/PIDN;
NIDN=zeros(D);
for j=1:7
    NIDN[j]=AIDN[j]/7
end
for j=8:D
    NIDN[j]=(AIDN[j]-AIDN[j-7])/7
end
BIDN=parse.(Float64,Array(Bcsv[150,5:qw]))/PIDN;
NDIDN=zeros(D);
for j=1:7
    NDIDN[j]=BIDN[j]/7
end
for j=8:D
    NDIDN[j]=(BIDN[j]-BIDN[j-7])/7
end

# Japan
PJPN=126.032481;
AJPN=parse.(Float64,Array(Acsv[157,5:qw]))/PJPN;
NJPN=zeros(D);
for j=1:7
    NJPN[j]=AJPN[j]/7
end
for j=8:D
    NJPN[j]=(AJPN[j]-AJPN[j-7])/7
end
BJPN=parse.(Float64,Array(Bcsv[157,5:qw]))/PJPN;
NDJPN=zeros(D);
for j=1:7
    NDJPN[j]=BJPN[j]/7
end
for j=8:D
    NDJPN[j]=(BJPN[j]-BJPN[j-7])/7
end

# South Korea
PKOR=51.319753;
AKOR=parse.(Float64,Array(Acsv[162,5:qw]))/PKOR;
NKOR=zeros(D);
for j=1:7
    NKOR[j]=AKOR[j]/7
end
for j=8:D
    NKOR[j]=(AKOR[j]-AKOR[j-7])/7
end
BKOR=parse.(Float64,Array(Bcsv[162,5:qw]))/PKOR;
NDKOR=zeros(D);
for j=1:7
    NDKOR[j]=BKOR[j]/7
end
for j=8:D
    NDKOR[j]=(BKOR[j]-BKOR[j-7])/7
end

# Malaysia
PMYS=32.7;
AMYS=parse.(Float64,Array(Acsv[178,5:qw]))/PMYS;
NMYS=zeros(D);
for j=1:7
    NMYS[j]=AMYS[j]/7
end
for j=8:D
    NMYS[j]=(AMYS[j]-AMYS[j-7])/7
end
BMYS=parse.(Float64,Array(Bcsv[178,5:qw]))/PMYS;
NDMYS=zeros(D);
for j=1:7
    NDMYS[j]=BMYS[j]/7
end
for j=8:D
    NDMYS[j]=(BMYS[j]-BMYS[j-7])/7
end

# Mexico
PMEX=130.482814;
AMEX=parse.(Float64,Array(Acsv[185,5:qw]))/PMEX;
NMEX=zeros(D);
for j=1:7
    NMEX[j]=AMEX[j]/7
end
for j=8:D
    NMEX[j]=(AMEX[j]-AMEX[j-7])/7
end
BMEX=parse.(Float64,Array(Bcsv[185,5:qw]))/PMEX;
NDMEX=zeros(D);
for j=1:7
    NDMEX[j]=BMEX[j]/7
end
for j=8:D
    NDMEX[j]=(BMEX[j]-BMEX[j-7])/7
end

# Nepal
PNPL=29.735589;
ANPL=parse.(Float64,Array(Acsv[194,5:qw]))/PNPL;
NNPL=zeros(D);
for j=1:7
    NNPL[j]=ANPL[j]/7
end
for j=8:D
    NNPL[j]=(ANPL[j]-ANPL[j-7])/7
end
BNPL=parse.(Float64,Array(Bcsv[194,5:qw]))/PNPL;
NDNPL=zeros(D);
for j=1:7
    NDNPL[j]=BNPL[j]/7
end
for j=8:D
    NDNPL[j]=(BNPL[j]-BNPL[j-7])/7
end

# Pakistan 205 190
PPAK=225.791619;
APAK=parse.(Float64,Array(Acsv[208,5:qw]))/PPAK;
NPAK=zeros(D);
for j=1:7
    NPAK[j]=APAK[j]/7
end
for j=8:D
    NPAK[j]=(APAK[j]-APAK[j-7])/7
end
BPAK=parse.(Float64,Array(Bcsv[208,5:qw]))/PPAK;
NDPAK=zeros(D);
for j=1:7
    NDPAK[j]=BPAK[j]/7
end
for j=8:D
    NDPAK[j]=(BPAK[j]-BPAK[j-7])/7
end

# Philippines
PPHI=111.249116;
APHI=parse.(Float64,Array(Acsv[214,5:qw]))/PPHI;
NPHI=zeros(D);
for j=1:7
    NPHI[j]=APHI[j]/7
end
for j=8:D
    NPHI[j]=(APHI[j]-APHI[j-7])/7
end
BPHI=parse.(Float64,Array(Bcsv[214,5:qw]))/PPHI;
NDPHI=zeros(D);
for j=1:7
    NDPHI[j]=BPHI[j]/7
end
for j=8:D
    NDPHI[j]=(BPHI[j]-BPHI[j-7])/7
end

# Sri Lanka
PLKA=21.516097;
ALKA=parse.(Float64,Array(Acsv[240,5:qw]))/PLKA;
NLKA=zeros(D);
for j=1:7
    NLKA[j]=ALKA[j]/7
end
for j=8:D
    NLKA[j]=(ALKA[j]-ALKA[j-7])/7
end
BLKA=parse.(Float64,Array(Bcsv[240,5:qw]))/PLKA;
NDLKA=zeros(D);
for j=1:7
    NDLKA[j]=BLKA[j]/7
end
for j=8:D
    NDLKA[j]=(BLKA[j]-BLKA[j-7])/7
end

# Thailand
PTHA=70.000662;
ATHA=parse.(Float64,Array(Acsv[250,5:qw]))/PTHA;
NTHA=zeros(D);
for j=1:7
    NTHA[j]=ATHA[j]/7
end
for j=8:D
    NTHA[j]=(ATHA[j]-ATHA[j-7])/7
end
BTHA=parse.(Float64,Array(Bcsv[250,5:qw]))/PTHA;
NDTHA=zeros(D);
for j=1:7
    NDTHA[j]=BTHA[j]/7
end
for j=8:D
    NDTHA[j]=(BTHA[j]-BTHA[j-7])/7
end

# United Staes
PUSA=333.225477;
AUSA=parse.(Float64,Array(Acsv[256,5:qw]))/PUSA;
NUSA=zeros(D);
for j=1:7
    NUSA[j]=AUSA[j]/7
end
for j=8:D
    NUSA[j]=(AUSA[j]-AUSA[j-7])/7
end
BUSA=parse.(Float64,Array(Bcsv[256,5:qw]))/PUSA;
NDUSA=zeros(D);
for j=1:7
    NDUSA[j]=BUSA[j]/7
end
for j=8:D
    NDUSA[j]=(BUSA[j]-BUSA[j-7])/7
end

# United Kingdom
PGBR=68.294438;
AGBR=parse.(Float64,Array(Acsv[271,5:qw]))/PGBR;
NGBR=zeros(D);
for j=1:7
    NGBR[j]=AGBR[j]/7
end
for j=8:D
    NGBR[j]=(AGBR[j]-AGBR[j-7])/7
end
BGBR=parse.(Float64,Array(Bcsv[271,5:qw]))/PGBR;
NDGBR=zeros(D);
for j=1:7
    NDGBR[j]=BGBR[j]/7
end
for j=8:D
    NDGBR[j]=(BGBR[j]-BGBR[j-7])/7
end

p1=plot([AJPN APHI AMYS AIDN AIND ANPL APAK ABGD ALKA ATHA],
    grid=false,
    linewidth=2,
    title="COVID-19: cases per 1M \n data sourced by JHU Coronavirus Resource Center",
    right_margin=Plots.Measures.Length(:mm, 10.0),
    left_margin=Plots.Measures.Length(:mm, 5.0),
    xlabel="date",
    xticks = ([1 floor(D/4)  floor(D/2) floor(3*D/4) D;], [l0 l1 l2 l3 l4]),
    yaxis="cases/1M",
    legendfont=font(10),
    label=["Japan" "Philippines" "Malaysia" "Indonesia" "India" "Nepal" "Pakistan" "Bangladesh" "Sri Lanka" "Thailand"],
    palette = :seaborn_bright,
    legend = :topleft)
savefig("./jhu/jhu_asia_cases.png")

p2=plot([BJPN BPHI BMYS BIDN BIND BNPL BPAK BBGD BLKA BTHA],
    grid=false,
    linewidth=2,
    title="COVID-19: deaths per 1M \n data sourced by JHU Coronavirus Resource Center",
    right_margin=Plots.Measures.Length(:mm, 10.0),
    left_margin=Plots.Measures.Length(:mm, 5.0),
    xticks = ([1 floor(D/4)  floor(D/2) floor(3*D/4) D;], [l0 l1 l2 l3 l4]),
    xlabel="date",
    yaxis="deaths/1M",
    legendfont=font(10),
    label=["Japan" "Philippines" "Malaysia" "Indonesia" "India" "Nepal" "Pakistan" "Bangladesh" "Sri Lanka" "Thailand"],
   palette = :seaborn_bright,
     legend = :topleft)
savefig("./jhu/jhu_asia_deaths.png")

p3=plot([NJPN NPHI NMYS NIDN NIND NNPL NPAK NBGD NLKA NTHA],
    grid=false,
    linewidth=2,
    title="COVID-19: 7-day average of new cases per 1M \n data sourced by JHU Coronavirus Resource Center",
    right_margin=Plots.Measures.Length(:mm, 10.0),
    left_margin=Plots.Measures.Length(:mm, 5.0),
    xlabel="date",
    xticks = ([1 floor(D/4)  floor(D/2) floor(3*D/4) D;], [l0 l1 l2 l3 l4]),
    yaxis="cases/1M",
    legendfont=font(10),
    label=["Japan" "Philippines" "Malaysia" "Indonesia" "India" "Nepal" "Pakistan" "Bangladesh" "Sri Lanka" "Thailand"],
    palette = :seaborn_bright,
    legend = :topleft)
savefig("./jhu/jhu_asia_new_cases.png")

p4=plot([NDJPN NDPHI NDMYS NDIDN NDIND NDNPL NDPAK NDBGD NDLKA NDTHA],
    grid=false,
    linewidth=2,
    title="COVID-19: 7-day average deaths per 1M \n data sourced by JHU Coronavirus Resource Center",
    right_margin=Plots.Measures.Length(:mm, 10.0),
    left_margin=Plots.Measures.Length(:mm, 5.0),
    xticks = ([1 floor(D/4)  floor(D/2) floor(3*D/4) D;], [l0 l1 l2 l3 l4]),
    xlabel="date",
    yaxis="deaths/1M",
    legendfont=font(10),
    label=["Japan" "Philippines" "Malaysia" "Indonesia" "India" "Nepal" "Pakistan" "Bangladesh" "Sri Lanka" "Thailand"],
   palette = :seaborn_bright,
     legend = :topleft)
savefig("./jhu/jhu_asia_recent_deaths.png")

p5=plot([AUSA AGBR AARG ABRA ACOL AMEX AFIE],
    grid=false,
    linewidth=2,
    title="COVID-19: cases per 1M \n data sourced by JHU Coronavirus Resource Center",
    right_margin=Plots.Measures.Length(:mm, 15.0),
    left_margin=Plots.Measures.Length(:mm, 5.0),
    xticks = ([1 floor(D/4)  floor(D/2) floor(3*D/4) D;], [l0 l1 l2 l3 l4]),
    xlabel="date",
    yaxis="cases/1M",
    legendfont=font(10),
    label=["United States" "United Kingdom" "Argentina" "Brazil" "Colombia" "Mexico" "FRA+ITA+ESP"],
    palette = :seaborn_bright,
    legend = :topleft)

p6=plot([BUSA BGBR BARG BBRA BCOL BMEX BFIE],
    grid=false,
    linewidth=2,
    title="COVID-19: deaths per 1M \n data sourced by JHU Coronavirus Resource Center",
    right_margin=Plots.Measures.Length(:mm, 10.0),
    left_margin=Plots.Measures.Length(:mm, 5.0),
    xticks = ([1 floor(D/4)  floor(D/2) floor(3*D/4) D;], [l0 l1 l2 l3 l4]),
    xlabel="date",
    yaxis="deaths/1M",
    legendfont=font(10),
    label=["United States" "United Kingdom" "Argentina" "Brazil" "Colombia" "Mexico" "FRA+ITA+ESP"],
    palette = :seaborn_bright,
    legend = :topleft)

p7=plot([NUSA NGBR NARG NBRA NCOL NMEX NFIE],
    grid=false,
    linewidth=2,
    title="COVID-19: 7-day-average of new cases per 1M \n data sourced by JHU Coronavirus Resource Center",
    right_margin=Plots.Measures.Length(:mm, 10.0),
    left_margin=Plots.Measures.Length(:mm, 5.0),
    xticks = ([1 floor(D/4)  floor(D/2) floor(3*D/4) D;], [l0 l1 l2 l3 l4]),
    xlabel="date",
    yaxis="cases/1M",
    legendfont=font(10),
    label=["United States" "United Kingdom" "Argentina" "Brazil" "Colombia" "Mexico" "FRA+ITA+ESP"],
    palette = :seaborn_bright,
    legend = :topleft)

p8=plot([NDUSA NDGBR NDARG NDBRA NDCOL NDMEX NDFIE],
    grid=false,
    linewidth=2,
    title="COVID-19: 7-day average deaths per 1M \n data sourced by JHU Coronavirus Resource Center",
    right_margin=Plots.Measures.Length(:mm, 10.0),
    left_margin=Plots.Measures.Length(:mm, 5.0),
    xticks = ([1 floor(D/4)  floor(D/2) floor(3*D/4) D;], [l0 l1 l2 l3 l4]),
    xlabel="date",
    yaxis="deaths/1M",
    legendfont=font(10),
    label=["United States" "United Kingdom" "Argentina" "Brazil" "Colombia" "Mexico" "FRA+ITA+ESP"],
   palette = :seaborn_bright,
     legend = :topleft)

plot(p1, p2, p3, p4,
     layout=(2,2),
     size=(1200,800),
     left_margin=Plots.Measures.Length(:mm, 5.0),
     right_margin=Plots.Measures.Length(:mm, 15.0),
     top_margin=Plots.Measures.Length(:mm, 5.0),
     bottom_margin=Plots.Measures.Length(:mm, 5.0))
savefig("./jhu/jhu_asia.png")

plot(p5, p6, p7, p8,
     layout=(2,2),
     size=(1200,800),
     left_margin=Plots.Measures.Length(:mm, 5.0),
     right_margin=Plots.Measures.Length(:mm, 15.0),
     top_margin=Plots.Measures.Length(:mm, 5.0),
     bottom_margin=Plots.Measures.Length(:mm, 5.0))
savefig("./jhu/jhu_world.png")


text:
using JLD
using BenchmarkTools
using UnicodePlots

res = JLD.load(ARGS[1])
serial_results = res["results"]["Serial"]
dagger_results = res["results"]["Dagger"]
pf = res["peakflops"]
renders = res["renders"]

theory_flops(nrow, ncol, nfeatures) = 11 * ncol * nrow * nfeatures + 2 * (ncol + nrow) * nfeatures

scale_set = sort([key=>parse(Int, lstrip(last(split(key, ':')), ' ')) for key in keys(serial_results)]; by=x->x[2])

expected_seconds = [theory_flops(2001*scale, 10002, 12)/pf for scale in last.(scale_set)]
plt = lineplot(last.(scale_set), expected_seconds; name="Speed of Light", title="NNMF Execution Times", ylabel="Best Time (s)", xlabel="Problem size", width=60, height=25, xlim=[minimum(last.(scale_set)), maximum(last.(scale_set))], ylim=[0, 20])

serial_times = [minimum(serial_results[scale]).time/(10^9) for scale in first.(scale_set)]
lineplot!(plt, last.(scale_set), serial_times; name="Serial")

nw_set = sort([key=>parse(Int, lstrip(last(split(key, ':')), ' ')) for key in keys(dagger_results[first(first(scale_set))])]; by=x->x[2])
for (nw,nw_val) in nw_set
    dagger_times = [minimum(dagger_results[scale][nw]).time/(10^9) for scale in first.(scale_set)]
    lineplot!(plt, last.(scale_set), dagger_times; name="Dagger $nw")
end
display(plt)
println()


text:
using Plots; pyplot();
using Distributions
using LinearAlgebra
using KernelFunctions
# Translational invariants kernels

default(lw=3.0,titlefontsize=28,tickfontsize=18)

x₀ = 0.0; l = 0.1
n_grid = 101
fill(x₀,n_grid,1)
xrange = reshape(collect(range(-3,3,length=n_grid)),:,1)

k = transform(SqExponentialKernel(),1.0)
K1 = kernelmatrix(k,xrange,obsdim=1)
p = heatmap(K1,yflip=true,colorbar=false,framestyle=:none,background_color=RGBA(0.0,0.0,0.0,0.0))
savefig(joinpath(@__DIR__,"src","assets","heatmap_sqexp.png"))

k = @kernel Matern32Kernel FunctionTransform(x->(sin.(x)).^2)
K2 = kernelmatrix(k,xrange,obsdim=1)
p = heatmap(K2,yflip=true,colorbar=false,framestyle=:none,background_color=RGBA(0.0,0.0,0.0,0.0))
savefig(joinpath(@__DIR__,"src","assets","heatmap_matern.png"))

k = transform(PolynomialKernel(c=0.0,d=2.0), LinearTransform(randn(3,1)))
K3 = kernelmatrix(k,xrange,obsdim=1)
p = heatmap(K3,yflip=true,colorbar=false,framestyle=:none,background_color=RGBA(0.0,0.0,0.0,0.0))
savefig(joinpath(@__DIR__,"src","assets","heatmap_poly.png"))

k = 0.5*SqExponentialKernel()*transform(LinearKernel(),0.5) + 0.4*(@kernel Matern32Kernel() FunctionTransform(x->sin.(x)))
K4 = kernelmatrix(k,xrange,obsdim=1)
p = heatmap(K4,yflip=true,colorbar=false,framestyle=:none,background_color=RGBA(0.0,0.0,0.0,0.0))
savefig(joinpath(@__DIR__,"src","assets","heatmap_prodsum.png"))

plot(heatmap.([K1,K2,K3,K4],yflip=true,colorbar=false)...,layout=(2,2))
savefig(joinpath(@__DIR__,"src","assets","heatmap_combination.png"))

##

for k in [SqExponentialKernel,ExponentialKernel]
    K = kernelmatrix(k(),xrange,obsdim=1)
    v = rand(MvNormal(K+1e-7I))
    plot(xrange,v,lab="",title="f(x)",framestyle=:none) |> display
    savefig(joinpath(@__DIR__,"src","assets","GP_sample_$(k).png"))
    plot(xrange,kernel.(k(),x₀,xrange),lab="",ylims=(0,1.1),title="k(0,x)") |> display
    savefig(joinpath(@__DIR__,"src","assets","kappa_function_$(k).png"))
end

for k in [GammaExponentialKernel(1.0,1.5)]
    sparse =1
    while !isposdef(kernelmatrix(k,xrange*sparse,obsdim=1) + 1e-5I); sparse += 1; end
    v = rand(MvNormal(kernelmatrix(k,xrange*sparse,obsdim=1)+1e-7I))
    plot(xrange,v,lab="",title="f(x)",framestyle=:none) |> display
    savefig(joinpath(@__DIR__,"src","assets","GP_sample_GammaExponentialKernel.png"))
    plot(xrange,kernel.(k,x₀,xrange),lab="",ylims=(0,1.1),title="k(0,x)") |> display
    savefig(joinpath(@__DIR__,"src","assets","kappa_function_GammaExponentialKernel.png"))
end

for k in [MaternKernel,Matern32Kernel,Matern52Kernel]
    K = kernelmatrix(k(),xrange,obsdim=1)
    v = rand(MvNormal(K+1e-7I))
    plot(xrange,v,lab="",title="f(x)",framestyle=:none) |> display
    savefig(joinpath(@__DIR__,"src","assets","GP_sample_$(k).png"))
    plot(xrange,kernel.(k(),x₀,xrange),lab="",ylims=(0,1.1),title="k(0,x)") |> display
    savefig(joinpath(@__DIR__,"src","assets","kappa_function_$(k).png"))
end

for k in [RationalQuadraticKernel]
    K = kernelmatrix(k(),xrange,obsdim=1)
    v = rand(MvNormal(K+1e-7I))
    plot(xrange,v,lab="",title="f(x)",framestyle=:none) |> display
    savefig(joinpath(@__DIR__,"src","assets","GP_sample_$(k).png"))
    plot(xrange,kernel.(k(),x₀,xrange),lab="",ylims=(0,1.1),title="k(0,x)") |> display
    savefig(joinpath(@__DIR__,"src","assets","kappa_function_$(k).png"))
end


text:
using MPIMeasurements
using Unitful

hr = MPIScanner("IselRobot.toml")
initRefZYX(hr)
moveAbs(hr,10.0Unitful.mm,0.0Unitful.mm,0.0Unitful.mm)
# Basic Test
# using LibSerialPort
# serialports = list_ports()
# s = open("/dev/ttyS0", 19200)
# baurate 19200, 8 Datenbit, 1 Stoppbit, keine Parität
# r1=write(s, "@07\r")
# r2=write(s, "@0R7\r")
# close(s)