text:
#---------------------------------------------------------
# # [LazyGrids ndgrid](@id 1-ndgrid)
#---------------------------------------------------------

# This page explains the `ndgrid` method(s) in the Julia package
# [`LazyGrids`](https://github.com/JuliaArrays/LazyGrids.jl).

# ### Setup

# Packages needed here.

using LazyGrids: ndgrid, ndgrid_array
using LazyGrids: btime, @timeo # not exported; just for timing tests here
using BenchmarkTools: @benchmark
using InteractiveUtils: versioninfo

# ### Overview

# We begin with simple illustrations.

# The Julia method `ndgrid_array` in this package
# is comparable to Matlab's `ndgrid` function.
# It is given a long name here
# to discourage its use,
# because the lazy `ndgrid` version is preferable.
# The package provides `ndgrid_array`
# mainly for testing and timing comparisons.

(xa, ya) = ndgrid_array(1.0:3.0, 1:2)

#
xa

#
ya

# This package provides a "lazy" version of `ndgrid` that appears to the user
# to be the same, but under the hood it is not storing huge arrays.

(xl, yl) = ndgrid(1.0:3.0, 1:2)

#
xl

#
yl

# The following example illustrates the memory savings
# (thanks to Julia's powerful `AbstractArray` type):

(xl, yl) = ndgrid(1:100, 1:200)
(xa, ya) = ndgrid_array(1:100, 1:200)
sizeof(xl), sizeof(xa)

# One can do everything with a lazy array that one would expect
# from a "normal" array, e.g., multiplication and summation:

sum(xl * xl'), sum(xa * xa')

# ### Using lazy `ndgrid`

# Many applications with multiple variables
# involve evaluating functions over a *grid* of values.

# As a simple example (for illustration),
# one can numerically approximate the area of the unit circle
# by sampling that circle over a grid of x,y values,
# corresponding to numerical evaluation of the double integral
# ``∫ ∫ 1_{\{x^2 + y^2 < 1\}} \, dx \, dy``.
# There are many ways to implement this approximation in Julia,
# given a vector of `x` and `y` samples.

Δ = 1/2^10
x = range(-1, stop=1, step=Δ)
y = copy(x)

@inline circle(x::Real, y::Real) = abs2(x) + abs2(y) < 1
@inline circle(xy::NTuple{2}) = circle(xy...)

# The documentation below has many timing comparisons.
# The times in the Julia comments are on a 2017 iMac with Julia 1.6.1;
# the times printed out are whatever server GitHub actions uses.
# Using a trick to [capture output](https://fredrikekre.github.io/Literate.jl/v2/generated/example/#Output-Capturing),
# let's find out:

io = IOBuffer()
versioninfo(io)
split(String(take!(io)), '\n')

# A basic double loop is the C/Fortran way.
# It uses minimal memory (only 48 bytes) but is somewhat slow.

function method0(x,y) # basic double loop
    sum = 0.0
    for x in x, y in y
        sum += circle(x,y)
    end
    return sum * Δ^2
end

area0 = method0(x,y)
t = @benchmark method0($x,$y) # 10.5 ms (3 allocations: 48 bytes)
btime(t)

# The loop version does not look much like the math.
# It often seems natural to think of a grid of x,y values
# and simply sum over that grid, accounting for the grid spacing,
# using a function like this:

area(xx,yy) = sum(circle.(xx,yy)) * Δ^2

# Users coming from Matlab who are unfamiliar with its newer broadcast
# capabilities might use an `ndgrid` of arrays, like in the following code,
# to compute the area.
# But this array approach is much slower and uses much more memory,
# so it does not scale well to higher dimensions.

function area_array(x, y)
    (xa, ya) = ndgrid_array(x, y)
    return area(xa, ya)
end
@assert area_array(x, y) ≈ area0
t = @benchmark area_array($x, $y) # 21.4 ms (11 allocations: 64.57 MiB)
btime(t)

# To be fair, one might have multiple uses of the grids `xa,ya`
# so perhaps they should be excluded from the timing.
# Separating that allocation makes the timing look faster,
# but it still uses a lot of memory,
# both for allocating the grids, and for the `circle.` broadcast
# in the `area` function above:

(xa, ya) = ndgrid_array(x, y)
@assert area(xa, ya) ≈ area0
t = @benchmark area($xa, $ya) # 5.2 ms (7 allocations: 516.92 KiB)
btime(t)

# The main point of this package is to provide
# a lazy version of `ndgrid` that uses minimal memory.

(xl, yl) = ndgrid(x, y)
@assert xl == xa
@assert yl == ya
sizeof(xa), sizeof(xl)

# Now there is essentially no memory overhead for the grids,
# but memory is still used for the `circle.` broadcast.

@assert area(xl,yl) ≈ area0
t = @benchmark area($xl,$yl) # 3.7 ms (7 allocations: 516.92 KiB)
btime(t)

# Furthermore, creating this lazy ndgrid is so efficient
# that we can include its construction time
# and still have performance comparable to the array version
# that had pre-allocated arrays.

function area_lazy(x, y)
    (xl, yl) = ndgrid(x, y)
    return area(xl, yl)
end
@assert area_lazy(x, y) ≈ area0
t = @benchmark area_lazy($x, $y) # 3.7 ms (7 allocations: 516.92 KiB)
btime(t)

# ### More details

# The comparisons below here might be more
# for the curiosity of the package developers
# than for most users...

# One can preallocate memory to store the `circle.` array,
# to avoid additional memory during the area calculation:

out = Array{Float64}(undef, length(x), length(y))
function area!(xx, yy)
    global out .= circle.(xx,yy)
    return sum(out) * Δ^2
end
@assert area!(xl,yl) ≈ area0
t = @benchmark area!(xl,yl) # 4.8 ms (4 allocations: 128 bytes)
btime(t)

# Interestingly, the lazy version is *faster* than the array version,
# presumably because of the overheard of moving data from RAM to CPU:

@assert area!(xa,ya) ≈ area0
t = @benchmark area!(xa,ya) # 6.2 ms (4 allocations: 80 bytes)
btime(t)

# One can avoid allocating the output array by using a loop
# with [CartesianIndices](https://julialang.org/blog/2016/02/iteration):

function area_ci(xx, yy)
    size(xx) == size(yy) || throw("size")
    sum = 0.0
    @inbounds for c in CartesianIndices(xx)
        sum += circle(xx[c], yy[c])
    end
    return sum * Δ^2
end

# With this approach the lazy version is a bit faster
# than the array version:

@assert area_ci(xl,yl) ≈ area0
t = @benchmark area_ci(xl,yl) # 5.2 ms (3 allocations: 48 bytes)
btime(t)

#
@assert area_ci(xa,ya) ≈ area0
t = @benchmark area_ci(xa,ya) # 5.9 ms (3 allocations: 48 bytes)
btime(t)

# Alternatively one can use a linear index for loop,
# that also avoids the extra memory of `circle.` above,
# but is slower, especially for the lazy arrays
# that are optimized for Cartesian indexing:

function area_for2(xx,yy)
    size(xx) == size(yy) || throw("size")
    sum = 0.0
    @inbounds for i in 1:length(xx)
        sum += circle(xx[i], yy[i])
    end
    return sum * Δ^2
end
@assert area_for2(xa, ya) ≈ area0
t = @benchmark area_for2($xa, $ya) # 5.9 ms (3 allocations: 48 bytes)
btime(t)

#
@assert area_for2(xl, yl) ≈ area0
t = @benchmark area_for2($xl, $yl) # 15.4 ms (3 allocations: 48 bytes)
btime(t)

# Some Julia users would
# [recommend using broadcast](https://discourse.julialang.org/t/meshgrid-function-in-julia/48679/25).
# In this case, broadcast is reasonably fast, but still uses a lot of memory
# for the `circle.` output in the simplest implementation.

areab(x,y) = sum(circle.(x,y')) * Δ^2
@assert areab(x,y) ≈ area0
t = @benchmark areab($x,$y) # 11.6 ms (7 allocations: 516.92 KiB)
btime(t)

# Using `zip` can avoid the "extra" memory beyond the grids,
# but seems to have some undesirable overhead,
# presumably because `zip` uses linear indexing:

areaz(xa, ya) = sum(circle, zip(xa,ya)) * Δ^2
@assert areaz(xa, ya) ≈ area0
t = @benchmark areaz($xa, $ya) # 3.9 ms (3 allocations: 48 bytes)
btime(t)

#
@assert areaz(xl, yl) ≈ area0
t = @benchmark areaz($xl, $yl) # 12.2 ms (3 allocations: 48 bytes)
btime(t)

# One can also ensure low memory by using a product iterator,
# but the code starts to look pretty different from the math at this point
# and it is not much faster than broadcast here.

areap(x,y) = sum(circle, Iterators.product(x, y)) * Δ^2
@assert areap(x,y) ≈ area0
t = @benchmark areap($x, $y) # 9.9 ms (3 allocations: 48 bytes)
btime(t)

# ### 3D case

# A 3D example is finding (verifying) the volume of a unit sphere.

sphere(x::Real,y::Real,z::Real) = abs2(x) + abs2(y) + abs2(z) < 1
sphere(r::NTuple) = sum(abs2, r) < 1

# Storing three 3D arrays of size 2049^3 Float64 would take 192GB,
# so already we must greatly reduce the sampling to use either
# `broadcast` or `ndgrid_array`.
# Furthermore, the `broadcast` requires annoying `reshape` steps:

Δc = 1/2^8 # coarse grid
xc = range(-1, stop=1, step=Δc)
yc = xc
zc = xc
nc = length(zc)
3 * nc^3 * 8 / 1024^3 # GB prediction

# Here is broadcast in 3D (yuch!):

vol_br(x,y,z,Δ) = sum(sphere.(
        repeat(x, 1, length(y), length(z)),
        repeat(reshape(y, (1, :, 1)), length(x), 1, length(z)),
        repeat(reshape(z, (1, 1, :)), length(x), length(y), 1),
    )) * Δ^3
vol_br([0.],[0.],[0.],Δc) # warm-up
@timeo vol0 = vol_br(xc,yc,zc,Δc) # 2.7 sec, 3.0 GiB, roughly (4/3)π

function vol_ci(xx, yy, zz, Δ)
    size(xx) == size(yy) == size(zz) || throw("size")
    sum = 0.0
    @inbounds for c in CartesianIndices(xx)
        sum += sphere(xx[c], yy[c], zz[c])
    end
    return sum * Δ^3
end

# Here is the lazy version:
(xlc, ylc, zlc) = ndgrid(xc, yc, zc) # warm-up
@timeo (xlc, ylc, zlc) = ndgrid(xc, yc, zc); # 0.000022 sec (1.8 KiB)

#
vol_ci([0.], [0.], [0.], Δc) # warm-up
@timeo vol_ci(xlc, ylc, zlc, Δc) # 0.2 sec, 1.4 MiB

# Creating the grid of arrays itself is quite slow, even for the coarse grid:
(xac, yac, zac) = ndgrid_array(xc, yc, zc) # warm-up
@timeo (xac, yac, zac) = ndgrid_array(xc, yc, zc) # 1.8 sec 3.0GiB

# Once created, the array version is no faster than the lazy version:
@timeo vol_ci(xac, yac, zac, Δc) # 0.2 seconds (1.1 MiB)

# Using `zip` is more concise (but slower):
vol_zip(xx, yy, zz, Δ) = sum(sphere, zip(xx,yy,zz)) * Δ^3

#
vol_zip([0.], [0.], [0.], Δc) # warm-up
@timeo vol_zip(xlc, ylc, zlc, Δc) # 1.0 sec, 26 MiB

# Using zip for the array version seems to have less overhead
# so that is a potential for future improvement:

@assert vol_zip(xac, yac, zac, Δc) ≈ vol0
@timeo vol_zip(xac, yac, zac, Δc) # 0.19 sec, 16 byte

# Importantly, with the lazy ndgrid now we can return to the fine scale;
# it takes a few seconds, but it is feasible because of the low memory.
z = copy(x)
@timeo (xlf, ylf, zlf) = ndgrid(x, y, z) # 0.000023 sec

#
@timeo vol_ci(xlf, ylf, zlf, Δ) # 12.7 sec, 16 bytes

# I was hoping that with a lazy grid, now we could explore
# higher-dimensional spheres.  But with the current `zip` overhead
# it was too slow, even with coarse grid.
# @timeo (π^2/2, sum(sphere, zip(ndgrid(xc,xc,xc,xc)...)) * Δc^4)

# Probably I need to learn more about stuff like `pairs(IndexCartesian(), A)`
# [e.g., this PR](https://github.com/JuliaLang/julia/pull/38150).
# Another day...


text:
using DD.BDD

@testset "BDD1" begin
    b = bdd(Symbol)
    header!(b, :x)
    header!(b, :y)
    header!(b, :z)
    println(b.headers)
end

@testset "BDD2" begin
    b = bdd(Symbol)
    header!(b, :x)
    header!(b, :y)
    x = var!(b, :x)
    y = var!(b, :y)
    z = and(b, x, y)
    println(todot(b, z))

    tmp = node!(b, header!(b, :x), b.zero, b.one)
    tmp = node!(b, header!(b, :y), b.zero, tmp)
    @test z.id == tmp.id
end

@testset "BDD3" begin
    b = bdd(Symbol)
    x = var!(b, :x)
    y = var!(b, :y)
    z = or(b, x, y)
    println(todot(b, z))

    tmp = node!(b, header!(b, :x), b.zero, b.one)
    tmp = node!(b, header!(b, :y), tmp, b.one)
    @test z.id == tmp.id
end

@testset "BDD4" begin
    b = bdd(Symbol)
    x = var!(b, :x)
    y = var!(b, :y)
    z = var!(b, :z)
    result1 = and(b, or(b, x, y), z)
    println(todot(b, result1))

    result2 = or(b, and(b, not(b, x), not(b, y)), not(b, z))
    println(todot(b, result2))

    result3 = not(b, result1)
    @test result3.id == result2.id
end

ext: jl


text:
import Base: merge, getindex, keys, values, setindex!

export Options, merge, getindex

available_keywords = Set{Symbol}([])

"""
    Options

Type that wraps a dictionary used for options.

### Fields

- `dict` -- the wrapped dictionary
"""
struct Options
    dict::Dict{Symbol,Any}
    Options(args::Pair{Symbol,<:Any}...) = new(Dict{Symbol,Any}(args))
    Options(dict::Dict{Symbol,Any}) = new(dict)
end

"""
    keys(op::Options)

Return the keys of the given options object.

### Input

- `op`    -- options object

### Examples

Obtain the keys of some options with one element:

```jldoctest options_setindex
julia> op = Options(:T=>1.0)
Reachability.Options(Dict{Symbol,Any}(Pair{Symbol,Any}(:T, 1.0)))

julia> collect(keys(op))
1-element Array{Symbol,1}:
 :T
```
"""
keys(op::Options) = keys(op.dict)

"""
    values(op::Options)

Return the values of the given options object.

### Input

- `op`    -- options object

### Examples

Obtain the values of some options with one element:

```jldoctest options_setindex
julia> op = Options(:T=>1.0)
Reachability.Options(Dict{Symbol,Any}(Pair{Symbol,Any}(:T, 1.0)))

julia> collect(values(op))
1-element Array{Any,1}:
 1.0
```
"""
values(op::Options) = values(op.dict)

"""
    merge(op1, opn)

Merges two `Options` objects by just falling back to the wrapped `Dict` fields.
Values are inserted in the order in which the function arguments occur, i.e.,
for conflicting keys a later object overrides a previous value.

### Input

- `op1` -- first options object
- `opn` -- list of options objects
"""
function merge(op1::Options, opn::Options...)::Options
    dict = Dict{Symbol,Any}(op1.dict)
    for i in 1 : length(opn)
        merge!(dict, opn[i].dict)
    end
    return Options(dict)
end

"""
    getindex(op::Options, sym::Symbol)

Returns the value stored for key `sym`.

### Input

- `op`  -- options object
- `sym` -- key
"""
function getindex(op::Options, sym::Symbol)
    return getindex(op.dict, sym)
end

"""
    setindex!(op, value key)

Store the given value at the given key in the options.

### Input

- `op`    -- options object
- `value` -- value
- `key`   -- key

### Examples

Create an empty options object and add an input:

```jldoctest options_setindex
julia> Options()
Reachability.Options(Dict{Symbol,Any}())

julia> op[:T] = 1.0
1.0
```
"""
setindex!(op::Options, value, key) = setindex!(op.dict, value, key)

"""
    validate_solver_options_and_add_default_values!(options)

Validates that the given solver options are supported and adds default values for most
unspecified options.

### Input

- `options` -- an `Options` object, a dictionary of options

Supported options:

- `:verbosity`     -- controls logging output
- `:logfile`       -- name of a log file
- `:mode`          -- main analysis mode
- `:approx_model`  -- model for bloating/continuous time analysis
- `:property`      -- a safety property
- `:δ`             -- time step; alias: `:sampling_time`
- `:N`             -- number of time steps
- `:T`             -- time horizon; alias `:time_horizon`
- `:algorithm`     -- algorithm backend
- `:vars`          -- variables of interest
- `:partition`     -- block partition; elements are 2D vectors containing the
                      start and end of a block
- `:block_types`   -- short hand to set `:block_types_init` and
                      `:block_types_iter`
- `:block_types_init` -- set type for the approximation of the initial states
                         for each block
- `:block_types_iter` -- set type for the approximation of the states ``X_k``,
                         ``k>0``, for each block
- `:ε`             -- short hand to set `:ε_init` and `:ε_iter`
- `:set_type`      -- short hand to set `:set_type_init` and `:set_type_iter`
- `:ε_init`        -- error bound for the approximation of the initial states
                      (during decomposition)
- `:set_type_init` -- set type for the approximation of the initial states
                      (during decomposition)
- `:ε_iter`        -- error bound for the approximation of the states ``X_k``,
                      ``k>0``
- `:set_type_iter` -- set type for the approximation of the states ``X_k``,
                      ``k>0``
- `:ε_proj`        -- error bound for the approximation of the states during
                      projection
- `:set_type_proj` -- set type for the approximation of the states during
                      projection
- `:lazy_expm`     -- switch for using lazy matrix exponential all the time
- `:assume_sparse` -- switch for sparse matrices
- `:pade_expm`     -- switch for using Pade approximant method
- `:lazy_X0`       -- switch for keeping the initial states a lazy set
- `:lazy_sih`      -- switch for using a lazy symmetric interval hull during the
                      discretization
- `:template_directions`       -- short hand to set `template_directions_init`
                                  and `template_directions_iter`
- `:template_directions_init`  -- directions to use for the approximation of the
                                  initial states (during decomposition)
- `:template_directions_iter`  -- directions to use for the approximation of the
                                  states ``X_k``, ``k>0``, for each block
- `:coordinate_transformation` -- coordinate transformation method
- `:assume_homogeneous`        -- switch for ignoring inputs
- `:projection_matrix`         -- projection matrix
- `:project_reachset`          -- switch for applying projection
- `:eager_checking`            -- switch for early terminating property checks
- `:lazy_inputs_interval`      -- length of interval in which the inputs are
                                  handled as a lazy set (``-1`` for 'never');
                                  generally may also be a predicate over indices
- `:lazy_expm_discretize`      -- switch to use lazy matrix exponential in the
                                  discretization phase (see also `:lazy_expm`)
- `:plot_vars`     -- variables for projection and plotting;
                      alias: `:output_variables`
- `:n`             -- system's dimension

Internal options (inputs are ignored or even illegal):

- `:blocks`        -- list of all interesting block indices in the partition

We add default values for almost all undefined options, i.e., modify the input
options. The benefit is that the user can print the options that were actually
used. For supporting aliases, we create another copy that is actually used where
we only keep those option names that are used internally.
"""
function validate_solver_options_and_add_default_values!(options::Options)::Options
    global LOGGER

    dict = options.dict
    options_copy = Options()
    dict_copy = options_copy.dict

    # first read the verbosity option and set global log level accordingly
    LOGGER = haskey(dict, :verbosity) ?
        configure_logger(dict[:verbosity]) :
        configure_logger()
    check_aliases_and_add_default_value!(dict, dict_copy, [:verbosity], getlevel(LOGGER))
    if haskey(dict, :logfile)
        add_file_logger(dict[:logfile])
        check_aliases!(dict, dict_copy, [:logfile])
    end

    # check for aliases and use default values for unspecified options
    check_aliases_and_add_default_value!(dict, dict_copy, [:mode], "reach")
    check_aliases_and_add_default_value!(dict, dict_copy, [:approx_model], "forward")
    check_aliases_and_add_default_value!(dict, dict_copy, [:property], nothing)
    check_aliases_and_add_default_value!(dict, dict_copy, [:algorithm], "explicit")
    check_aliases_and_add_default_value!(dict, dict_copy, [:vars], 1:options.dict[:n])
    check_aliases_and_add_default_value!(dict, dict_copy, [:lazy_expm], false)
    check_aliases_and_add_default_value!(dict, dict_copy, [:assume_sparse], false)
    check_aliases_and_add_default_value!(dict, dict_copy, [:pade_expm], false)
    check_aliases_and_add_default_value!(dict, dict_copy, [:lazy_X0], false)
    check_aliases_and_add_default_value!(dict, dict_copy, [:lazy_sih], false)
    check_aliases_and_add_default_value!(dict, dict_copy, [:coordinate_transformation], "")
    check_aliases_and_add_default_value!(dict, dict_copy, [:assume_homogeneous], false)
    check_aliases_and_add_default_value!(dict, dict_copy, [:projection_matrix], nothing)
    check_aliases_and_add_default_value!(dict, dict_copy, [:project_reachset], dict_copy[:projection_matrix] == nothing)
    check_aliases_and_add_default_value!(dict, dict_copy, [:eager_checking], true)
    check_aliases_and_add_default_value!(dict, dict_copy, [:lazy_expm_discretize],
                                         dict_copy[:lazy_expm])
    check_aliases_and_add_default_value!(dict, dict_copy, [:n], nothing)

    # special options: δ, N, T
    check_and_add_δ_N_T!(dict, dict_copy)

    # special option: plot_vars
    check_and_add_plot_vars!(dict, dict_copy)

    # special options: ε, ε_init, ε_iter, ε_proj,
    #                  set_type, set_type_init, set_type_iter, set_type_proj,
    #                  template_directions, template_directions_init,
    #                  template_directions_iter
    check_and_add_approximation!(dict, dict_copy)

    # special options: partition, block_types, block_types_init, block_types_iter
    check_and_add_partition_block_types!(dict, dict_copy)

    # special option: lazy_inputs_interval
    check_and_add_lazy_inputs_interval!(dict, dict_copy)

    # validate that all input keywords are recognized
    check_valid_option_keywords(dict)

    # validation
    for kv_pair in dict_copy
        key::Symbol = kv_pair.first
        value = kv_pair.second

        # define type/domain constraints for each known key
        domain_constraints = (v  ->  true)
        if key == :verbosity
            expected_type = Union{String, Int}
        elseif key == :logfile
            expected_type = String
        elseif key == :mode
            expected_type = String
            domain_constraints = (v::String  ->  v in ["reach", "check"])
            if value == "check" && dict_copy[:property] == nothing
                error("No property has been defined.")
            end
        elseif key == :approx_model
            expected_type = String
            domain_constraints = (v::String  ->  v in ["forward", "backward",
                                                       "firstorder", "nobloating"])
        elseif key == :property
            expected_type = Union{Property, Void}
        elseif key == :δ
            expected_type = Float64
            domain_constraints = (v::Float64  ->  v > 0.)
        elseif key == :N
            expected_type = Int
            domain_constraints = (v::Int  ->  v > 0)
        elseif key == :T
            expected_type = Float64
            domain_constraints = (v::Float64  ->  v > 0.)
        elseif key == :algorithm
            expected_type = String
            domain_constraints = (v::String  ->  v in ["explicit", "wrap"])
        elseif key == :vars
            expected_type = AbstractVector{Int}
            domain_constraints =
                (v::AbstractVector{Int}  ->  length(v) > 0 && all(x -> x > 0, v))
        elseif key == :partition
            expected_type = AbstractVector{<:AbstractVector{Int}}
        elseif key == :block_types
            expected_type = Union{Void,
                Dict{Type{<:LazySet}, AbstractVector{<:AbstractVector{Int}}}}
        elseif key == :block_types_init
            expected_type =
                Dict{Type{<:LazySet}, AbstractVector{<:AbstractVector{Int}}}
        elseif key == :block_types_iter
            expected_type =
                Dict{Type{<:LazySet}, AbstractVector{<:AbstractVector{Int}}}
        elseif key == :blocks
            expected_type = AbstractVector{Int}
            domain_constraints =
                (v::AbstractVector{Int}  ->  length(v) > 0 && all(x -> x > 0, v))
        elseif key == :ε
            expected_type = Float64
            domain_constraints = (v  ->  v > 0.)
        elseif key == :ε_init
            expected_type = Float64
            domain_constraints = (v  ->  v > 0.)
        elseif key == :ε_iter
            expected_type = Float64
            domain_constraints = (v  ->  v > 0.)
        elseif key == :ε_proj
            expected_type = Float64
            domain_constraints = (v  ->  v > 0.)
        elseif key == :set_type
            expected_type = Union{Type{HPolygon}, Type{Hyperrectangle},
                                  Type{LazySets.Interval}}
        elseif key == :set_type_init
            expected_type = Union{Type{HPolygon}, Type{Hyperrectangle},
                                  Type{LazySets.Interval}}
        elseif key == :set_type_iter
            expected_type = Union{Type{HPolygon}, Type{Hyperrectangle},
                                  Type{LazySets.Interval}}
        elseif key == :set_type_proj
            expected_type = Union{Type{HPolygon}, Type{Hyperrectangle},
                                  Type{LazySets.Interval}}
        elseif key == :lazy_expm
            expected_type = Bool
        elseif key == :assume_sparse
            expected_type = Bool
        elseif key == :pade_expm
            expected_type = Bool
        elseif key == :lazy_X0
            expected_type = Bool
        elseif key == :lazy_sih
            expected_type = Bool
        elseif key == :template_directions
            expected_type = Symbol
            domain_constraints = (v::Symbol  ->  v in [:box, :oct, :boxdiag,
                                                       :nothing])
        elseif key == :template_directions_init
            expected_type = Symbol
            domain_constraints = (v::Symbol  ->  v in [:box, :oct, :boxdiag,
                                                       :nothing])
        elseif key == :template_directions_iter
            expected_type = Symbol
            domain_constraints = (v::Symbol  ->  v in [:box, :oct, :boxdiag,
                                                       :nothing])
        elseif key == :coordinate_transformation
            expected_type = String
            domain_constraints = (v::String  ->  v in ["", "schur"])
        elseif key == :assume_homogeneous
            expected_type = Bool
        elseif key == :projection_matrix
            expected_type = Union{AbstractMatrix, Void}
        elseif key == :project_reachset
            expected_type = Bool
        elseif key == :eager_checking
            expected_type = Bool
        elseif key == :lazy_inputs_interval
            expected_type = Union{Int, Function, Void}
            domain_constraints = (v  ->  !(v isa Int) || v >= -1)
        elseif key == :lazy_expm_discretize
            expected_type = Bool
            if !value && dict_copy[:lazy_expm]
                error("cannot use option $(:lazy_expm) with deactivated " *
                      "option $(:lazy_expm_discretize)")
            end
        elseif key == :plot_vars
            expected_type = Vector{Int}
            domain_constraints = (v::Vector{Int}  ->  length(v) == 2)
        elseif key == :n
            expected_type = Int
            domain_constraints = (v::Int  ->  v > 0)
        else
            error(get_unrecognized_key_message(key))
        end

        # check value type
        if !(value isa expected_type)
            error("Option :$key must be of '$expected_type' type.")
        end
        # check value domain constraints
        if !domain_constraints(value)
            error("$value is not a valid value for option :$key.")
        end
    end

    return options_copy
end

"""
Handling of the special trio `:δ`, `:N`, `:T`.
Usually two of them should be defined and the third one is automatically inferred.
If three of them are defined, they must be consistent.
If only `:T` is defined, we use `:N = 100`.

INPUT:

- `dict`      -- a dictionary of options
- `dict_copy` -- a copy of the dictionary of options for internal names

Supported options:
- `:δ` (time step)
- alias: `:sampling_time`
- `:N` (number of time steps)
- `:T` (time horizon)
- alias: `:time_horizon`
"""
function check_and_add_δ_N_T!(dict::Dict{Symbol,Any}, dict_copy::Dict{Symbol,Any})
    δ_aliases = [:δ, :sampling_time]
    T_aliases = [:T, :time_horizon]
    check_aliases!(dict, dict_copy, δ_aliases)
    check_aliases!(dict, dict_copy, [:N])
    check_aliases!(dict, dict_copy, T_aliases)

    defined = 0
    if haskey(dict_copy, :δ)
        value = dict_copy[:δ]
        if !(value isa Float64 && value > 0.)
            error("$value is not a valid value for option $δ_aliases.")
        end
        defined += 1
    end
    if haskey(dict_copy, :N)
        value = dict_copy[:N]
        if !(value isa Int && value > 0)
            error("$value is not a valid value for option :N.")
        end
        defined += 1
    end
    if haskey(dict_copy, :T)
        value = dict_copy[:T]
        if !(value isa Float64 && value > 0.)
            error("$value is not a valid value for option $T_aliases.")
        end
        defined += 1
    end

    if defined == 1 && haskey(dict_copy, :T)
        N = 100
        dict[:N] = N
        dict_copy[:N] = N
        δ = dict_copy[:T] / dict_copy[:N]
        dict[:δ] = δ
        dict_copy[:δ] = δ
    elseif defined == 2
        if !haskey(dict_copy, :δ)
            δ = dict_copy[:T] / dict_copy[:N]
            dict[:δ] = δ
            dict_copy[:δ] = δ
        end
        if !haskey(dict_copy, :N)
            N = (Int)(ceil(dict_copy[:T] / dict_copy[:δ]))
            dict[:N] = N
            dict_copy[:N] = N
        end
        if !haskey(dict_copy, :T)
            T = dict_copy[:N] * dict_copy[:δ]
            dict[:T] = T
            dict_copy[:T] = T
        end
    elseif defined == 3
        N_computed = (Int)(ceil(dict_copy[:T] / dict_copy[:δ]))
        if N_computed != dict_copy[:N]
            error("Values for :δ ($(dict_copy[:δ])), :N ($(dict_copy[:N])), and :T ($(dict_copy[:T])) were defined, but they are inconsistent.")
        end
    else
        error("Need two of the following options: :δ, :N, :T (or at least :T and using default values)")
    end
end

"""
Handling of the special option `:plot_vars`.

### Input

- `dict`      -- dictionary of options
- `dict_copy` -- copy of the dictionary of options for internal names

### Notes:

If no value is given, we take the first two dimensions from `:vars`.
If `:vars` has only one element, we use time for the other.
"""
function check_and_add_plot_vars!(dict::Dict{Symbol,Any},
                                 dict_copy::Dict{Symbol,Any})
    check_aliases!(dict, dict_copy, [:plot_vars, :output_variables])
    if !haskey(dict_copy, :plot_vars)
        vars = dict_copy[:vars]
        if length(vars) == 1
            plot_vars = [0, vars[1]]
        else
            plot_vars = [vars[1], vars[2]]
        end
        dict_copy[:plot_vars] = plot_vars
    else
        # sanity check
        vars = dict_copy[:vars]
        for i in dict_copy[:plot_vars]
            if i != 0 && i ∉ vars
                error("$(dict_copy[:plot_vars]) is not a subset of the variables $vars")
            end
        end
    end
end

"""
    check_and_add_approximation!(dict::Dict{Symbol,Any},
                                 dict_copy::Dict{Symbol,Any})

Handling of the special options `:ε` and `:set_type` resp. the subcases
`:ε_init`, `:ε_iter`, `:ε_proj`, `:set_type_init`, `:set_type_iter`,
`:set_type_proj`, `:template_directions`, `:template_directions_init`, and
`:template_directions_iter`.

### Input

- `dict`      -- dictionary of options
- `dict_copy` -- copy of the dictionary of options for internal names

### Notes:

`:ε` and `:set_type`, if defined, are used as fallbacks (if the more specific
options are undefined).
If both `:ε_init` and `:set_type_init` are defined, they must be consistent.
"""
function check_and_add_approximation!(dict::Dict{Symbol,Any},
                                      dict_copy::Dict{Symbol,Any})
    # fallback options
    check_aliases!(dict, dict_copy, [:ε])
    check_aliases!(dict, dict_copy, [:set_type])
    check_aliases!(dict, dict_copy, [:template_directions])
    ε = haskey(dict, :ε) ? dict[:ε] : Inf

    # fallback set type
    if haskey(dict, :set_type)
        # use the provided set type
        set_type = dict[:set_type]
    elseif ε < Inf
        # use polygons
        set_type = HPolygon
    elseif !haskey(dict, :partition)
        # use intervals
        set_type = Interval
        set_type_proj = Interval
    else
        # use hyperrectangles
        set_type = Hyperrectangle
    end

    ε_init = (haskey(dict, :set_type_init) && dict[:set_type_init] == HPolygon) ||
             (!haskey(dict, :set_type_init) && set_type == HPolygon) ? ε : Inf
    check_aliases_and_add_default_value!(dict, dict_copy, [:ε_init], ε_init)

    set_type_init = dict_copy[:ε_init] < Inf ? HPolygon : set_type
    check_aliases_and_add_default_value!(dict, dict_copy, [:set_type_init], set_type_init)

    template_directions_init = haskey(dict, :template_directions_init) ?
        dict[:template_directions_init] : haskey(dict, :template_directions) ?
        dict[:template_directions] : :nothing
    check_aliases_and_add_default_value!(dict, dict_copy,
        [:template_directions_init], template_directions_init)

    ε_iter = (haskey(dict, :set_type_iter) && dict[:set_type_iter] == HPolygon) ||
             (!haskey(dict, :set_type_iter) && set_type == HPolygon) ? ε : Inf
    check_aliases_and_add_default_value!(dict, dict_copy, [:ε_iter], ε_iter)

    set_type_iter = dict_copy[:ε_iter] < Inf ? HPolygon : set_type
    check_aliases_and_add_default_value!(dict, dict_copy, [:set_type_iter], set_type_iter)

    template_directions_iter = haskey(dict, :template_directions_iter) ?
        dict[:template_directions_iter] : haskey(dict, :template_directions) ?
        dict[:template_directions] : :nothing
    check_aliases_and_add_default_value!(dict, dict_copy,
        [:template_directions_iter], template_directions_iter)

    ε_proj = (haskey(dict, :set_type_proj) && dict[:set_type_proj] == HPolygon) ||
             (!haskey(dict, :set_type_proj) && set_type == HPolygon) ? ε :
             min(dict_copy[:ε_init], dict_copy[:ε_iter], ε)
    check_aliases_and_add_default_value!(dict, dict_copy, [:ε_proj], ε_proj)

    set_type_proj = dict_copy[:ε_proj] < Inf ? HPolygon :
        set_type == Interval && length(dict_copy[:plot_vars]) > 1 ?
        Hyperrectangle : set_type
    check_aliases_and_add_default_value!(dict, dict_copy, [:set_type_proj], set_type_proj)

    @assert (dict_copy[:ε_init] == Inf || dict_copy[:set_type_init] == HPolygon) &&
        (dict_copy[:ε_iter] == Inf || dict_copy[:set_type_iter] == HPolygon) &&
        (dict_copy[:ε_proj] == Inf || dict_copy[:set_type_proj] == HPolygon) (
            "ε-close approximation is only supported with the HPolygon set type")
end

"""
    check_and_add_partition_block_types!(dict::Dict{Symbol,Any},
                                         dict_copy::Dict{Symbol,Any})

Handling of the special options `:partition` and `:block_types` resp. the
subcases `:block_types_init` and `:block_types_iter`.

### Input

- `dict`      -- dictionary of options
- `dict_copy` -- copy of the dictionary of options for internal names
"""
function check_and_add_partition_block_types!(dict::Dict{Symbol,Any},
                                              dict_copy::Dict{Symbol,Any})
    check_aliases!(dict, dict_copy, [:partition])
    if !haskey(dict_copy, :partition)
        partition = [[i] for i in 1:dict[:n]]
        dict_copy[:partition] =  partition
    end

    block_types = nothing
    if haskey(dict, :block_types)
        for (key, value) in dict[:block_types]
            @assert key <: LazySet "the keys of the `partition` dictionary should be lazy sets"
            @assert typeof(value) <: AbstractVector{<:AbstractVector{Int}} "the keys of the `partition` dictionary should be vectors of vectors"
        end
        block_types = convert(Dict{Type{<:LazySet}, AbstractVector{<:AbstractVector{Int}}}, dict[:block_types])
    elseif haskey(dict, :set_type)
        block_types = Dict{Type{<:LazySet}, AbstractVector{<:AbstractVector{Int}}}(dict[:set_type] => copy(dict_copy[:partition]))
    end
    check_aliases_and_add_default_value!(dict, dict_copy, [:block_types], block_types)
    dict_copy[:block_types] = block_types

    block_types_init = haskey(dict, :block_types_init) ?
        dict[:block_types_init] :
        block_types != nothing ? block_types :
            Dict{Type{<:LazySet}, AbstractVector{<:AbstractVector{Int}}}(
                dict_copy[:set_type_init] => copy(dict_copy[:partition])
            )
    check_aliases_and_add_default_value!(dict, dict_copy, [:block_types_init],
                                         block_types_init)

    block_types_iter = haskey(dict, :block_types_iter) ?
        dict[:block_types_iter] :
        block_types != nothing ? block_types :
            Dict{Type{<:LazySet}, AbstractVector{<:AbstractVector{Int}}}(
                dict_copy[:set_type_iter] => copy(dict_copy[:partition])
            )
    check_aliases_and_add_default_value!(dict, dict_copy, [:block_types_iter],
                                         block_types_iter)

    # TODO add check that arguments are consistent

    # compute :blocks
    if haskey(dict, :blocks)
        error("option :blocks is not allowed")
    end
    dict_copy[:blocks] = compute_blocks(dict_copy[:vars], dict_copy[:partition])
end

function compute_blocks(vars, partition)
    blocks = Vector{Int}()
    sizehint!(blocks, length(vars))
    next = 0
    var_idx = 1
    for (i, block) in enumerate(partition)
        next += length(block)
        if vars[var_idx] <= next
            push!(blocks, i)
            var_idx += 1
            while var_idx <= length(vars) && vars[var_idx] <= next
                var_idx += 1
            end
            if var_idx > length(vars)
                break
            end
        end
    end
    @assert var_idx == length(vars) + 1
    sizehint!(blocks, length(blocks))
    return blocks
end

"""
    check_and_add_lazy_inputs_interval!(dict::Dict{Symbol,Any},
                                        dict_copy::Dict{Symbol,Any})

Handling of the special option `:lazy_inputs_interval`.

### Input

- `dict`      -- dictionary of options
- `dict_copy` -- copy of the dictionary of options for internal names

### Notes

Originally, this option was used to overapproximate a lazy set of inputs at
every multiple of ``k`` steps, where ``k`` was the controllable option; hence
the term `interval` in the name.
Meanwhile, this option was generalized to a predicate over integers.
The predicate returns `true` iff the lazy set should be overapproximated.
However, we still support index numbers as input, which are then translated into
a predicate.

The default value for this option is `nothing` or the index ``0``, which results
in overapproximation in each iteration.
Note that internally this has a different effect than passing the predicate that
always returns `true` because functions cannot be checked for equality.
Thus this option should simply not be set when not wanting this behavior (or the
value ``0`` should be used).

The input ``-1`` is interpreted as the predicate that always returns `false`.
"""
function check_and_add_lazy_inputs_interval!(dict::Dict{Symbol,Any},
                                             dict_copy::Dict{Symbol,Any})
    check_aliases!(dict, dict_copy, [:lazy_inputs_interval])
    if haskey(dict_copy, :lazy_inputs_interval)
        if dict_copy[:lazy_inputs_interval] == -1
            dict_copy[:lazy_inputs_interval] = (k -> false)
        elseif dict_copy[:lazy_inputs_interval] == 0
            dict_copy[:lazy_inputs_interval] = nothing
        elseif dict_copy[:lazy_inputs_interval] isa Int
            m = dict_copy[:lazy_inputs_interval]
            dict_copy[:lazy_inputs_interval] = (k -> k % m == 0)
        end
    else
        dict_copy[:lazy_inputs_interval] = nothing
    end
end

"""
    check_aliases!(dict, dict_copy, aliases)

This function has several purposes:

- translate aliases to the option that is used internally
- check that not several aliases were used at the same time

INPUT:

- `dict` -- a dictionary of options
- `dict_copy` -- a copy of the dictionary of options for internal names
- `aliases` -- option aliases; the first name is the one we use internally
"""
function check_aliases!(dict::Dict{Symbol,Any}, dict_copy::Dict{Symbol,Any}, aliases::Vector{Symbol})
    # find aliases and check consistency in case several aliases have been used
    print_warning::Bool = false
    for alias in aliases
        # add alias to global keyword set
        push!(available_keywords, alias)

        if haskey(dict, alias)
            # alias was used
            if haskey(dict_copy, aliases[1])
                # several aliases were used
                if (dict_copy[aliases[1]] != dict[alias])
                    error("Several option aliases were used with different values for aliases $aliases.")
                else
                    print_warning = true
                end
            else
                dict_copy[aliases[1]] = dict[alias]
            end
        end
    end
    if print_warning
        warn("Several option aliases were used for aliases $aliases.")
    end
end

"""
    check_aliases_and_add_default_value!(dict, dict_copy, aliases, default_value, [modify_dict])

This function has several purposes:
- translate aliases to the option that is used internally
- check that not several aliases were used at the same time
- assign default values for undefined options

INPUT:

- `dict`          -- a dictionary of options
- `dict_copy`     -- a copy of the dictionary of options for internal names
- `aliases`       -- option aliases; the first name is the one we use internally
- `default_value` -- the default value for the option
- `modify_dict`   -- (optional, default: `false`) indicates if `dict` should be
                     modified
"""
function check_aliases_and_add_default_value!(dict::Dict{Symbol,Any}, dict_copy::Dict{Symbol,Any}, aliases::Vector{Symbol}, default_value::Any, modify_dict::Bool=false)
    check_aliases!(dict, dict_copy, aliases)

    if !haskey(dict_copy, aliases[1])
        # no alias and no value
        # TODO<notification> print message to user for each default option, depending on verbosity level
        if (modify_dict)
            dict[aliases[1]] = default_value
        end
        dict_copy[aliases[1]] = default_value
    end
end

"""
    check_valid_option_keywords(dict)

Validate that all input keywords are recognized.

### Input

- `dict` -- input dictionary
"""
function check_valid_option_keywords(dict)
    for kv_pair in dict
        key::Symbol = kv_pair.first
        if !in(key, available_keywords)
            error(get_unrecognized_key_message(key))
        end
    end
end

"""
    get_unrecognized_key_message(key)

Create an error message for an unrecognized option key.

### Input

- `key` -- unrecognized option key

### Output

The error message.
"""
function get_unrecognized_key_message(key)
    return "Unrecognized option '$key' found. See " *
        "`keys(Reachability.available_keywords.dict)` for all valid keywords."
end


text:
function euler_residuals(model, s::ListOfPoints, x::Vector{<:ListOfPoints}, dr, dprocess, parms::SVector; keep_J_S=false, set_dr=true) #, jres=nothing, S_ij=nothing)
    #
    if set_dr ==true
      set_values!(dr,x)
    end

    N_s = length(s) # Number of gris points for endo_var
    n_s = length(s[1]) # Number of states
    n_x = length(x[1][1]) # Number of controls

    n_ms = length(x)  # number of exo states today
    n_mst = n_inodes(dprocess,1)  # number of exo states tomorrow
    d = length(s[1])

    # TODO: allocate properly...
    res = deepcopy(x)
    for i_m=1:length(res)
        res[i_m][:] *= 0.0
    end

    if keep_J_S
        jres = zeros((n_ms,n_mst,N_s,n_x,n_x))
        fut_S = zeros((n_ms,n_mst,N_s,n_s))
        J_ij = Vector{SMatrix{n_x,n_x,Float64,n_x*n_x}}[to_LOJ(jres[i,j,:,:,:]) for i=1:size(jres,1), j=1:size(jres,2)]
        S_ij = Vector{Point{d}}[to_LOP(fut_S[i,j,:,:]) for i=1:size(fut_S,1), j=1:size(fut_S,2)]
    end

    for i_ms in 1:n_ms
        m = node(Point, dprocess,i_ms)
        for I_ms in 1:n_mst
            M = inode(Point, dprocess, i_ms, I_ms)
            w = iweight(dprocess, i_ms, I_ms)
            S = transition(model, m, s, x[i_ms], M, parms)
            X = dr(i_ms, I_ms, S)
            if keep_J_S==true
                rr, rr_XM = arbitrage(model,Val{(0,6)},m,s,x[i_ms],M,S,X,parms)
                J_ij[i_ms,I_ms][:] = w*rr_XM
                S_ij[i_ms,I_ms][:] = S
            else
                rr = arbitrage(model, m,s,x[i_ms],M,S,X,parms)
            end
            res[i_ms][:] += w*rr
        end
    end
    if keep_J_S
        return res,J_ij,S_ij
    else
        return res
    end
end

####################################
# one in place filtering step: M.r #
####################################

function d_filt_dx!(Π_i::Vector{ListOfPoints{n_x}}, π_i::Vector{ListOfPoints{n_x}}, M_ij, S_ij::Matrix{ListOfPoints{d}}, dumdr::Dolo.CachedDecisionRule) where n_x where d
    n_m,n_mt = size(M_ij)
    Dolo.set_values!(dumdr,π_i)
    for i in 1:n_m
        Π_i[i][:] *= 0
        for j in 1:n_mt
            A = M_ij[i,j]
            B = dumdr(i,j,S_ij[i,j])
            addmul!(Π_i[i],A,B)
        end
    end
end

function preinvert!(π_i, D_i, J_ij, S_ij)

    for i=1:length(D_i)
        invert!(D_i[i])
    end
    for i=1:size(J_ij,1)
        for j=1:size(J_ij,2)
            premult!(D_i[i],J_ij[i,j])
        end
    end
    for i=1:length(D_i)
        premult!(D_i[i], π_i[i])
    end

    return π_i, J_ij, S_ij

end

function invert_jac(R_i, M_ij, S_ij, dumdr; tol::Float64=1e-10,
            maxit::Int=1000, verbose::Bool=false)

    dprocess = dumdr.process

    lam = -1.0
    lam_max = -1.0

    if verbose==true
    print("Starting inversion")
    end

    verbose && println(repeat("-", 35))
    verbose && @printf "%-6s%-12s%-5s\n" "err" "gain" "gain_max"
    verbose && println(repeat("-", 35))
    precomputed=false

    π_i = deepcopy(R_i)
    Π_i = deepcopy(π_i)
    tot = deepcopy(π_i) # should premult by D_i

    err_0 = maxabs(π_i)
    err = 10000

    errors = [err_0]

    it = 0
    while it<maxit && err>tol
        it +=1
        d_filt_dx!(Π_i, π_i, M_ij, S_ij, dumdr)
        π_i = Π_i
        err = maxabs(Π_i)
        push!(errors, err)
        lam = err/err_0
        lam_max = max(lam_max, lam)
        err_0 = err
        verbose && @printf "%-6f%-12.2e%-5.2e\n" err lam lam_max
        for i=1:length(tot)
            tot[i] += Π_i[i]
        end
    end

    return tot, it, lam, errors

end

####
# The LinearThing represents
# - (I-M) with (`L*...``)
# - M with (`mul(L,...)``)
####

mutable struct LinearThing
   M_ij
   S_ij
   I
   counter::Int
end

import LinearAlgebra
import LinearAlgebra: *, mul!
using LinearAlgebra

LinearThing(M,S,I) = LinearThing(M,S,I,0)

eltype(L::LinearThing) = Float64
function shape(L::LinearThing)
   n_m = size(L.M_ij,1)
   N = size(L.M_ij[1,1],1)
   n_x = size(L.M_ij[1,1][1],1)
   return (n_x, N, n_m)
end
size(L::LinearThing,d) = prod(shape(L))

function *(L::LinearThing,x::AbstractVector{<:AbstractVector{Point{n_x}}}) where n_x
   xx = deepcopy(x)
   Dolo.d_filt_dx!(xx, x, L.M_ij, L.S_ij, L.I)
   L.counter += 1
   return xx
end

function *(L::LinearThing,m::AbstractArray{Float64, 3})
   n_x,N,n_m = size(m)
   # TODO remove copy there
   x = [copy(reshape(reinterpret(SVector{n_x,Float64},vec(m[:,:,i])),(N,))) for i=1:n_m]
   y = deepcopy(x)
   xx = L*y
   rr = [reshape(reinterpret(Float64, vec(xx[i])), (n_x,N)) for i=1:length(xx)]
   rrr = cat(rr...; dims=3)
   return reshape(rrr, n_x,N,n_m)
end

function *(L::LinearThing,v::AbstractVector{Float64})
   m = copy(v)
   sh = shape(L)
   n_x = sh[1]
   vv = reshape(reinterpret(Point{n_x},vec(m)),(sh[2],sh[3]))
   # x = [view(vv,:,i) for i=1:sh[3]]
   x = [vv[:,i] for i=1:sh[3]]
   y = x-L*x
   mm = reshape(m, sh...)
   mmm = mm-L*mm
   return mmm[:]
end

function mul(L::LinearThing,v::AbstractVector{Float64})
   m = copy(v)
   sh = shape(L)
   mm = reshape(m, sh...)
   mmm = L*mm
   return mmm[:]
end

function A_mul_B!(w::AbstractVector{Float64},L::LinearThing,v::AbstractVector{Float64})
   w[:] = L*v
   return w
end

function mul!(w,L::LinearThing,v)
   w[:] = L*v
   return w
end

####

mutable struct ImprovedTimeIterationResult
  dr::AbstractDecisionRule
  N::Int
  f_x::Float64
  d_x::Float64
  x_converged::Bool
  complementarities::Bool
  # Time_search::
  tol::Float64
  Lambda::Float64
  N_invert::Float64
  N_search::Float64
  radius::Float64
  trace_data
end

converged(r::ImprovedTimeIterationResult) = r.x_converged

function Base.show(io::IO, r::ImprovedTimeIterationResult)
  @printf io "Results of Improved Time Iteration Algorithm\n"
  @printf io " * Number of iterations: %s\n" string(r.N)
  @printf io " * Complementarities: %s\n" string(r.complementarities)
  @printf io " * Decision Rule type: %s\n" string(typeof(r.dr))
  @printf io " * Convergence: %s\n" converged(r)
  @printf io " * Contractivity: %s\n" string(r.Lambda)
  @printf io "   * |x - x'| < %.1e: %s\n" r.tol r.x_converged
end


text:
f() = g()
g() = println.(stacktrace())

f()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module CS2_1
using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test
using Optim
#using Evolutionary
using JLD2

## CS2_1  - Case-study example: Monitoring air data sensor faults
#           Robust least order LTI synthesis
println("Case study CS2_1")

## Part 1 - Model setup
# load matrices of the aircraft multiple-model SYSACM, actuator model ACT,
# output-feedback gain K, percentages of mass variations massi

cd(joinpath(pkgdir(FaultDetectionTools), "Examples"))
@load "cs2data.jld2"

N = size(SYSACSM,1)
p, m = size(SYSACSM[1].D)
T = Float64

# build minimal realizations of AC-models with actuator faults
# set dimensions
nom = 7;             # nominal system
# set sensor indices and set dimensions of inputs
#     [ AoA VCAS ]
sen = [  2,    4  ];
md = 2; mu = m-md; mf = length(sen); n = maximum(order.(SYSACSM));

idm = eye(p);

syssen = similar(Vector{DescriptorStateSpace},N)
for i = 1:N
    syssen[i] = [ SYSACSM[i] idm[:,sen]]
end

syssenf = fdimodset(syssen; c = Vector(1:mu), d = mu .+ Vector(1:md), f = (mu+md) .+ (1:mf));

## Part 2 - Setup of the synthesis specifications

# compute the achievable strong specifications for constant faults
S_strong = fdigenspec(syssenf[nom], atol = 1.e-6, atol3 = 1.e-6, FDtol = 1.e-5, FDGainTol = 0.001, FDfreq = 0, sdeg = -0.05);

@test size(S_strong,1) == 2^mf-1

# define SFDI, the signatures for isolation of single faults
SFDI = eye(mf) .> 0;

## Part 3 - Multiple filter synthesis using Procedure EFDI

Q = similar(Vector{DescriptorStateSpace},N)
for i = 1:N
    Qi = efdisyn( syssenf[i], SFDI;  atol = 1.e-6, sdeg = -5, smarg = -0.05, simple = false, FDfreq = 0, FDGainTol = 0.0001, rdim = 1)[1];
    Q[i] = vcat(Qi.sys...)
end
@test order.(Q) == zeros(Int,N)

##  Part 4 - Assesment of synthesis results
# form extended system [ Gu Gd Gf; I 0 0] with outputs extended with u
syse = similar(Vector{DescriptorStateSpace},N)
Nref = similar(Vector{Float64},N)
for i = 1:N
    syse[i] = [syssenf[i].sys;eye(mu,mu+md+mf)]
    rezi = Q[i]*syse[i] - [zeros(mf,mu+md) eye(2)]
    Nref[i]  = glinfnorm(rezi)[1]
end

##  Part 5 - Assesment of synthesis results for nominal synthesis
# form extended system [ Gu Gd Gf; I 0 0] with outputs extended with u
Nnom = similar(Vector{Float64},N)
for i = 1:N
    rezi = Q[nom]*syse[i] - [zeros(mf,mu+md) eye(2)]
    Nnom[i]  = glinfnorm(rezi)[1]
end

##  Part 6 - Assesment of synthesis results for constant interpolated gain
D0 = zeros(mf,p+mu)
D = zeros(T,mf,p+mu,N)
for i = 1:N
    D[:,:,i] = Q[i].D
end

for i = 1:mf
    for j = 1:p+mu
        D0[i,j] = fit(view(massi,:), view(D,i,j,:), 0)(0)
    end
end

N0 = similar(Vector{Float64},N)
for i = 1:N
    rezi = D0*syse[i] - [zeros(mf,mu+md) eye(2)]
    N0[i]  = glinfnorm(rezi)[1]
end

## Part 7 - Multiple model synthesis of a constant gain
# to be implemented using optimization-based tools

# optimal solution computed with MATLAB's SYSTUNE
Dopt0 = [
    -0.7998    1.0000   -0.0315   -0.0000    0.0069   -0.2488    4.1934    0.3858   -0.1468   -0.0773    0.0347
    0.0567         0   -0.0207    1.0000   -1.4180   -0.0067    4.5920    0.2673   -0.0255   -0.1018   -0.0228
]

# define parameterized constant filter gain: those set to false are fixed
mask = trues(mf,mu+p)
mask[1,2] = false
mask[1,4] = false
mask[2,2] = false
mask[2,4] = false
#Dinit = Dopt0  # this is the optimal solution
#Dinit = D0  # this  the normal initalization, unfortunately not converging
#Dres = Dinit[.!mask]
# use a 10% perturbed initialization
#Dinit = Dopt0 .* (1 .+ .1*randn(mf,mu+p))

function maxerror(x)
    err = 0.0
    D = similar(Matrix{Float64},mf,mu+p)
    D[mask] = x
    D[.!mask] = Dres
    for i in 1:N
        rezi = D*syse[i] - [zeros(mf,mu+md) eye(2)]
        err = max(err, glinfnorm(rezi)[1])
        #err += glinfnorm(rezi)[1]^2
        #err += norm(dcgain(rezi))^2
    end
    return err
end

# initialize with the optimal solution computed with MATLAB's SYSTUNE
Dinit = Dopt0  # this is the optimal solution
Dres = Dinit[.!mask]

println("Ninit = $(maxerror(Dinit[mask]))")
println("Please wait a few seconds - this optimization will finish soon")

res = optimize(maxerror, Dinit[mask], NelderMead(), Optim.Options(iterations = 20000))

println("Nfinal = $(maxerror(Optim.minimizer(res)))")

# comment out the next text for a more realistic try
# Dinit = Dopt0 .* (1 .+ .1*randn(mf,mu+p))
# Dres = Dopt0[.!mask]
# println("Ninit = $(maxerror(Dinit[mask]))")
# println("Please wait - this optimization take some time")

# res = optimize(maxerror, Dinit[mask], NelderMead(), Optim.Options(iterations = 20000))

# println("Nfinal = $(maxerror(Optim.minimizer(res)))")

Dopt = similar(Matrix{Float64},mf,mu+p)
Dopt[mask] = Optim.minimizer(res)
Dopt[.!mask] = Dres

Nopt = similar(Vector{Float64},N)
for i = 1:N
    rezi = Dopt*syse[i] - [zeros(mf,mu+md) eye(2)]
    Nopt[i]  = glinfnorm(rezi)[1]
end

norms = [Nref        Nnom  N0    Nopt]
println(" Table 8.4  Robustness analysis results for constant approximations")
println(" Nref         Nnom          N0        Nopt")
display(norms)
@test norm(Nref,Inf) < 1.e-7 &&  norm(Nnom,Inf) > 643 &&   norm(N0,Inf) > 504 &&   norm(Nopt,Inf) > 0.76

end  # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
#Here we will optimize the learning rate and sigmoid parameter for the XOR neural network

include("XOR_MD.jl")
include("Kernals.jl")
include("gaussian_process.jl")

#Set the random seed:
srand(1234) #Seed for stalzer srand(1234)

"""
We will compare the effect of randomly selecting a learning rate and sigmoid
hyperparameter vs the use of Bayesian Optimization for finding the optimial values wrt the MSE.
Suppose we have limited computing time of 100000 epochs and that we have N tries to
minimise the MSE. Let us say that the learning rate is between a and b

"""

#Initialise Layers and params ==========================================

Layer_1=uniform(0,1,2,2)

Layer_2=uniform(0,1,2,1)

epochs=1000   # was 1000 below, please change back
a=0.001  #Change to 0.001
b=1

c=0.001 #Was the same as above
d=1

N=20

#Curry the sigmoid functions:

function hyper_curry(h)
    return (x->sigmoid(x,h))
end

function hyper_curry_deriv(h)
    return (x->sigmoid_deriv(x,h))
end

# Random Learning Rates Examples ========================================

@time begin
    Random_Learning_Rates=uniform(a,b,N,1)
    Random_Hyperparameters=uniform(c,d,N,1)
    Random_Mat=cat(2,Random_Learning_Rates,Random_Hyperparameters)
    Random_MSE=zeros(N)

    #Random_Mat conjoins Random_Learning_Rates and Random_Hyperparameters
    # Random_Mat is a Nx2 matrix where Random_Mat[1,:] is the first entry
    #with LR_1 and hyperparemeter 1.

    for i=1:length(Random_Learning_Rates)
        node_function=hyper_curry(Random_Mat[i,2])
        node_deriv=hyper_curry_deriv(Random_Mat[i,2])
        learning_rate=Random_Mat[i,1]
        Random_MSE[i]=Train_Neural_Net_Loop(epochs,Layer_1,Layer_2,learning_rate,node_function,node_deriv)[3]
        println("Epoch Complete")
    end

    println("Random Learning Rates Training Completed")

end

# Random Test Plotting=========================================================

using PyPlot
# fig = figure("pyplot_subplot_mixed",figsize=(7,7))
ax=axes()

surf(reshape(Random_Learning_Rates,size(Random_MSE)),reshape(Random_Hyperparameters,size(Random_MSE)),Random_MSE,alpha=0.65,color="#40d5bb")
title("MSE for 20 Randomly Selected Parameter Values")
xlabel("Learning Rate")
ylabel("Hyper-Parameter")
zlabel("Mean Square Error")
grid("on")
show()
#12a3b4 aqua

#Bayesian Optimization Examples===================================================

#Here are the points we can pick from in the Optimization
@time begin

    #Initialise Layers and params ==========================================

    LR_Test=linspace(a,b,50)
    HP_Test=linspace(c,d,50)

    #Here is the carteisan product of these written as a vector
    Test=gen_points([LR_Test,HP_Test])[1]

    #We first have to pick a random point to begin bayesian optimization:

    #currently starts with the midpoint, possibly randomise this:
    Bayesian_Points=[Test[Int(round(length(Test)/2))]]

    #Bayesian_Points is an vector of arrays where in each array first entry is LR second entry is Hyper-Parameters:

    #Define hyperparemeter functions:
    node_function=hyper_curry(Bayesian_Points[1][2])
    node_deriv=hyper_curry_deriv(Bayesian_Points[1][2])

    #Define Learning Rate:
    learning_rate=Bayesian_Points[1][1]

    #Run first train before Bayesian Optimization:
    Bayesian_MSE=[Train_Neural_Net_Loop(epochs,Layer_1,Layer_2,learning_rate,node_function,node_deriv)[3]]

    #Begin Bayesian Optimization:

    for k=2:N
        D=[(Bayesian_Points[i],Bayesian_MSE[i]) for i=1:length(Bayesian_Points)]
        mu, sigma, D=gaussian_process_chol(std_exp_square_ker,D,1e-6,Test)
        # println("Gaussian Process Complete","\r")
        mu=reshape(mu,length(mu));
        sigma=reshape(sigma,length(sigma))

        new_point=findmin(mu-sigma)[2]

        #Here we will need to change the number 2 to k
        Bayesian_Points=cat(1,Bayesian_Points,[Test[new_point]])

        learning_rate=Bayesian_Points[k][1]

        node_function=hyper_curry(Bayesian_Points[k][2])

        node_deriv=hyper_curry_deriv(Bayesian_Points[k][2])

        value_to_be_appended=Train_Neural_Net_Loop(epochs,Layer_1,Layer_2,learning_rate,node_function,node_deriv)[3]

        if value_to_be_appended !=Bayesian_MSE[k-1]
            Bayesian_MSE=cat(1,Bayesian_MSE,[value_to_be_appended])
            println("Epoch Complete")
        else
            println("Found Optimum on the ", k-1, " iteration of ", N, " iterations")
            Bayesian_Points=Bayesian_Points[1:length(Bayesian_Points)-1]

            break
        end

    end

end

# Bayesian Plotting =========================================================

println(" The optimium is located at ",Bayesian_Points[end])

LR=[Bayesian_Points[i][1] for i=1:length(Bayesian_Points)]
HP=[Bayesian_Points[i][2] for i=1:length(Bayesian_Points)]

planex=[0,0,1,1]
planey=[0,1,0,1]
planem=[minimum(Random_MSE),minimum(Random_MSE),minimum(Random_MSE),minimum(Random_MSE)]

#Move this to the bottom
println("Bayesian_Learning_Rates Training Complete")
println("The minimum MSE by Bayesian Optimization was", minimum(Bayesian_MSE))
println("The minimum MSE by Random Selection was", minimum(Random_MSE))

using PyPlot
# fig = figure("pyplot_subplot_mixed",figsize=(7,7))
# ax=axes()
surf(LR,HP,Bayesian_MSE,alpha=0.65,color="#40d5bb")
surf(planex,planey,planem,alpha=0.3,color="#aa231f")

title("MSE BO - (6 point convergance)")
xlabel("Learning Rate")
ylabel("Hyper-Parameter")
zlabel("Mean Square Error")
grid("off")
show()


text:
module ModificationCaches

using LinearAlgebra
using SparseArrays

export
    VectorModificationCache,
    MatrixModificationCache,
    ProblemModificationCache,
    WarmStartCache,
    processupdates!

import OSQP

mutable struct VectorModificationCache{T}
    data::Vector{T}
    dirty::Bool
    VectorModificationCache(data::Vector{T}) where {T} = new{T}(copy(data), false)
end
Base.setindex!(cache::VectorModificationCache, x, i::Integer) = (cache.dirty = true; cache.data[i] = x)
Base.setindex!(cache::VectorModificationCache, x, ::Colon) = (cache.dirty = true; cache.data .= x)
Base.getindex(cache::VectorModificationCache, i::Integer) = cache.data[i]

function processupdates!(model, cache::VectorModificationCache, updatefun::Function)
    if cache.dirty
        updatefun(model, cache.data)
        cache.dirty = false
    end
end

struct MatrixModificationCache{T}
    cartesian_indices::Vector{CartesianIndex{2}}
    cartesian_indices_set::Set{CartesianIndex{2}} # to speed up checking whether indices are out of bounds in setindex!
    cartesian_indices_per_row::Dict{Int, Vector{CartesianIndex{2}}}
    modifications::Dict{CartesianIndex{2}, T}
    vals::Vector{T}
    inds::Vector{Int}

    function MatrixModificationCache(S::SparseMatrixCSC{T}) where T
        cartesian_indices = Vector{CartesianIndex{2}}(undef, nnz(S))
        cartesian_indices_per_row = Dict{Int, Vector{CartesianIndex{2}}}()
        sizehint!(cartesian_indices_per_row, nnz(S))
        @inbounds for col = 1 : S.n, k = S.colptr[col] : (S.colptr[col+1]-1) # from sparse findn
            row = S.rowval[k]
            I = CartesianIndex(row, col)
            cartesian_indices[k] = I
            push!(get!(() -> CartesianIndex{2}[], cartesian_indices_per_row, row), I)
        end
        modifications = Dict{CartesianIndex{2}, Int}()
        sizehint!(modifications, nnz(S))
        new{T}(cartesian_indices, Set(cartesian_indices), cartesian_indices_per_row, modifications, T[], Int[])
    end
end

function Base.setindex!(cache::MatrixModificationCache, x, row::Integer, col::Integer)
    I = CartesianIndex(row, col)
    @boundscheck I ∈ cache.cartesian_indices_set || throw(ArgumentError("Changing the sparsity pattern is not allowed."))
    cache.modifications[I] = x
end

function Base.setindex!(cache::MatrixModificationCache, x::Real, ::Colon)
    # used to zero out the entire matrix
    @boundscheck x == 0 || throw(ArgumentError("Changing the sparsity pattern is not allowed."))
    for I in cache.cartesian_indices
        cache.modifications[I] = 0
    end
end

function Base.setindex!(cache::MatrixModificationCache, x::Real, row::Integer, ::Colon)
    # used to zero out a row
    @boundscheck x == 0 || throw(ArgumentError("Changing the sparsity pattern is not allowed."))
    for I in cache.cartesian_indices_per_row[row]
        cache.modifications[I] = 0
    end
end

function Base.getindex(cache::MatrixModificationCache, row::Integer, col::Integer)
    cache.modifications[CartesianIndex(row, col)]
end

function processupdates!(model, cache::MatrixModificationCache, updatefun::Function)
    dirty = !isempty(cache.modifications)
    if dirty
        nmods = length(cache.modifications)
        resize!(cache.vals, nmods)
        resize!(cache.inds, nmods)
        count = 1
        for i in eachindex(cache.cartesian_indices)
            I = cache.cartesian_indices[i]
            if haskey(cache.modifications, I)
                cache.vals[count] = cache.modifications[I]
                cache.inds[count] = i
                count += 1
            end
        end
        updatefun(model, cache.vals, cache.inds)
        empty!(cache.modifications)
    end
end

# More OSQP-specific from here on:
struct ProblemModificationCache{T}
    P::MatrixModificationCache{T}
    q::VectorModificationCache{T}
    A::MatrixModificationCache{T}
    l::VectorModificationCache{T}
    u::VectorModificationCache{T}

    ProblemModificationCache{T}() where {T} = new{T}()
    function ProblemModificationCache(P::SparseMatrixCSC, q::Vector{T}, A::SparseMatrixCSC, l::Vector{T}, u::Vector{T}) where T
        MC = MatrixModificationCache
        VC = VectorModificationCache
        new{T}(MC(triu(P)), VC(q), MC(A), VC(l), VC(u))
    end
end

function processupdates!(model::OSQP.Model, cache::ProblemModificationCache)
    if cache.l.dirty && cache.u.dirty
        # Special case because setting just l or u may cause an 'upper bound must be greater than or equal to lower bound' error
        OSQP.update_bounds!(model, cache.l.data, cache.u.data)
        cache.l.dirty = false
        cache.u.dirty = false
    end

    processupdates!(model, cache.P, OSQP.update_P!)
    processupdates!(model, cache.q, OSQP.update_q!)
    processupdates!(model, cache.A, OSQP.update_A!)
    processupdates!(model, cache.l, OSQP.update_l!)
    processupdates!(model, cache.u, OSQP.update_u!)
end

struct WarmStartCache{T}
    x::VectorModificationCache{T}
    y::VectorModificationCache{T}

    WarmStartCache{T}() where {T} = new{T}()
    function WarmStartCache{T}(n::Integer, m::Integer) where T
        new{T}(VectorModificationCache(zeros(T, n)), VectorModificationCache(zeros(T, m)))
    end
end

function processupdates!(model::OSQP.Model, cache::WarmStartCache)
    if cache.x.dirty && cache.y.dirty
        # Special case because setting warm start for x only zeroes the stored warm start for y and vice versa.
        OSQP.warm_start_x_y!(model, cache.x.data, cache.y.data)
        cache.x.dirty = false
        cache.y.dirty = false
    end
    processupdates!(model, cache.x, OSQP.warm_start_x!)
    processupdates!(model, cache.y, OSQP.warm_start_y!)
end

end


text:
using Elly
using Test
using Distributed

function test_yarn(host="localhost", rmport=8032, schedport=8030)
    limitedtestenv = (get(ENV, "CI", "false") == "true")

    yarnclnt = YarnClient(host, rmport)
    nnodes = nodecount(yarnclnt)
    @test nnodes > 0
    println("number of yarn nodes: $nnodes")

    nlist = nodes(yarnclnt)
    show(stdout, nlist)

    yarncm = YarnManager(yarnhost=host, rmport=rmport, schedport=schedport, launch_timeout=60);

    env = Dict{String,String}()
    for envname in ("USER", "LD_LIBRARY_PATH", "USERNAME", "HOME", "PATH", "LOGNAME", "JULIA_LOAD_PATH", "LIBDIR")
        if envname in keys(ENV)
            env[envname] = ENV[envname]
        end
    end
    if !("JULIA_LOAD_PATH" in keys(env))
        home = ENV["HOME"]
        env["JULIA_LOAD_PATH"] = join(Base.LOAD_PATH, ':') * ":$(home)/.julia/dev:$(home)/.julia/packages"
    end
    if !("JULIA_DEPOT_PATH" in keys(env))
        env["JULIA_DEPOT_PATH"] = join(Base.DEPOT_PATH, ':')
    end
    println("starting workers with environment:")
    for (n,v) in env
        println("    - $n => $v")
    end

    addprocs(yarncm; np=1, env=env);
    @test nprocs() == 2

    if !limitedtestenv
        addprocs(yarncm; np=2, env=env);
        @test nprocs() == 4
    end

    @everywhere println("hi")
    rmprocs(workers())
    @test nprocs() == 1

    Elly.disconnect(yarncm)
end


text:
######################
# Duffing Oscillator #
######################

Adaptive=1
LONG=Int(1e4)
δt=1e-3

#Kinematic Parameter
Λ=3

#Nominal Trajectory
ω=5
Ampl=0.4

#Controller Parameters
K=1e5
B=-1
A=1e-5

# Exact model parameters
αₑ=1
δₑ=0.2
βₑ=1

# Approx model parameters
αₐ=0.8
δₐ=0.1
βₐ=0.9

###############
# Exact Model #
###############

function exact(q,q_p,u)
  out=αₑ*q+δₑ*q_p+βₑ*q^3+u
  return out #q_pp
end

################
# Approx Model #
################

function approx(q,q_p,q_pp)
  out=q_pp-αₐ*q-δₐ*q_p-βₐ*q^3
  return out #u
end

########################
# Deformation Function #
########################
function sigmoid(x)
	s=x/(1+abs(x))
  return s
end
function G(past_input,past_response,desired)
  out=(K+past_input)*(1 + B * tanh(A * (past_response - desired))) - K
  return out
end #qDef_pp

######################################
# Arrays for Plotting and simulation #
######################################

qN_mem=zeros(LONG)
qN_p_mem=zeros(LONG)
qN_pp_mem=zeros(LONG)

qDes_pp_mem=zeros(LONG)

qDef_pp_mem=zeros(LONG)
u_mem=zeros(LONG)

q_mem=zeros(LONG)
q_p_mem=zeros(LONG)
q_pp_mem=zeros(LONG)

time_mem=zeros(LONG)

error_int=0
past_input=0
past_response=0

#simulation
for t=1:LONG-1
  time_mem[t]=t*δt
  #Nominal trajectory
  qN_mem[t]=Ampl*sin(ω*δt*t)#sin(ω*time_mem[t])
  qN_p_mem[t]=ω*Ampl*cos(ω*δt*t)
  qN_pp_mem[t]=-ω^2*Ampl*sin(ω*δt*t)
  #Kinematic Block
  qDes_pp_mem[t]=Λ^3*error_int+3*Λ^2*(qN_mem[t]-q_mem[t])+3*Λ*(qN_p_mem[t]-q_p_mem[t])+qN_pp_mem[t]
  # Deformation
  if Adaptive==1 && t>10
    qDef_pp_mem[t]=G(past_input,past_response,qDes_pp_mem[t])
  else
    qDef_pp_mem[t]=qDes_pp_mem[t]
  end
  #Approx model
  u=approx(q_mem[t],q_p_mem[t],qDef_pp_mem[t])
  u_mem[t]=u
  past_input=qDef_pp_mem[t]
  #Exact model
  q_pp_mem[t]=exact(q_mem[t],q_p_mem[t],u)
  past_response=q_pp_mem[t]
  #Integrals
  q_p_mem[t+1]=q_p_mem[t]+δt*q_pp_mem[t]
  q_mem[t+1]=q_mem[t]+δt*q_p_mem[t]
  error_int=error_int+δt*(qN_mem[t]-q_mem[t])
end #for (simulation)

using PyPlot

figure(1)
title("Tracking error")
plot(time_mem[3:LONG-1],qN_mem[3:LONG-1]-q_mem[3:LONG-1])
grid("on")
xlabel("Time [s]")
ylabel("error [m]")
#xlabel(L"\frac{x}{y}")   ##(L"") latex-es felirathoz

figure(2)
title("Nominal and Simulated")
plot(time_mem[3:LONG-1],qN_mem[3:LONG-1],color="red")
plot(time_mem[3:LONG-1],q_mem[3:LONG-1],color="green")
grid("on")
xlabel("Time [s]")
ylabel("q [m]")

max_stars_repo_name: krisztiankosi/Control

max_issues_repo_name: krisztiankosi/Control

max_forks_repo_name: krisztiankosi/Control


text:
module TestManager
  using FactCheck
  import Helper.getArg
  import Helper
  import Manager
  import Config
  import Creature
  import ManagerTypes
  import Mutator
  import Code
  import World
  include("Helper.jl")

  facts("Checking if mutations mechanism works") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}(:orgRainMutationPeriod=>2, :orgRainMutationPercent=>0.1, :plugIncluded=>["Status"]))
    local mutations = d.man.plugins["Status"].mps

    consume(d.task)
    @fact d.man.plugins["Status"].mps - mutations --> 0
    consume(d.task)
    @fact d.man.plugins["Status"].mps - mutations --> 1

    Manager.destroy(d.man)
  end
  facts("Checking if mutations mechanism works with specified amount") do
    local d = create([Helper.Point(1,1)], Dict(:orgRainMutationPeriod=>2, :orgRainMutationPercent=>1.0, :plugIncluded=>["Status"]))
    local mutations = d.man.plugins["Status"].mps

    code(d, :plus, d.orgs[1])
    code(d, :plus, d.orgs[1])

    consume(d.task)
    @fact d.man.plugins["Status"].mps - mutations --> 0
    consume(d.task)
    @fact d.man.plugins["Status"].mps - mutations --> 2

    Manager.destroy(d.man)
  end
  facts("Checking if correct amount of organisms are created on start") do
    local d = create(Helper.Point[], Dict{Symbol, Any}(:orgStartAmount=>5))

    @fact length(d.man.positions) === length(d.man.organisms) === 0 --> true
    consume(d.task)
    @fact length(d.man.positions) === length(d.man.organisms) === 5 --> true

    Manager.destroy(d.man)
  end
  # orgStartAmount, orgStartEnergy
  facts("Checking if organisms, which are created on start contain correct amount of energy") do
    local d = create(Helper.Point[], Dict{Symbol, Any}(:orgStartAmount=>3,:orgStartEnergy=>7))

    @fact length(d.man.positions) === length(d.man.organisms) === 0 --> true
    consume(d.task)
    @fact length(d.man.positions) === length(d.man.organisms) === 3 --> true
    for org in d.man.organisms @fact org[2].energy --> 7 end

    Manager.destroy(d.man)
  end
  # orgEnergySpendPeriod
  facts("Checking amount of energy grabbing from organisms per period") do
    local d = create([Helper.Point(1,1), Helper.Point(2,2), Helper.Point(3,3)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>2, :orgGarbagePeriod=>1, :plugIncluded=>["Status"]))

    code(d, :plus, d.orgs[1])
    code(d, :minus, d.orgs[2])
    code(d, :minus, d.orgs[2])
    code(d, :not, d.orgs[3])
    code(d, :not, d.orgs[3])
    code(d, :not, d.orgs[3])
    @fact d.orgs[1].energy --> 100
    @fact d.orgs[2].energy --> 100
    @fact d.orgs[3].energy --> 100
    # orgEnergySpendPeriod === 2
    # so we need run 4 iterations to decrease energy on 6 points
    consume(d.task)
    consume(d.task)
    consume(d.task)
    consume(d.task)
    @fact d.orgs[1].energy --> 98
    @fact d.orgs[2].energy --> 96
    @fact d.orgs[3].energy --> 94

    Manager.destroy(d.man)
  end
  # orgClonePeriod
  facts("Checking organisms clonning ability") do
    local d = create([Helper.Point(5,5)], Dict{Symbol, Any}(:orgClonePeriod=>3))
    local orgAmount = length(d.man.organisms)

    consume(d.task)
    @fact length(d.man.organisms) - orgAmount --> 0
    consume(d.task)
    @fact length(d.man.organisms) - orgAmount --> 0
    consume(d.task)
    @fact length(d.man.organisms) - orgAmount --> 1
    @fact length(d.man.positions) - orgAmount --> 1

    Manager.destroy(d.man)
  end
  # orgStartEnergy, orgEnergySpendPeriod
  facts("Checking if organisms dying if no energy works") do
    local d = create([Helper.Point(5,5)], Dict{Symbol, Any}(:orgStartEnergy=>2, :orgEnergySpendPeriod=>2))
    local org = d.orgs[1]

    code(d, :plus, org)
    consume(d.task)
    @fact org.energy --> 2
    consume(d.task)
    @fact org.energy --> 1
    consume(d.task)
    @fact org.energy --> 1
    consume(d.task)
    @fact org.energy --> 0
    @fact length(d.man.organisms) --> 0

    Manager.destroy(d.man)
  end
  # orgStartEnergy, orgAlivePeriod, worldMinOrgs
  facts("Checking if organisms dying, because of age works") do
    local d = create([Helper.Point(5,5)], Dict{Symbol, Any}(:orgStartEnergy=>10, :orgAlivePeriod=>2, :worldMinOrgs=>0))
    local org = d.orgs[1]

    code(d, :plus, org)
    consume(d.task)
    @fact length(d.man.organisms) --> 1
    consume(d.task)
    @fact length(d.man.organisms) --> 1
    consume(d.task)
    @fact length(d.man.organisms) --> 0
    @fact org.age                 --> 2

    Manager.destroy(d.man)
  end
  #orgClonePeriod, orgCloneMutation
  facts("Checking organisms mutations on clone") do
    local d         = create([Helper.Point(5,5)], Dict{Symbol, Any}(:orgClonePeriod=>2, :orgCloneMutation=>1.0, :plugIncluded=>["Status"]))
    local org       = d.orgs[1]
    local mutations = d.man.plugins["Status"].mps

    code(d, :plus, org)
    code(d, :plus, org)

    @fact org.mutationsFromStart        --> 1
    @fact length(d.man.organisms)       --> 1

    consume(d.task)
    @fact length(d.man.organisms)       --> 1
    @fact org.mutationsFromStart        --> 1

    consume(d.task)
    consume(d.task)
    @fact length(d.man.organisms)       --> 2
    @fact org.mutationsFromStart        --> 1
    @fact d.man.organisms[org.id+1].mutationsFromStart  --> 3

    Manager.destroy(d.man)
  end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
abstract type ConvergenceMetric end

description(m::ConvergenceMetric) =
    error("`description` is not implemented for $(typeof(m)).")

"""
    converged(metric)

判断是否已收敛
"""
converged(m::ConvergenceMetric) = diff(m) <= tolerance(m)

diff(m::ConvergenceMetric) = m.Δ #error("`diff` is not implemented for $(typeof(m)).")

tolerance(m::ConvergenceMetric) = m.tol #error("`tolerance` is not implemented for $(typeof(m)).")

"""
    assess!(metric, state)

根据状态计算并判断是否已收敛
"""
assess!(m::ConvergenceMetric, s::AbstractOptimizerState) =
    error("`assess!` is not implemented for $(typeof(m)).")

const ConvergenceMetrics = Vector{ConvergenceMetric}

copy(cm::ConvergenceMetrics) = map(copy, cm)

# Utilities

abschange(curr, prev) = Float64(abs(curr - prev))
relchange(curr, prev) = Float64(abs(curr - prev) / abs(curr))

maxdiff(x::AbstractArray, y::AbstractArray) = mapreduce((a, b) -> abs(a - b), max, x, y)
abschange(curr::T, prev) where {T<:AbstractArray} = maxdiff(curr, prev)
relchange(curr::T, prev) where {T<:AbstractArray} = maxdiff(curr, prev) / maximum(abs, curr)

# Single value

"""
Absolute difference convergence for single objective optimization.

This convergence metric allows to estimate an absolute difference between consecutive
states of the optimization algorithm, and triggers convergence when,

- `|f(x) - f(x')| < ε`

where `ε` is a tolerance value, `x` and `x'` are previous and current minimizers
found by the optimization algorithm.
"""
mutable struct AbsDiff{T} <: ConvergenceMetric
    tol::T
    Δ::Float64
    value::T
end
AbsDiff(tol::T) where {T<:AbstractFloat} = AbsDiff(tol, Inf, zero(T))
AbsDiff() = AbsDiff(1e-12)
description(m::AbsDiff) = "|f(x) - f(x')|"
function assess!(m::AbsDiff, state::AbstractOptimizerState)
    val = value(state)
    m.Δ = abschange(val, m.value)
    m.value = val
    converged(m)
end
copy(m::AbsDiff) = AbsDiff(copy(m.tol), copy(m.Δ), copy(m.value))

"""
Relative difference convergence metric for single objective optimization.

This convergence metric allows to estimate a relative difference between consecutive
states of the optimization algorithm, and triggers convergence when,

- `|f(x) - f(x')|/|f(x')| < ε`

where `ε` is a tolerance value, `x` and `x'` are previous and current minimizers
found by the optimization algorithm.
"""
mutable struct RelDiff{T} <: ConvergenceMetric
    tol::T
    Δ::Float64
    value::T
end
RelDiff(tol::T) where {T<:AbstractFloat} = RelDiff(tol, Inf, zero(T))
RelDiff() = RelDiff(1e-12)
description(m::RelDiff) = "|f(x) - f(x')|/|f(x')|"
function assess!(m::RelDiff, state::AbstractOptimizerState)
    val = value(state)
    m.Δ = relchange(val, m.value)
    m.value = val
    converged(m)
end
copy(m::RelDiff) = RelDiff(copy(m.tol), copy(m.Δ), copy(m.value))

"""
    gd(A,R)

Calculate a generational distance between set `A` and the reference set `R`.
This metric measures the convergence, i.e. closeness of the non-dominated solutions
to the Pareto front, of a population.

*Note:* Parameters are column-major matrices.
"""
function gd(A::AbstractMatrix, R::AbstractMatrix)
    na, da = size(A)
    nr, dr = size(R)
    (na == 0 || nr == 0) && return Inf
    sum = 0
    for a in eachrow(A)
        sum += minimum(norm(a - r) for r in eachrow(R))
    end
    sum / na
end
gd(A::T, R::T) where {T<:Population} = gd(objectives(A), objectives(R))

"""
    igd(S,R)

Calculate an inverted generational distance, [`gd`](@ref), between set `S` and the reference set `R`.
Parameters are column-major matrices.
"""
igd(S, R) = gd(R, S)

"""
(Inverse) generational distance convergence metric for multi-objective optimization.

This convergence metric allows to estimate a (inverse) generational distance between
consecutive states of the multi-objective optimization algorithm, and triggers convergence when,

- ``GD(S,R) = \\frac{1}{|S|}\\left(\\sum_{s \\in S} \\min_{r \\in R}||f(s) - f(r)||^2\\right)^{1/2} < \\varepsilon``

where `ε` is a tolerance value, `S` and `R` are discrete Pareto front approximations
for previous and current state of the optimization algorithm.

The inverse generational distance calculated as `IGD(S,R) = GR(R,S)`.
To use, `IGD` metric pass `true` value to the metric constructor argument, e.g.

- `GD()` create a generational distance metric
- `GD(true)` create an inverse generational metric
"""
mutable struct GD{T} <: ConvergenceMetric
    tol::T
    Δ::Float64
    value::Matrix{T}
    inverted::Bool
end
GD(tol::T, inv = false) where {T<:AbstractFloat} = GD(tol, Inf, zeros(T, 0, 0), inv)
GD(inv = false) = GD(1e-5, inv)
function description(m::GD)
    prefix = m.inverted ? "I" : ""
    "|$(prefix)GD(P) - $(prefix)GD(P')|"
end
function assess!(m::GD, state::AbstractOptimizerState)
    val = value(state)
    m.Δ = m.inverted ? igd(val, m.value) : gd(val, m.value)
    m.value = val
    converged(m)
end
copy(m::GD) = GD(copy(m.tol), copy(m.Δ), copy(m.value), copy(m.inverted))

"""
    spread(S,R)

Returns a diversity metric of a population of set `S` to the reference set `R`.
"""
function spread(S::AbstractMatrix, R::AbstractMatrix)
    n = size(S, 2)
    m = size(R, 2)
    Δₖ = [
        minimum(
            norm(view(S, :, i) - view(R, :, j)) for
            i in 1:n if view(S, :, i) != view(R, :, j)
        ) for j = 1:m
    ]
    Δ = mean(Δₖ)
    sum(abs.(Δₖ .- Δ)) / m * Δ
end

function spread(S::AbstractMatrix)
    n = size(S, 2)
    n == 1 && return NaN
    Δₖ = [
        minimum(
            norm(view(S, :, i) - view(S, :, j)) for
            j in 1:n if view(S, :, i) != view(S, :, j)
        ) for i = 1:n
    ]
    Δ = mean(Δₖ)
    sum(abs.(Δₖ .- Δ)) / n * Δ
end

##########################
# Convergence Assessment #
##########################

function assess_convergence!(state)
    for metric in state.metrics
        converged = assess!(metric, state)
        if converged
            return true
        end
    end
    return false
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -6.482349417201714 <= q <= 6.592803073996197)

                     add_NL_constraint(m, :(sigmoid(-0.7946398993851154 + 0.38713964342476537*sigmoid(0.024809189775459206 + 0.8902991643012754*$(x[1]) + 0.015536230415071905*$(x[2]) + -0.8172113269588737*$(x[3]) + 0.9746036973952403*$(x[4])) + -0.9438792568465639*sigmoid(-0.017899887293421735 + -0.33137455646820024*$(x[1]) + 0.5594894653083848*$(x[2]) + -0.34768764708208755*$(x[3]) + 0.43788143230656695*$(x[4]))) + sigmoid(0.8346867556655706 + -0.9900146410670576*sigmoid(0.024809189775459206 + 0.8902991643012754*$(x[1]) + 0.015536230415071905*$(x[2]) + -0.8172113269588737*$(x[3]) + 0.9746036973952403*$(x[4])) + -0.7397521454164151*sigmoid(-0.017899887293421735 + -0.33137455646820024*$(x[1]) + 0.5594894653083848*$(x[2]) + -0.34768764708208755*$(x[3]) + 0.43788143230656695*$(x[4]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module DFUtil

using DataFrames
using Dates
using Pipe
using CSV

export sum_columns, group_data_into_periods, match_row, to_json, to_json_var, to_js_value
export de_miss_rows, include_or_exclude, firstrow, lastrow, nthrow
export pQuarter, pWeek, pYear, pMonth
export to_csv_text, from_csv_text
export concat!, leftjoiner, tryRename, tryRename!, dropByName!, colsToGrid

function force_vector(v)
	if isa(v, String) || isa(v, Symbol)
		return [v]
	end
	v = collect(Iterators.flatten([v]))
	return isa(v, Vector) ? v : [v]
end

"""
	include_or_exclude(df, includes, excludes)

Given a DataFrame either restrict it to the includes list, or remove the exludes list

Sounds a bit daft but it is to provide KW options in functions

If both are given, the `includes` take priority

# Arguments
- `df` the DataFrame
- `includes` the String, Symbol or list of those to include
- `excludes` the String, Symbol or list of those to exclude

# Examples

	df = DataFrame([[1,2,3], [10,20,30]], ["a", "b"])

	julia> include_or_exclude(df, includes=["a"])
	3×1 DataFrame
	Row │ a
		│ Int64
	─────┼───────
	1 │     1
	2 │     2
	3 │     3

	julia> include_or_exclude(df, excludes="a")
	3×1 DataFrame
	Row │ b
		│ Int64
	─────┼───────
	1 │    10
	2 │    20
	3 │    30

	julia> include_or_exclude(df, includes=["a"], excludes=["a", "b"])
	3×1 DataFrame
	Row │ a
		│ Int64
	─────┼───────
	1 │     1
	2 │     2
	3 │     3

"""
function include_or_exclude(df; includes=nothing, excludes=nothing)
	if excludes === nothing && includes === nothing
		return df
	end

	if excludes === nothing
		# force vector to makes sure it returns a DataFrame
		return df[!, force_vector(includes)]
	end

	if includes === nothing
		return df[!, Not(excludes)]
	end
	# we have includes *and* excludes.
	# Includes take priority

	# makes sure they are both lists, and the same datatype
	includes = map(i->string(i), force_vector(includes))
	excludes = map(e->string(e), force_vector(excludes))

	return df[!, Not(filter(e->! (e in includes), excludes))]
end
"""
	sum_columns(df, group_by::Vector{String}=Vector{String}())

Sum the columns of the dataframe, optionally grouping by group_by

All the columns not in the grouping *must* have a `+` method

# Arguments

- `df` DataFrame to sum
- `group_by` Vector of String names of the columns to group by
- `replace_with` by default DataFrames takes the column names and appends _function. Instead use this string.
"""
sum_columns(df; group_by=Vector{String}(), replace_with="s") = @pipe map(String, group_by) |>
		groupby(df, _) |>
		combine(_,  [ c => c->sum(skipmissing(c)) for c in filter(n->!(n in group_by), names(df)) ]) |>
		rename(_, map(c-> c=>replace(c, "_function"=>replace_with), names(_)))

"""
	group_data_into_periods(df::DataFrame, date_column::Union{Symbol, AbstractString}, period::Union{Symbol, AbstractString}; andgrpby = Vector{String}())

Using the given data_column, group the dataframe into periods.

All columns except the grouping ones *must* have methods for `+``

# Arguments
- `df` DataFrame to group
- `date_column` the column to use as the grouping
- `period` Period to use, options are: :Quarter, :Month, :Year, :Week (or as strings)
- `andgrpby` optional keyword to use additional groupings

# Examples
	group_data_into_periods(df, :SaleDate, :Week)
	group_data_into_periods(df, "SaleDate", :Year, "BranchId")
	group_data_into_periods(df, :SaleDate, "Quarter", [:Area, :BranchId)]
"""
function group_data_into_periods(df, date_column, period; andgrpby=Vector{String}())
	if isa(period, String)
		period = Symbol(period)
	end
	if isa(andgrpby, String) || isa(andgrpby, Symbol)
		andgrpby = [andgrpby]
	end

	period_fns = Dict(:Quarter => pQuarter, :Month => pMonth, :Year => pYear, :Week => pWeek)

	@pipe transform(df, date_column => ByRow(dt -> period_fns[period](dt)) => "$(date_column)_function") |>
		select(_, Not(date_column)) |>
		sum_columns(_, group_by=vcat(map(String, andgrpby), ["$(date_column)_function"]))
end

# period functions
"""
	pQuarter(dt)

turn a Date / Datetime into its eqivalent YearQn representation e.g. 2001Q1
"""
pQuarter(dt)   = ismissing(dt) ? missing : string(Dates.year(dt)) * "Q" * string(Dates.quarterofyear(dt))
"""
pMonth(dt)

turn a Date / Datetime into its eqivalent YearMonth representation e.g. 2001-01
"""
pMonth(dt) = ismissing(dt) ? missing : string(Dates.year(dt)) * "-" * ("00" * string(Dates.month(dt)))[end-1:end]
"""
	pYear(dt)

turn a Date / Datetime into its eqivalent Year representation e.g. 2001
"""
pYear(dt)  = ismissing(dt) ? missing : Dates.year(dt)
"""
	pWeek(dt)

Turn a Date / Datetime into its eqivalent YearWeek representation e.g. 2001-01
Preserves 2000-01-01 becoming 1999-52
"""
function pWeek(dt) # this is the tricky one because 2000-01-01 is 1999-52
	if ismissing(dt)
		return  missing
	end
	m = Dates.month(dt)
	w = Dates.week(dt)
	y = Dates.year(dt)
	if m == 1 && w > 10
		y = y - 1
	end
	"$y-" * ("00" * string(w))[end-1:end] # right pad with 0
end
"""
	match_row(df, col, val)

Just a shortcut for `filter(row -> row[col] == val, df)`

"""
match_row(df, col, val) = filter(row -> row[col] == val, df)

"""
	to_json_var(io::IO, df::DataFrame, keys::Union{AbstractString, Vector{AbstractString}}, var="object")

Output the dataframe to json using to_json but wrap it inside a Javascript JSON.parse() for direct use in Javascript

# Arguments
- `io` IO handle to write to
- `df` DataFrame to write
- `keys` Keys to lift outside the object - see [`to_json`](@ref) for explanation of that
- `var` JavaScript object name it will be assign to with `var \$var = ....`

# Examples

	df = DataFrame([[1,2],[3,4],[4,5]],["a", "b", "c"]), "a")

	to_json_var(stdout, df, "a", "object")

	"var object = JSON.parse('{ \"1\" : {\"a\" : \"1\", \"b\" : \"3\", \"c\" : \"4\" } , \"2\" : {\"a\" : \"2\", \"b\" : \"4\", \"c\" : \"5\" }  }');"
"""
function to_json_var(io::IO, df::DataFrame, keys, var="object")
	print(io, "var $var = JSON.parse('")
	to_json(io, df, keys)
	println(io, "');");
end

"""
	to_js_value(io::IO, df::DataFrame, keys::Union{AbstractString, Vector{AbstractString}}, jsdeclaration="var object")

Output the dataframe to json using to_json but wrap it inside a Javascript JSON.parse() for direct use in Javascript

# Arguments
- `io` IO handle to write to
- `df` DataFrame to write
- `keys` Keys to lift outside the object - see [`to_json`](@ref) for explanation of that
- `var` JavaScript object name it will be assign to with `var \$var = ....`

# Examples

	df = DataFrame([[1,2],[3,4],[4,5]],["a", "b", "c"]), "a")

	to_js_value(stdout, df, "a", "export const vals")

	"export const vals = JSON.parse('{ \"1\" : {\"a\" : \"1\", \"b\" : \"3\", \"c\" : \"4\" } , \"2\" : {\"a\" : \"2\", \"b\" : \"4\", \"c\" : \"5\" }  }');"
"""
function to_js_value(io::IO, df::DataFrame, keys, jsdeclaration="var object")
	print(io, "$jsdeclaration = JSON.parse('")
	to_json(io, df, keys)
	println(io, "');");
end

"""
	to_json(io::IO, data::DataFrame, keys::Union{AbstractString, Vector{AbstractString}})

Output the dataframe to Json, grouping by the given keys

	{ "row[keys[1]]" : { "row[keys[2]]" : { "names(row)[1]" : "row[names(row[1])]", "names(row)[2]" : "row[names(row[2])]", ..., "names(row)[end]" : "row[names(row[end])]"}}}

The motivation for this is to enable referencing the objects via their nested keys in JavaScript

	julia> df = DataFrame([[10, 10, 11], [1, 2, 3], ["Fred", "Luke", "Alice"]], ["Area", "EmpId", "Name"])
	3×3 DataFrame
	 Row │ Area   EmpId  Name
		 │ Int64  Int64  String
	─────┼──────────────────────
	   1 │    10      1  Fred
	   2 │    10      2  Luke
	   3 │    11      1  Alice

	   julia> to_json(stdout, df, ["Area", "EmpId"])

	   { "10" : {"1" : {"Area" : "10", "EmpId" : "1", "Name" : "Fred" } , "2" : {"Area" : "10", "EmpId" : "2", "Name" : "Luke" }  }, "11" : {"3" : {"Area" : "11", "EmpId" : "3", "Name" : "Alice" }  } }

in Javascript

	let fred = object[10][1]
	let luke = object[10][2]
	let alice = object[11][1]

# Arguments
- `io` IO handle to write to
- `df` DataFrame to write
- `keys` Keys to lift outside the object

# Examples

	df = DataFrame([[1,2],[3,4],[4,5]],["a", "b", "c"]), "a")

	to_json(stdout, df, "a")

	"{ \"1\" : {\"a\" : \"1\", \"b\" : \"3\", \"c\" : \"4\" } , \"2\" : {\"a\" : \"2\", \"b\" : \"4\", \"c\" : \"5\" }  }"

	to_json(stdout, df, ["a", "b"])

	"{ \"1\" : {\"3\" : {\"a\" : \"1\", \"b\" : \"3\", \"c\" : \"4\" }  }, \"2\" : {\"4\" : {\"a\" : \"2\", \"b\" : \"4\", \"c\" : \"5\" }  } }")

"""
function to_json(io, df, keys)
	print(io, "{ ")
	if isa(keys, Array)
		_to_json(io, groupby(df, keys[1]), keys, 1)
	else
		_to_json(io, groupby(df, keys), [keys], 1)
	end
	print(io, " }")
end

# turn the first row of a DF into a DataFrameRow, a helper for _to_json
row(grp) = eachrow(grp)[1]

# helper for to_json
function _to_json(io, data, keys, depth)
	prefix = ""
	if depth == length(keys)
		print_data_row(io, row(data[1]), keys[depth], "")
		for grp in data[2:end]
			print_data_row(io, row(grp), keys[depth], ", ")
		end
	else
		for grp in data
			k = kesc(grp[!, keys[depth]][1])
			print(io, prefix, "\"$k\" : {")
			_to_json(io, groupby(grp, keys[depth+1]), keys, depth+1)
			print(io, " }");
			prefix = ", "
		end
	end
end

# key / value escapers
jesc(v) = replace(string(v), "'"=>"\\'")
kesc(k) = replace(jesc(k), " "=>"")

# print a key value pair
print_term(io, row, key, prefix="") = print(io, prefix, "\"$(kesc(key))\" : ", "\"$(jesc(row[key]))\"");

# print all key value pairs, with appropriate comma-ing
function print_terms(io, row)
	print_term(io, row, names(row)[1], "")
	foreach(term -> print_term(io, row, term, ", "), names(row)[2:end])
end

# print a compete row, with single pkey
function print_data_row(io, row, pkey, prefix)
	print(io, prefix, "\"$(kesc(row[pkey]))\" : {")
	print_terms(io, row)
	print(io, " } ");
end
"""
	nthrow(df)

Return the nth row of a dataframe as a DataFrameRow, if n is out of bounds, return nothing
"""
function nthrow(df, n)
	if size(df)[1] < n || n < 1
		return
	end
	eachrow(df)[n]
end

"""
	firstrow(df)

Return the first row as a DataFrameRow
"""
firstrow(df) = nthrow(df, 1)

"""
	lastrow(df)

	Return the last row of df as a DataFrameRow
"""
lastrow(df) = nthrow(df, size(df)[1])

"""
	exclude_rows(df, fn)

Exclude any row which has any column return true for the function

# Examples

	exclude_rows(df, ismissing)
	exclude_rows(df, isnan)

"""

function exclude_rows(df, fn)

	function pass(row)
		for n in names(df)
			if fn(row[n])
				return false
			end
		end
		return true
	end

	filter(pass, df)
end

"""
	de_miss_rows(df)

For a given DataFrame, remove any rows in which any column has a missing value

# Examples

	df = DataFrame([[1,missing,3], [10,20,30]], ["a", "b"])
	3×2 DataFrame
	 Row │ a        b
		 │ Int64?   Int64?
	─────┼─────────────────
	   1 │       1      10
	   2 │ missing      20
	   3 │       3      30

	julia> dm = de_miss_rows(df)
	2×2 DataFrame
	 Row │ a       b
		 │ Int64?  Int64?
	─────┼────────────────
	   1 │      1      10
	   2 │      3      30

"""
de_miss_rows(df) = exclude_rows(df, ismissing)

"""
	to_csv_text(df)

turn the dataframe into the string representation of the CSV

# Examples

	df = DataFrame(
	[[53.685335, 53.785335, 53.885335],
	[-0.416387,-0.426387,-0.436387],
	[53.479544, 53.479544,53.479544],
	[-2.954101,-2.954101,-2.954101]], ["lat1", "lon1", "lat2", "lon2"])

	to_csv_text(df)
	"lat1,lon1,lat2,lon2\\n53.685335,-0.416387,53.479544,-2.954101\\n53.785335,-0.426387,53.479544,-2.954101\\n53.885335,-0.436387,53.479544,-2.954101\\n"
"""
function to_csv_text(df)
	out = IOBuffer()
	CSV.write(out, df)
	seek(out, 0)
	read(out, String)
end

"""
	from_csv_text(csv_text)

Create a dataframe from the string represenation of a CSV

# Examples

	from_csv_text("lat1,lon1,lat2,lon2\\n53.685335,-0.416387,53.479544,-2.954101\\n53.785335,-0.426387,53.479544,-2.954101\\n53.885335,-0.436387,53.479544,-2.954101\\n")

	3×4 DataFrame
	Row │ lat1     lon1       lat2     lon2
		│ Float64  Float64    Float64  Float64
   ─────┼──────────────────────────────────────
	  1 │ 53.6853  -0.416387  53.4795  -2.9541
	  2 │ 53.7853  -0.426387  53.4795  -2.9541
	  3 │ 53.8853  -0.436387  53.4795  -2.9541

"""
from_csv_text(csv_text) = CSV.read(IOBuffer(csv_text), DataFrame)

"""
	concat!(df, args...)

All of the dataframes supplied as args, appended to df
"""
concat!(df, args...) = reduce((adf, df)->append!(adf, df), args, init=df)

"""
	leftjoiner(df, on, args...; kw...)

# Arguments
- `df` The initial DataFrame
- `on` The field(s) to join with
- `args...` The DataFrames to join to the df
- `kw...` KW args passed to leftjoin
"""
leftjoiner(df, on, args...; kw...) = reduce((accumdf, nextdf)->leftjoin(accumdf, nextdf, on=on, kw...), args, init=df)

"""
	tryRename(df, renamelist)
	tryRename!(df, renamelist)

Attempts to renames the given fields, however, unlike DataFrames.rename, doesn't complain if the column doesn't exist

"""
tryRename(df, renamelist::Pair) = tryRename(df, [renamelist])
tryRename(df, renamelist::Vector{Pair{T,T}}) where T <: Any = tryRename(df, map(p->Symbol(p[1])=>Symbol(p[2]), renamelist))
tryRename(df, renamelist::Vector{Pair{Symbol, Symbol}}) = tryRenameFn(df, filter(p->p[1] in map(Symbol, names(df)), renamelist), rename)

tryRename!(df, renamelist::Pair) = tryRename(df, [renamelist])
tryRename!(df, renamelist::Vector{Pair{T,T}}) where T <: Any = tryRename(df, map(p->Symbol(p[1])=>Symbol(p[2]), renamelist))
tryRename!(df, renamelist::Vector{Pair{Symbol, Symbol}}) = tryRenameFn(df, filter(p->p[1] in map(Symbol, names(df)), renamelist), rename!)
tryRenameFn(df, renamelist, fn) = fn(df, renamelist)

"""
	dropByName!(df, fn)
	dropByName(df, fn)

Drop the columns whose names return true for the given function
# Example
	dropByName!(orders, n->startswith(n, "Sale") || startswith(n, "Order"))
"""
dropByName!(df, fn) = select!(df, Not(map(fn, names(df))))
dropByName(df, fn) = select(df, Not(map(fn, names(df))))

"""
	colsToGrid(df, x, y, v, col1name)

For a given DataFrame, turn the values of one of the columns into the name of the output columns
and use a given column as the values. That's harder to describe than illustrate!

# Examples

	julia> df = DataFrame([["a", "b", "c"], ["d", "e", "f"], [1,2,3], [4,5,6]], ["c1", "c2", "c3", "v"])
	3×4 DataFrame
	Row │ c1      c2      c3     v
		│ String  String  Int64  Int64
   ─────┼──────────────────────────────
	  1 │ a       d           1      4
	  2 │ b       e           2      5
	  3 │ c       f           3      6

   julia> colsToGrid(df, :c1, :c2, :v, "see")

   3×4 DataFrame
   Row │ see     d        e        f
	   │ String  Int64?   Int64?   Int64?
  ─────┼───────────────────────────────────
	 1 │ a             4  missing  missing
	 2 │ b       missing        5  missing
	 3 │ c       missing  missing        6

	 #  result = DataFrame([["a", "b", "c"], [4, missing, missing], [missing, 5, missing], [missing, missing, 6]], ["see", "d", "e", "f"])
"""
function colsToGrid(df, x, y, v, col1name)
	xs = unique(df[!, x]) |> sort
	ys = unique(df[!, y]) |> sort
	xref = Dict(xx=>Dict([r[y]=>r[v] for r in eachrow(filter(row->row[x] == xx, df))]) for xx in xs)
	grid = DataFrame()
	grid[!, col1name] = xs
	for yy in ys
		grid[!, yy] = map(k->get(xref[k], yy, missing), xs)
	end
	grid
end

###
end