text:
"""
    Griddly.load!(gdy_reader::GDYReader,path_yaml_file::String)
The output of this function is the corresponding grid of your GDY file. The
GDYReader will also link the images and shaders ressources to the objects in your
grid.

```jldoctest
julia> using Griddly

julia> image_path = joinpath(@__DIR__,"..","resources","images");

julia> shader_path = joinpath(@__DIR__,"..","resources","shaders");

julia> gdy_path = joinpath(@__DIR__,"..","resources","games");

julia> gdy_reader = Griddly.GDYReader(image_path,shader_path);

julia> grid = Griddly.load!(gdy_reader,joinpath(gdy_path,"Single-Player/GVGAI/sokoban.yaml"))
```
"""
Griddly.load!

"""
    Griddly.load_string!(gdy_reader::GDYReader,gdy_string::String)
Same as Griddly.load!() except that instead of providing the path to your game
GDY file you can give as input the corresponding string
"""
Griddly.load_string!

"""
    Griddly.get_player_count(grid)
Will return the number of player required to play any game from the grid
"""
Griddly.get_player_count

"""
    Griddly.get_all_available_actions(grid)
This will return all the actions name enable in this grid. Actions available for
the player or for other objects in the grid (example: ennemies random walk)
"""
Griddly.get_all_available_actions

"""
    Griddly.get_player_available_actions(grid)
Get action names that the player can execute
"""
Griddly.get_player_available_actions

"""
    Griddly.get_non_player_available_actions(grid)
Get action names that only other objects than the agent in the grid can do
"""
Griddly.get_non_player_available_actions

"""
    Griddly.get_input_ids(grid,action_name)
Each action can have several sub-actions. For instance, your action name can be
"move" and you can go LEFT,UP,RIGHT or DOWN. With this example, the output will
be [1,2,3,4] corresponding to each directions
"""
Griddly.get_input_ids

"""
    Griddly.get_player_actions_dict(grid)
This will get you a dictionnary where the keys are the action name available for
an agent and where corresponding value is the list of possible ids for this action.
"""
Griddly.get_player_actions_dict

"""
    Griddly.get_map_object_ids_char(grid)
Will get you the mapping between an object id and the corresponding character used
to represent it in the GDY file
"""
Griddly.get_map_object_ids_char

"""
    Griddly.get_avatar_object(grid)
If the avatar's name is defined in your GDY file it will return it
"""
Griddly.get_avatar_object

"""
    Griddly.load_level!(grid,lvl_number)
This will set the current level in the grid to the level index corresponding to
lvl number in your GDY file
"""
Griddly.load_level!

"""
    Griddly.load_level_string!(grid,lvl_string)
This will set the current level in the grid define by the string you provided. For
instance if we say we have load the sokoban grid as in the Griddly.load! example
we can do:
```jldoctest
julia> lvl_str ="       wwwww
                        w...w
                        wh.hw
                        wh.bw
                        wb.Aw
                        w...w
                        wwwww
                "
julia> Griddly.load_level_string!(grid,lvl_string)
```
"""
Griddly.load_level_string!

"""
    Griddly.vector_obs(grid)
This will return a Vector observation of the current state of the grid. If you
want to know more about Vector observation see the Observer section.
"""
Griddly.vector_obs

"""
    Griddly.create_game(grid,observer::ObserverType)
Will instantiate a game from a grid and define the observer type.
```jldoctest
julia> game = Griddly.create_game(grid,Griddly.SPRITE_2D)
```
"""
Griddly.create_game

"""
    Griddly.register_player!(game,player_name,observer_type)
This function will create and register a player in the current game. Especially,
a player can get his own __ObserverType__. Moreover, the player only got a partial
observation of the grid.
```jldoctest
julia> player = Griddly.register_player!(game,"Tux",Griddly.BLOCK_2D)
```
"""
Griddly.register_player!

"""
    Griddly.get_num_players(game)
This will get you the number of player registered for your game instance
"""
Griddly.get_num_players

"""
    Griddly.init!(game)
Will initialize the game, connect the dot with the grid which you created
your game from, it also configure both game observer(global observer) and
all players observer (partial observer).
"""
Griddly.init!

"""
    Griddly.reset!(game)
Will reset the grid to the initial configuration of the last level you loaded in.
Will also reset positions and parameters of all the observers.
"""
Griddly.reset!

"""
    Griddly.observe(game)
Will return a NumpyWrapper object of, depending on your __ObserverType__, representing the current state of the grid.
If you want to get the data of those observation you can use: `Griddly.get_data(observation)`. If you only want to get the shape
of your observations you can use: `Griddly.get_shape(observation)`. The shape will depends on your __ObserverType__.

Example:
```jldoctest
julia> using Griddly

julia> image_path = joinpath(@__DIR__,"..","resources","images");

julia> shader_path = joinpath(@__DIR__,"..","resources","shaders");

julia> gdy_path = joinpath(@__DIR__,"..","resources","games");

julia> gdy_reader = Griddly.GDYReader(image_path,shader_path);

julia> grid = Griddly.load!(gdy_reader,joinpath(gdy_path,"Single-Player/GVGAI/sokoban.yaml"))

julia> game = Griddly.create_game(grid,Griddly.SPRITE_2D)

julia> player1 = Griddly.register_player!(game,"Tux", Griddly.BLOCK_2D)

julia> Griddly.init!(game)

julia> Griddly.load_level!(grid,1)

julia> Griddly.reset!(game)

julia> initial_sprites = convert(Array{Int,3},Griddly.get_data(Griddly.observe(game))))
```
"""
Griddly.observe(::Griddly.GameProcess)

"""
    Griddly.step_player!(player,action_name,[actions_array])
This will make your player perform the action defined by action name with the ids
of actions array. For instance, let's say you want your player to move to the LEFT
and then UP you will do:
```jldoctest
julia> Griddly.step_player!(player1,"move",[1 2])
```
Indeed, the 1st id of "move" is LEFT and 2 is UP.
"""
Griddly.step_player!

"""
    Griddly.observe(player)
Will return a NumpyWrapper object of, depending on the __ObserverType__ you
attached to your player, representing partially (from player point of view) the
current state of the grid. If you want to get the data of those observation you
can use: `Griddly.get_data(observation)`.If you only want to get the shape of your
observations you can use: `Griddly.get_shape(observation)`. The shape will depends
on your __ObserverType__.

Example:
```jldoctest
julia> using Griddly

julia> image_path = joinpath(@__DIR__,"..","resources","images");

julia> shader_path = joinpath(@__DIR__,"..","resources","shaders");

julia> gdy_path = joinpath(@__DIR__,"..","resources","games");

julia> gdy_reader = Griddly.GDYReader(image_path,shader_path);

julia> grid = Griddly.load!(gdy_reader,joinpath(gdy_path,"Single-Player/GVGAI/sokoban.yaml"))

julia> game = Griddly.create_game(grid,Griddly.SPRITE_2D)

julia> player1 = Griddly.register_player!(game,"Tux", Griddly.BLOCK_2D)

julia> Griddly.init!(game)

julia> Griddly.load_level!(grid,1)

julia> Griddly.reset!(game)

julia> initial_player_tiles = convert(Array{Int,3},Griddly.get_data(Griddly.observe(player1))))
```
"""
Griddly.observe(::Griddly.Player)


text:
using GtkBuilderAid
using Gtk
using Cairo
using Base.Test

# Force loading a few things into memory
# Need to do this to make the builder work in early versions
dummy = GtkWindow("DUMMY")
dummy_box = GtkBox(:v)
push!(dummy_box, GtkLabel("dummy"))
dummy_button_box = GtkButtonBox(:h)
push!(dummy_box, dummy_button_box)
push!(dummy_button_box, GtkButton("dummy"))
destroy(dummy)

include("aidbuild.jl")
include("misc.jl")

ext: jl


text:
path = dirname(@__FILE__)

# Set up
m = AnSchorfheide(testing = true)

kal1, kal2 = JLD2.jldopen("$path/../reference/kalman_cat_args.jld2", "r") do file
    read(file, "kal1"), read(file, "kal2")
end

# Concatenate Kalmans
kal12 = cat(m, kal1, kal2)

# Test equality
exp_kal12 = JLD2.jldopen("$path/../reference/kalman_cat_out.jld2", "r") do file
    read(file, "kal12")
end

@testset "Testing Kalman output concatenation" begin
    for arg in fieldnames(typeof(kal1))
        @test exp_kal12[arg] ≈ kal12[arg]
    end
end

nothing


text:
using JuLIP
using JuLIP.Potentials

println("-------------------------------------------------")
println("   Variable Cell Test")
println("-------------------------------------------------")
calc = lennardjones(r0=rnn(:Al))
at = set_pbc!(bulk(:Al) * 2, true)
set_calculator!(at, calc)

println("Check atoms deform correctly with the cell")
# perturb the cell shape
set_defm!(at, defm(at) + 0.1*rand(JMatF), updatepositions=true)
# check that under this deformation the atom positions are still
# in a homogeneous lattice (i.e. they are deformed correctly with the cell)
@test maxnorm(forces(at)) < 1e-12
# now perturb the atom positions as well
rattle!(at, 0.1)

# set the constraint >>> this means the deformed cell defines F0 and X0
set_constraint!(at, VariableCell(at))

print("check that energy, forces, virial, stress, dofs, gradient evaluate ... ")
energy(at)
forces(at)
gradient(at)
# JuLIP.hessian_pos(at)
# JuLIP.hessian(at)     # this has no implementation for variable cells yet
virial(at)
@test stress(at) == - virial(at) / det(defm(at))
println("ok.")

print("Check that setting and getting dofs is consistent ... ")
x = dofs(at)
@assert length(x) == 3 * length(at) + 9
y = x + 0.01 * rand(size(x))
set_dofs!(at, y)
z = dofs(at)
@test vecnorm(y - z, Inf) < 1e-14
println("ok")

# now perform the finite-difference test
set_dofs!(at, x)
@test JuLIP.Testing.fdtest(calc, at, verbose=true, rattle=0.1)

println("Check virial for SW potential")
si = bulk(:Si, cubic=true) * 2
rattle!(si, 0.1)
set_defm!(si, defm(si) + 0.03 * rand(JMatF))
sw = StillingerWeber()
set_constraint!(si, VariableCell(si))
@test JuLIP.Testing.fdtest(calc, si, verbose=true, rattle=0.1)

println("-------------------------------------------------")
println("Test optimisation with VariableCell")
# start with a clean `at`
at = bulk(:Al) * 2   # cubic=true,
set_calculator!(at, calc)
set_constraint!(at, VariableCell(at))

println("For the initial state, stress/virial is far from 0:")
@show vecnorm(virial(at), Inf)
JuLIP.Solve.minimise!(at, verbose=2)
println("After optimisation, stress/virial should be 0:")
@show vecnorm(virial(at), Inf)
@test vecnorm(virial(at), Inf) < 1e-4

println("-------------------------------------------------")
println("And now with pressure . . .")
set_constraint!(at, VariableCell(at, pressure=10.0123))
JuLIP.Testing.fdtest(calc, at, verbose=true, rattle=0.1)
at = bulk(:Al) * 2
set_calculator!(at, calc)
set_constraint!(at, VariableCell(at, pressure=0.01))
JuLIP.Solve.minimise!(at, verbose = 2)
@show vecnorm(virial(at), Inf)
@show vecnorm(gradient(at), Inf)
@test vecnorm(gradient(at), Inf) < 1e-4
println("note it is correct that virial is O(1) since we applied pressure")


text:
##############################################################################
##
## @col
##
##############################################################################

"""
    @col(kw)

`@col` transforms an expression of the form `:z = :x + :y` into it's equivalent in
DataFrames's `source => fun => destination` syntax.

### Details

Parsing follows the same convention as other DataFramesMeta.jl macros, such as `@with`. All
terms in the expression that are `Symbol`s are treated as columns in the data frame, except
`Symbol`s wrapped in `^`. To use a variable representing a column name, wrap the variable
in `cols`.

`@col` constructs an anonymous function `fun` based on the given expression. It then creates
a `source => fun => destination` pair that is suitable for the `select`, `transform`, and
`combine` functions in DataFrames.jl.

### Examples

```julia
julia> @col :z = :x + :y
[:x, :y] => (##595 => :z)
```

In the above example, `##595` is an anonymous function equivalent to the following

```julia
(_x, _y) -> _x + _y
```

```jldoctest
julia> using DataFramesMeta;

julia> df = DataFrame(x = [1, 2], y = [3, 4]);

julia> import DataFramesMeta: @col;

julia> DataFrames.transform(df, @col :z = :x .* :y)
2×3 DataFrame
│ Row │ x     │ y     │ z     │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 1     │ 3     │ 3     │
│ 2   │ 2     │ 4     │ 8     │

```
"""
macro col(kw)
    esc(fun_to_vec(kw))
end

##############################################################################
##
## @byrow
##
##############################################################################
"""
    @byrow

Broadcast operations within DataFramesMeta.jl macros.

`@byrow` is not a "real" Julia macro but rather serves as a "flag"
to indicate that the anonymous function created by DataFramesMeta
to represent an operation should be applied "by-row".

If an expression starts with `@byrow`, either of the form `@byrow :y = f(:x)`
in transformations or `@byrow f(:x)` in `@orderby`, `@where`, and `@with`,
then the anonymous function created by DataFramesMeta is wrapped in the
`DataFrames.ByRow` function wrapper, which broadcasts the function so that it run on each row.

### Examples

```julia
julia> df = DataFrame(a = [1, 2, 3, 4], b = [5, 6, 7, 8]);

julia> @transform(df, @byrow :c = :a * :b)
4×3 DataFrame
 Row │ a      b      c
     │ Int64  Int64  Int64
─────┼─────────────────────
   1 │     1      5      5
   2 │     2      6     12
   3 │     3      7     21
   4 │     4      8     32

julia> @where(df, @byrow :a == 1 ? true : false)
1×2 DataFrame
 Row │ a      b
     │ Int64  Int64
─────┼──────────────
   1 │     1      5
```

To avoid writing `@byrow` multiple times when performing multiple
operations, it is allowed to use`@byrow` at the beginning of a block of
operations. All transformations in the block will operate by row.

```julia
julia> @where df @byrow begin
           :a > 1
           :b < 5
       end
1×2 DataFrame
 Row │ a      b
     │ Int64  Int64
─────┼──────────────
   1 │     2      4
```

### Comparison with `@eachrow`

To re-cap, the `@eachrow` macro roughly transforms

```julia
@eachrow df begin
    :a * :b
end
```

to

```julia
begin
    function tempfun(a, b)
        for i in eachindex(a)
            a[i] * b[i]
        end
    end
    tempfun(df.a, df.b)
    df
end
```

The function `*` is applied by-row. But the result of those operations
is not stored anywhere, as with `for`-loops in Base Julia.
Rather, `@eachrow` and `@eachrow!` return data frames.

Now consider `@byrow`. `@byrow` transforms

```julia
@with df @byrow begin
    :a * :b
end
```

to

```julia
tempfun(a, b) = a * b
tempfun.(df.a, df.b)
```

In contrast to `@eachrow`, `@with` combined with `@byrow` returns a vector of the
broadcasted multiplication and not a data frame.

Additionally, transformations applied using `@eachrow!` modify the input
data frame. On the contrary, `@byrow` does not update columns.

```julia
julia> df = DataFrame(a = [1, 2], b = [3, 4]);

julia> @with df @byrow begin
           :a = 500
       end
2-element Vector{Int64}:
 500
 500

julia> df
2×2 DataFrame
 Row │ a      b
     │ Int64  Int64
─────┼──────────────
   1 │     1      3
   2 │     2      4
```

### Comparison with `@.` and Base broadcasting

Base Julia provides the broadasting macro `@.` and in many cases `@.`
and `@byrow` will give equivalent results. But there are important
deviations in behavior. Consider the setup

```julia
df = DataFrame(a = [1, 2], b = [3, 4])
```

* Control flow. `@byrow` allows for operations of the form `if ... else`
  and `a ? b : c` to be applied by row. These expressions cannot be
  broadcasted in Base Julia. `@byrow` also allows for expressions of
  the form `a && b` and `a || b` to be applied by row, something that
  is not possible in Julia versions below 1.7.

```
julia> @with df @byrow begin
           if :a == 1
               5
           else
               10
           end
       end
2-element Vector{Int64}:
  5
 10

julia> @with df @. begin
           if :a == 1
               5
           else
               10
           end
       end # will error
```

* Broadcasting objects that are not columns. `@byrow` constructs an
  anonymous function which accepts only the columns of the input data frame
  and broadcasts that function. Consequently, it does not broadcast
  referenced objects which are not columns.

```julia
julia> df = DataFrame(a = [1, 2], b = [3, 4]);
julia> @with df @byrow :x + [5, 6]
```

  will error, because the `:x` in the above expression refers
  to a scalar `Int`, and you cannot do `1 + [5, 6]`.

  On the other hand

```julia
@with df @. :x + [5, 6]
```

  will succeed, as `df.x` is a 2-element vector as is `[5, 6]`.

  Because `ByRow` inside `transform` blocks does not internally
  use broadcasting in all circumstances, in the rare instance
  that a column in a data frame is a custom vector type that
  implements custom broadcasting, this custom behavior will
  not be called with `@byrow`.

* Broadcasting expensive calls. In Base Julia, broadcasting
  evaluates calls first and then broadcasts the result. Because
  `@byrow` constructs an anonymous function and evaluates
  that function for every row in the data frame, expensive functions
  will be evaluated many times.

```julia
julia> function expensive()
           sleep(.5)
           return 1
       end;

julia> @time @with df @byrow :a + expensive();
  1.037073 seconds (51.67 k allocations: 3.035 MiB, 3.19% compilation time)

julia> @time @with df :a .+ expensive();
  0.539900 seconds (110.67 k allocations: 6.525 MiB, 7.05% compilation time)

```

  This problem comes up when using the `@.` macro as well,
  but can easily be fixed with `\$`.

```julia
julia> @time @with df @. :a + expensive();
  1.036888 seconds (97.55 k allocations: 5.617 MiB, 3.20% compilation time)

julia> @time @with df @. :a + \$expensive();
  0.537961 seconds (110.68 k allocations: 6.525 MiB, 6.73% compilation time)
```

  No such solution currently exists with `@byrow`.
"""
macro byrow(args...)
    throw(ArgumentError("@byrow is deprecated outside of DataFramesMeta macros."))
end

##############################################################################
##
## @with
##
##############################################################################

function exec(df, p::Pair)
    cols = first(p)
    fun = last(p)
    fun(map(c -> DataFramesMeta.getsinglecolumn(df, c), cols)...)
end
exec(df, s::Union{Symbol, AbstractString}) = df[!, s]

getsinglecolumn(df, s::DataFrames.ColumnIndex) = df[!, s]
getsinglecolumn(df, s) = throw(ArgumentError("Only indexing with Symbols, strings and integers " *
    "is currently allowed with cols"))

function with_helper(d, body)
    # Make body an expression to force the
    # complicated method of fun_to_vec
    # in the case of QuoteNode
    t = fun_to_vec(Expr(:block, body); no_dest=true)
    :(DataFramesMeta.exec($d, $t))
end

"""
    @with(d, expr)

`@with` allows DataFrame columns keys to be referenced as symbols.

### Arguments

* `d` : an AbstractDataFrame type
* `expr` : the expression to evaluate in `d`

### Details

`@with` works by parsing the expression body for all columns indicated
by symbols (e.g. `:colA`). Then, a function is created that wraps the
body and passes the columns as function arguments. This function is
then called. Operations are efficient because:

- A pseudo-anonymous function is defined, so types are stable.
- Columns are passed as references, eliminating DataFrame indexing.

The following

```julia
@with(d, :a .+ :b .+ 1)
```

becomes

```julia
tempfun(a, b) = a .+ b .+ 1
tempfun(d[!, :a], d[!, :b])
```

If an expression is wrapped in `^(expr)`, `expr` gets passed through untouched.
If an expression is wrapped in  `cols(expr)`, the column is referenced by the
variable `expr` rather than a symbol.

If the expression provide to `@with` begins with `@byrow`, the function
created by the `@with` block is broadcasted along the columns of the
data frame.

### Examples

```jldoctest
julia> using DataFramesMeta

julia> y = 3;

julia> df = DataFrame(x = 1:3, y = [2, 1, 2]);

julia> x = [2, 1, 0];

julia> @with(df, :y .+ 1)
3-element Array{Int64,1}:
 3
 2
 3

julia> @with(df, :x + x)
3-element Array{Int64,1}:
 3
 3
 3

julia> @with df begin
            res = 0.0
            for i in 1:length(:x)
                res += :x[i] * :y[i]
            end
            res
        end
10.0

julia> @with(df, df[:x .> 1, ^(:y)]) # The ^ means leave the :y alone
2-element Array{Int64,1}:
 1
 2

julia> colref = :x;

julia> @with(df, :y + cols(colref)) # Equivalent to df[!, :y] + df[!, colref]
3-element Array{Int64,1}:
 3
 3
 5

julia> @with df @byrow :x * :y
3-element Vector{Int64}:
 2
 2
 6

```

!!! note
    `@with` creates a function, so the scope within `@with` is a local scope.
    Variables in the parent can be read. Writing to variables in the parent scope
    differs depending on the type of scope of the parent. If the parent scope is a
    global scope, then a variable cannot be assigned without using the `global` keyword.
    If the parent scope is a local scope (inside a function or let block for example),
    the `global` keyword is not needed to assign to that parent scope.
"""
macro with(d, body)
    esc(with_helper(d, body))
end

##############################################################################
##
## @subset and subset! - select row subsets
##
##############################################################################

function subset_helper(x, args...)
    exprs, outer_flags = create_args_vector(args...)
    t = (fun_to_vec(ex; no_dest=true, outer_flags=outer_flags) for ex in exprs)
    quote
        $subset($x, $(t...); skipmissing=true)
    end
end

function where_helper(x, args...)
    exprs, outer_flags = create_args_vector(args...)
    t = (fun_to_vec(ex; no_dest=true, outer_flags=outer_flags) for ex in exprs)
    quote
        $subset($x, $(t...); skipmissing=true)
    end
end

"""
    @subset(d, i...)

Select row subsets in `AbstractDataFrame`s and `GroupedDataFrame`s.

### Arguments

* `d` : an AbstractDataFrame or GroupedDataFrame
* `i...` : expression for selecting rows

Multiple `i` expressions are "and-ed" together.

If given a `GroupedDataFrame`, `@subset` applies transformations by
group, and returns a fresh `DataFrame` containing the rows
for which the generated values are all `true`.

Inputs to `@subset` can come in two formats: a `begin ... end` block, in which case each
line is a separate selector, or as multiple arguments.
For example the following two statements are equivalent:

```julia
@subset df begin
    :x .> 1
    :y .< 2
end
```

and

```
@subset(df, :x .> 1, :y .< 2)
```

!!! note
    `@subset` treats `missing` values as `false` when filtering rows.
    Unlike `DataFrames.subset` and other Boolean operations with
    `missing`, `@subset` will *not* error on missing values, and
    will only keep `true` values.

If an expression provided to `@subset` begins with `@byrow`, operations
are applied "by row" along the data frame. To avoid writing `@byrow` multiple
times, `@orderby` also allows `@byrow`to be placed at the beginning of a block of
operations. For example, the following two statements are equivalent.

```
@subset df @byrow begin
    :x > 1
    :y < 2
end
```

and

```
@subset df
    @byrow :x > 1
    @byrow :y < 2
end
```

### Examples

```jldoctest
julia> using DataFramesMeta, Statistics

julia> df = DataFrame(x = 1:3, y = [2, 1, 2]);

julia> globalvar = [2, 1, 0];

julia> @subset(df, :x .> 1)
2×2 DataFrame
 Row │ x      y
     │ Int64  Int64
─────┼──────────────
   1 │     2      1
   2 │     3      2

julia> @subset(df, :x .> globalvar)
2×2 DataFrame
 Row │ x      y
     │ Int64  Int64
─────┼──────────────
   1 │     2      1
   2 │     3      2

julia> @subset df begin
    :x .> globalvar
    :y .== 3
end
0×2 DataFrame

julia> d = DataFrame(n = 1:20, x = [3, 3, 3, 3, 1, 1, 1, 2, 1, 1,
                                    2, 1, 1, 2, 2, 2, 3, 1, 1, 2]);

julia> g = groupby(d, :x);

julia> @subset(g, :n .> mean(:n))
8×2 DataFrame
 Row │ n      x
     │ Int64  Int64
─────┼──────────────
   1 │    12      1
   2 │    13      1
   3 │    15      2
   4 │    16      2
   5 │    17      3
   6 │    18      1
   7 │    19      1
   8 │    20      2

julia> @subset g begin
           :n .> mean(:n)
           :n .< 20
       end
7×2 DataFrame
 Row │ n      x
     │ Int64  Int64
─────┼──────────────
   1 │    12      1
   2 │    13      1
   3 │    15      2
   4 │    16      2
   5 │    17      3
   6 │    18      1
   7 │    19      1

julia> d = DataFrame(a = [1, 2, missing], b = ["x", "y", missing]);

julia> @subset(d, :a .== 1)
1×2 DataFrame
│ Row │ a      │ b       │
│     │ Int64? │ String? │
├─────┼────────┼─────────┤
│ 1   │ 1      │ x       │
```
"""
macro subset(x, args...)
    esc(subset_helper(x, args...))
end

"""
    @where(x, args...)

Deprecated version of `@subset`, see `?@subset` for details.
"""
macro where(x, args...)
    @warn "`@where is deprecated, use `@subset`  with `@skipmissing` instead."
    esc(where_helper(x, args...))
end

function subset!_helper(x, args...)
    exprs, outer_flags = create_args_vector(args...)
    t = (fun_to_vec(ex; no_dest=true, outer_flags=outer_flags) for ex in exprs)
    quote
        $subset!($x, $(t...); skipmissing=true)
    end
end

"""
    @subset!(d, i...)

Select row subsets in `AbstractDataFrame`s and `GroupedDataFrame`s,
mutating the underlying data-frame in-place.

### Arguments

* `d` : an AbstractDataFrame or GroupedDataFrame
* `i...` : expression for selecting rows

Multiple `i` expressions are "and-ed" together.

If given a `GroupedDataFrame`, `@subset!` applies transformations by
group, and returns a fresh `DataFrame` containing the rows
for which the generated values are all `true`.

Inputs to `@subset!` can come in two formats: a `begin ... end` block, in which case each
line is a separate selector, or as multiple arguments.
For example the following two statements are equivalent:

```julia
@subset! df begin
    :x .> 1
    :y .< 2
end
```

and

```
@subset!(df, :x .> 1, :y .< 2)
```

!!! note
    `@subset!` treats `missing` values as `false` when filtering rows.
    Unlike `DataFrames.subset!` and other Boolean operations with
    `missing`, `@subset!` will *not* error on missing values, and
    will only keep `true` values.

If an expression provided to `@subset!` begins with `@byrow`, operations
are applied "by row" along the data frame. To avoid writing `@byrow` multiple
times, `@orderby` also allows `@byrow`to be placed at the beginning of a block of
operations. For example, the following two statements are equivalent.

```
@subset! df @byrow begin
    :x > 1
    :y < 2
end
```

and

```
@subset! df
    @byrow :x > 1
    @byrow :y < 2
end
```

### Examples

```jldoctest
julia> using DataFramesMeta, Statistics

julia> df = DataFrame(x = 1:3, y = [2, 1, 2]);

julia> globalvar = [2, 1, 0];

julia> @subset!(df, :x .> 1)
2×2 DataFrame
 Row │ x      y
     │ Int64  Int64
─────┼──────────────
   1 │     2      1
   2 │     3      2

julia> @subset!(df, :x .> globalvar)
2×2 DataFrame
 Row │ x      y
     │ Int64  Int64
─────┼──────────────
   1 │     2      1
   2 │     3      2

julia> @subset! df begin
    :x .> globalvar
    :y .== 3
end
0×2 DataFrame

julia> d = DataFrame(n = 1:20, x = [3, 3, 3, 3, 1, 1, 1, 2, 1, 1,
                                    2, 1, 1, 2, 2, 2, 3, 1, 1, 2]);

julia> g = groupby(d, :x);

julia> @subset!(g, :n .> mean(:n))
8×2 DataFrame
 Row │ n      x
     │ Int64  Int64
─────┼──────────────
   1 │    12      1
   2 │    13      1
   3 │    15      2
   4 │    16      2
   5 │    17      3
   6 │    18      1
   7 │    19      1
   8 │    20      2

julia> @subset! g begin
           :n .> mean(:n)
           :n .< 20
       end
7×2 DataFrame
 Row │ n      x
     │ Int64  Int64
─────┼──────────────
   1 │    12      1
   2 │    13      1
   3 │    15      2
   4 │    16      2
   5 │    17      3
   6 │    18      1
   7 │    19      1

julia> d = DataFrame(a = [1, 2, missing], b = ["x", "y", missing]);

julia> @subset!(d, :a .== 1)
1×2 DataFrame
│ Row │ a      │ b       │
│     │ Int64? │ String? │
├─────┼────────┼─────────┤
│ 1   │ 1      │ x       │
```
"""
macro subset!(x, args...)
    esc(subset!_helper(x, args...))
end

##############################################################################
##
## @orderby
##
##############################################################################

function orderby_helper(x, args...)
    exprs, outer_flags = create_args_vector(args...)
    t = (fun_to_vec(ex; gensym_names = true, outer_flags = outer_flags) for ex in exprs)
    quote
        $DataFramesMeta.orderby($x, $(t...))
    end
end

function orderby(x::AbstractDataFrame, @nospecialize(args...))
    t = DataFrames.select(x, args...; copycols = false)
    x[sortperm(t), :]
end

function orderby(x::GroupedDataFrame, @nospecialize(args...))
    throw(ArgumentError("@orderby with a GroupedDataFrame is reserved"))
end

function orderby(x::SubDataFrame, @nospecialize(args...))
    t = DataFrames.select(x, args...)
    x[sortperm(t), :]
end

"""
    @orderby(d, i...)

Sort rows by values in one of several columns or a transformation of columns.
Always returns a fresh `DataFrame`. Does not accept a `GroupedDataFrame`.

When given a `DataFrame`, `@orderby` applies the transformation
given by its arguments (but does not create new columns) and sorts
the given `DataFrame` on the result, returning a new `DataFrame`.

Inputs to `@orderby` can come in two formats: a `begin ... end` block, in which case each
line in the block is a separate ordering operation, and as mulitple
arguments. For example, the following two statements are equivalent:

```julia
@orderby df begin
    :x
    -:y
end
```

and

```
@orderby(df, :x, -:y)
```

### Arguments

* `d` : an AbstractDataFrame
* `i...` : expression for sorting

If an expression provided to `@orderby` begins with `@byrow`, operations
are applied "by row" along the data frame. To avoid writing `@byrow` multiple
times, `@orderby` also allows `@byrow`to be placed at the beginning of a block of
operations. For example, the following two statements are equivalent.

```
@orderby df @byrow begin
    :x^2
    :x^3
end
```

and

```
@orderby df
    @byrow :x^2
    @byrow :x^3
end
```

### Examples

```jldoctest
julia> using DataFramesMeta, Statistics

julia> d = DataFrame(x = [3, 3, 3, 2, 1, 1, 1, 2, 1, 1], n = 1:10,
                     c = ["a", "b", "c", "d", "e", "f", "g", "h", "i", "j"]);

julia> @orderby(d, -1 .* :n)
10×3 DataFrame
 Row │ x      n      c
     │ Int64  Int64  String
─────┼──────────────────────
   1 │     1     10  j
   2 │     1      9  i
   3 │     2      8  h
   4 │     1      7  g
   5 │     1      6  f
   6 │     1      5  e
   7 │     2      4  d
   8 │     3      3  c
   9 │     3      2  b
  10 │     3      1  a

julia> @orderby(d, sortperm(:c, rev = true))
10×3 DataFrame
 Row │ x      n      c
     │ Int64  Int64  String
─────┼──────────────────────
   1 │     1     10  j
   2 │     1      9  i
   3 │     2      8  h
   4 │     1      7  g
   5 │     1      6  f
   6 │     1      5  e
   7 │     2      4  d
   8 │     3      3  c
   9 │     3      2  b
  10 │     3      1  a

julia> @orderby d begin
    :x
    abs.(:n .- mean(:n))
end
10×3 DataFrame
 Row │ x      n      c
     │ Int64  Int64  String
─────┼──────────────────────
   1 │     1      5  e
   2 │     1      6  f
   3 │     1      7  g
   4 │     1      9  i
   5 │     1     10  j
   6 │     2      4  d
   7 │     2      8  h
   8 │     3      3  c
   9 │     3      2  b
  10 │     3      1  a

julia> @orderby d @byrow :x^2
10×3 DataFrame
 Row │ x      n      c
     │ Int64  Int64  String
─────┼──────────────────────
   1 │     1      5  e
   2 │     1      6  f
   3 │     1      7  g
   4 │     1      9  i
   5 │     1     10  j
   6 │     2      4  d
   7 │     2      8  h
   8 │     3      1  a
   9 │     3      2  b
  10 │     3      3  c
```
"""
macro orderby(d, args...)
    esc(orderby_helper(d, args...))
end

##############################################################################
##
## transform & @transform
##
##############################################################################

function transform_helper(x, args...)
    exprs, outer_flags = create_args_vector(args...)
    t = (fun_to_vec(ex; gensym_names = false, outer_flags = outer_flags) for ex in exprs)
    quote
        $DataFrames.transform($x, $(t...))
    end
end

"""
    @transform(d, i...)

Add additional columns or keys based on keyword arguments.

### Arguments

* `d` : an `AbstractDataFrame`, or `GroupedDataFrame`
* `i...` : keyword arguments defining new columns or keys

### Returns

* `::AbstractDataFrame` or `::GroupedDataFrame`

Inputs to `@transform` can come in two formats: a `begin ... end` block,
in which case each line in the block is a separate
transformation, (`:y = f(:x)`), or as a series of
keyword arguments. For example, the following are
equivalent:

```julia
@transform df begin
    :a = :x
    :b = :y
end
```

and

```
@transform(df, :a = :x, :b = :y)
```

`@transform` uses the syntax `@byrow` to wrap transformations in
the `ByRow` function wrapper from DataFrames, apply a function row-wise,
similar to broadcasting. For example, the call

```
@transform(df, @byrow :y = :x == 1 ? true : false)
```

becomes

```
transform(df, :x => ByRow(x -> x == 1 ? true : false) => :y)
```

a transformation which cannot be conveniently expressed
using broadcasting.

To avoid writing `@byrow` multiple times when performing multiple
transformations by row, `@transform` allows `@byrow` at the
beginning of a block of transformations (i.e. `@byrow begin... end`).
All transformations in the block will operate by row.

### Examples

```jldoctest
julia> using DataFramesMeta

julia> df = DataFrame(A = 1:3, B = [2, 1, 2]);

julia> @transform df begin
           :a = 2 * :A
           :x = :A .+ :B
       end
3×4 DataFrame
 Row │ A      B      a      x
     │ Int64  Int64  Int64  Int64
─────┼────────────────────────────
   1 │     1      2      2      3
   2 │     2      1      4      3
   3 │     3      2      6      5

julia> @transform df @byrow :z = :A * :B
3×3 DataFrame
 Row │ A      B      z
     │ Int64  Int64  Int64
─────┼─────────────────────
   1 │     1      2      2
   2 │     2      1      2
   3 │     3      2      6

julia> @transform df @byrow begin
           :x = :A * :B
           :y = :A == 1 ? 100 : 200
       end

3×4 DataFrame
 Row │ A      B      x      y
     │ Int64  Int64  Int64  Int64
─────┼────────────────────────────
   1 │     1      2      2    100
   2 │     2      1      2    200
   3 │     3      2      6    200

```
"""
macro transform(x, args...)
    esc(transform_helper(x, args...))
end

##############################################################################
##
## transform! & @transform!
##
##############################################################################

function transform!_helper(x, args...)
    exprs, outer_flags = create_args_vector(args...)
    t = (fun_to_vec(ex; gensym_names = false, outer_flags = outer_flags) for ex in exprs)
    quote
        $DataFrames.transform!($x, $(t...))
    end
end

"""
    @transform!(d, i...)

Mutate `d` inplace to add additional columns or keys based on keyword arguments and return it.
No copies of existing columns are made.

### Arguments

* `d` : an `AbstractDataFrame`, or `GroupedDataFrame`
* `i...` : keyword arguments defining new columns or keys

### Returns

* `::DataFrame`

Inputs to `@transform!` can come in two formats: a `begin ... end` block,
in which case each line in the block is a separate
transformation, (`:y = f(:x)`), or as a series of
keyword arguments. For example, the following are
equivalent:

```julia
@transform! df begin
    :a = :x
    :b = :y
end
```

and

```
@transform!(df, :a = :x, :b = :y)
```

`@transform!` uses the syntax `@byrow` to wrap transform!ations in
the `ByRow` function wrapper from DataFrames, apply a function row-wise,
similar to broadcasting. For example, the call

```
@transform!(df, @byrow :y = :x == 1 ? true : false)
```

becomes

```
transform!(df, :x => ByRow(x -> x == 1 ? true : false) => :y)
```

a transformation which cannot be conveniently expressed
using broadcasting.

To avoid writing `@byrow` multiple times when performing multiple
transform!ations by row, `@transform!` allows `@byrow` at the
beginning of a block of transform!ations (i.e. `@byrow begin... end`).
All transform!ations in the block will operate by row.

### Examples

```jldoctest
julia> using DataFramesMeta

julia> df = DataFrame(A = 1:3, B = [2, 1, 2]);

julia> df2 = @transform!(df, :a = 2 * :A, :x = :A .+ :B)
3×4 DataFrame
│ Row │ A     │ B     │ a     │ x     │
│     │ Int64 │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┼───────┤
│ 1   │ 1     │ 2     │ 2     │ 3     │
│ 2   │ 2     │ 1     │ 4     │ 3     │
│ 3   │ 3     │ 2     │ 6     │ 5     │

julia> df === df2
true
```
"""
macro transform!(x, args...)
    esc(transform!_helper(x, args...))
end

##############################################################################
##
## @select - select and transform columns
##
##############################################################################

function select_helper(x, args...)
    exprs, outer_flags = create_args_vector(args...)
    t = (fun_to_vec(ex; gensym_names = false, outer_flags = outer_flags) for ex in exprs)
    quote
        $DataFrames.select($x, $(t...))
    end
end

"""
    @select(d, e...)

Select and transform columns.

### Arguments

* `d` : an `AbstractDataFrame` or `GroupedDataFrame`
* `e` :  keyword arguments specifying new columns in terms of existing columns
  or symbols to specify existing columns

### Returns

* `::AbstractDataFrame`

Inputs to `@select` can come in two formats: a `begin ... end` block,
in which case each line in the block is a separate
transformation or selector, or as a series of
arguments and keyword arguments. For example, the following are
equivalent:

```julia
@select df begin
    :x
    :y = :a .+ :b
end
```

and

```
@select(df, :x, :y = :a .+ :b)
```

`@select` uses the syntax `@byrow` to wrap transformations in
the `ByRow` function wrapper from DataFrames, apply a function row-wise,
similar to broadcasting. For example, the call

```
@select(df, @byrow :y = :x == 1 ? true : false)
```

becomes

```
select(df, :x => ByRow(x -> x == 1 ? true : false) => :y)
```

a transformation which cannot be conveniently expressed
using broadcasting.

To avoid writing `@byrow` multiple times when performing multiple
transformations by row, `@select` allows `@byrow` at the
beginning of a block of selectations (i.e. `@byrow begin... end`).
All transformations in the block will operate by row.

### Examples

```jldoctest
julia> using DataFramesMeta

julia> df = DataFrame(a = repeat(1:4, outer = 2), b = repeat(2:-1:1, outer = 4), c = 1:8);

julia> @select(df, :c, :a)
8×2 DataFrame
 Row │ c      a
     │ Int64  Int64
─────┼──────────────
   1 │     1      1
   2 │     2      2
   3 │     3      3
   4 │     4      4
   5 │     5      1
   6 │     6      2
   7 │     7      3
   8 │     8      4

julia> @select df begin
           :c
           :x = :b + :c
       end
8×2 DataFrame
 Row │ c      x
     │ Int64  Int64
─────┼──────────────
   1 │     1      3
   2 │     2      3
   3 │     3      5
   4 │     4      5
   5 │     5      7
   6 │     6      7
   7 │     7      9
   8 │     8      9
```
"""
macro select(x, args...)
    esc(select_helper(x, args...))
end

##############################################################################
##
## @select! - in-place select and transform columns
##
##############################################################################

function select!_helper(x, args...)
    exprs, outer_flags = create_args_vector(args...)
    t = (fun_to_vec(ex; gensym_names = false, outer_flags = outer_flags) for ex in exprs)
    quote
        $DataFrames.select!($x, $(t...))
    end
end

"""
    @select!(d, e...)

Mutate `d` in-place to retain only columns or transformations specified by `e` and return it. No copies of existing columns are made.

### Arguments

* `d` : an AbstractDataFrame
* `e` :  keyword arguments specifying new columns in terms of existing columns
  or symbols to specify existing columns

### Returns

* `::DataFrame`

Inputs to `@select!` can come in two formats: a `begin ... end` block,
in which case each line in the block is a separate
transformation or selector, or as a series of
arguments and keyword arguments. For example, the following are
equivalent:

`@select!` uses the syntax `@byrow` to wrap transformations in
the `ByRow` function wrapper from DataFrames, apply a function row-wise,
similar to broadcasting. For example, the call

```
@select!(df, @byrow :y = :x == 1 ? true : false)
```

becomes

```
select!(df, :x => ByRow(x -> x == 1 ? true : false) => :y)
```

a transformation which cannot be conveniently expressed
using broadcasting.

To avoid writing `@byrow` multiple times when performing multiple
transformations by row, `@select!` allows `@byrow` at the
beginning of a block of select!ations (i.e. `@byrow begin... end`).
All transformations in the block will operate by row.

### Examples

```jldoctest
julia> using DataFrames, DataFramesMeta

julia> df = DataFrame(a = repeat(1:4, outer = 2), b = repeat(2:-1:1, outer = 4), c = 1:8);

julia> df2 = @select!(df, :c, :a)
8×2 DataFrame
 Row │ c      a
     │ Int64  Int64
─────┼──────────────
   1 │     1      1
   2 │     2      2
   3 │     3      3
   4 │     4      4
   5 │     5      1
   6 │     6      2
   7 │     7      3
   8 │     8      4

julia> df === df2
true

julia> df = DataFrame(a = repeat(1:4, outer = 2), b = repeat(2:-1:1, outer = 4), c = 1:8);

julia> df2 = @select! df begin
           :c
           :x = :b + :c
       end
8×2 DataFrame
 Row │ c      x
     │ Int64  Int64
─────┼──────────────
   1 │     1      3
   2 │     2      3
   3 │     3      5
   4 │     4      5
   5 │     5      7
   6 │     6      7
   7 │     7      9
   8 │     8      9

julia> df === df2
true
```
"""
macro select!(x, args...)
    esc(select!_helper(x, args...))
end

##############################################################################
##
## @combine - summarize a grouping operation
##
##############################################################################

function combine_helper(x, args...; deprecation_warning = false)
    deprecation_warning && @warn "`@based_on` is deprecated. Use `@combine` instead."

    exprs, outer_flags = create_args_vector(args...)

    fe = first(exprs)
    if length(exprs) == 1 &&
        !(fe isa QuoteNode || onearg(fe, :cols)) &&
        !(fe.head == :(=) || fe.head == :kw)

        @warn "Returning a Table object from @by and @combine now requires `cols(AsTable)` on the LHS."

        exprs = ((:(cols(AsTable) = $fe)),)
    end

    t = (fun_to_vec(ex; gensym_names = false, outer_flags = outer_flags) for ex in exprs)

    quote
        $DataFrames.combine($x, $(t...))
    end
end

"""
    @combine(x, args...)

Summarize a grouping operation

### Arguments

* `x` : a `GroupedDataFrame` or `AbstractDataFrame`
* `args...` : keyword arguments defining new columns

Inputs to `@combine` can come in two formats: a `begin ... end` block,
in which case each line in the block is a separate
transformation, or as a series of keyword arguments.
For example, the following are equivalent:

```
@combine df begin
    :mx = mean(:x)
    :sx = std(:x)
end
```

and

```
@combine(df, :mx = mean(:x), :sx = std(:x))
```

### Examples

```jldoctest
julia> using DataFramesMeta

julia> d = DataFrame(
            n = 1:20,
            x = [3, 3, 3, 3, 1, 1, 1, 2, 1, 1, 2, 1, 1, 2, 2, 2, 3, 1, 1, 2]);

julia> g = groupby(d, :x);

julia> @combine(g, :nsum = sum(:n))
3×2 DataFrame
 Row │ x      nsum
     │ Int64  Int64
─────┼──────────────
   1 │     1     99
   2 │     2     84
   3 │     3     27

julia> @combine g begin
           :x2 = 2 * :x
           :nsum = sum(:n)
       end
20×3 DataFrame
 Row │ x      x2     nsum
     │ Int64  Int64  Int64
─────┼─────────────────────
   1 │     1      2     99
   2 │     1      2     99
   3 │     1      2     99
   4 │     1      2     99
   5 │     1      2     99
   6 │     1      2     99
   7 │     1      2     99
   8 │     1      2     99
   9 │     1      2     99
  10 │     2      4     84
  11 │     2      4     84
  12 │     2      4     84
  13 │     2      4     84
  14 │     2      4     84
  15 │     2      4     84
  16 │     3      6     27
  17 │     3      6     27
  18 │     3      6     27
  19 │     3      6     27
  20 │     3      6     27

```
"""
macro combine(x, args...)
    esc(combine_helper(x, args...))
end

"""
    @based_on(d, args...)

Deprecated version of `@combine`, see: [`@combine`](@ref)
"""
macro based_on(x, args...)
    esc(combine_helper(x, args...; deprecation_warning = true))
end

##############################################################################
##
## @by - grouping
##
##############################################################################

function by_helper(x, what, args...)
    # Only allow one argument when returning a Table object
    exprs, outer_flags = create_args_vector(args...)
    fe = first(exprs)
    if length(exprs) == 1 &&
        !(fe isa QuoteNode || onearg(fe, :cols)) &&
        !(fe.head == :(=) || fe.head == :kw)

        @warn "Returning a Table object from @by and @combine now requires `cols(AsTable)` on the LHS."

        exprs = ((:(cols(AsTable) = $fe)),)
    end

    t = (fun_to_vec(ex; gensym_names = false, outer_flags = outer_flags) for ex in exprs)

    quote
        $DataFrames.combine($groupby($x, $what), $(t...))
    end
end

"""
    @by(d::AbstractDataFrame, cols, e...)

Split-apply-combine in one step.

### Arguments

* `d` : an AbstractDataFrame
* `cols` : a column indicator (Symbol, Int, Vector{Symbol}, etc.)
* `e` :  keyword arguments specifying new columns in terms of column groupings

### Returns

* `::DataFrame`

Transformation inputs to `@by` can come in two formats: a `begin ... end` block,
in which case each line in the block is a separate
transformation, or as a series of keyword arguments.
For example, the following are equivalent:

```
@by df :g begin
    :mx = mean(:x)
    :sx = std(:x)
end
```

and

```
@by(df, :g, mx = mean(:x), sx = std(:x))
```

### Examples

```jldoctest
julia> using DataFramesMeta, Statistics

julia> df = DataFrame(
            a = repeat(1:4, outer = 2),
            b = repeat(2:-1:1, outer = 4),
            c = 1:8);

julia> @by(df, :a, :d = sum(:c))
4×2 DataFrame
 Row │ a      d
     │ Int64  Int64
─────┼──────────────
   1 │     1      6
   2 │     2      8
   3 │     3     10
   4 │     4     12

julia> @by df :a begin
           :d = 2 * :c
       end
8×2 DataFrame
 Row │ a      d
     │ Int64  Int64
─────┼──────────────
   1 │     1      2
   2 │     1     10
   3 │     2      4
   4 │     2     12
   5 │     3      6
   6 │     3     14
   7 │     4      8
   8 │     4     16

julia> @by(df, :a, :c_sum = sum(:c), :c_mean = mean(:c))
4×3 DataFrame
 Row │ a      c_sum  c_mean
     │ Int64  Int64  Float64
─────┼───────────────────────
   1 │     1      6      3.0
   2 │     2      8      4.0
   3 │     3     10      5.0
   4 │     4     12      6.0

julia> @by df :a begin
           :c = :c
           :c_mean = mean(:c)
       end
8×3 DataFrame
 Row │ a      c      c_mean
     │ Int64  Int64  Float64
─────┼───────────────────────
   1 │     1      1      3.0
   2 │     1      5      3.0
   3 │     2      2      4.0
   4 │     2      6      4.0
   5 │     3      3      5.0
   6 │     3      7      5.0
   7 │     4      4      6.0
   8 │     4      8      6.0

```
"""
macro by(x, what, args...)
    esc(by_helper(x, what, args...))
end


text:
"""
  Implements a rate model corresponding to LIF-integrator
"""

struct RateParameters
    tau::Float64
    reset::Float64
    thresholdmean::Float64
    thresholds::Array{Array{Float64, 1}, 1}
    refractory_period::Float64
    spiketrace_timeconstant::Float64
    activation_function::Function
    plasticityrule
    in_amplitude::Float64
    target_amplitude::Float64
    learningrate::Float64
end
export RateParameters
RateParameters(; tau = 10., reset = 0., thresholdmean = 20., thresholds = Array{Array{Float64, 1}, 1}(),
             refractory_period = 2., spiketrace_timeconstant = 20., #refractory_period = 2.
             activation_function = LIFactivationfunction!, plasticityrule = deltarule,
               in_amplitude = 100., target_amplitude = 0.2, learningrate = 1e1) =
RateParameters(tau, reset, thresholdmean, thresholds, refractory_period,
           spiketrace_timeconstant, activation_function, plasticityrule,
           in_amplitude, target_amplitude, learningrate)

mutable struct RateNet
  layer_sizes::Array{Int64, 1} #numnber of neurons per layer
  v::Array{Array{Float64, 1}, 1} #membrane potential/current
  input::Array{Array{Float64, 1}, 1} #external input/bias
  rates::Array{Array{Float64, 1}, 1} #firing rate (nonlinear function of v)
  targets::Array{Float64, 1} #only for last layer
  weights::Array{Array{Float64, 2}, 1}
  parameters::RateParameters
end
export RateNet
function RateNet(layer_sizes, parameters = RateParametes())
  nl = length(layer_sizes)
  RateNet(layer_sizes,
  [zeros(layer_sizes[i]) for i in 1:nl],
  [zeros(layer_sizes[i]) for i in 1:nl],
  [zeros(layer_sizes[i]) for i in 1:nl],
  zeros(layer_sizes[end]),
  [randn(layer_sizes[i+1], layer_sizes[i]) ./ (sqrt(layer_sizes[i])) for i in 1:nl - 1],
  parameters)
end

function Linear!(input, output, parameters, layerindex)
  for i in 1:length(input)
    output[i] = input[i]
  end
end
export Linear!
function ReLU!(input, output, parameters, layerindex)
  for i in 1:length(input)
    output[i] = max(0, input[i] - parameters.thresholds[layerindex][i])
  end
end
export ReLU!
function LIFactivationfunction!(input, output, parameters, layerindex)
  for i in 1:length(input)
    if input[i] <= parameters.thresholds[layerindex][i]
      output[i] = 0.
    else #factor 1000 for unit Hz
      output[i] = 1000. / (parameters.refractory_period - parameters.tau *
                  log(1 - parameters.thresholds[layerindex][i] / input[i]))
    end
  end
end
export LIFactivationfunction!

function forwardprop!(net::RateNet)
  net.v[1] = net.input[1]
  net.parameters.activation_function(net.v[1], net.rates[1], net.parameters, 1)
	for i in 2:length(net.layer_sizes)
		BLAS.gemv!('N', 1., net.weights[i-1], net.rates[i-1], 0., net.v[i])
		BLAS.axpy!(1., net.input[i], net.v[i])
		net.parameters.activation_function(net.v[i], net.rates[i], net.parameters, i)
	end
end

struct SGD
    learningrate::Float64
end
export SGD
struct ADAM
    learningrate::Float64
    β1::Float64
    β2::Float64
    ϵ::Float64
    βs::Array{Float64, 1}
    m::Array{Float64, 2}
    v::Array{Float64, 2}
end
export ADAM
function ADAM(lr, nin, nout; β1 = .9, β2 = .999, ϵ = 1e-8)
    ADAM(lr, β1, β2, ϵ, [β1; β2], zeros(nout, nin), zeros(nout, nin))
end

function update!(net::RateNet, opt::SGD)
  BLAS.ger!(opt.learningrate, net.targets-net.rates[end], net.rates[end-1], net.weights[end])
end
function update!(net::RateNet, opt::ADAM)
    Δw = (net.targets - net.rates[end]) * net.rates[end-1]'
    @. opt.m = opt.β1 * opt.m + (1 - opt.β1) * Δw
    @. opt.v = opt.β2 * opt.v + (1 - opt.β2) * Δw^2
    @. net.weights[end] += opt.learningrate * opt.m/(1 - opt.βs[1]) /
                            √(opt.v / (1 - opt.βs[2]) + opt.ϵ)
    opt.βs .*= [opt.β1; opt.β2]
end
using ProgressMeter

function learnMNIST!(net::RateNet, opt,	images::Array{Float64, 2}, labels, iterations::Int64)
    for i in 2:length(net.input) net.input[i] = ones(length(net.input[i])) .* net.parameters.thresholdmean end
    	@showprogress for i in 1:iterations
        patternindex = rand(1:size(images)[2])
        net.input[1] = net.parameters.in_amplitude .* images[:, patternindex] .+ net.parameters.thresholdmean
        net.targets = [Int(j == labels[patternindex] + 1) for j in 1:10] .* net.parameters.target_amplitude
    	forwardprop!(net)
    	update!(net, opt)
    end
end
export learnMNIST!

function geterrors(net::RateNet, images::Array{Float64, 2}, labels)
  for i in 2:length(net.input) net.input[i] = ones(length(net.input[i])) .* net.parameters.thresholdmean end
	error = 0
	noftest = size(images)[2]
	for i in 1:noftest
		net.input[1] = net.parameters.in_amplitude .* images[:, i] .+ net.parameters.thresholdmean
		forwardprop!(net)
		error += findmax(net.rates[end])[2] != labels[i] + 1
	end
	error/noftest
end
export geterrors


text:
import DoubleFloats: Double64

issquare(A::AbstractMatrix) = isequal(size(A)...)

"""
    row_scaling!(A, D::AbstractVector, tol=maximum(size(A))^2 * eps(real(eltype(A))))

Divide each row of `A` by its ∞-norm. This minimizes the ∞ condition number
for all possible scalings. In order to avoid scaling near zero rows back to
1 we only scale a row if its ∞-norm is larger than `tol`.
"""
function row_scaling!(A::AbstractMatrix{T}, D::AbstractVector, tol=max(size(A)...)^2 * eps(real(T))) where {T}
    D .= 0.0
    d_max = convert(real(T), -Inf)
    @inbounds for j=1:size(A, 2), i=1:size(A, 1)
        D[i] = @fastmath max(D[i], abs(A[i,j]))
        d_max = @fastmath max(d_max, D[i])
    end

    # In the case that the row norms are all less than 1 we have a lower bound
    # the maximal row norm
    # To avoid that we scale a zero matrix, we need a tol for wich we don't scale any more
    bound = clamp(d_max, tol, 1.0)
    D .= max.(D, bound)

    @inbounds for j=1:size(A, 2), i=1:size(A, 1)
        A[i, j] /= D[i]
    end
    A
end

@enum ActiveFactorization begin
    QR_FACTORIZATION
    LU_FACTORIZATION
end

"""
    Jacobian{T, <:Factorization}

Data structure holding a Jacobian matrix `J` with it's factorization.
Additional contains a vector `D` of scalings which row equilibarate `J`.

    Jacobian(J)

Setup a Jacobian with Jacobian matrices of the form `J`.
"""
mutable struct Jacobian{T}
    J::Matrix{T} # Jacobian matrix
    D::Vector{Float64} # Inverse of scaling factors D
    # lu is only there for square systems
    lu::Union{Nothing, LinearAlgebra.LU{T,Matrix{T}}} # LU Factorization of D * J
    qr::LinearAlgebra.QRPivoted{T,Matrix{T}} # Pivoted QR Factorization of D * J
    active_factorization::ActiveFactorization
    b::Vector{T} # Vector for overdetermined systems
    r::Vector{T} # Vector for iterative refinement
    perm::Vector{Int} # work copy of the permutation in the qr
    qr_work::Vector{ComplexF64} # preallocated memory for pivoted qr
    qr_rwork::Vector{Float64} # preallocated memory for pivoted qr
    ormqr_work::Vector{ComplexF64}
    # Numerical Informations
    corank::Int # The numerical corank
    corank_proposal::Int
    cond::Float64 # Estimate of the Jacobian
    # The relative number of digits lost during the solution of the linear systems
    # in Newton's method. See `solve_with_digits_lost!` in utilities/linear_algebra.jl
    # for how this is computed.
    digits_lost::Union{Nothing,Float64}
end

function Jacobian(A::AbstractMatrix, corank::Int=0)
    m, n = size(A)
    lu = m ≠ n ? nothing : LinearAlgebra.lu!(Random.randn!(complex.(float.(A))))

    work = Vector{ComplexF64}(undef, 1)
    ormqr_work = Vector{ComplexF64}(undef, 1)
    rwork = Vector{Float64}(undef, 2n)
    jpvt = zeros(BlasInt, n)
    B = ComplexF64.(A)
    tau = similar(B, min(m, n))
    qr = LinearAlgebra.QRPivoted(geqp3!(B, jpvt, tau, work, rwork)...)
    active_factorization = lu === nothing ? QR_FACTORIZATION : LU_FACTORIZATION
    J = copy(B)
    D = ones(m)
    r = similar(J, size(J, 1))
    b = copy(r)
    perm = zeros(Int, length(qr.p))
    cond = 1.0
    digits_lost = nothing
    Jacobian(J, D, lu, qr, active_factorization, b, r, perm, work, rwork, ormqr_work,
        corank, 0, cond, digits_lost)
end

function update_rank!(Jac::Jacobian)
    Jac.corank = Jac.corank_proposal
    Jac
end

function reset!(jacobian::Jacobian)
    jacobian.cond = 1.0
    jacobian.corank = 0
    jacobian.corank_proposal = 0
    jacobian.digits_lost = nothing
    jacobian
end

"""
    update!(jacobian, J, corank::Int=jacobian.corank)

Update the `jacobian` struct with the new Matrix `J`.
"""
function update!(Jac::Jacobian, J::AbstractMatrix; update_infos::Bool=false)
    copyto!(Jac.J, J)
    updated_jacobian!(Jac; update_infos=update_infos)
    Jac
end

"""
    update_jacobian!(jacobian)

This computes a factorization of the stored Jacobian matrix.
Call this instead of `update!` if `jacobian.J` got updated.
"""
function updated_jacobian!(Jac::Jacobian{ComplexF64}; update_infos::Bool=false) where {T}
    if !update_infos && issquare(Jac.J) && Jac.corank == 0
        Jac.lu !== nothing || return Jac
        copyto!(Jac.lu.factors, Jac.J)
        lu_factorization!(Jac.lu.factors, nothing, Jac.lu.ipiv)
        Jac.active_factorization = LU_FACTORIZATION
    else
        copyto!(Jac.qr.factors, Jac.J)

        m, n = size(Jac.J)
        if m == n
            # only apply row scaling for square matrices since
            # otherwise we change the problem
            row_scaling!(Jac.qr.factors, Jac.D, 1)
        end
        # this computes a pivoted qr factorization, i.e.,
        # qr!(Jac.qr.factors, Val(true)) but without allocating new memory
        geqp3!(Jac.qr.factors, Jac.qr.jpvt, Jac.qr.τ, Jac.qr_work, Jac.qr_rwork)

        ε = max(n,m) * eps() * 10
        # check rank 0
        rnm = min(n,m)
        r₁ = abs(real(Jac.qr.factors[1,1]))
        corank = 0
        if r₁ < ε
            corank = rnm
            Jac.cond = inv(ε)
        else
            for i in 2:rnm
                rᵢ = abs(real(Jac.qr.factors[i,i]))
                if ε * r₁ > rᵢ
                    corank = rnm - i + 1
                    break
                end
            end
            Jac.cond = max(1.0, r₁) / abs(real(Jac.qr.factors[rnm, rnm]))
        end
        Jac.corank_proposal = corank

        Jac.active_factorization = QR_FACTORIZATION
    end
    Jac
end

"""
    solve!([x, ], Jac::Jacobian, b)

Solve `Jac.J * x = b` inplace. Assumes `Jac` contains already the correct factorization.
This stores the result in `x`.
"""
function solve!(x, Jac::Jacobian, b::AbstractVector; update_digits_lost::Bool=false, refinement_step::Bool=false)
    if Jac.active_factorization == LU_FACTORIZATION
        lu = Jac.lu
        if lu !== nothing
            custom_ldiv!(x, lu, b)
        end
    elseif Jac.active_factorization == QR_FACTORIZATION
        if issquare(Jac.J)
            # we applied row scaling for square matrices
            @inbounds for i in eachindex(Jac.b)
                Jac.b[i] = b[i] / Jac.D[i]
            end
            custom_ldiv!(x, Jac.qr, Jac.b, Jac.corank, Jac.perm, Jac.ormqr_work )
        else
            Jac.b .= b
            custom_ldiv!(x, Jac.qr, Jac.b, Jac.corank, Jac.perm, Jac.ormqr_work )
        end
    end
    if issquare(Jac.J) && Jac.corank == 0 && update_digits_lost && !refinement_step
        norm_x = euclidean_norm(x)
        norm_δx = iterative_refinement_step!(x, Jac, b)
        Jac.digits_lost = log₁₀(norm_δx / eps(norm_x))
    elseif !issquare(Jac.J)
        Jac.digits_lost = nothing
    end
    x
end
solve!(A::Jacobian, b) = solve!(b, A::Jacobian, b)

function custom_ldiv!(x, LU::LinearAlgebra.LU, b::AbstractVector)
    if x !== b
        copyto!(x, b)
    end
     _ipiv!(LU, x)
     ldiv_unit_lower!(LU.factors, x)
     ldiv_upper!(LU.factors, x)
     x
end

const LA = LinearAlgebra

# This is an adoption of the Julia implementation of ldiv! for QRPivoted
function custom_ldiv!(x, A::LinearAlgebra.QRPivoted{ComplexF64},
                      b::Vector{ComplexF64}, corank::Int,
                      perm::Vector{Int}, ormqr_work::Vector{ComplexF64})
    mA, nA = size(A.factors)
    nr = min(mA,nA)
    nrhs = length(b)
    if nr == 0
        return x
    end
    if corank == nr
        fill!(x, zero(eltype(x)))
        return x
    end
    rank = nr - corank

    if rank == nA == nr
        # The following is equivalent to LA.lmul!(LA.adjoint(A.Q), b) but
        # we can also pass the preallocated work vector
        ormqr!('L', 'C', A.factors, A.τ, b, ormqr_work)
        ldiv_upper!(A.factors, b; singular_exception=false)
        @inbounds for i in 1:nr
            x[i] = b[i]
            perm[i] = A.p[i]
        end
        Base.invpermute!!(x, perm)
    else
        C, τ = LA.LAPACK.tzrzf!(A.factors[1:rank,:])
        LA.ldiv!(LA.UpperTriangular(C[1:rank,1:rank]),view(LA.lmul!(LA.adjoint(A.Q), b), 1:rank))
        for i in 1:rank
            x[i] = b[i]
        end
        for i in rank+1:nA
            x[i] = zero(eltype(x))
        end
        LA.LAPACK.ormrz!('L', eltype(x)<:Complex ? 'C' : 'T', C, τ, reshape(x, length(x), 1))
        invpermute!(x, A.p)
    end
    return x
end

"""
    iterative_refinement_step!(x, A, b, fac, T, r)

Apply one step of iterative refinement where the residual is computed with precision `T`.
Returns the euclidean norm of the update.
"""
function iterative_refinement_step!(x, Jac::Jacobian, b, ::Type{T}=eltype(x)) where T
    residual!(Jac.r, Jac.J, x, b, T)
    δx = solve!(Jac.r, Jac, Jac.r; refinement_step=true)
    norm_δx = euclidean_norm(δx)
    @inbounds for i in eachindex(x)
        x[i] = convert(T, x[i]) - convert(T, δx[i])
    end

    return norm_δx
end

"""
    residual!(u, A, x, b, ::Type{T}=eltype(u))

Compute the residual `Ax-b` in precision `T` and store in `u`.
"""
function residual!(u::AbstractVector, A, x, b, ::Type{T}=eltype(u)) where {T}
    @boundscheck size(A, 1) == length(b) && size(A,2) == length(x)
    m, n = size(A)
    @inbounds for i in 1:m
        dot = zero(T)
        for j in 1:n
            dot = muladd(convert(T, A[i,j]), convert(T, x[j]), dot)
        end
        u[i] = dot - convert(T, b[i])
    end
    u
end

###########################
# CUSTOM LU Factorization #
###########################

# This is an adoption of LinearAlgebra.generic_lufact!
# with 3 changes:
# 1) For choosing the pivot we use abs2 instead of abs
# 2) Instead of using the robust complex division we
#    just use the naive division. This shouldn't
#    lead to problems since numerical diffuculties only
#    arise for very large or small exponents
# 3) We fold lu! and ldiv! into one routine
#    this has the effect that we do not need to allocate
#    the pivot vector anymore and also avoid the allocations
#    coming from the LU wrapper
@inline function lu_factorization!(A::AbstractMatrix{T},
                           b::Union{AbstractVector, Nothing}=nothing,
                           ipiv::Union{Vector{I}, Nothing}=nothing,
                           ::Val{Pivot} = Val(true)) where {T,I,Pivot}
    m, n = size(A)
    minmn = min(m,n)
    # LU Factorization
    @inbounds begin
        for k = 1:minmn
            # find index max
            kp = k
            if Pivot
                amax = zero(real(T))
                for i = k:m
                    if T <: Complex
                        absi = abs2(A[i,k])
                    else
                        absi = abs(A[i,k])
                    end
                    if absi > amax
                        kp = i
                        amax = absi
                    end
                end
            end
            if !(ipiv isa Nothing)
                ipiv[k] = kp
            end
            if !iszero(A[kp,k])
                if k != kp
                    # Interchange
                    for i = 1:n
                        tmp = A[k,i]
                        A[k,i] = A[kp,i]
                        A[kp,i] = tmp
                    end
                    if !(b isa Nothing)
                        b[k], b[kp] = b[kp], b[k]
                    end
                end
                # Scale first column
                Akkinv = @fastmath inv(A[k,k])
                for i = k+1:m
                    A[i,k] *= Akkinv
                end
            end
            # Update the rest
            for j = k+1:n
                A_kj = A[k,j]
                for i = k+1:m
                    A[i,j] -= A[i,k] * A_kj
                end
            end
        end
    end
    A
end

_ipiv!(A::LinearAlgebra.LU, b::AbstractVector) = apply_ipiv!(b, A.ipiv)

function apply_ipiv!(b::AbstractVector, ipiv)
    @inbounds for i = 1:length(ipiv)
        if i != ipiv[i]
            _swap_rows!(b, i, ipiv[i])
        end
    end
    b
end

function _swap_rows!(B::StridedVector, i::Integer, j::Integer)
    B[i], B[j] = B[j], B[i]
    B
end

@inline function ldiv_upper!(A::AbstractMatrix, b::AbstractVector, x::AbstractVector = b; singular_exception::Bool=true)
    n = size(A, 2)
    for j in n:-1:1
        @inbounds singular_exception && iszero(A[j,j]) && throw(LinearAlgebra.SingularException(j))
        @inbounds xj = x[j] = (@fastmath A[j,j] \ b[j])
        for i in 1:(j-1)
            @inbounds b[i] -= A[i,j] * xj
        end
    end
    b
end
@inline function ldiv_unit_lower!(A::AbstractMatrix, b::AbstractVector, x::AbstractVector = b)
    n = size(A, 2)
    @inbounds for j in 1:n
        xj = x[j] = b[j]
        for i in j+1:n
            b[i] -= A[i,j] * xj
        end
    end
    x
end

import LinearAlgebra: BlasInt

"""
    geqp3!(A::AbstractMatrix{ComplexF64},
        jpvt::AbstractVector{BlasInt},
        tau::AbstractVector{ComplexF64},
        work::Vector{ComplexF64},
        rwork=Vector{Float64}(undef, 2*size(A,2)))

Version of `LAPACK.geqp3!` which also accepts the work buffers.
"""
function geqp3!(A::AbstractMatrix{ComplexF64},
		jpvt::AbstractVector{BlasInt},
		tau::AbstractVector{ComplexF64},
		work::Vector{ComplexF64},
		rwork=Vector{Float64}(undef, 2*size(A,2)))
	m,n = size(A)
	lda = stride(A,2)
	if lda == 0
		return return A, tau, jpvt
	end # Early exit
    jpvt .= BlasInt(0)
	lwork = BlasInt(-1)
	info = Ref{BlasInt}()
	for i = 1:2
		ccall((LinearAlgebra.LAPACK.@blasfunc(zgeqp3_), LinearAlgebra.LAPACK.liblapack), Cvoid,
			  (Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF64}, Ref{BlasInt},
			   Ptr{BlasInt}, Ptr{ComplexF64}, Ptr{ComplexF64}, Ref{BlasInt},
			   Ptr{Float64}, Ptr{BlasInt}),
			  m, n, A, lda,
			  jpvt, tau, work, lwork,
			  rwork, info)
		LinearAlgebra.LAPACK.chklapackerror(info[])
		if i == 1
			lwork = BlasInt(real(work[1]))
			resize!(work, lwork)
		end
	end
	return A, tau, jpvt
end

"""
    ormqr!(side, trans, A, tau, C, work)

Version of `LAPACK.ormqr!` which also accepts the work buffer.
Computes `Q * C` (`trans = N`), `transpose(Q) * C` (`trans = T`), `adjoint(Q) * C`
(`trans = C`) for `side = L` or the equivalent right-sided multiplication
for `side = R` using `Q` from a `QR` factorization of `A` computed using
`geqrf!`. `C` is overwritten.
"""
function ormqr!(side::AbstractChar, trans::AbstractChar, A::Matrix{ComplexF64},
                tau::Vector{ComplexF64}, C::Vector{ComplexF64},
                work::Vector{ComplexF64})
    m,n = (size(C, 1), 1)
    mA  = size(A, 1)
    k   = length(tau)
    if side == 'L' && m != mA
        throw(DimensionMismatch("for a left-sided multiplication, the first dimension of C, $m, must equal the second dimension of A, $mA"))
    end
    if side == 'R' && n != mA
        throw(DimensionMismatch("for a right-sided multiplication, the second dimension of C, $m, must equal the second dimension of A, $mA"))
    end
    if side == 'L' && k > m
        throw(DimensionMismatch("invalid number of reflectors: k = $k should be <= m = $m"))
    end
    if side == 'R' && k > n
        throw(DimensionMismatch("invalid number of reflectors: k = $k should be <= n = $n"))
    end
    lwork = BlasInt(-1)
    info  = Ref{BlasInt}()
    for i = 1:2  # first call returns lwork as work[1]
        ccall((LA.LAPACK.@blasfunc(zunmqr_), LA.LAPACK.liblapack), Cvoid,
              (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
               Ref{BlasInt}, Ptr{ComplexF64}, Ref{BlasInt}, Ptr{ComplexF64},
               Ptr{ComplexF64}, Ref{BlasInt}, Ptr{ComplexF64}, Ref{BlasInt},
               Ptr{BlasInt}),
              side, trans, m, n,
              k, A, max(1,stride(A,2)), tau,
              C, max(1, stride(C,2)), work, lwork,
              info)
        LA.LAPACK.chklapackerror(info[])
        if i == 1
            lwork = BlasInt(real(work[1]))
            resize!(work, lwork)
        end
    end
    C
end

"""
    hnf(A, T=elytpe(A))

    Compute the hermite normal form `H` of `A` by overwriting `A` and the corresponding transformation
    matrix `U` using precision T. This is `A*U == H` and `H` is a lower triangular matrix.

    The implementation follows the algorithm described in [1].

    [1] Kannan, Ravindran, and Achim Bachem. "Polynomial algorithms for computing the Smith and Hermite normal forms of an integer matrix." SIAM Journal on Computing 8.4 (1979): 499-507.
"""
function hnf(A, T=eltype(A))
    H = similar(A, T)
    H .= A
    U = similar(A, T)
    hnf!(H, U)
    H, U
end

# use checked arithmethic
⊡(x,y) = Base.checked_mul(x, y)
⊞(x,y) = Base.checked_add(x, y)

"""
    hnf!(H, U, A)

Inplace version of [hnf](@ref).

    hnf!(A, U)

Inplace version of [hnf](@ref) overwriting `A` with `H`.
"""
hnf!(H, U, A) = hnf!(copyto!(H, A), U)
function hnf!(A, U)
    n = size(A, 1)
    U .= 0
    @inbounds for i in 1:n
        U[i,i] = one(eltype(U))
    end
    @inbounds for i in 1:(n-1)
        ii = i ⊞ 1
        for j in 1:i
            if !iszero(A[j, j]) || !iszero(A[j, ii])
                # 4.1
                r, p, q = gcdx(A[j, j], A[j, ii])
                # 4.2
                d_j = -A[j, ii] ÷ r
                d_ii = A[j, j] ÷ r
                for k = 1:n
                    a_kj, a_kii = A[k,j], A[k,ii]
                    A[k, j] = a_kj ⊡ p ⊞ a_kii ⊡ q
                    A[k, ii] = a_kj ⊡ d_j ⊞ a_kii ⊡ d_ii

                    u_kj, u_kii = U[k,j], U[k,ii]
                    U[k, j] = u_kj ⊡ p ⊞ u_kii ⊡ q
                    U[k, ii] = u_kj ⊡ d_j ⊞ u_kii ⊡ d_ii
                end
            end
            # 4.3
            if j > 1
                reduce_off_diagonal!(A, U, j)
            end
        end
        #5
        reduce_off_diagonal!(A, U, ii)
    end
    # Special case for 1 × 1 matrices to guarantee that the diagonal is positive
    # This comes up in polyhedral homotopy for univariate polynomials
    if n == 1
        U[1,1] = flipsign(U[1,1], A[1,1])
        A[1,1] = flipsign(A[1,1], A[1,1])
    end

    nothing
end

@inline function reduce_off_diagonal!(A, U, k)
    n = size(A, 1)
    @inbounds if A[k,k] < 0
        for i in 1:n
            A[i, k] = -A[i, k]
            U[i, k] = -U[i, k]
        end
    end
    @inbounds for z in 1:(k-1)
        if !iszero(A[z,z])
            r = -ceil(eltype(A), A[k,z] / A[z,z])
            for i in 1:n
                U[i,z] = U[i,z] ⊞ r ⊡ U[i,k]
                A[i,z] = A[i,z] ⊞ r ⊡ A[i,k]
            end
        end
    end
    nothing
end


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module QT

using ExtractMacro
using ..Interface
using ..Common

using ...RRRMC # this is silly but it's required for correct cross-linking in docstrings, apparently

export GraphQT, Qenergy, transverse_mag,
       GraphQuant, Renergies, overlaps

import ..Interface: energy, delta_energy, neighbors, allΔE,
                    update_cache!, delta_energy_residual

import Base: iterate, length

"""
    Qenergy(X::DoubleGraph, C::Config)

When using the Suzuki-Trotter transformation to simulate quantum
systems in a transverse magnetic field with a replicated classical
system, this function should be used to obtain the average value of the
Hamiltonian observable (divided by the number of spins).
"""
Qenergy(::AbstractGraph, C::Config) = error("not implemented")

"""
    transverse_mag(X::DoubleGraph, C::Config, β::Float64)

When using the Suzuki-Trotter transformation to simulate quantum
systems in a transverse magnetic field with a replicated classical
system, this function should be used to obtain the average value of the
transverse magnetization observable.
"""
transverse_mag(X::DoubleGraph, C::Config, β::Float64) = transverse_mag(inner_graph(X), C, β)

const MAXDIGITS = 8

mutable struct GraphQT{fourK} <: DiscrGraph{Float64}
    N::Int
    M::Int
    Nk::Int
    function GraphQT{fourK}(N::Integer, M::Integer) where {fourK}
        M > 2 || throw(ArgumentError("M must be greater than 2, given: $M"))
        isa(fourK, Float64) || throw(ArgumentError("invalid parameter fourK, expected Float64, given: $(typeof(fourK))"))
        round(fourK, digits=MAXDIGITS) ≠ fourK && throw(ArgumentError("up to $MAXDIGITS decimal digits supported in fourK, given: $(fourK)"))
        N % M == 0 || throw(ArgumentError("N must be divisible by M, given: N=$N M=$M"))
        Nk = N ÷ M
        return new{fourK}(N, M, Nk)
    end
end

@doc """
    GraphQT{fourK}(N::Integer, M::Integer) <: DiscrGraph

An auxiliary `DiscrGraph` used to implement the interactions in the
Suzuki-Trotter dimension when simulating quantum spin systems in a
transverse field.

It is only useful when implementing other graph types; see [`GraphQuant`](@ref).
""" -> GraphQT{fourK}(N::Integer, M::Integer)

GraphQT(X::GraphQT{oldK}, newK::Float64) where {oldK} = GraphQT{newK}(X.N, X.M)

function energy0(X::GraphQT, C::Config)
    @assert X.N == C.N
    @extract X : M Nk
    @extract C : s
    n = 0
    for i = 1:Nk
        sj = s[i + (M-1) * Nk]
        for k = 1:M
            sk = s[i + (k-1) * Nk]
            n -= 1 - 2 * (sk ⊻ sj)
            sj = sk
        end
    end
    return n
end

energy(X::GraphQT{fourK}, C::Config) where {fourK} = energy0(X, C) * fourK / 4

function delta_energy(X::GraphQT{fourK}, C::Config, move::Int) where {fourK}
    @assert X.N == C.N
    @assert 1 ≤ move ≤ C.N
    @extract C : s

    k1, k2 = neighbors(X, move)
    sk = s[move]
    s1 = s[k1]
    s2 = s[k2]

    #Δ = ((2sk - 1) * ((2s1 - 1) + (2s2 - 1))) / 2
    #Δ = ((1 - 2 * (sk ⊻ s1)) + (1 - 2 * (sk ⊻ s2))) / 2
    Δ = ((sk ⊻ ~s1) - (sk ⊻ s2))
    #Δ = ((s1==sk) - (s2≠sk))
    #Δ = (s1 == s2) * (2 * (s1 == sk) - 1)

    return Δ * fourK
end

function neighbors(X::GraphQT, i::Int)
    @extract X : N Nk
    return (i - Nk + N * (i ≤ Nk), i + Nk - N * (i + Nk > N))
end
#neighbors(X::GraphQT, i::Int) = (mod1(i-X.Nk, X.N), mod1(i+X.Nk, X.N))

@generated allΔE(::Type{GraphQT{fourK}}) where {fourK} = (0.0, fourK)

function transverse_mag(X::GraphQT{fourK}, C::Config, β::Float64) where {fourK}
    @extract X : N M
    p = -energy0(X, C) / N

    # x = log(coth(β * Γ / M))
    x = β * fourK / 2

    return cosh(x) - p * sinh(x)
end

# Add Transverse field to (almost) any AbstractGraph

mutable struct GraphQuant{fourK,G<:AbstractGraph} <: DoubleGraph{DiscrGraph{Float64},Float64}
    N::Int
    M::Int
    Nk::Int
    X0::GraphQT{fourK}
    X1::Vector{G}
    C1::Vector{Config}
    β::Float64
    Γ::Float64
    function GraphQuant{fourK,G}(N::Integer, M::Integer, β::Float64, Γ::Float64, g0::G, Gconstr, args...) where {fourK,G}
        X0 = GraphQT{fourK}(N, M)
        Nk = X0.Nk
        #J = gen_J(Nk)
        X1 = Array{G}(undef, M)
        X1[1] = g0
        for k = 2:M
            X1[k] = Gconstr(args...)
        end
        C1 = [Config(Nk, init=false) for k = 1:M]
        return new{fourK,G}(N, M, Nk, X0, X1, C1, β, Γ)
    end
end

"""
    GraphQuant(N::Integer, M::Integer, Γ::Float64, β::Float64, Gconstr, args...) <: DoubleGraph

A `DoubleGraph` which implements a quantum Ising spin model in a transverse magnetic field,
using the Suzuki-Trotter transformation. This allows to model the quantum transverse field
case for any classical Ising model previously defined. This kind of graph can be simulated
efficiently with [`rrrMC`](@ref).

`N` is the number of spins, `M` the number of Suzuki-Trotter replicas, `Γ` the transverse
field, `β` the inverse temperature. `GConstr` is the (classical) graph constructor, and `args` the
arguments to the contructor.

See also [`Qenergy`](@ref) and [`transverse_mag`](@ref).
"""
function GraphQuant(Nk::Integer, M::Integer, Γ::Float64, β::Float64, Gconstr, args...)
    @assert Γ ≥ 0
    fourK = round(2/β * log(coth(β * Γ / M)), digits=MAXDIGITS)
    # H0 = Nk * M / 2β * log(sinh(2β * Γ / M) / 2) # useless!!!
    g0 = Gconstr(args...)
    G = typeof(g0)
    return GraphQuant{fourK,G}(Nk * M, M, β, Γ, g0, Gconstr, args...)
end

function update_cache!(X::GraphQuant, C::Config, move::Int)
    @extract X : X1 Nk C1
    #@extract C : s
    k = (move - 1) ÷ Nk + 1
    i = mod1(move, Nk)

    spinflip!(X1[k], C1[k], i)

    #s1 = C1[k].s
    #copyto!(s1, 1, s, (k-1)*Nk + 1, Nk)
    #@assert C1[k].s == s[((k-1)*Nk + 1):k*Nk]
end

function energy(X::GraphQuant, C::Config)
    @assert X.N == C.N
    @extract X : M Nk X0 X1 C1
    @extract C : s

    E = energy(X0, C)

    for k = 1:M
        s1 = C1[k].s
        copyto!(s1, 1, s, (k-1)*Nk + 1, Nk)
        E += energy(X1[k], C1[k]) / M
    end

    return E
end

function Renergies(X::GraphQuant)
    @extract X : M X1 C1

    Es = zeros(M)

    for k = 1:M
        Es[k] = energy(X1[k], C1[k])
    end

    return Es
end

function overlaps(X::GraphQuant)
    @extract X : M Nk C1

    aux = BitArray(undef, Nk)
    ovs = zeros(M ÷ 2)

    for k1 = 1:(M-1)
        @extract C1[k1] : s1=s
        for k2 = (k1+1):M
            @extract C1[k2] : s2=s
            map!(⊻, aux, s1, s2)
            δ = min(k2 - k1, M + k1 - k2)
            ovs[δ] += Nk - 2 * sum(aux)
        end
    end

    # for i = 1:Nk
    #     for k1 = 1:(M-1)
    #         s1 = s[i + (k1-1) * Nk]
    #         for k2 = (k1+1):M
    #             δ = min(k2 - k1, M + k1 - k2)
    #             s2 = s[i + (k2-1) * Nk]
    #             ovs[δ] += 1 - 2 * (sk1 ⊻ sk2)
    #         end
    #     end
    # end

    # note: for each node in the loop, there are
    # 2 edges linking to other nodes at a given distance δ,
    # except if M is even and δ=M÷2, in which case there is
    # only 1 edge. So the multiplicity is M for all overlaps
    # except for the last entry for even M, then it's M/2.
    for δ = 1:((M-1)÷2)
        ovs[δ] /= M * Nk
    end
    iseven(M) && (ovs[M÷2] /= M * Nk / 2)

    return ovs
end

function Qenergy(X::GraphQuant, C::Config)
    @assert X.N == C.N
    @extract X : M N X0 X1 C1 β Γ
    #@extract C : s

    E = -Γ * transverse_mag(X0, C, β)

    for k = 1:M
        #s1 = C1[k].s
        #copyto!(s1, 1, s, (k-1)*Nk + 1, Nk)
        #@assert C1[k].s == s[((k-1)*Nk + 1):k*Nk]
        E += energy(X1[k], C1[k]) / N
    end

    return E
end

function delta_energy_residual(X::GraphQuant, C::Config, move::Int)
    @extract X : M Nk X1 C1
    #@extract C : s

    k = (move - 1) ÷ Nk + 1
    #s1 = C1[k].s
    #copyto!(s1, 1, s, (k-1)*Nk + 1, Nk)
    #@assert C1[k].s == s[((k-1)*Nk + 1):k*Nk]

    i = mod1(move, Nk)
    return delta_energy(X1[k], C1[k], i) / M
end

function delta_energy(X::GraphQuant, C::Config, move::Int)
    return delta_energy(X.X0, C, move) +
           delta_energy_residual(X, C, move)
end

mutable struct QNeighbIter{T}
    tn1::Int
    tn2::Int
    off::Int
    r::T
    QNeighbIter(tn1, tn2, off, r::T) where {T} = new{T}(tn1, tn2, off, r)
end

iterate(qn::QNeighbIter) = (qn.tn1, (nothing, Val(1)))
iterate(qn::QNeighbIter, st::Tuple{Nothing,Val{1}}) = (qn.tn2, (nothing, Val(2)))
function iterate(qn::QNeighbIter, st::Tuple{Nothing,Val{2}})
    rr = iterate(qn.r)
    rr ≡ nothing && return nothing
    v, i = rr
    return (v + qn.off, (i, nothing))
end
function iterate(qn::QNeighbIter, st::Tuple{<:Any,Nothing})
    rr = iterate(qn.r, st[1])
    rr ≡ nothing && return nothing
    v, i = rr
    return (v + qn.off, (i, nothing))
end

length(qn::QNeighbIter) = 2 + length(qn.r)

function neighbors(X::GraphQuant, i::Int)
    @extract X : N Nk X0 X1
    j1, j2 = neighbors(X0, i)

    k = (i - 1) ÷ Nk + 1
    j = mod1(i, Nk)

    return QNeighbIter(j1, j2, (k-1)*Nk, neighbors(X1[k], j))
end

end # module


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMI
using Flux
using DifferentialEquations: Tsit5

FMUPath = joinpath(dirname(@__FILE__), "..", "model", "SpringPendulumExtForce1D.fmu")

t_start = 0.0
t_step = 0.01
t_stop = 5.0
tData = t_start:t_step:t_stop

# generate traing data
myFMU = fmiLoad(FMUPath)
fmiInstantiate!(myFMU; loggingOn=false)
fmiSetupExperiment(myFMU, t_start, t_stop)
fmiSetReal(myFMU, "mass_s0", 1.3)   # increase amplitude, invert phase
fmiEnterInitializationMode(myFMU)
fmiExitInitializationMode(myFMU)
realSimData = fmiSimulateCS(myFMU, t_start, t_stop; recordValues=["mass.a"], setup=false, saveat=tData)

# reset FMU for use as NeuralFMU
fmiReset(myFMU)
fmiSetupExperiment(myFMU, t_start, t_stop)
fmiEnterInitializationMode(myFMU)
fmiExitInitializationMode(myFMU)

# setup traing data
extF = zeros(length(tData)) # no external force
accData = fmi2SimulationResultGetValues(realSimData, "mass.a")

# loss function for training
function losssum()
    solution = problem(t_step; inputs=extF)

    accNet = collect(data[2] for data in solution)

    Flux.Losses.mse(accNet, accData)
end

# callback function for training
global iterCB = 0
global lastLoss = 0.0
function callb()
    global iterCB += 1
    global lastLoss

    if iterCB == 1
        lastLoss = losssum()
    end

    if iterCB % 50 == 0
        loss = losssum()
        @info "Loss: $loss"
        @test loss < lastLoss
        lastLoss = loss
    end
end

# NeuralFMU setup
numInputs = length(myFMU.modelDescription.inputValueReferences)
numOutputs = length(myFMU.modelDescription.outputValueReferences)

net = Chain(inputs -> fmiInputDoStepCSOutput(myFMU, t_step, inputs),
            Dense(numOutputs, 16, tanh),
            Dense(16, 16, tanh),
            Dense(16, numOutputs))

problem = CS_NeuralFMU(myFMU, net, (t_start, t_stop); saveat=tData)
@test problem != nothing

solutionBefore = problem(t_step; inputs=extF)
ts = collect(data[1] for data in solutionBefore)
@test length(ts) == length(tData)
@test abs(ts[1] - (t_start+t_step)) < 1e-10
@test abs(ts[end] - (t_stop+t_step)) < 1e-10

# train it ...
p_net = Flux.params(problem)

optim = ADAM()
Flux.train!(losssum, p_net, Iterators.repeated((), 300), optim; cb=callb)

# check results
solutionAfter = problem(t_step; inputs=extF)
ts = collect(data[1] for data in solutionAfter)
@test length(ts) == length(tData)
@test abs(ts[1] - (t_start+t_step)) < 1e-10
@test abs(ts[end] - (t_stop+t_step)) < 1e-10

fmiUnload(myFMU)


text:
## BOUNDEDNESS TRAIT

# For random search and perhaps elsewhere, we need a variation on the
# built-in boundedness notions:
abstract type PositiveUnbounded <: Unbounded end
abstract type Other <: Unbounded end

boundedness(::NumericRange{<:Any,<:Bounded}) = Bounded
boundedness(::NumericRange{<:Any,<:LeftUnbounded}) = Other
boundedness(::NumericRange{<:Any,<:DoublyUnbounded}) = Other
function boundedness(r::NumericRange{<:Any,<:RightUnbounded})
    if r.lower >= 0
        return PositiveUnbounded
    end
    return Other
end

"""
    MLJTuning.grid([rng, ] prototype, fields, iterators)

Given a vector of (nested) field names `fields`, and a vector of
`iterators` (each element being a vector of values for the
corresponding field) return a vector of models generated by cloning
and mutating the fields of `prototype` through all possible
combinations of values given by the iterator values.

Specification of an `AbstractRNG` object `rng` implies shuffling of
the results. Otherwise models are ordered, with the first
hyperparameter referenced cycling fastest.

"""
grid(rng::AbstractRNG, prototype::Model, fields, iterators) =
    shuffle(rng, grid(prototype, fields, iterators))

function grid(prototype::Model, fields, iterators)

    A = MLJBase.unwind(iterators...)

    N = size(A, 1)
    map(1:N) do i
        clone = deepcopy(prototype)
        for k in eachindex(fields)
            recursive_setproperty!(clone, fields[k], A[i,k])
        end
        clone
    end
end

## PRE-PROCESSING OF USER-SPECIFIED CARTESIAN RANGE OBJECTS

"""
    process_grid_range(user_specified_range, resolution, verbosity)

Convert a user-specified range (see [`Grid`](@ref)) into a tuple  of
tuples `(ranges, resolutions)`.

For example, if `r1`, `r2` are `NumericRange`s and `s` is a
NominalRange` with 5 values, then we have:

    julia> MLJTuning.process_grid_range([(r1, 3), r2, s], 42, 1) ==
                            ((r1, r2, s), (3, 42, 5))
    true

If `verbosity` > 0, then a warning is issued if a `Nominal` range is
paired with a resolution.

"""
process_grid_range(user_specified_range, args...) =
    throw(ArgumentError("Unsupported range. "))

process_grid_range(usr::Union{ParamRange,Tuple{ParamRange,Int}}, args...) =
    process_grid_range([usr, ], args...)

function process_grid_range(user_specified_range::AbstractVector,
                    resolution, verbosity)
    # r unpaired:
    stand(r) = throw(ArgumentError("Unsupported range. "))
    stand(r::NumericRange) = (r, resolution)
    stand(r::NominalRange) = (r, length(r.values))

    # (r, res):
    stand(t::Tuple{NumericRange,Integer}) = t
    function stand(t::Tuple{NominalRange,Integer})
        verbosity < 0 ||
            @warn  "Ignoring a resolution specified for a `NominalRange`. "
        return (first(t), length(first(t).values))
    end

    ret = zip(stand.(user_specified_range)...) |> collect
    return first(ret), last(ret)
end

"""
    process_random_range(user_specified_range,
                         bounded,
                         positive_unbounded,
                         other)

Utility to convert a user-specified range (see [`RandomSearch`](@ref))
into an n-tuple of `(field, sampler)` pairs.

"""
process_random_range(user_specified_range, args...) =
    throw(ArgumentError("Unsupported range #1. "))

const DIST = Distributions.Distribution

process_random_range(user_specified_range::Union{ParamRange, Tuple{Any,Any}},
                     args...) =
    process_random_range([user_specified_range, ], args...)

function process_random_range(user_specified_range::AbstractVector,
                              bounded,
                              positive_unbounded,
                              other)
    # r not paired:
    stand(r) = throw(ArgumentError("Unsupported range #2. "))
    stand(r::NumericRange) = stand(r, boundedness(r))
    stand(r::NumericRange, ::Type{<:Bounded}) = (r.field, sampler(r, bounded))
    stand(r::NumericRange, ::Type{<:Other}) = (r.field, sampler(r, other))
    stand(r::NumericRange, ::Type{<:PositiveUnbounded}) =
        (r.field, sampler(r, positive_unbounded))
    stand(r::NominalRange) = (n = length(r.values);
                              (r.field, sampler(r, fill(1/n, n))))
    # (r, d):
    stand(t::Tuple{ParamRange,Any}) = stand(t...)
    stand(r, d) = throw(ArgumentError("Unsupported range #3. "))
    stand(r::NominalRange, d::AbstractVector{Float64}) =  _stand(r, d)
    stand(r::NumericRange, d:: Union{DIST, Type{<:DIST}}) = _stand(r, d)
    _stand(r, d) = (r.field, sampler(r, d))

    # (field, s):
    stand(t::Tuple{Union{Symbol,Expr},Any}) = t

    return  Tuple(stand.(user_specified_range))

    # ret = zip(stand.(user_specified_range)...) |> collect
    # return first(ret), last(ret)
end