text:
@info "Running interactive tests"

using PkgTemplates: @with_kw_noshow

const CR = "\r"
const LF = "\n"
const UP = "\eOA"
const DOWN = "\eOB"
const ALL = "a"
const NONE = "n"
const DONE = "d"
const SIGINT = "\x03"

# Because the plugin selection dialog prints directly to stdin in the same way
# as we do here, and our input prints happen first, we're going to need to insert
# the plugin selection prints ourselves, and then "undo" the extra ones at the end
# by consuming whatever is left in stdin.
const NDEFAULTS = length(PT.default_plugins())
const SELECT_DEFAULTS = (CR * DOWN)^NDEFAULTS

struct FromString
    s::String
end

@testset "Interactive mode" begin
    @testset "Input conversion" begin
        generic(T, x) = PT.convert_input(PT.Plugin, T, x)
        @test generic(String, "foo") == "foo"
        @test generic(Float64, "1.23") == 1.23
        @test generic(Int, "01") == 1
        @test generic(Bool, "yes") === true
        @test generic(Bool, "True") === true
        @test generic(Bool, "No") === false
        @test generic(Bool, "false") === false
        @test generic(Vector{Int}, "1, 2, 3") == [1, 2, 3]
        @test generic(Vector{String}, "a, b,c") == ["a", "b", "c"]
        @test generic(FromString, "hello") == FromString("hello")
        if VERSION < v"1.1"
            @test_broken generic(Union{String, Nothing}, "nothing") === nothing
        else
            @test generic(Union{String, Nothing}, "nothing") === nothing
        end

        @test_throws ArgumentError generic(Int, "hello")
        @test_throws ArgumentError generic(Float64, "hello")
        @test_throws ArgumentError generic(Bool, "hello")
    end

    @testset "input_tips" begin
        @test isempty(PT.input_tips(Int))
        @test PT.input_tips(Vector{String}) == ["comma-delimited"]
        @test PT.input_tips(Union{Vector{String}, Nothing}) ==
            ["'nothing' for nothing", "comma-delimited"]
        @test PT.input_tips(Union{String, Nothing}) == ["'nothing' for nothing"]
        @test PT.input_tips(Union{Vector{Secret}, Nothing}) ==
            ["'nothing' for nothing", "comma-delimited", "name only"]
    end

    @testset "Interactive name/type pair collection" begin
        name = gensym()
        @eval begin
            PT.@plugin struct $name <: PT.Plugin
                x::Int = 0
                y::String = ""
            end

            @test PT.interactive_pairs($name) == [:x => Int, :y => String]
            PT.customizable(::Type{$name}) = (:x => PT.NotCustomizable, :y => Float64, :z => Int)
            @test PT.interactive_pairs($name) == [:y => Float64, :z => Int]
        end
    end

    @testset "Simulated inputs" begin
        @testset "Default template" begin
            print(
                stdin.buffer,
                CR,        # Select user
                DONE,      # Finish menu
                USER, LF,  # Enter user
            )
            @test Template(; interactive=true) == Template(; user=USER)
        end

        @testset "Custom options, accept defaults" begin
            print(
                stdin.buffer,
                ALL, DONE,        # Customize all fields
                "user", LF,       # Enter user (don't assume we have default for this one).
                LF,               # Enter authors
                LF,               # Enter dir
                CR,               # Enter host
                CR,               # Enter julia
                SELECT_DEFAULTS,  # Pre-select default plugins
                DONE,             # Select no additional plugins
                DONE^NDEFAULTS,   # Don't customize plugins
            )
            @test Template(; interactive=true) == Template(; user="user")
            readavailable(stdin.buffer)
        end

        @testset "Custom options, custom values" begin
            nversions = VERSION.minor + 1
            print(
                stdin.buffer,
                ALL, DONE,           # Customize all fields
                "user", LF,          # Enter user
                "a, b", LF,          # Enter authors
                "~", LF,             # Enter dir
                DOWN^3, CR,          # Choose "Other" for host
                "x.com", LF,         # Enter host
                DOWN^nversions, CR,  # Choose "Other" for julia
                "0.7", LF,           # Enter Julia version
                SELECT_DEFAULTS,     # Pre-select default plugins
                DONE,                # Select no additional plugins
                DONE^NDEFAULTS,      # Don't customize plugins
            )
            @test Template(; interactive=true) == Template(;
                user="user",
                authors=["a", "b"],
                dir="~",
                host="x.com",
                julia=v"0.7",
            )
            readavailable(stdin.buffer)
        end

        @testset "Disabling default plugins" begin
            print(
                stdin.buffer,
                CR, DOWN^5, CR, DONE,  # Customize user and plugins
                USER, LF,                # Enter user
                SELECT_DEFAULTS,         # Pre-select default plugins
                UP, CR, UP^2, CR, DONE,  # Disable TagBot and Readme
                DONE^(NDEFAULTS - 2),    # Don't customize plugins
            )
            @test Template(; interactive=true) == Template(;
                user=USER,
                plugins=[!Readme, !TagBot],
            )
            readavailable(stdin.buffer)
        end

        @testset "Plugins" begin
            print(
                stdin.buffer,
                ALL, DONE,       # Customize all fields
                "true", LF,      # Enable ARM64
                "no", LF,        # Disable coverage
                "1.1,v1.2", LF,  # Enter extra versions
                "x.txt", LF,     # Enter file
                "Yes", LF,       # Enable Linux
                "false", LF,     # Disable OSX
                "TRUE", LF,      # Enable Windows
                "YES",  LF,      # Enable x64
                "NO", LF,        # Disable x86
            )
            @test PT.interactive(TravisCI) == TravisCI(;
                arm64=true,
                coverage=false,
                extra_versions=["1.1", "v1.2"],
                file="x.txt",
                linux=true,
                osx=false,
                windows=true,
                x64=true,
                x86=false,
            )

            print(
                stdin.buffer,
                DOWN^2, CR,      # Select GitLabCI
                DOWN^2, CR, DONE,  # Customize index_md
                "x.txt", LF,     # Enter index file
            )
            @test PT.interactive(Documenter) == Documenter{GitLabCI}(; index_md="x.txt")

            print(
                stdin.buffer,
                CR, DOWN, CR, DONE,  # Customize name and destination
                "COPYING", LF,       # Enter destination
                CR,                  # Choose MIT for name (it's at the top)
            )
            @test PT.interactive(License) == License(; destination="COPYING", name="MIT")
        end

        @testset "Quotes" begin
            print(
                stdin.buffer,
                CR, DOWN^2, CR, DONE,  # Customize user and dir
                "\"me\"", LF,          # Enter user with quotes
                "\"~\"", LF,           # Enter home dir with quotes
            )
            result = Template(; interactive=true) == Template(; user="me", dir="~")
            if get(ENV, "CI", "false") == "true"
                @test_broken result
            else
                @test result
            end

            print(
                stdin.buffer,
                DOWN^2, CR, DONE,                        # Customize extra_versions
                "\"1.1.1\", \"^1.5\", \"nightly\"", LF,  # Enter versions with quotes
            )
            @test PT.interactive(TravisCI) == TravisCI(;
                extra_versions=["1.1.1", "^1.5", "nightly"],
            )
        end

        @testset "Union{T, Nothing} weirdness" begin
            print(
                stdin.buffer,
                DOWN, CR, DONE,  # Customize changelog
                "hello", LF,     # Enter changelog
            )
            @test PT.interactive(TagBot) == TagBot(; changelog="hello")

            print(
                stdin.buffer,
                DOWN, CR, DONE,  # Customize changelog
                "nothing", LF,   # Set to null
            )
            @test PT.interactive(TagBot) == TagBot(; changelog=nothing)
        end

        @testset "Only one field" begin
            print(
                stdin.buffer,
                DOWN, CR, DONE,  # Select "None" option
            )
            @test PT.interactive(Codecov) == Codecov()
        end

        @testset "Missing user" begin
            print(
                stdin.buffer,
                DONE,            # Customize nothing
                "username", LF,  # Enter user after it's prompted
            )

            patch = @patch PkgTemplates.default_user() = ""
            apply(patch) do
                @test Template(; interactive=true) == Template(; user="username")
            end
        end

        @testset "Interrupts" begin
            print(
                stdin.buffer,
                SIGINT,        # Send keyboard interrupt
            )
            @test Template(; interactive=true) === nothing
        end

        println()
    end
end


text:
import Base: size, length, show, iterate, lastindex, isempty, convert, ndims, float, eltype, copy, ==
if VERSION >= v"0.7-"
    using Dates
else
    using Base.Dates
end

################################################################################
# TYPE DEFINITION ##############################################################
################################################################################

abstract type AbstractTS end

"""
Time series type aimed at efficiency and simplicity.

Motivated by the `xts` package in R and the `pandas` package in Python.
"""
mutable struct TS{V<:Real,T<:TimeType}
    values::Matrix{V}
    index::Vector{T}
    fields::Vector{Symbol}
    function TS{V,T}(values::AbstractArray{V}, index::AbstractVector{T}, fields::Vector{Symbol}) where {V<:Real,T<:TimeType}
        @assert size(values,1)==length(index) "Length of index not equal to number of value rows."
        @assert size(values,2)==length(fields) "Length of fields not equal to number of columns in values."
        order = sortperm(index)
        return new(values[order,:], index[order], SANITIZE_NAMES ? namefix.(fields) : fields)
    end
end

TS(v::AbstractArray{V}, t::AbstractVector{T}, f::Union{Symbol,String,Char}) where {V,T} = TS{V,T}(v, t, [Symbol(f)])
TS(v::AbstractArray{V}, t::AbstractVector{T}, f::Union{AbstractVector{Symbol},AbstractVector{String},AbstractVector{Char}}) where {V,T} = TS{V,T}(v, t, Symbol.(f))
TS(v::AbstractArray{V}, t::AbstractVector{T}) where {V,T} = TS{V,T}(v, t, autocol(1:size(v,2)))
TS(v::AbstractArray{V}, t::T, f) where {V,T} = TS{V,T}(v, [t], f)
TS(v::V, t::AbstractVector{T}, f) where {V,T} = TS{V,T}([v], t, f)
TS(v::V, t::T, f) where {V,T} = TS{V,T}([v][:,:], [t], f)
TS(v::V, t::T) where {V,T} = TS{V,T}([v], [t], [:A])
TS(v::AbstractArray{V}) where {V} = TS{V,Date}(v, autoidx(size(v,1)), autocol(1:size(v,2)))
TS() = TS{Float64,Date}(Matrix{Float64}(UndefInitializer(),0,0), Date[], Symbol[])

# Conversions ------------------------------------------------------------------
convert(::Type{TS{Float64}}, x::TS{Bool}) = TS{Float64}(map(Float64, x.values), x.index, x.fields)
convert(::Type{TS{Int}}, x::TS{Bool}) = TS{Int}(map(Int, x.values), x.index, x.fields)
convert(::Type{TS{Bool}}, x::TS{V}) where {V<:Real} = TS{Bool}(map(V, x.values), x.index, x.fields)
convert(x::TS{Bool}) = convert(TS{Int}, x::TS{Bool})
# convert{V}(::Type{TS}, x::Array{V}) = TS{V,Date}(x, [Dates.Date() for i in 1:size(x,1)])
# convert(x::TS{Bool}) = convert(TS{Float64}, x::TS{Bool})
const ts = TS

################################################################################
# BASIC UTILITIES ##############################################################
################################################################################
size(x::TS) = size(x.values)
size(x::TS, dim::Int) = size(x.values, dim)
length(x::TS) = prod(size(x))::Int
lastindex(x::TS) = lastindex(x.values)
lastindex(x::TS, d) = lastindex(x.values, d)
iterate(x::TS) = size(x,1) == 0 ? nothing : (x.index[1], x.values[1,:]), 2
iterate(x::TS, i::Int) = i == lastindex(x, 1) ? nothing : ((x.index[i], x.values[i,:]), i+1)
isempty(x::TS) = (isempty(x.index) && isempty(x.values))
first(x::TS) = x[1]
last(x::TS) = x[end]
ndims(::TS) = 2
float(x::TS) = ts(float(x.values), x.index, x.fields)
eltype(x::TS) = eltype(x.values)
copy(x::TS) = TS(x.values, x.index, x.fields)

ext: jl


text:
# This file is part of the TaylorIntegration.jl package; MIT licensed

# jetcoeffs!
@doc doc"""
    jetcoeffs!(eqsdiff::Function, t, x, params)

Returns an updated `x` using the recursion relation of the
derivatives obtained from the differential equations
``\dot{x}=dx/dt=f(x, p, t)``.

`eqsdiff` is the function defining the RHS of the ODE,
`x` contains the Taylor1 expansion of the dependent variable(s) and
`t` is the independent variable, and `params` are the parameters appearing on the
function defining the differential equation. See [`taylorinteg`](@ref) for examples
and convention for `eqsdiff`. Note that `x` is of type `Taylor1{U}` where
`U<:Number`; `t` is of type `Taylor1{T}` where `T<:Real`.

Initially, `x` contains only the 0-th order Taylor coefficient of
the current system state (the initial conditions), and `jetcoeffs!`
computes recursively the high-order derivates back into `x`.

"""
function jetcoeffs!(eqsdiff::Function, t::Taylor1{T}, x::Taylor1{U}, params) where
        {T<:Real, U<:Number}
    order = x.order
    for ord in 0:order-1
        ordnext = ord+1

        # Set `taux`, `xaux`, auxiliary Taylor1 variables to order `ord`
        @inbounds taux = Taylor1( t.coeffs[1:ordnext] )
        @inbounds xaux = Taylor1( x.coeffs[1:ordnext] )

        # Equations of motion
        dx = eqsdiff(xaux, params, taux)

        # Recursion relation
        @inbounds x[ordnext] = dx[ord]/ordnext
    end
    nothing
end

@doc doc"""
    jetcoeffs!(eqsdiff!::Function, t, x, dx, xaux, params)

Mutates `x` in-place using the recursion relation of the
derivatives obtained from the differential equations
``\dot{x}=dx/dt=f(x, p, t)``.

`eqsdiff!` is the function defining the RHS of the ODE,
`x` contains the Taylor1 expansion of the dependent variables and
`t` is the independent variable, and `params` are the parameters appearing on the
function defining the differential equation. See [`taylorinteg`](@ref) for examples
and convention for `eqsdiff`. Note that `x` is of type `Vector{Taylor1{U}}`
where `U<:Number`; `t` is of type `Taylor1{T}` where `T<:Real`. In this case,
two auxiliary containers `dx` and `xaux` (both of the same type as `x`) are
needed to avoid allocations.

Initially, `x` contains only the 0-th order Taylor coefficient of
the current system state (the initial conditions), and `jetcoeffs!`
computes recursively the high-order derivates back into `x`.

"""
function jetcoeffs!(eqsdiff!::Function, t::Taylor1{T},
        x::AbstractArray{Taylor1{U}, N}, dx::AbstractArray{Taylor1{U}, N},
        xaux::AbstractArray{Taylor1{U}, N}, params) where {T<:Real, U<:Number, N}
    order = x[1].order
    for ord in 0:order-1
        ordnext = ord+1

        # Set `taux`, auxiliary Taylor1 variable to order `ord`
        @inbounds taux = Taylor1( t.coeffs[1:ordnext] )
        # Set `xaux`, auxiliary vector of Taylor1 to order `ord`
        for j in eachindex(x)
            @inbounds xaux[j] = Taylor1( x[j].coeffs[1:ordnext] )
        end

        # Equations of motion
        eqsdiff!(dx, xaux, params, taux)

        # Recursion relations
        for j in eachindex(x)
            @inbounds x[j][ordnext] = dx[j][ord]/ordnext
        end
    end
    nothing
end

# __jetcoeffs
"""
    __jetcoeffs!(::Val{bool}, f, t, x, params, tmpTaylor, arrTaylor)
    __jetcoeffs!(::Val{bool}, f, t, x, dx, xaux, params, tmpTaylor, arrTaylor)

Chooses a method of [`jetcoeffs!`](@ref) (hard-coded) or the generated by
[`@taylorize`](@ref)) depending on `Val{bool}` (`bool::Bool`).
"""
@inline __jetcoeffs!(::Val{false}, f, t, x, params, tmpTaylor, arrTaylor) =
    jetcoeffs!(f, t, x, params)
@inline __jetcoeffs!(::Val{true},  f, t, x, params, tmpTaylor, arrTaylor) =
    jetcoeffs!(Val(f), t, x, params, tmpTaylor, arrTaylor)
@inline __jetcoeffs!(::Val{false}, f, t, x, dx, xaux, params, tmpTaylor, arrTaylor) =
    jetcoeffs!(f, t, x, dx, xaux, params)
@inline __jetcoeffs!(::Val{true},  f, t, x, dx, xaux, params, tmpTaylor, arrTaylor) =
    jetcoeffs!(Val(f), t, x, dx, params, tmpTaylor, arrTaylor)

# _allocate_jetcoeffs! fallbacks; needed to define the corresponding parsed method
_allocate_jetcoeffs!(::Taylor1{T}, x::Taylor1{S}, params) where {T,S} =
    (typeof(x)[], Vector{typeof(x)}[])
_allocate_jetcoeffs!(::Taylor1{T}, x::AbstractArray{Taylor1{S}, N},
        ::AbstractArray{Taylor1{S}, N}, params) where {T,S,N} =
    (eltype(x)[], Vector{eltype(x)}[])

# stepsize
"""
    stepsize(x, epsilon) -> h

Returns a maximum time-step for a the Taylor expansion `x`
using a prescribed absolute tolerance `epsilon` and the last two
Taylor coefficients of (each component of) `x`.

Note that `x` is of type `Taylor1{U}` or `Vector{Taylor1{U}}`, including
also the cases `Taylor1{TaylorN{U}}` and `Vector{Taylor1{TaylorN{U}}}`.

Depending of `eltype(x)`, i.e., `U<:Number`, it may be necessary to overload
`stepsize`, specializing it on the type `U`, to avoid type instabilities.
"""
function stepsize(x::Taylor1{U}, epsilon::T) where {T<:Real, U<:Number}
    R = promote_type(typeof(norm(constant_term(x), Inf)), T)
    ord = x.order
    h = convert(R, Inf)
    z = zero(R)
    for k in (ord-1, ord)
        @inbounds aux = norm( x[k], Inf)
        aux == z && continue
        aux1 = _stepsize(aux, epsilon, k)
        h = min(h, aux1)
    end
    return h::R
end

function stepsize(q::AbstractArray{Taylor1{U}, N}, epsilon::T) where
        {T<:Real, U<:Number, N}
    R = promote_type(typeof(norm(constant_term(q[1]), Inf)), T)
    h = convert(R, Inf)
    for i in eachindex(q)
        @inbounds hi = stepsize( q[i], epsilon )
        h = min( h, hi )
    end

    # If `isinf(h)==true`, we use the maximum (finite)
    # step-size obtained from all coefficients as above.
    # Note that the time step is independent from `epsilon`.
    if isinf(h)
        h = zero(R)
        for i in eachindex(q)
            @inbounds hi = _second_stepsize(q[i], epsilon)
            h = max( h, hi )
        end
    end
    return h::R
end

"""
    _stepsize(aux1, epsilon, k)

Helper function to avoid code repetition.
Returns ``(epsilon/aux1)^(1/k)``.
"""
@inline function _stepsize(aux1::U, epsilon::T, k::Int) where {T<:Real, U<:Number}
    aux = epsilon / aux1
    kinv = 1 / k
    return aux^kinv
end

"""
    _second_stepsize(x, epsilon)

Corresponds to the "second stepsize control" in Jorba and Zou
(2005) paper. We use it if [`stepsize`](@ref) returns `Inf`.
"""
function _second_stepsize(x::Taylor1{U}, epsilon::T) where {T<:Real, U<:Number}
    R = promote_type(typeof(norm(constant_term(x), Inf)), T)
    x == zero(x) && return convert(R, Inf)
    ord = x.order
    z = zero(R)
    u = one(R)
    h = z
    for k in 1:ord-2
        @inbounds aux = norm( x[k], Inf)
        aux == z && continue
        aux1 = _stepsize(aux, u, k)
        h = max(h, aux1)
    end
    return h::R
end

#taylorstep
@doc doc"""
    taylorstep!(f, t, x, t0, order, abstol, params, tmpTaylor, arrTaylor, parse_eqs=true) -> δt
    taylorstep!(f!, t, x, dx, xaux, t0, order, abstol, params, tmpTaylor, arrTaylor, parse_eqs=true) -> δt

One-step Taylor integration for the one-dependent variable ODE ``\dot{x}=dx/dt=f(x, p, t)``
with initial conditions ``x(t_0)=x_0``.
Returns the time-step `δt` of the actual integration carried out (δt is positive).

Here, `f` is the function defining the RHS of the ODE (see
[`taylorinteg`](@ref)), `t` is the
independent variable, `x` contains the Taylor expansion of the dependent
variable, `order` is
the degree  used for the `Taylor1` polynomials during the integration
`abstol` is the absolute tolerance used to determine the time step
of the integration, and `params` are the parameters entering the ODE
functions.
For several variables, `dx` and `xaux`, both of the same type as `x`,
are needed to save allocations. Finally, `parse_eqs` is a switch
to force *not* using (`parse_eqs=false`) the specialized method of `jetcoeffs!`
created with [`@taylorize`](@ref); the default is `true` (parse the equations).
Finally, `parse_eqs` is a switch
to force *not* using (`parse_eqs=false`) the specialized method of `jetcoeffs!`
created with [`@taylorize`](@ref); the default is `true` (parse the equations).

"""
function taylorstep!(f, t::Taylor1{T}, x::Taylor1{U}, abstol::T, params,
        tmpTaylor, arrTaylor, parse_eqs::Bool=true) where {T<:Real, U<:Number}

    # Compute the Taylor coefficients
    __jetcoeffs!(Val(parse_eqs), f, t, x, params, tmpTaylor, arrTaylor)

    # Compute the step-size of the integration using `abstol`
    δt = stepsize(x, abstol)
    if isinf(δt)
        δt = _second_stepsize(x, abstol)
    end

    return δt
end

function taylorstep!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}},
        dx::Vector{Taylor1{U}}, xaux::Vector{Taylor1{U}}, abstol::T, params,
        tmpTaylor, arrTaylor, parse_eqs::Bool=true) where {T<:Real, U<:Number}

    # Compute the Taylor coefficients
    __jetcoeffs!(Val(parse_eqs), f!, t, x, dx, xaux, params, tmpTaylor, arrTaylor)

    # Compute the step-size of the integration using `abstol`
    δt = stepsize(x, abstol)

    return δt
end

"""
    _determine_parsing!(parse_eqs::Bool, f, t, x, params)
    _determine_parsing!(parse_eqs::Bool, f, t, x, dx, params)

Check if the parsed method of `jetcoeffs!` exists and check it
runs without error.
"""
function _determine_parsing!(parse_eqs::Bool, f, t, x, params)
    parse_eqs = parse_eqs && !isempty(methodswith(Val{f}, jetcoeffs!)) &&
                             !isempty(methodswith(Val{f}, _allocate_jetcoeffs!))

    tmpTaylor, arrTaylor = _allocate_jetcoeffs!(t, x, params)

    if parse_eqs
        try
            tmpTaylor, arrTaylor = _allocate_jetcoeffs!(Val(f), t, x, params)
            jetcoeffs!(Val(f), t, x, params, tmpTaylor, arrTaylor)
        catch
            @warn("""Unable to use the parsed method of `jetcoeffs!` for `$f`,
            despite of having `parse_eqs=true`, due to some internal error.
            Using `parse_eqs = false`.""")
            parse_eqs = false
        end
    end

    return parse_eqs, tmpTaylor, arrTaylor
end
function _determine_parsing!(parse_eqs::Bool, f, t, x, dx, params)
    parse_eqs = parse_eqs && !isempty(methodswith(Val{f}, jetcoeffs!)) &&
                             !isempty(methodswith(Val{f}, _allocate_jetcoeffs!))

    tmpTaylor, arrTaylor = _allocate_jetcoeffs!(t, x, dx, params)

    if parse_eqs
        try
            tmpTaylor, arrTaylor = _allocate_jetcoeffs!(Val(f), t, x, dx, params)
            jetcoeffs!(Val(f), t, x, dx, params, tmpTaylor, arrTaylor)
        catch
            @warn("""Unable to use the parsed method of `jetcoeffs!` for `$f`,
            despite of having `parse_eqs=true`, due to some internal error.
            Using `parse_eqs = false`.""")
            parse_eqs = false
        end
    end

    return parse_eqs, tmpTaylor, arrTaylor
end

# taylorinteg
@doc doc"""
    taylorinteg(f, x0, t0, tmax, order, abstol, params[=nothing]; kwargs... )
    taylorinteg(f, x0, t0, tmax, order, abstol, Val(true), params[=nothing]; kwargs... )

General-purpose Taylor integrator for the explicit ODE ``\dot{x}=f(x, p, t)``.
The initial condition are specified by `x0` at time `t0`, and any parameters
encoded in `params`. The initial condition `x0` may be of type `T<:Number`
or `Vector{T}`, with `T` including `TaylorN{T}`; the latter case
is of interest for [jet transport applications](@ref jettransport).

The equations of motion are specified by the function `f`; we follow the same
convention of `DifferentialEquations.jl` to define this function, i.e.,
`f(x, p, t)` or `f!(dx, x, p, t)`; see the examples below.

It returns a vector with the values of time (independent variable),
and a vector with the computed values of
the dependent variable(s), and if the method used involves `Val(true)` it also
outputs the Taylor polynomial solutions obtained at each time step.
The integration stops when time is larger than `tmax`, in which case the last returned
value(s) correspond to that time, or when the number of saved steps is larger
than `maxsteps`.

The integration method uses polynomial expansions on the independent variable
of order `order`; the parameter `abstol` serves to define the
time step using the last two Taylor coefficients of the expansions.
Make sure you use a *large enough* `order` to assure convergence.

Currently, the recognized keyword arguments are:
- `maxsteps[=500]`: maximum number of integration steps.
- `parse_eqs[=true]`: use the specialized method of `jetcoeffs!` created
    with [`@taylorize`](@ref).

## Examples

For one dependent variable the function `f` defines the RHS of the equation of
motion, returning the value of ``\dot{x}``. The arguments of
this function are `(x, p, t)`, where `x` are the dependent variables, `p` are
the paremeters and `t` is the independent variable.

For several (two or more) dependent variables, the function `f!` defines
the RHS of the equations of motion, mutating (in-place) the (preallocated) vector
with components of ``\dot{x}``. The arguments of this function are `(dx, x, p, t)`,
where `dx` is the preallocated vector of ``\dot{x}``, `x` are the dependent
variables, `p` are the paremeters entering the ODEs and `t` is the independent
variable. The function may return this vector or simply `nothing`.

```julia
using TaylorIntegration

f(x, p, t) = x^2

tv, xv = taylorinteg(f, 3, 0.0, 0.3, 25, 1.0e-20, maxsteps=100 )

function f!(dx, x, p, t)
    for i in eachindex(x)
        dx[i] = x[i]^2
    end
    return nothing
end

tv, xv = taylorinteg(f!, [3, 3], 0.0, 0.3, 25, 1.0e-20, maxsteps=100 )

tv, xv, polynv = taylorinteg(f!, [3, 3], 0.0, 0.3, 25, 1.0e-20, maxsteps=100, Val(true) )
```

"""
function taylorinteg(f, x0::U, t0::T, tmax::T, order::Int, abstol::T,
    params = nothing; maxsteps::Int=500, parse_eqs::Bool=true) where {T<:Real, U<:Number}

    # Allocation
    tv = Array{T}(undef, maxsteps+1)
    xv = Array{U}(undef, maxsteps+1)

    # Initialize the Taylor1 expansions
    t = Taylor1( T, order )
    x = Taylor1( x0, order )

    # Initial conditions
    nsteps = 1
    @inbounds t[0] = t0
    @inbounds tv[1] = t0
    @inbounds xv[1] = x0
    sign_tstep = copysign(1, tmax-t0)

    # Determine if specialized jetcoeffs! method exists
    parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f, t, x, params)

    # Integration
    while sign_tstep*t0 < sign_tstep*tmax
        δt = taylorstep!(f, t, x, abstol, params, tmpTaylor, arrTaylor, parse_eqs) # δt is positive!
        # Below, δt has the proper sign according to the direction of the integration
        δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
        x0 = evaluate(x, δt) # new initial condition
        @inbounds x[0] = x0
        t0 += δt
        @inbounds t[0] = t0
        nsteps += 1
        @inbounds tv[nsteps] = t0
        @inbounds xv[nsteps] = x0
        if nsteps > maxsteps
            @warn("""
            Maximum number of integration steps reached; exiting.
            """)
            break
        end
    end

    #return tv, xv
    return view(tv,1:nsteps), view(xv,1:nsteps)
end

function taylorinteg(f, x0::U, t0::T, tmax::T, order::Int, abstol::T,
        ::Val{true}, params = nothing;
        maxsteps::Int=500, parse_eqs::Bool=true) where {T<:Real, U<:Number}

    # Allocation
    tv = Array{T}(undef, maxsteps+1)
    xv = Array{U}(undef, maxsteps+1)
    polynV = Array{Taylor1{U}}(undef, maxsteps+1)

    # Initialize the Taylor1 expansions
    t = Taylor1( T, order )
    x = Taylor1( x0, order )

    # Initial conditions
    nsteps = 1
    @inbounds t[0] = t0
    @inbounds tv[1] = t0
    @inbounds xv[1] = x0
    @inbounds polynV[1] = deepcopy(x)
    sign_tstep = copysign(1, tmax-t0)

    # Determine if specialized jetcoeffs! method exists
    parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f, t, x, params)

    # Integration
    while sign_tstep*t0 < sign_tstep*tmax
        δt = taylorstep!(f, t, x, abstol, params, tmpTaylor, arrTaylor, parse_eqs) # δt is positive!
        # Below, δt has the proper sign according to the direction of the integration
        δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
        x0 = evaluate(x, δt) # new initial condition
        @inbounds polynV[nsteps+1] = deepcopy(x)
        @inbounds x[0] = x0
        t0 += δt
        @inbounds t[0] = t0
        nsteps += 1
        @inbounds tv[nsteps] = t0
        @inbounds xv[nsteps] = x0
        if nsteps > maxsteps
            @warn("""
            Maximum number of integration steps reached; exiting.
            """)
            break
        end
    end

    return view(tv,1:nsteps), view(xv,1:nsteps), view(polynV, 1:nsteps)
end

function taylorinteg(f!, q0::Array{U,1}, t0::T, tmax::T, order::Int, abstol::T,
        params = nothing; maxsteps::Int=500, parse_eqs::Bool=true) where {T<:Real, U<:Number}

    # Allocation
    tv = Array{T}(undef, maxsteps+1)
    dof = length(q0)
    xv = Array{U}(undef, dof, maxsteps+1)

    # Initialize the vector of Taylor1 expansions
    t = Taylor1(T, order)
    x = Array{Taylor1{U}}(undef, dof)
    dx = Array{Taylor1{U}}(undef, dof)
    xaux = Array{Taylor1{U}}(undef, dof)
    for i in eachindex(q0)
        @inbounds x[i] = Taylor1( q0[i], order )
        @inbounds dx[i] = Taylor1( zero(q0[i]), order )
    end

    # Initial conditions
    @inbounds t[0] = t0
    x .= Taylor1.(q0, order)
    x0 = deepcopy(q0)
    @inbounds tv[1] = t0
    @inbounds xv[:,1] .= q0
    sign_tstep = copysign(1, tmax-t0)

    # Determine if specialized jetcoeffs! method exists
    parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f!, t, x, dx, params)

    # Integration
    nsteps = 1
    while sign_tstep*t0 < sign_tstep*tmax
        δt = taylorstep!(f!, t, x, dx, xaux, abstol, params, tmpTaylor, arrTaylor, parse_eqs) # δt is positive!
        # Below, δt has the proper sign according to the direction of the integration
        δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
        evaluate!(x, δt, x0) # new initial condition
        for i in eachindex(x0)
            @inbounds x[i][0] = x0[i]
            @inbounds dx[i] = Taylor1( zero(x0[i]), order )
        end
        t0 += δt
        @inbounds t[0] = t0
        nsteps += 1
        @inbounds tv[nsteps] = t0
        @inbounds xv[:,nsteps] .= x0
        if nsteps > maxsteps
            @warn("""
            Maximum number of integration steps reached; exiting.
            """)
            break
        end
    end

    return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:)
end

# Method to also return (output) the Taylor polynomials of the integration
function taylorinteg(f!, q0::Array{U,1}, t0::T, tmax::T, order::Int, abstol::T,
        ::Val{true}, params = nothing;
        maxsteps::Int=500, parse_eqs::Bool=true) where {T<:Real, U<:Number}

    # Allocation
    tv = Array{T}(undef, maxsteps+1)
    dof = length(q0)
    xv = Array{U}(undef, dof, maxsteps+1)
    polynV = Array{Taylor1{U}}(undef, dof, maxsteps+1)

    # Initialize the vector of Taylor1 expansions
    t = Taylor1(T, order)
    x = Array{Taylor1{U}}(undef, dof)
    dx = Array{Taylor1{U}}(undef, dof)
    xaux = Array{Taylor1{U}}(undef, dof)
    for i in eachindex(q0)
        @inbounds x[i] = Taylor1( q0[i], order )
        @inbounds dx[i] = Taylor1( zero(q0[i]), order )
    end

    # Initial conditions
    @inbounds t[0] = t0
    x .= Taylor1.(q0, order)
    x0 = deepcopy(q0)
    @inbounds tv[1] = t0
    @inbounds xv[:,1] .= q0
    @inbounds polynV[:,1] .= deepcopy.(x)
    sign_tstep = copysign(1, tmax-t0)

    # Determine if specialized jetcoeffs! method exists
    parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f!, t, x, dx, params)

    # Integration
    nsteps = 1
    while sign_tstep*t0 < sign_tstep*tmax
        δt = taylorstep!(f!, t, x, dx, xaux, abstol, params, tmpTaylor, arrTaylor, parse_eqs) # δt is positive!
        # Below, δt has the proper sign according to the direction of the integration
        δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
        evaluate!(x, δt, x0) # new initial condition
        # Store the Taylor polynomial solution
        @inbounds polynV[:,nsteps+1] .= deepcopy.(x)
        for i in eachindex(x0)
            @inbounds x[i][0] = x0[i]
            @inbounds dx[i] = Taylor1( zero(x0[i]), order )
        end
        t0 += δt
        @inbounds t[0] = t0
        nsteps += 1
        @inbounds tv[nsteps] = t0
        @inbounds xv[:,nsteps] .= x0
        if nsteps > maxsteps
            @warn("""
            Maximum number of integration steps reached; exiting.
            """)
            break
        end
    end

    return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:),
        view(transpose(view(polynV, :, 1:nsteps)), 1:nsteps, :)
end

# Integrate and return results evaluated at given time
@doc doc"""
    taylorinteg(f, x0, trange, order, abstol, params[=nothing]; keyword... )

General-purpose Taylor integrator for the explicit ODE
``\dot{x}=f(t,x)`` with initial condition specified by `x0::{T<:Number}`
or `x0::Vector{T}` at time `t0`.

The method returns a vector with the computed values of
the dependent variable(s), evaluated *only* at the times specified by
the range `trange`.

## Examples

```julia
xv = taylorinteg(f, 3, 0.0:0.001:0.3, 25, 1.0e-20, maxsteps=100 )

xv = taylorinteg(f!, [3, 3], 0.0:0.001:0.3, 25, 1.0e-20, maxsteps=100 );

```

"""
function taylorinteg(f, x0::U, trange::AbstractVector{T},
        order::Int, abstol::T, params = nothing;
        maxsteps::Int=500, parse_eqs::Bool=true) where {T<:Real, U<:Number}

    # Check if trange is increasingly or decreasingly sorted
    @assert (issorted(trange) ||
        issorted(reverse(trange))) "`trange` or `reverse(trange)` must be sorted"

    # Allocation
    nn = length(trange)
    xv = Array{U}(undef, nn)
    fill!(xv, T(NaN))

    # Initialize the Taylor1 expansions
    t = Taylor1( T, order )
    x = Taylor1( x0, order )

    # Initial conditions
    @inbounds t0, t1, tmax = trange[1], trange[2], trange[end]
    sign_tstep = copysign(1, tmax-t0)
    @inbounds t[0] = t0
    @inbounds xv[1] = x0

    # Determine if specialized jetcoeffs! method exists
    parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f, t, x, params)

    # Integration
    iter = 2
    nsteps = 1
    while sign_tstep*t0 < sign_tstep*tmax
        δt = taylorstep!(f, t, x, abstol, params, tmpTaylor, arrTaylor, parse_eqs)# δt is positive!
        # Below, δt has the proper sign according to the direction of the integration
        δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
        x0 = evaluate(x, δt) # new initial condition
        tnext = t0+δt
        # Evaluate solution at times within convergence radius
        while sign_tstep*t1 < sign_tstep*tnext
            x1 = evaluate(x, t1-t0)
            @inbounds xv[iter] = x1
            iter += 1
            @inbounds t1 = trange[iter]
        end
        if δt == tmax-t0
            @inbounds xv[iter] = x0
            break
        end
        @inbounds x[0] = x0
        t0 = tnext
        @inbounds t[0] = t0
        nsteps += 1
        if nsteps > maxsteps
            @warn("""
            Maximum number of integration steps reached; exiting.
            """)
            break
        end
    end
    return xv
end

function taylorinteg(f!, q0::Array{U,1}, trange::AbstractVector{T},
        order::Int, abstol::T, params = nothing; maxsteps::Int=500,
        parse_eqs::Bool=true) where {T<:Real, U<:Number}

    # Check if trange is increasingly or decreasingly sorted
    @assert (issorted(trange) ||
        issorted(reverse(trange))) "`trange` or `reverse(trange)` must be sorted"

    # Allocation
    nn = length(trange)
    dof = length(q0)
    x0 = similar(q0, eltype(q0), dof)
    x1 = similar(x0)
    fill!(x0, T(NaN))
    xv = Array{eltype(q0)}(undef, dof, nn)
    for ind in 1:nn
        @inbounds xv[:,ind] .= x0
    end

    # Initialize the vector of Taylor1 expansions
    t = Taylor1( T, order )
    x = Array{Taylor1{U}}(undef, dof)
    dx = Array{Taylor1{U}}(undef, dof)
    xaux = Array{Taylor1{U}}(undef, dof)
    for i in eachindex(q0)
        @inbounds x[i] = Taylor1( q0[i], order )
        @inbounds dx[i] = Taylor1( zero(q0[i]), order )
    end

    # Initial conditions
    @inbounds t[0] = trange[1]
    @inbounds t0, t1, tmax = trange[1], trange[2], trange[end]
    sign_tstep = copysign(1, tmax-t0)
    x .= Taylor1.(q0, order)
    @inbounds x0 .= q0
    @inbounds xv[:,1] .= q0

    # Determine if specialized jetcoeffs! method exists
    parse_eqs, tmpTaylor, arrTaylor = _determine_parsing!(parse_eqs, f!, t, x, dx, params)

    # Integration
    iter = 2
    nsteps = 1
    while sign_tstep*t0 < sign_tstep*tmax
        δt = taylorstep!(f!, t, x, dx, xaux, abstol, params, tmpTaylor, arrTaylor, parse_eqs) # δt is positive!
        # Below, δt has the proper sign according to the direction of the integration
        δt = sign_tstep * min(δt, sign_tstep*(tmax-t0))
        evaluate!(x, δt, x0) # new initial condition
        tnext = t0+δt
        # Evaluate solution at times within convergence radius
        while sign_tstep*t1 < sign_tstep*tnext
            evaluate!(x, t1-t0, x1)
            @inbounds xv[:,iter] .= x1
            iter += 1
            @inbounds t1 = trange[iter]
        end
        if δt == tmax-t0
            @inbounds xv[:,iter] .= x0
            break
        end
        for i in eachindex(x0)
            @inbounds x[i][0] = x0[i]
            @inbounds dx[i] = Taylor1( zero(x0[i]), order )
        end
        t0 = tnext
        @inbounds t[0] = t0
        nsteps += 1
        if nsteps > maxsteps
            @warn("""
            Maximum number of integration steps reached; exiting.
            """)
            break
        end
    end

    return transpose(xv)
end

# Generic functions
for R in (:Number, :Integer)
    @eval begin
        function taylorinteg(f, xx0::S, tt0::T, ttmax::U, order::Int, aabstol::V,
                params = nothing; maxsteps::Int=500, parse_eqs::Bool=true) where
                    {S<:$R, T<:Real, U<:Real, V<:Real}

            # In order to handle mixed input types, we promote types before integrating:
            t0, tmax, abstol, bfloat = promote(tt0, ttmax, aabstol, one(Float64))
            x0, tfloat1 = promote(xx0, t0)

            taylorinteg(f, x0, t0, tmax, order, abstol, params, maxsteps=maxsteps,
                parse_eqs=parse_eqs)
        end

        function taylorinteg(f, q0::Array{S,1}, tt0::T, ttmax::U, order::Int, aabstol::V,
                params = nothing; maxsteps::Int=500, parse_eqs::Bool=true) where
                    {S<:$R, T<:Real, U<:Real, V<:Real}

            #promote to common type before integrating:
            t0, tmax, abstol, afloat = promote(tt0, ttmax, aabstol, one(Float64))
            elq0, tt0 = promote(q0[1], t0)
            #convert the elements of q0 to the common, promoted type:
            q0_ = convert(Array{typeof(elq0)}, q0)

            taylorinteg(f, q0_, t0, tmax, order, abstol, params, maxsteps=maxsteps,
                parse_eqs=parse_eqs)
        end

        function taylorinteg(f, xx0::S, trange::AbstractVector{T},
                order::Int, aabstol::U, params = nothing;
                maxsteps::Int=500, parse_eqs::Bool=true) where
                    {S<:$R, T<:Real, U<:Real}
                    #
                t0, abstol, bfloat = promote(trange[1], aabstol, one(Float64))
                x0, tfloat1 = promote(xx0, t0)
                taylorinteg(f, x0, vec(trange.*one(t0)), order, abstol,
                    params, maxsteps=maxsteps, parse_eqs=parse_eqs)
        end

        function taylorinteg(f, q0::Array{S,1}, trange::AbstractVector{T},
                order::Int, aabstol::U, params = nothing;
                maxsteps::Int=500, parse_eqs::Bool=true) where
                    {S<:$R, T<:Real, U<:Real}
                    #
                t0, abstol, bfloat = promote(trange[1], aabstol, one(Float64))
                elq0, tt0 = promote(q0[1], t0)
                q0_ = convert(Array{typeof(elq0)}, q0)

                taylorinteg(f, q0_, vec(trange.*one(t0)), order, abstol,
                    params, maxsteps=maxsteps, parse_eqs=parse_eqs)
        end

        function taylorinteg(f, xx0::S, tt0::T, ttmax::U, order::Int, aabstol::V,
                ::Val{true}, params = nothing; maxsteps::Int=500, parse_eqs::Bool=true) where
                    {S<:$R, T<:Real, U<:Real, V<:Real}

            # In order to handle mixed input types, we promote types before integrating:
            t0, tmax, abstol, bfloat = promote(tt0, ttmax, aabstol, one(Float64))
            x0, tfloat1 = promote(xx0, t0)

            return taylorinteg(f, x0, t0, tmax, order, abstol, Val(true), params,
                maxsteps=maxsteps, parse_eqs=parse_eqs)
        end

        function taylorinteg(f, q0::Array{S,1}, tt0::T, ttmax::U, order::Int, aabstol::V,
                ::Val{true}, params = nothing; maxsteps::Int=500, parse_eqs::Bool=true) where
                    {S<:$R, T<:Real, U<:Real, V<:Real}

            #promote to common type before integrating:
            t0, tmax, abstol, afloat = promote(tt0, ttmax, aabstol, one(Float64))
            elq0, tt0 = promote(q0[1], t0)
            #convert the elements of q0 to the common, promoted type:
            q0_ = convert(Array{typeof(elq0)}, q0)

            return taylorinteg(f, q0_, t0, tmax, order, abstol, Val(true), params,
                maxsteps=maxsteps, parse_eqs=parse_eqs)
        end
    end
end


text:
module Go

export fix, exists,
       black, empty, white, other, point, Position, IllegalMove, move!,
       evaluate, random_expression, moves,
       evolve, evolve2, dump_expression, undump,
       replay_direct, replay_pair, replay_latest,
       longest_lived, final_gen,
       parse_log, plot_tramlines, ancestors, multi, multig, children

using AutoHashEquals
using StatefulIterators
using SHA
import Drawing; D = Drawing
import ColorTypes; C = ColorTypes
using ColorVectorSpace
import Base.show, Base.delete!, Base.isless, Base.hash, Base.==
import FixedPointNumbers; N = FixedPointNumbers

# todo
# start on net / pre-processing
# graphs, evolution, tournaments etc

include("utils.jl")
include("emulation.jl")
include("position.jl")
include("expression.jl")
include("play.jl")
include("evolution.jl")
include("log.jl")

end


text:
using GNSSBenchmarks
using Test

@testset "GNSSBenchmarks.jl" begin
    # Write your tests here.
end


text:
using BattleshipSimulator
using Base.Test

@test typeof(layoutBoard(10,10)) == Array{Int64,2}

max_stars_repo_name: mstobb/BattleshipSimulator

max_issues_repo_name: mstobb/BattleshipSimulator

max_forks_repo_name: mstobb/BattleshipSimulator


text:
using LinearMultistep
using Test

@testset "LinearMultistep.jl" begin
    # Write your tests here.
end


text:
# This file contains some code that is chared across several tests

using StaticArrays, Test, ForwardDiff
using NeighbourLists: SMat, SVec
using Test, StaticArrays, ForwardDiff

# ------ generate random configurations -------

function rand_config(N)
   C = SMat( diagm(0 => 2.0 .+ 0.2 * rand(3)) * N )
   X = [ C' * rand(SVec)   for i = 1:ceil(Int, abs(det(C))) ÷ 4 + 2 ]
   pbc = SVec(rand(Bool, 3))
   return X, C, pbc
end

# --------- MANY BODY CODE THAT IS SHARED ACROSS TESTS ------------

sqrt1p(x) = expm1( 0.5 * log1p( x ) )   # sqrt(1 + x) - 1
sqrt1p_d(x) = 0.5 / sqrt(1+x)

function fnbody(r, r0, rcut)
   E = sum( exp.(1.0 .- r./r0) )
   F = prod( (r./rcut .- 1.0) .* (r .< rcut) )
   return sqrt1p(E) * F^2
end

function fnbody_d(r, r0, rcut)
   E = sum( exp.(1.0 .- r./r0) )
   F = prod( (r./rcut .- 1.0) .* (r .< rcut) )
   ∇f = (- sqrt1p_d(E) * F^2 / r0) .* (exp.(1.0 .- r./r0))
   ∇f += (sqrt1p(E) * 2.0 / rcut * F^2) .* (r./rcut .- 1.0 .- 1e-14).^(-1)
   return ∇f
end

# Generate a MODEL N-Body function
function gen_fnbody(rcut, r0=1.0)
   return r->fnbody(r, r0, rcut), r -> fnbody_d(r, r0, rcut)
end

println("Checking that `fnbody` is correct...")
fnbody_ad(rs, r0, rcut) = ForwardDiff.gradient( t -> fnbody(t, r0, rcut), rs )
for n = 1:5
   r = rand(SVector{3, Float64})
   @test fnbody_d(r, 1.0, 2.0) ≈ fnbody_ad(r, 1.0, 2.0)
end

# global assembly of n-body energies and forces

n_body(X, f, M, rcut, C, pbc = (false, false, false),
            nlist = PairList(X, rcut, C, pbc, sorted = true)) =
   n_body(f, M, nlist)

n_body(f, M, nlist::PairList) =
   NeighbourLists.maptosites!(f, zeros(nsites(nlist)),
                                 NeighbourLists.nbodies(M, nlist)) |> sum

grad_n_body(X, df, M, rcut, C, pbc = (false, false, false),
            nlist = PairList(X, rcut, C, pbc, sorted = true)) =
   grad_n_body(df, M, nlist)

grad_n_body(df, M, nlist::PairList) =
   NeighbourLists.maptosites_d!(df, zeros(SVec{Float64}, nsites(nlist)),
                                  NeighbourLists.nbodies(M, nlist))


text:
"""
    trace_map(p::Plot, axis::Symbol=:x)

Return an array of `length(p.data)` that maps each element of `p.data` into an
integer for which number axis of kind `axis` that trace belogs to. `axis` can
either be `x` or `y`. If `x` is given, return the integer for which x-axis the
trace belongs to. Similar for `y`.
"""
function trace_map(p::Plot, axis=:x)
    out = fill(1, length(p.data))
    ax_key = axis == :x ? :xaxis : :yaxis
    for (i, t) in enumerate(p.data)
        if haskey(t, ax_key)
            out[i] = parse(Int, t[ax_key][2:end])
        end
    end
    out
end

trace_map(p::SyncPlot, axis=:x) = trace_map(p.plot, axis)


text:
using EvolutionaryStrategies
using Cambrian
using Test

cfg = get_config("test.yaml")

function test_evo(T::DataType, fitness::Function)
    println("Evolving ", T, " on ", fitness)
    e = T(cfg, fitness)
    @test length(e.population) == e.config.n_population

    step!(e)
    best = e.elites[end]
    @test best.fitness[1] <= 0.0
    @test e.gen == 1

    step!(e)
    new_best = e.elites[end]
    @test new_best.fitness[1] >= best.fitness[1]

    @timev run!(e)
    final_best = e.elites[end]
    @test final_best.fitness[1] >= new_best.fitness[1]
    return final_best.fitness[1]
end

@testset "Evolution xNES" begin
    @test abs(test_evo(xNES, sphere)) < 1e-5
    @test abs(test_evo(xNES, rosenbrock)) < 1.0
end

@testset "Evolution sNES" begin
    @test abs(test_evo(sNES, sphere)) < 1e-5
    @test abs(test_evo(sNES, rosenbrock)) < 1.0
end

# @testset "CMA Evolution Strategy" begin
#     test_evo(CMAES, sphere)
#     test_evo(CMAES, rosenbrock)
# end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export col

"""
    col(;kwargs...)
    col(children::Any;kwargs...)
    col(children_maker::Function;kwargs...)

A Col component.
Utility component to pass to Grid. For more information, see: https://mantine.dev/core/grid/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): Col content
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `className` (String; optional): Often used with CSS to style elements with common properties
- `lg` (Real; optional): Col span at (min-width: theme.breakpoints.lg)
- `md` (Real; optional): Col span at (min-width: theme.breakpoints.md)
- `offset` (Real; optional): Column left offset
- `sm` (Real; optional): Col span at (min-width: theme.breakpoints.sm)
- `span` (Real; optional): Default col span
- `style` (Dict; optional): Inline style override
- `xl` (Real; optional): Col span at (min-width: theme.breakpoints.xl)
- `xs` (Real; optional): Col span at (min-width: theme.breakpoints.xs)
"""
function col(; kwargs...)
        available_props = Symbol[:children, :id, :className, :lg, :md, :offset, :sm, :span, :style, :xl, :xs]
        wild_props = Symbol[]
        return Component("col", "Col", "dash_mantine_components", available_props, wild_props; kwargs...)
end

col(children::Any; kwargs...) = col(;kwargs..., children = children)
col(children_maker::Function; kwargs...) = col(children_maker(); kwargs...)


text:
"""
    pywith(f, o, d=nothing)

Equivalent to `with o as x: f(x)` in Python, where `x` is a `Py`.

On success, the value of `f(x)` is returned.

If an exception occurs but is suppressed then `d` is returned.
"""
function pywith(f, o, d = nothing)
    o = Py(o)
    t = pytype(o)
    exit = t.__exit__
    value = t.__enter__(o)
    exited = false
    try
        f(value)
    catch exc
        if exc isa PyException
            exited = true
            if pytruth(exit(o, exc.t, exc.v, exc.b))
                rethrow()
            else
                d
            end
        else
            rethrow()
        end
    finally
        exited || exit(o, pybuiltins.None, pybuiltins.None, pybuiltins.None)
    end
end
export pywith


text:
####################################################
#              Library OptimizedNeuronalNetwork
####################################################
module OptimizedNeuralNetwork
include("activation_functions.jl")
include("one_layer_neuronal_network.jl")
include("forward_propagation.jl")
include("metric_estimation.jl")
include("weight-initializer-algorithm/weight-initializer-algorithm.jl")
include("backpropagation.jl")
end


text:
using MultiScaleArrays
using OrdinaryDiffEq, DiffEqBase, Test, StochasticDiffEq

#=
struct Cell{B} <: AbstractMultiScaleArrayLeaf{B}
    values::Vector{B}
end
struct Population{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArray{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end
struct Tissue{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArray{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end
struct Embryo{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArrayHead{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end
=#

cell1 = Cell([1.0;2.0;3.0])
cell2 = Cell([4.0;5])
population = construct(Population,deepcopy([cell1,cell2])) # Make a Population from cells
cell3 = Cell([3.0;2.0;5.0])
cell4 = Cell([4.0;6])
population2 = construct(Population,deepcopy([cell3,cell4]))
tissue1 = construct(Tissue,deepcopy([population,population2])) # Make a Tissue from Populations
tissue2 = construct(Tissue,deepcopy([population2,population]))
embryo = construct(Embryo,deepcopy([tissue1,tissue2])) # Make an embryo from Tissues

@test getindices(embryo) == 1:20
@test getindices(embryo,1) == 1:10
@test getindices(embryo,2) == 11:20
@test getindices(embryo,1,1) == 1:5
@test getindices(embryo,1,2) == 6:10
@test getindices(embryo,1,1,1) == 1:3
@test getindices(embryo,1,1,2) == 4:5
@test getindices(embryo,1,2,1) == 6:8
@test getindices(embryo,1,2,2) == 9:10
@test getindices(embryo,1,2,2,1) == 9:9
@test getindices(embryo,1,2,2,2) == 10:10
@test getindices(embryo,2,1) == 11:15
@test getindices(embryo,2,2) == 16:20
@test getindices(embryo,2,1,1) == 11:13
@test getindices(embryo,2,1,2) == 14:15
@test getindices(embryo,2,2,1) == 16:18
@test getindices(embryo,2,2,2) == 19:20
@test getindices(embryo,2,2,2,1) == 19:19
@test getindices(embryo,2,2,2,2) == 20:20


text:
# distance to center of the simulation cell
function gen_dr_center( r, LatVecs )
    #
    cx = 0.5*sum(LatVecs[:,1])
    cy = 0.5*sum(LatVecs[:,2])
    cz = 0.5*sum(LatVecs[:,3])
    center = [cx,cy,cz]
    #
    Npoints = size(r)[2]
    dr = Array{Float64}(undef,Npoints)
    for ip=1:Npoints
        dr[ip] = norm( r[:,ip] - center )
    end
    return dr
end


text:
abstract type Material end

struct CheckerMaterial <: Material
    scale::Float64
    reflectiveness::Float64
end
CheckerMaterial(scale) = CheckerMaterial(scale, 0)

function sample(cm::CheckerMaterial, ray::Ray3, position::Vector3, normal::Vector3)
    abs(Int64(floor(position.x * cm.scale) + floor(position.z * cm.scale))) & 1 == 0 ? Black : White
end

struct PhongMaterial <: Material
    diffuse::ColorRGB
    specular::ColorRGB
    shininess::Int64
    reflectiveness::Float64
end
PhongMaterial(diffuse, specular, shininess) = PhongMaterial(diffuse, specular, shininess, 0)

function sample(cm::PhongMaterial, ray::Ray3, position::Vector3, normal::Vector3)
    ndotl = normal ⋅ LightDirection()
    𝐡 = LightDirection() - ray.direction |> normalize
    ndoth = normal ⋅ 𝐡
    diffuse_term = cm.diffuse * max(ndotl, 0)
    specular_term = cm.specular * max(ndoth, 0)^cm.shininess
    modulate(LightColor(), diffuse_term + specular_term)
end


text:
using TaglessTrees

add(args...) = call(+, args...)

t = add(term(2), term(3))
evaluate(t)
trace(t)


text:
# Callbacks

function simple_trace(θ::AbstractArray, t::Int, n::Int, v::Real)
    @printf("Iter %d: avg.loss = %.4e\n",
        t, v / n)
end

function gtcompare_trace(θg::AbstractArray)
    function _trace(θ::AbstractArray, t::Int, n::Int, v::Real)
        dev = vecnorm(θ - θg)
        @printf("Iter %d: avg.loss = %.4e, deviation = %.4e\n",
            t, v / n, dev)
    end
    return _trace
end


text:
"""
    eval_interp(m::Module, ex)

evaluate the interpolation operator in `ex` inside given module `m`.
"""
function eval_interp(m::Module, ex)
    ex isa Expr || return ex
    if ex.head === :$
        x = ex.args[1]
        if x isa Symbol && isdefined(m, x)
            return Base.eval(m, x)
        else
            return ex
        end
    end
    return Expr(ex.head, map(x->eval_interp(m, x), ex.args)...)
end

function eval_literal(m::Module, ex)
    ex isa Expr || return ex
    if ex.head === :call && all(is_literal, ex.args[2:end])
        return Base.eval(m, ex)
    end
    return Expr(ex.head, map(x->eval_literal(m, x), ex.args)...)
end

replace_symbol(x::Symbol, name::Symbol, value) = x === name ? value : x
replace_symbol(x, ::Symbol, value) = x # other expressions

function replace_symbol(ex::Expr, name::Symbol, value)
    Expr(ex.head, map(x->replace_symbol(x, name, value), ex.args)...)
end

function subtitute(ex::Expr, replace::Pair)
    name, value = replace
    return replace_symbol(ex, name, value)
end

"""
    name_only(ex)

Remove everything else leaving just names, currently supports
function calls, type with type variables, subtype operator `<:`
and type annotation `::`.

# Example

```julia
julia> using Expronicon

julia> name_only(:(sin(2)))
:sin

julia> name_only(:(Foo{Int}))
:Foo

julia> name_only(:(Foo{Int} <: Real))
:Foo

julia> name_only(:(x::Int))
:x
```
"""
function name_only(@nospecialize(ex))
    ex isa Symbol && return ex
    ex isa Expr || error("unsupported expression $ex")
    ex.head in [:call, :curly, :(<:), :(::), :where, :function, :(=), :(->)] && return name_only(ex.args[1])
    error("unsupported expression $ex")
end

"""
    rm_lineinfo(ex)

Remove `LineNumberNode` in a given expression.
"""
function rm_lineinfo(ex)
    @match ex begin
        Expr(:macrocall, name, line, args...) => Expr(:macrocall, name, line, map(rm_lineinfo, args)...)
        Expr(head, args...) => Expr(head, map(rm_lineinfo, filter(x->!(x isa LineNumberNode), args))...)
        _ => ex
    end
end

"""
    prettify(ex)

Prettify given expression, remove all `LineNumberNode` and
extra code blocks.
"""
function prettify(ex)
    ex isa Expr || return ex
    ex = rm_lineinfo(ex)
    ex = flatten_blocks(ex)
    ex = rm_nothing(ex)
    return ex
end

"""
    flatten_blocks(ex)

Remove hierachical expression blocks.
"""
function flatten_blocks(ex)
    ex isa Expr || return ex
    ex.head === :block || return Expr(ex.head, map(_flatten_blocks, ex.args)...)
    has_block = any(ex.args) do x
        x isa Expr && x.head === :block
    end

    if has_block
        return flatten_blocks(_flatten_blocks(ex))
    end
    return Expr(ex.head, map(flatten_blocks, ex.args)...)
end

function _flatten_blocks(ex)
    ex isa Expr || return ex
    ex.head === :block || return Expr(ex.head, map(flatten_blocks, ex.args)...)

    args = []
    for stmt in ex.args
        if stmt isa Expr && stmt.head === :block
            for each in stmt.args
                push!(args, flatten_blocks(each))
            end
        else
            push!(args, flatten_blocks(stmt))
        end
    end
    return Expr(:block, args...)
end

function rm_nothing(ex)
    @match ex begin
        Expr(:block, args...) => Expr(:block, filter(x->x!==nothing, args)...)
        Expr(head, args...) => Expr(head, map(rm_nothing, args)...)
        _ => ex
    end
end

"""
    rm_annotations(x)

Remove type annotation of given expression.
"""
function rm_annotations(x)
    x isa Expr || return x
    if x.head == :(::)
        if length(x.args) == 1 # anonymous
            return gensym("::$(x.args[1])")
        else
            return x.args[1]
        end
    elseif x.head in [:(=), :kw] # default values
        return rm_annotations(x.args[1])
    else
        return Expr(x.head, map(rm_annotations, x.args)...)
    end
end


text: println("An example was run!")

max_stars_repo_path: Chapter01/08. Useful options for interaction with Julia/example.jl


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Deprecated functions.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

#                       Introduced in TerminalPager v0.1
# ==============================================================================

export @dpr
macro dpr(expr)
    Base.depwarn("@dpr is deprecated, use @help instead.", :dpr)
    return :(@help $(esc(expr)))
end


text:
@testset "src/analysis/validation.jl" begin
    time_step_h = network_mn["time_step"]/3600

    @testset "compare nodes" begin
	    ## Check 'compare nodes'
	    for (i, node_dict) in network_mn["nw"]["1"]["node"]
            node_id = string(node_dict["name"])
            node_info = network_mn["nw"]["1"]["node"][node_id]
            elevation = node_info["elevation"]
            node_df = get_node_dataframe(network_mn, solution, wntr_network, wntr_result,
                                         node_id)

            ## I encountered this error and hence changed the syntax for all dataframe usage
            ## (JJS 5/10/21):
            #ERROR: ArgumentError: syntax df[column] is not supported use df[!, column] instead
            for t in 1:length(node_df[!, "time"])
                @test check_difference(node_df[!, "time"][t], time_step_h*t, tol)
                @test check_difference(node_df[!, "elevation"][t], elevation, tol)
                wntr_head = getproperty(wntr_result.node["head"], node_id).values[t]
                @test check_difference(node_df[!, "head_wntr"][t], wntr_head, tol)
                @test check_difference(node_df[!, "head_watermodels"][t],
                                       solution["nw"][string(t)]["node"][node_id]["h"], tol)
                wntr_pres = getproperty(wntr_result.node["pressure"], node_id).values[t]
                @test check_difference(node_df[!, "pressure_wntr"][t], wntr_pres, tol)
                @test check_difference(node_df[!, "pressure_watermodels"][t],
                                       solution["nw"][string(t)]["node"][node_id]["p"], tol)
            end
        end
    end

    @testset "compare tanks" begin
        ## Check 'compare tanks'
        for (i, tank_dict) in network_mn["nw"]["1"]["tank"]
            tank_id = i
            tank_node_id = string(network_mn["nw"]["1"]["tank"][i]["node"])
            tank_name = tank_node_id
            diameter = wntr_network.nodes._data[tank_name].diameter
            tank_df = get_tank_dataframe(network_mn,solution,wntr_network,wntr_result,
                                         tank_id)

            for t in 1:length(tank_df[!, "time"])
                @test check_difference(tank_df[!, "time"][t], time_step_h*t, tol)
                wntr_vol = getproperty(wntr_result.node["pressure"],
                                       tank_name).values[t]*(0.25*pi*diameter^2)
                @test check_difference(tank_df[!, "volume_wntr"][t], wntr_vol, tol)
                @test check_difference(tank_df[!, "volume_watermodels"][t],
                     solution["nw"][string(t)]["node"][tank_name]["p"]*(0.25*pi*diameter^2),
                                       tol)
                wntr_pres = getproperty(wntr_result.node["pressure"], tank_name).values[t]
                @test check_difference(tank_df[!, "level_wntr"][t], wntr_pres, tol)
                @test check_difference(tank_df[!, "level_watermodels"][t],
                                       solution["nw"][string(t)]["node"][tank_name]["p"],
                                       tol)
            end
        end
    end

    @testset "compare pipes" begin
        ## check 'compare pipes'
        for (i,pipe_dict) in network_mn["nw"]["1"]["pipe"]
            pipe_id = i
            pipe_name = "pipe"*pipe_id
            node_fr = string(network_mn["nw"]["1"]["pipe"][pipe_id]["node_fr"])
            node_to = string(network_mn["nw"]["1"]["pipe"][pipe_id]["node_to"])
            pipe_df = get_pipe_dataframe(network_mn,solution,wntr_network,wntr_result,
                                         pipe_id)

            for t in 1:length(pipe_df[!, "time"])
                @test check_difference(pipe_df[!, "time"][t], time_step_h*t, tol)
                wntr_flow = getproperty(wntr_result.link["flowrate"], pipe_name).values[t]
                @test check_difference(pipe_df[!, "flow_wntr"][t], wntr_flow, tol)
                @test check_difference(pipe_df[!, "flow_watermodels"][t],
                                       solution["nw"][string(t)]["pipe"][pipe_id]["q"], tol)
                wntr_head_nf = getproperty(wntr_result.node["head"], node_fr).values[t]
                wntr_head_nt = getproperty(wntr_result.node["head"], node_to).values[t]
                @test check_difference(pipe_df[!, "head_loss_wntr"][t],
                                       wntr_head_nf - wntr_head_nt, tol)
                @test check_difference(pipe_df[!, "head_loss_watermodels"][t],
                                       solution["nw"][string(t)]["node"][node_fr]["h"] -
                                       solution["nw"][string(t)]["node"][node_to]["h"], tol)
            end
        end

        ## check 'compare shortpipes'
        for (i,short_pipe_dict) in network_mn["nw"]["1"]["short_pipe"]
            short_pipe_id = i
            short_pipe_name = "short_pipe"*short_pipe_id
            node_fr = string(network_mn["nw"]["1"]["short_pipe"][short_pipe_id]["node_fr"])
            node_to = string(network_mn["nw"]["1"]["short_pipe"][short_pipe_id]["node_to"])
            short_pipe_df = get_short_pipe_dataframe(network_mn,solution,wntr_network,
                                                     wntr_result, short_pipe_id)

            for t in 1:length(short_pipe_df[!, "time"])
                @test check_difference(short_pipe_df[!, "time"][t], time_step_h*t, tol)
                wntr_flow = getproperty(wntr_result.link["flowrate"],
                                        short_pipe_name).values[t]
                @test check_difference(short_pipe_df[!, "flow_wntr"][t], wntr_flow, tol)
                @test check_difference(short_pipe_df[!, "flow_watermodels"][t],
                              solution["nw"][string(t)]["short_pipe"][short_pipe_id]["q"],
                                       tol)
                wntr_head_nf = getproperty(wntr_result.node["head"], node_fr).values[t]
                wntr_head_nt = getproperty(wntr_result.node["head"], node_to).values[t]
                @test check_difference(short_pipe_df[!, "head_loss_wntr"][t],
                                       wntr_head_nf - wntr_head_nt, tol)
                @test check_difference(short_pipe_df[!, "head_loss_watermodels"][t],
                                       solution["nw"][string(t)]["node"][node_fr]["h"] -
                                       solution["nw"][string(t)]["node"][node_to]["h"], tol)
            end
        end
    end

    @testset "compare valves" begin
        ## check 'compare valves'
        for (i,valve_dict) in network_mn["nw"]["1"]["valve"]
            valve_id = i
            valve_name = "valve"*valve_id
            node_fr = string(network_mn["nw"]["1"]["valve"][valve_id]["node_fr"])
            node_to = string(network_mn["nw"]["1"]["valve"][valve_id]["node_to"])
            valve_df = get_valve_dataframe(network_mn,solution,wntr_network,wntr_result,
                                           valve_id)

            for t in 1:length(valve_df[!, "time"])
                @test check_difference(valve_df[!, "time"][t], time_step_h*t, tol)
                wntr_st = getproperty(wntr_result.link["status"], valve_name).values[t]
                @test check_difference(valve_df[!, "status_wntr"][t], wntr_st, tol)
                @test check_difference(valve_df[!, "status_watermodels"][t],
                                     solution["nw"][string(t)]["valve"][valve_id]["status"],
                                       tol)
                wntr_flow = getproperty(wntr_result.link["flowrate"], valve_name).values[t]
                @test check_difference(valve_df[!, "flow_wntr"][t], wntr_flow, tol)
                @test check_difference(valve_df[!, "flow_watermodels"][t],
                                       solution["nw"][string(t)]["valve"][valve_id]["q"],
                                       tol)
                wntr_head_nf = getproperty(wntr_result.node["head"], node_fr).values[t]
                wntr_head_nt = getproperty(wntr_result.node["head"], node_to).values[t]
                @test check_difference(valve_df[!, "head_loss_wntr"][t],
                                       wntr_head_nf - wntr_head_nt, tol)
                @test check_difference(valve_df[!, "head_loss_watermodels"][t],
                                       solution["nw"][string(t)]["node"][node_fr]["h"] -
                                       solution["nw"][string(t)]["node"][node_to]["h"], tol)
            end
        end
    end

    @testset "compare pumps" begin
        ## check 'compare pumps'
        function compute_pump_power(pump_obj,q,dh,wntr_data)
            if pump_obj.efficiency == nothing
                eff_curve = nothing
            else
                eff_curve = wntr_data.get_curve(pump_obj.efficiency).points
            end
            if eff_curve != nothing
                for j in 1:length(eff_curve)
                    if (j == 1) & (q < eff_curve[j][1])
                        # use the first Q-eff point when q < q_min on the eff curve
                        eff = eff_curve[j][2]/100
                        break
                    elseif j == length(eff_curve)
                        # use the last Q-eff point when q > q_max on the eff curve
                        eff = eff_curve[j][2]/100
                    elseif (q >= eff_curve[j][1]) & (q < eff_curve[j+1][1])
                        # linear interpolation of the efficiency curve
                        # why not use Interpolations package? see `plots/pumps.jl`, JJS
                        # 5/13/21
                        a1 = (eff_curve[j+1][1]-q)/(eff_curve[j+1][1]-eff_curve[j][1])
                        a2 = (q-eff_curve[j][1])/(eff_curve[j+1][1]-eff_curve[j][1])
                        eff = (a1*eff_curve[j][2]+a2*eff_curve[j+1][2])/100
                        break
                    end
                end
            elseif wntr_data.options.energy.global_efficiency != nothing
                 # e.g. convert 60 to 60% (0.6)
                eff = wntr_data.options.energy.global_efficiency/100
            else
                Memento.info(_LOGGER, "No pump efficiency provided, 75% is used")
                eff = 0.75
            end

            power = _WM._GRAVITY * _WM._DENSITY * dh * q / eff # Watt
            return power
        end

        for (i,pump_dict) in network_mn["nw"]["1"]["pump"]
            pump_id = i
            pump_name = "pump"*pump_id
            pump_obj = wntr_network.get_link(pump_name)
            node_fr = string(network_mn["nw"]["1"]["pump"][pump_id]["node_fr"])
            node_to = string(network_mn["nw"]["1"]["pump"][pump_id]["node_to"])
            pump_df = get_pump_dataframe(network_mn,solution,wntr_network,wntr_result,
                                         pump_id)

            for t in 1:length(pump_df[!, "time"])
                @test check_difference(pump_df[!, "time"][t], time_step_h*t, tol)
                wntr_st = getproperty(wntr_result.link["status"], pump_name).values[t]
                @test check_difference(pump_df[!, "status_wntr"][t], wntr_st, tol)
                @test check_difference(pump_df[!, "status_watermodels"][t],
                                       solution["nw"][string(t)]["pump"][pump_id]["status"],
                                       tol)
                wntr_flow = getproperty(wntr_result.link["flowrate"], pump_name).values[t]
                @test check_difference(pump_df[!, "flow_wntr"][t], wntr_flow, tol)
                @test check_difference(pump_df[!, "flow_watermodels"][t],
                                       solution["nw"][string(t)]["pump"][pump_id]["q"], tol)
                wntr_head_nf = getproperty(wntr_result.node["head"], node_fr).values[t]
                wntr_head_nt = getproperty(wntr_result.node["head"], node_to).values[t]
                @test check_difference(pump_df[!, "head_gain_wntr"][t],
                                       wntr_st*(wntr_head_nt - wntr_head_nf), tol)
                @test check_difference(pump_df[!, "head_gain_watermodels"][t],
                                       solution["nw"][string(t)]["pump"][pump_id]["status"]*
                                       (solution["nw"][string(t)]["node"][node_to]["h"]-
                                        solution["nw"][string(t)]["node"][node_fr]["h"]),
                                       tol)
                power_wntr = wntr_st*compute_pump_power(pump_obj, wntr_flow,
                                   wntr_st*(wntr_head_nt - wntr_head_nf), wntr_network)
                power_watermodels = solution["nw"][string(t)]["pump"][pump_id]["status"]*
                compute_pump_power(pump_obj,solution["nw"][string(t)]["pump"][pump_id]["q"],
                                   solution["nw"][string(t)]["pump"][pump_id]["status"]*
                                   (solution["nw"][string(t)]["node"][node_to]["h"] -
                                    solution["nw"][string(t)]["node"][node_fr]["h"]),
                                   wntr_network)
                @test check_difference(pump_df[!, "power_wntr"][t], power_wntr, tol)
                @test check_difference(pump_df[!, "power_watermodels"][t],
                                       power_watermodels, tol)
                energy_price =
                    network_mn["nw"][string(t)]["pump"][pump_id]["energy_price"]*3600 # $/Wh
                @test check_difference(pump_df[!, "cost_wntr"][t],
                                       energy_price*time_step_h*power_wntr,tol)
                @test check_difference(pump_df[!, "cost_watermodels"][t],
                                       energy_price*time_step_h*power_watermodels,tol)
            end
        end
    end
end


text:
"""
This holds all the scripts that are relevant to the downloading of ECMWF reanalysis
data, which includes the following functionalities:
    - Creation of python download scripts
    - Creation of respective folders to put data in within the ClimateERA directory

"""

# Creation of ClimateERA Download Scripts

function eradscriptcreate(dataID::Integer,emod::Dict,epar::Dict,ereg::Dict)

    if !(epar["level"] == "sfc")
          fname = "$(emod["prefix"])-$(ereg["fol"])-$(epar["ID"])-$(epar["level"])hPa";
    else; fname = "$(emod["prefix"])-$(ereg["fol"])-$(epar["ID"])-$(epar["level"])";
    end

    fID = open("$(fname).py","w");

    write(fID,"#!/usr/bin/env python\n");

    if dataID == 1
        write(fID,"import cdsapi\n");
        write(fID,"c = cdsapi.Client()\n\n");
    else
        write(fID,"from ecmwfapi import ECMWFDataServer\n");
        write(fID,"server = ECMWFDataServer()\n\n");
    end

    return fname,fID

end

function eradscriptheader(fID,dataID::Integer,emod::Dict,epar::Dict)

    parID = epar["ID"];

    if dataID == 1
        write(fID,"c.retrieve(\"$(emod["moduleprint"])\",\n");
        write(fID,"    {\n");
        write(fID,"        \"product_type\": \"reanalysis\",\n");
    else
        write(fID,"server.retrieve({\n");
        write(fID,"    \"class\": \"ei\",\n");
        write(fID,"    \"dataset\": \"interim\",\n");
        write(fID,"    \"stream\": \"oper\",\n");

        if !all(occursin.(["cape","prcp"],parID))
              write(fID,"    \"type\": \"an\",\n");
        else; write(fID,"    \"type\": \"fc\",\n");
        end

        write(fID,"    \"type\": \"an\",\n");
        write(fID,"    \"expver\": \"1\",\n");
    end

end

function eradscriptpprint(fID,dataID::Integer,emod::Dict,epar::Dict)

    pre = epar["level"];

    if dataID == 1; var = epar["era5"];
        write(fID,"        \"variable\": \"$(var)\",\n");
        if !(pre == "sfc"); write(fID,"        \"pressure_level\": $(pre),\n"); end
    else; var = epar["erai"];
        write(fID,"    \"param\": $(var),\n");
        if pre == "sfc"
            write(fID,"    \"levtype\": \"sfc\",\n");
        else
            write(fID,"    \"levtype\": \"pl\",\n");
            write(fID,"    \"levelist\": \"$(pre)\",\n");
        end
    end

end

function eradscriptregion(fID,dataID::Integer,emod::Dict,ereg::Dict)

    estep = ereg["step"]; N,S,E,W = ereg["grid"]
    if W == E; E = W + 0.1 end
    if N == S; S = N - 0.1 end
    if dataID == 1
        if !ereg["isglobe"]
              write(fID,"        \"area\": [$(N),$(W),$(S),$(E)],\n");
        end
        write(fID,"        \"grid\": [$(estep),$(estep)],\n");
    else
        if !ereg["isglobe"]
              write(fID,"    \"area\": \"$(N)/$(W)/$(S)/$(E)\",\n");
        end
        write(fID,"    \"grid\": \"$(estep)/$(estep)\",\n");
    end

end

function eradscriptdprint(fID,dataID::Integer,epar::Dict,yr::Integer,mo::Integer)

    parID = epar["ID"]; ndy = daysinmonth(yr,mo); mo = mo2str(mo);
    dystr = "\"25\""; for dy = 26 : ndy; dystr = string(dystr,",\"$(dy2str(dy))\""); end

    if dataID == 1
        write(fID,"        \"year\": \"$(yr)\",\n");
        write(fID,"        \"month\": \"$(mo)\",\n");
        write(fID,"        \"day\":[\n");
        write(fID,"            \"01\",\"02\",\"03\",\"04\",\"05\",\"06\",\n");
        write(fID,"            \"07\",\"08\",\"09\",\"10\",\"11\",\"12\",\n");
        write(fID,"            \"13\",\"14\",\"15\",\"16\",\"17\",\"18\",\n");
        write(fID,"            \"19\",\"20\",\"21\",\"22\",\"23\",\"24\",\n");
        write(fID,"            $(dystr)\n");
        write(fID,"        ],\n");
        write(fID,"        \"time\":[\n");
        write(fID,"            \"00:00\",\"01:00\",\"02:00\",\"03:00\",\"04:00\",\n");
        write(fID,"            \"05:00\",\"06:00\",\"07:00\",\"08:00\",\"09:00\",\n");
        write(fID,"            \"10:00\",\"11:00\",\"12:00\",\"13:00\",\"14:00\",\n");
        write(fID,"            \"15:00\",\"16:00\",\"17:00\",\"18:00\",\"19:00\",\n");
        write(fID,"            \"20:00\",\"21:00\",\"22:00\",\"23:00\"\n");
        write(fID,"        ],\n");
    else
        write(fID,"    \"date\": \"$(yr)-$(mo)-01/to/$(yr)-$(mo)-$(ndy)\",\n");
        if !all(occursin.(["cape","prcp"],parID))
            write(fID,"    \"time\": \"00:00:00/06:00:00/12:00:00/18:00:00\",\n");
            write(fID,"    \"step\": \"0\",\n");
        else
            write(fID,"    \"time\": \"00:00:00/12:00:00\",\n");
            write(fID,"    \"step\": \"3/6/9/12\",\n");
        end
    end

end

function eradscripttarget(
    fID,dataID::Integer,fname::AbstractString,
    yr::Integer,mo::Integer
)

    fnc = joinpath("..","raw","$(yr)","$(fname)-$(yrmo2str(yr,mo)).nc")
    if dataID == 1
        write(fID,"        \"format\": \"netcdf\"\n");
        write(fID,"    },\n");
        write(fID,"    \"$(fnc)\")\n\n");
    else
        write(fID,"    \"format\": \"netcdf\",\n");
        write(fID,"    \"target\": \"$(fnc)\",\n");
        write(fID,"})\n\n");
    end

end

# Master ClimateERA Download Scripts.end

function eradscript(emod::Dict,epar::Dict,ereg::Dict,etime::Dict)

    dataID = emod["datasetID"];
    fname,fID = eradscriptcreate(dataID,emod,epar,ereg);

    for yr = etime["Begin"] : etime["End"]
        for mo = 1 : 12

            eradscriptheader(fID,dataID,emod,epar);
            eradscriptpprint(fID,dataID,emod,epar);
            eradscriptregion(fID,dataID,emod,ereg);
            eradscriptdprint(fID,dataID,epar,yr,mo);
            eradscripttarget(fID,dataID,fname,yr,mo);

        end
    end

    close(fID); return "$(fname).py"

end

function eradownload(
    init::Dict, eroot::Dict;
    modID::AbstractString, parID::AbstractString,
    regID::AbstractString="GLB", timeID::Union{Integer,Vector}=0,
    gres::Real=0
)

    emod,epar,ereg,etime = erainitialize(
        init;
        modID=modID,parID=parID,regID=regID,timeID=timeID,
        gres=gres
    );
    eradownload(emod,epar,ereg,etime,eroot)

end

function eradownload(emod::Dict,epar::Dict,ereg::Dict,etime::Dict,eroot::Dict)

    prelist = emod["levels"]; dataID = emod["datasetID"];

    jfol = joinpath(DEPOT_PATH[1],"files/ClimateERA/"); mkpath(jfol);
    if dataID == 1; dwnsh = joinpath(jfol,"erad5.sh");
    else;           dwnsh = joinpath(jfol,"eradi.sh");
    end

    if !isfile(dwnsh);
        if dataID == 1; template = joinpath(@__DIR__,"../extra/erad5_eg.sh");
        else;           template = joinpath(@__DIR__,"../extra/eradi_eg.sh");
        end
        cp(template,dwnsh,force=true);
    end

    @info "$(Dates.now()) - Creating download scripts, directories and subdirectories for data downloading, temporary storage, analysis and image creation..."

    for preii in prelist; epar["level"] = preii;

        @debug "$(Dates.now()) - Creating download scripts and directories at pressure level preii ..."
        fname = eradscript(emod,epar,ereg,etime);
        efol  = erafolder(emod,epar,ereg,etime,eroot);

        @debug "$(Dates.now()) - Moving download scripts to tmp folder $(efol["tmp"]) ..."
        mv(fname,joinpath(efol["tmp"],fname),force=true);
        cp(dwnsh,joinpath(efol["tmp"],"erad.sh"),force=true);

        @debug "$(Dates.now()) - Saving module and variable information into efol[\"var\"] directory ..."
        @save "info_par.jld2" emod epar;
        mv("info_par.jld2",joinpath(efol["var"],"info_par.jld2"),force=true)

    end

    efol = erafolder(emod,epar,ereg,etime,eroot);
    @info "$(Dates.now()) - Saving region information into efol[\"reg\"] directory ..."
    @save "info_reg.jld2" ereg;
    mv("info_reg.jld2",joinpath(efol["reg"],"info_reg.jld2"),force=true)

end

function eralvlsplit(
    init::Dict, eroot::Dict;
    modID::AbstractString, parID::AbstractString,
    regID::AbstractString="GLB", timeID::Union{Integer,Vector}=0,
    gres::Real=0, plist::Union{Integer,Vector{<:Real}}=0
)

    emod,epar,ereg,etime = erainitialize(
        init;
        modID=modID,parID=parID,regID=regID,timeID=timeID,
        gres=gres
    );
    eralvlsplit(emod,epar,ereg,etime,eroot,plist)

end

function eralvlsplit(
    emod::Dict,epar::Dict,ereg::Dict,etime::Dict,eroot::Dict,
    plist::Union{Integer,Vector{<:Real}}=0
)

    datevec = collect(Date(etime["Begin"],1):Month(1):Date(etime["End"],12));

    if typeof(plist) <: Array

        for dtii in datevec

            fol = joinpath(eroot["era"],ereg["fol"],epar["ID"],"raw",yr2str(dtii));
            fnc = "$(emod["prefix"])-$(ereg["fol"])-$(epar["ID"])-$(yrmo2str(dtii)).nc"
            eds = NCDataset(joinpath(fol,fnc))
            pre = eds["level"][:]

            for preii in plist
                ip = argmin(abs.(pre.-preii)); epar["level"] = pre[ip];
                efol = erafolder(emod,epar,ereg,etime,eroot,epar["level"]);
                dataii = eds[epar["IDnc"]][:,:,ip,:]*1
                erarawsave(dataii,emod,epar,ereg,dtii,eroot)
            end

            close(eds); rm(joinpath(fol,fnc),force=true)

        end

    elseif iszero(plist)

        for dtii in datevec

            fol = joinpath(eroot["era"],ereg["fol"],epar["ID"],"raw",yr2str(dtii));
            fnc = "$(emod["prefix"])-$(ereg["fol"])-$(epar["ID"])-$(yrmo2str(dtii)).nc"
            eds = NCDataset(joinpath(fol,fnc))
            pre = eds["level"][:]; np = length(pre)

            for ip in 1 : np; epar["level"] = pre[ip];

                efol = erafolder(emod,epar,ereg,etime,eroot,epar["level"]);
                dataii = eds[epar["IDnc"]][:,:,ip,:]*1
                erarawsave(dataii,emod,epar,ereg,dtii,eroot)

            end

            close(eds); rm(joinpath(fol,fnc),force=true)

        end

    else

        for dtii in datevec

            fol = joinpath(eroot["era"],ereg["fol"],epar["ID"],"raw",yr2str(dtii));
            fnc = "$(emod["prefix"])-$(ereg["fol"])-$(epar["ID"])-$(yrmo2str(dtii)).nc"
            eds = NCDataset(joinpath(fol,fnc))
            pre = eds["level"][:]; ip = argmin(abs.(pre.-plist)); epar["level"] = pre[ip];

            efol = erafolder(emod,epar,ereg,etime,eroot,epar["level"]);
            dataii = eds[epar["IDnc"]][:,:,ip,:]*1
            erarawsave(dataii,emod,epar,ereg,dtii,eroot)

            close(eds); rm(joinpath(fol,fnc),force=true)

        end

    end

end


text:
### A Pluto.jl notebook ###
# v0.12.6

using Markdown
using InteractiveUtils

# ╔═╡ 971fc5f0-1c52-11eb-27da-43003bbda193
begin
	using DifferentialEquations
	using Plots
end

# ╔═╡ 79e6b4e0-1c56-11eb-30a7-b78aa49dc552
function msd_sys(du,u,p,t)
  # Manual implementation of state transition matrix
  du[1] = u[2]
  du[2] = -p[1]*u[1] - p[2]*u[2]

  # Add acceleration due to gravity
  du[2] += p[3]

  # Add disturbance in the force
  du[2] +=  p[4]*sin(p[5]*t)

  # Add control input
  du[2] += p[6]
end

# ╔═╡ 943468fe-1c56-11eb-07ed-9d4f3578e7c6
begin
	g = 9.8 # gravity in metres per second per second
	mass = 0.1
	k_s = 2.5 # spring constant
	k_d = 0.1 # damping coefficient

	u0 = [-1., 0.]

	f_smc = 10.0

	p_simple = [k_s/mass, k_d/mass, g, 0, 0, 0]
	p_dist = [k_s/mass, k_d/mass, g, 5., 10., 0]
	p_smc = [k_s/mass, k_d/mass, g, 5., 10., 0.]

	tspan = (0.0,10.0)
end

# ╔═╡ a705ea40-1c56-11eb-1a72-136337ca099f
function smc_affect!(msd_ss_model)
    sigma = 20.0*msd_ss_model.u[1] + 1.0*msd_ss_model.u[2]
    if sigma < 0
        msd_ss_model.p[6] = +f_smc/mass
    else
        msd_ss_model.p[6] = -f_smc/mass
    end
end

# ╔═╡ b9fb9550-1c56-11eb-1080-b7ec4f1598fe
pcb = PeriodicCallback(smc_affect!, 10/10000)

# ╔═╡ b0c2973e-1c56-11eb-063b-cf48350c9806
begin
	prob_simple = ODEProblem(msd_sys,u0,tspan,p_simple)
	prob_dist = ODEProblem(msd_sys,u0,tspan,p_dist)
	prob_smc = ODEProblem(msd_sys,u0,tspan,p_smc)
end

# ╔═╡ dbb46550-1c56-11eb-32c5-3dd9bceda11a
begin
    raw_sol = solve(prob_simple,Tsit5());
    dist_sol = solve(prob_dist,Tsit5());
    smc_sol = solve(prob_smc,Tsit5(),callback=pcb);
end

# ╔═╡ fc2b3f70-1c56-11eb-1322-b9d6b5728032
begin
	plot(raw_sol,vars=(0,1),w=1,label="Raw", layout = (2,1), subplot=1)
	plot!(dist_sol,vars=(0,1),w=1,label="Disturbed")
	plot!(smc_sol,vars=(0,1),w=1,label="SMC")
	plot!(raw_sol,vars=(1,2),w=1,label="Raw",subplot=2)
	plot!(dist_sol,vars=(1,2),w=1,label="Disturbed",subplot=2)
	plot!(smc_sol,vars=(1,2),w=1,label="SMC",subplot=2,xlim=(-1.1,+1.5),ylim=(-5,11))
end

# ╔═╡ Cell order:
# ╠═971fc5f0-1c52-11eb-27da-43003bbda193
# ╠═79e6b4e0-1c56-11eb-30a7-b78aa49dc552
# ╠═943468fe-1c56-11eb-07ed-9d4f3578e7c6
# ╠═a705ea40-1c56-11eb-1a72-136337ca099f
# ╠═b9fb9550-1c56-11eb-1080-b7ec4f1598fe
# ╠═b0c2973e-1c56-11eb-063b-cf48350c9806
# ╠═dbb46550-1c56-11eb-32c5-3dd9bceda11a
# ╠═fc2b3f70-1c56-11eb-1322-b9d6b5728032


text:
"""
Function replaces the negative values in a datacube with missing.
"""
function replace_negative(datacube)
    for i in 1:prod(size(datacube))
        if ismissing(datacube[i])
            continue
        elseif datacube[i] < 0
                datacube[i] = missing
        else
            continue
        end
    end
    return datacube
end


text:
"""
    PararealCache <: AbstractTimeParallelCache

A composite type for the [`AbstractTimeParallelCache`](@ref) of [`Parareal`](@ref).

# Constructors
```julia
PararealCache(U, F, G, T)
PararealCache(problem::AbstractInitialValueProblem, parareal::Parareal)
```

# Arguments
- `U :: AbstractVector{<:Union{Number, AbstractVector{<:Number}}}` : chunks' initial conditions.
- `F :: AbstractVector{<:Union{Number, AbstractVector{<:Number}}}` : chunks' fine values.
- `G :: AbstractVector{<:Union{Number, AbstractVector{<:Number}}}` : chunks' coarse values.
- `T :: AbstractVector{<:Number}` : chunks' initial time.
"""
struct PararealCache{U_T, F_T, G_T, T_T} <: AbstractTimeParallelCache
    U::U_T
    F::F_T
    G::G_T
    T::T_T
end

function PararealCache(problem::AbstractInitialValueProblem, parareal::Parareal)
    @↓ P = parareal
    @↓ u0, (t0, tN) ← tspan = problem
    u0_T = typeof(u0)
    U = Vector{u0_T}(undef, P+1)
    U[1] = u0
    F = similar(U)
    G = similar(U)
    t0_T = typeof(t0)
    T = Vector{t0_T}(undef, P+1)
    T[1] = t0
    return PararealCache(U, F, G, T)
end


text:
using Documenter, MathProgBase

makedocs(
    format = :html,
    sitename = "MathProgBase",
    pages = [
        "index.md",
    ]
)

deploydocs(
    repo   = "github.com/JuliaOpt/MathProgBase.jl.git",
    target = "build",
    osname = "linux",
    julia  = "0.6",
    deps   = nothing,
    make   = nothing
)


text:
using Printf
using NPZ
include("spectre-defense/util.jl")

input_folder = ARGS[1]
output_path = ARGS[2]
target_label = ARGS[3]
eps = parse(Int16,ARGS[4])

reps = npzread(input_folder * "$(target_label).npy")'
n = size(reps)[2]
removed = round(Int, 1.5*eps)

println("Running PCA filter")
reps_pca, U = pca(reps, 1)
pca_poison_ind = k_lowest_ind(-abs.(mean(reps_pca[1, :]) .- reps_pca[1, :]), round(Int, 1.5*eps))
poison_removed = sum(pca_poison_ind[end-eps+1:end])
clean_removed = removed - poison_removed
@printf("%d poisons removed, %d cleans removed\n", poison_removed, clean_removed)
npzwrite(output_path, pca_poison_ind)

max_stars_repo_name: SewoongLab/backdoor-suite

max_issues_repo_name: SewoongLab/backdoor-suite

max_forks_repo_name: SewoongLab/backdoor-suite


text:
"""
Trait for `ODEOperator` that tells us whether the operator depends on the solution
(including its time derivatives), it is an affine operator that depends on time
or it is a constant operator (affine and time-indepedendent)
"""
abstract type OperatorType end
struct Nonlinear <: OperatorType end
struct Affine  <: OperatorType end
struct Constant  <: OperatorType end

"""
It represents the operator in an implicit ODE, i.e., A(t,u,∂tu) where the
implicit PDE reads A(t,u,∂tu) = 0, when ∂tu is the time derivative of u.
The trait `{C}` determines whether the operator is fully nonlinear, affine
or constant in time.
"""
abstract type ODEOperator{C<:OperatorType} <: GridapType end

"""
It represents an _affine_ operator in an implicit ODE, i.e., an ODE operator of
the form A(t,u,∂tu) = M(t)∂tu + K(t)u + f(t)
"""
const AffineODEOperator = ODEOperator{Affine}
# abstract type AffineODEOperator <: ODEOperator end

"""
It represents a constant operator in an implicit ODE, i.e., an ODE operator of
the form A(t,u,∂tu) = M∂tu + Ku + f
"""
const ConstantODEOperator = ODEOperator{Constant}

# @santiagobadia : I would consider in a future a more general case, in which
# the implicit ODE has an arbitrary order, i.e., A(t,u,u_t, u_nt) = 0.
# get_order(op::ODEOperator) = @notimplemented
# @santiagobadia :
# We probably want to consider second order time derivatives too, i.e.,
# A(t,u,u_t,u_tt) = 0 (wave propagation, elastodynamics) or even more...
"""
Returns the `OperatorType`, i.e., nonlinear, affine, or constant in time
"""
OperatorType(::ODEOperator{C}) where C = C
# OperatorType(::Type{<:ODEOperator{C}}) where C = C

"""
It provides A(t,u,∂tu) for a given (t,u,∂tu)
"""
function residual!(r::AbstractVector,op::ODEOperator,t::Real,u::AbstractVector,u_t::AbstractVector,ode_cache)
  @abstractmethod
end

"""
"""
function allocate_residual(op::ODEOperator,u::AbstractVector,ode_cache)
  @abstractmethod
end

"""
It adds [∂A/∂u](t,u,∂tu) for a given (t,u,∂tu) to a given matrix J
"""
function jacobian!(J::AbstractMatrix,op::ODEOperator,t::Real,u::AbstractVector,u_t::AbstractVector,ode_cache)
  @abstractmethod
  # Add values to J
end

"""
Add the contribution γ*[∂A/∂(∂tu)](u,∂tu) to the Jacobian matrix, where γ
is a scaling coefficient provided by the `ODESolver`, e.g., 1/Δt for Backward
Euler; It represents ∂(δt(u))/∂(u), in which δt(⋅) is the approximation of ∂t(⋅)
in the solver.
"""
function jacobian_t!(J::AbstractMatrix,op::ODEOperator,t::Real,u::AbstractVector,u_t::AbstractVector,dut_u::Real,ode_cache)
  @abstractmethod
  # Add values to J
end

"""
Add the contribution of both `jacobian!` and `jacobian_t!`,i.e.,
[∂A/∂u](t,u,∂tu)+γ*[∂A/∂(∂tu)](u,∂tu)
"""
function jacobian_and_jacobian_t!(J::AbstractMatrix,op::ODEOperator,t::Real,u::AbstractVector,u_t::AbstractVector,dut_u::Real,ode_cache)
  @abstractmethod
  # Add values to J
end

"""
"""
function allocate_jacobian(op::ODEOperator,u::AbstractVector,ode_cache)
  @abstractmethod
end

"""
Allocates the cache data required by the `ODESolution` for a given `ODEOperator`
"""
allocate_cache(op::ODEOperator) = @abstractmethod

#@fverdugo to be used as `cache = update_cache!(cache,op,t)`
update_cache!(cache,op::ODEOperator,t::Real) = @abstractmethod

"""
Tests the interface of `ODEOperator` specializations
"""
function test_ode_operator(op::ODEOperator,t::Real,u::AbstractVector,u_t::AbstractVector)
  cache = allocate_cache(op)
  cache = update_cache!(cache,op,0.0)
  r = allocate_residual(op,u,cache)
  residual!(r,op,t,u,u_t,cache)
  J = allocate_jacobian(op,u,cache)
  jacobian!(J,op,t,u,u_t,cache)
  jacobian_t!(J,op,t,u,u_t,1.0,cache)
  jacobian_and_jacobian_t!(J,op,t,u,u_t,1.0,cache)
  true
end


text:
using Distributions, Gadfly, Random
set_default_plot_size(6.6inch, 3.3inch)

Random.seed!(123)
d = rand(MvNormal([2, 2],[1.0 0.7; 0.7 1.0]), 50)'

pa= plot(x=d[:,1], y=d[:,2], Geom.point, layer(Stat.ellipse, Geom.polygon(preserve_order=true)))
pb= plot(x=d[:,1], y=d[:,2], Geom.point, Geom.ellipse(levels=[0.95, 0.99]))
hstack(pa,pb)


text:
export PseudoTerminal, pseudo_input

# +-----------------+
# | pseudo terminal |
# +-----------------+

mutable struct PseudoTerminal <: REPL.Terminals.UnixTerminal
    term_type::String
    in_stream::Base.IO
    out_stream::Base.IO
    err_stream::Base.IO
    raw::Bool
end

PseudoTerminal(in::Base.IO, out::Base.IO, err::Base.IO) = PseudoTerminal(
    get(ENV, "TERM", Sys.iswindows() ? "" : "dumb"),
    in, out, err,
    false
)

REPL.Terminals.raw!(t::PseudoTerminal, raw::Bool) = (t.raw = raw)

REPL.Terminals.displaysize(::PseudoTerminal) = (24, 80)

pseudo_input(key::String, t=term) = Base.print(t.in_stream, key)


text:
import Base: promote, neg_int, add_int, sub_int, mul_int, afoldl, @_inline_meta,
    power_by_squaring, throw_domerr_powbysq, to_power_type
import Base.Checked: checked_neg, checked_add, checked_sub, checked_mul, checked_abs,
    add_with_overflow, sub_with_overflow, mul_with_overflow, SignedInt,
    throw_overflowerr_binaryop

# convert multi-argument calls into nested two-argument calls
unchecked_add(a, b, c, xs...) = @unchecked (@_inline_meta; afoldl(+, (+)((+)(a, b), c), xs...))
unchecked_sub(a, b, c, xs...) = @unchecked (@_inline_meta; afoldl(-, (-)((-)(a, b), c), xs...))
unchecked_mul(a, b, c, xs...) = @unchecked (@_inline_meta; afoldl(*, (*)((*)(a, b), c), xs...))

checked_add(a, b, c, xs...) = @checked (@_inline_meta; afoldl(+, (+)((+)(a, b), c), xs...))
checked_sub(a, b, c, xs...) = @checked (@_inline_meta; afoldl(-, (-)((-)(a, b), c), xs...))
checked_mul(a, b, c, xs...) = @checked (@_inline_meta; afoldl(*, (*)((*)(a, b), c), xs...))

# passthrough for non-numbers
unchecked_neg(x) = Base.:-(x)
unchecked_add(x, y) = Base.:+(x, y)
unchecked_sub(x, y) = Base.:-(x, y)
unchecked_mul(x, y) = Base.:*(x, y)
unchecked_pow(x, y) = Base.:^(x, y)
unchecked_abs(x) = Base.abs(x)

checked_neg(x) = Base.:-(x)
checked_add(x, y) = Base.:+(x, y)
checked_sub(x, y) = Base.:-(x, y)
checked_mul(x, y) = Base.:*(x, y)
checked_pow(x, y) = Base.:^(x, y)
checked_abs(x) = Base.abs(x)

# promote unmatched number types to same type
unchecked_add(x::Number, y::Number) = unchecked_add(promote(x, y)...)
unchecked_sub(x::Number, y::Number) = unchecked_sub(promote(x, y)...)
unchecked_mul(x::Number, y::Number) = unchecked_mul(promote(x, y)...)
unchecked_mul(x::Number, y::Number) = unchecked_mul(promote(x, y)...)

checked_add(x::Number, y::Number) = checked_add(promote(x, y)...)
checked_sub(x::Number, y::Number) = checked_sub(promote(x, y)...)
checked_mul(x::Number, y::Number) = checked_mul(promote(x, y)...)
checked_pow(x::Number, y::Number) = checked_pow(promote(x, y)...)

# passthrough for same-type numbers that aren't integers
unchecked_add(x::T, y::T) where T <: Number = Base.:+(x, y)
unchecked_sub(x::T, y::T) where T <: Number = Base.:-(x, y)
unchecked_mul(x::T, y::T) where T <: Number = Base.:*(x, y)
unchecked_pow(x::T, y::T) where T <: Number = Base.:^(x, y)

checked_add(x::T, y::T) where T <: Number = Base.:+(x, y)
checked_sub(x::T, y::T) where T <: Number = Base.:-(x, y)
checked_mul(x::T, y::T) where T <: Number = Base.:*(x, y)
checked_pow(x::T, y::T) where T <: Number = Base.:^(x, y)

# core methods
unchecked_neg(x::T) where T <: BitInteger = neg_int(x)
unchecked_add(x::T, y::T) where T <: BitInteger = add_int(x, y)
unchecked_sub(x::T, y::T) where T <: BitInteger = sub_int(x, y)
unchecked_mul(x::T, y::T) where T <: BitInteger = mul_int(x, y)
unchecked_pow(x::T, y::S) where {T <: BitInteger, S <: BitInteger} = power_by_squaring(x, y)
unchecked_abs(x::T) where T <: SignedBitInteger = flipsign(x, x)

checked_neg(x::T) where T <: BitInteger = @checked T(0) - x
function checked_add(x::T, y::T) where T <: BitInteger
    @_inline_meta
    z, b = add_with_overflow(x, y)
    b && throw_overflowerr_binaryop(:+, x, y)
    z
end
function checked_sub(x::T, y::T) where T <: BitInteger
    @_inline_meta
    z, b = sub_with_overflow(x, y)
    b && throw_overflowerr_binaryop(:-, x, y)
    z
end
function checked_mul(x::T, y::T) where T <: BitInteger
    @_inline_meta
    z, b = mul_with_overflow(x, y)
    b && throw_overflowerr_binaryop(:*, x, y)
    z
end
function checked_pow(x::T, y::S) where {T <: BitInteger, S <: BitInteger}
    @_inline_meta
    z, b = pow_with_overflow(x, y)
    b && throw_overflowerr_binaryop(:^, x, y)
    z
end
function checked_abs(x::SignedBitInteger)
    @_inline_meta
    r = @unchecked ifelse(x < 0, -x, x)
    r < 0 && throw(OverflowError(string("checked arithmetic: cannot compute |x| for x = ", x, "::", typeof(x))))
    r
 end

 function pow_with_overflow(x_, p::Integer)
    x = to_power_type(x_)
    if p == 1
        return (copy(x), false)
    elseif p == 0
        return (one(x), false)
    elseif p == 2
        return mul_with_overflow(x, x)
    elseif p < 0
        isone(x) && return (copy(x), false)
        isone(-x) && return (iseven(p) ? one(x) : copy(x), false)
        throw_domerr_powbysq(x, p)
    end
    t = trailing_zeros(p) + 1
    p >>= t
    b = false
    while (t -= 1) > 0
        x, b1 = mul_with_overflow(x, x)
        b |= b1
    end
    y = x
    while p > 0
        t = trailing_zeros(p) + 1
        p >>= t
        while (t -= 1) >= 0
            x, b1 = mul_with_overflow(x, x)
            b |= b1
        end
        y, b1 = mul_with_overflow(y, x)
        b |= b1
    end
    return y, b
end
pow_with_overflow(x::Bool, p::Unsigned) = ((p==0) | x, false)
function pow_with_overflow(x::Bool, p::Integer)
    p < 0 && !x && throw_domerr_powbysq(x, p)
    return (p==0) | x, false
end


text:
const a = split("パトカー","")
const b = split("タクシー","")

for i in 1:length(a)
  print(a[i]*b[i])
end
println("")