text:
module FrankenTuples

export FrankenTuple, ftuple, @ftuple, ftcall

"""
    FrankenTuple{T<:Tuple, names, NT<:Tuple}

A `FrankenTuple` contains a `Tuple` of type `T` and a `NamedTuple` with names `names`
and types `NT`.
It acts like a cross between the two, like a partially-named tuple.

The named portion of a `FrankenTuple` can be accessed using `NamedTuple`, and the unnamed
portion can be accessed with `Tuple`.

# Examples
```jldoctest
julia> ft = FrankenTuple((1, 2), (a=1, b=2))
FrankenTuple((1, 2), (a = 1, b = 2))

julia> Tuple(ft)
(1, 2)

julia> NamedTuple(ft)
(a = 1, b = 2)
```
"""
struct FrankenTuple{T<:Tuple,names,NT<:Tuple}
    t::T
    nt::NamedTuple{names,NT}

    FrankenTuple{T,names,NT}(t, nt) where {T<:Tuple,names,NT<:Tuple} = new{T,names,NT}(t, nt)
    FrankenTuple{T}(t::Tuple) where {T<:Tuple} = new{T,(),Tuple{}}(t, NamedTuple())
    FrankenTuple{T,names,NT}(nt::NamedTuple) where {T<:Tuple,names,NT<:Tuple} =
        new{T,names,NT}((), nt)
    FrankenTuple{T,names,NT}(ft::FrankenTuple) where {T<:Tuple,names,NT<:Tuple} =
        new{T,names,NT}(convert(T, getfield(ft, :t)),
                        convert(NamedTuple{names,NT}, getfield(ft, :nt)))
end

FrankenTuple(t::T, nt::NamedTuple{names,NT}) where {T<:Tuple,names,NT<:Tuple} =
    FrankenTuple{T,names,NT}(t, nt)
FrankenTuple() = FrankenTuple{Tuple{},(),Tuple{}}((), NamedTuple())

FrankenTuple(t::Tuple) = FrankenTuple(t, NamedTuple())
FrankenTuple(nt::NamedTuple) = FrankenTuple((), nt)

FrankenTuple(ft::FrankenTuple) = ft

"""
    Tuple(ft::FrankenTuple)

Access the `Tuple` part of a `FrankenTuple`, i.e. the "plain," unnamed portion.
"""
Base.Tuple(ft::FrankenTuple) = getfield(ft, :t)

"""
    NamedTuple(ft::FrankenTuple)

Access the `NamedTuple` part of a `FrankenTuple`, i.e. the named portion.
"""
Base.NamedTuple(ft::FrankenTuple) = getfield(ft, :nt)

function Base.show(io::IO, ft::FrankenTuple)
    print(io, "FrankenTuple(")
    show(io, Tuple(ft))
    print(io, ", ")
    show(io, NamedTuple(ft))
    print(io, ')')
    nothing
end
Base.show(io::IO, ft::FrankenTuple{Tuple{},(),Tuple{}}) =
    print(io, "FrankenTuple()")

Base.convert(::Type{FrankenTuple{T,names,NT}}, t::Tuple) where {T<:Tuple,names,NT<:Tuple} =
    FrankenTuple{T,names,NT}(convert(T, t), NamedTuple())
Base.convert(::Type{FrankenTuple{T,names,NT}}, nt::NamedTuple) where {T<:Tuple,names,NT<:Tuple} =
    FrankenTuple{T,names,NT}((), convert(NamedTuple{names,NT}, nt))

"""
    isempty(ft::FrankenTuple)

Determine whether the given `FrankenTuple` is empty, i.e. has at least 1 element.
"""
Base.isempty(ft::FrankenTuple) = false
Base.isempty(ft::FrankenTuple{Tuple{},(),Tuple{}}) = true

"""
    length(ft::FrankenTuple)

Compute the number of elements in `ft`.
"""
Base.length(ft::FrankenTuple) = length(Tuple(ft)) + length(NamedTuple(ft))
Base.length(ft::FrankenTuple{Tuple{},(),Tuple{}}) = 0
Base.length(ft::FrankenTuple{<:Tuple,(),Tuple{}}) = length(Tuple(ft))
Base.length(ft::FrankenTuple{Tuple{},names}) where {names} = length(names)

"""
    getindex(ft::FrankenTuple, i)

Retrieve the value of `ft` at the given index `i`. When `i::Integer`, this gets the value
at index `i` in iteration order.
When `i::Symbol`, this gets the value from the named section with name `i`.
(`getproperty` can also be used for the `Symbol` case.)

# Examples
```jldoctest
julia> ftuple(1, 2; a=3, b=4)[3]
3

julia> ftuple(1, 2; a=3, b=4)[:a]
3
```
"""
function Base.getindex(ft::FrankenTuple, i::Integer)
    t = Tuple(ft)
    n = length(t)
    if i <= n
        getfield(t, i)
    else
        getfield(NamedTuple(ft), i - n)
    end
end
Base.getindex(ft::FrankenTuple, x::Symbol) = getfield(NamedTuple(ft), x)

Base.getproperty(ft::FrankenTuple, x::Symbol) = getfield(NamedTuple(ft), x)

"""
    firstindex(ft::FrankenTuple)

Retrieve the first index of `ft`, which is always 1.
"""
Base.firstindex(ft::FrankenTuple) = 1

"""
    lastindex(ft::FrankenTuple)

Retrieve the last index of `ft`, which is equivalent to its `length`.
"""
Base.lastindex(ft::FrankenTuple) = length(ft)

"""
    first(ft::FrankenTuple)

Get the first value in `ft` in iteration order.
`ft` must be non-empty.
"""
Base.first(ft::FrankenTuple) = @inbounds ft[1]
Base.first(ft::FrankenTuple{Tuple{},(),Tuple{}}) =
    throw(ArgumentError("FrankenTuple must be non-empty"))

"""
    Base.tail(ft::FrankenTuple)

Return the tail portion of `ft`: a new `FrankenTuple` with the first element of `ft`
removed.
`ft` must be non-empty.

# Examples
```jldoctest
julia> Base.tail(ftuple(a=4, b=5))
FrankenTuple((), (b = 5,))
```
"""
Base.tail(ft::FrankenTuple) = _tail(Tuple(ft), NamedTuple(ft))
# TODO: Should be able to get rid of the helper after VERSION >= v"1.1.0-DEV.553"
_tail(t::Tuple{}, nt::NamedTuple{(),Tuple{}}) =
    throw(ArgumentError("FrankenTuple must be non-empty"))
_tail(t::Tuple{}, nt::NamedTuple{names,<:Tuple}) where {names} =
    FrankenTuple(t, NamedTuple{Base.tail(names)}(nt))
_tail(t::Tuple, nt::NamedTuple) = FrankenTuple(Base.tail(t), nt)

"""
    iterate(ft::FrankenTuple[, state])

Iterate over `ft`.
This yields the values of the unnamed section first, then the values of the named section.

# Examples
```jldoctest
julia> ft = @ftuple (1, a=3, 2, b=4)
FrankenTuple((1, 2), (a = 3, b = 4))

julia> collect(ft)
4-element Array{Int64,1}:
 1
 2
 3
 4
```
"""
Base.iterate(ft::FrankenTuple, state...) =
    iterate(Iterators.flatten((Tuple(ft), NamedTuple(ft))), state...)

"""
    keys(ft::FrankenTuple)

Get the keys of the given `FrankenTuple`, i.e. the set of valid indices into `ft`.
The unnamed section of `ft` has 1-based integer keys and the named section is keyed by
name, given as `Symbol`s.

# Examples
```jldoctest
julia> keys(ftuple(1, 2; a=3, b=4))
(1, 2, :a, :b)
```
"""
Base.keys(ft::FrankenTuple) = (keys(Tuple(ft))..., keys(NamedTuple(ft))...)
Base.keys(ft::FrankenTuple{Tuple{},(),Tuple{}}) = ()
Base.keys(ft::FrankenTuple{Tuple{},N,<:Tuple}) where {N} = N
Base.keys(ft::FrankenTuple{<:Tuple,(),Tuple{}}) = keys(Tuple(ft))

"""
    values(ft::FrankenTuple)

Get the values of the given `FrankenTuple` in iteration order.
The values for the unnamed section appear before that of the named section.

# Examples
```jldoctest
julia> values(ftuple(1, 2; a=3, b=4))
(1, 2, 3, 4)
```
"""
Base.values(ft::FrankenTuple) = (Tuple(ft)..., NamedTuple(ft)...)
Base.values(ft::FrankenTuple{Tuple{},(),Tuple{}}) = ()
Base.values(ft::FrankenTuple{Tuple{},names,<:Tuple}) where {names} = values(NamedTuple(ft))
Base.values(ft::FrankenTuple{<:Tuple,(),Tuple{}}) = Tuple(ft)

"""
    pairs(ft::FrankenTuple)

Construct a `Pairs` iterator that associates the `keys` of `ft` with its `values`.

# Examples
```jldoctest
julia> collect(pairs(ftuple(1, 2; a=3, b=4)))
4-element Array{Pair{Any,Int64},1}:
  1 => 1
  2 => 2
 :a => 3
 :b => 4
```
"""
Base.pairs(ft::FrankenTuple) = Iterators.Pairs(ft, keys(ft))

"""
    eltype(ft::FrankenTuple)

Determine the element type of `ft`.
This is the immediate supertype of the elements in `ft` if they are not homogeneously typed.

# Examples
```jldoctest
julia> eltype(ftuple(1, 2; a=3, b=4))
Int64

julia> eltype(ftuple(0x0, 1))
Integer

julia> eltype(ftuple(a=2.0, b=0x1))
Real

julia> eltype(ftuple())
Union{}
```
"""
Base.eltype(::Type{FrankenTuple{T,names,V}}) where {T<:Tuple,names,V<:Tuple} =
    Base.promote_typejoin(eltype(T), eltype(V))

"""
    empty(ft::FrankenTuple)

Construct an empty `FrankenTuple`.
"""
Base.empty(@nospecialize ft::FrankenTuple) = FrankenTuple()

if VERSION < v"1.2.0-DEV.217"
    function Base.hasmethod(@nospecialize(f),
                            ::Type{FrankenTuple{T,names,NT}};
                            world=typemax(UInt)) where {T,names,NT}
        hasmethod(f, T; world=world) || return false
        m = which(f, T)
        kws = Base.kwarg_decl(m, Core.kwftype(typeof(f)))
        for kw in kws
            endswith(String(kw), "...") && return true
        end
        issubset(names, kws)
    end

else
    function Base.hasmethod(@nospecialize(f),
                            ::Type{FrankenTuple{T,names,NT}};
                            world=typemax(UInt)) where {T,names,NT}
        hasmethod(f, T, names; world=world)
    end
end

function Base.hasmethod(@nospecialize(f),
                        ::Type{FrankenTuple{T,(),Tuple{}}};
                        world=typemax(UInt)) where {T}
    hasmethod(f, T; world=world)
end

function Base.hasmethod(@nospecialize(f),
                        ::Type{FrankenTuple{T,names}};
                        world=typemax(UInt)) where {T,names}
    NT = Tuple{Iterators.repeated(Any, length(names))...}
    hasmethod(f, FrankenTuple{T,names,NT}; world=world)
end

"""
    hasmethod(f, ft::Type{<:FrankenTuple})

Determine whether the function `f` has a method with positional argument types matching
those in the unnamed portion of `ft` and with keyword arguments named in accordance with
those in the named portion of `ft`.

Note that the types in the named portion of `ft` do not factor into determining the
existence of a matching method because keyword arguments to not participate in dispatch.
Similarly, calling `hasmethod` with a `FrankenTuple` with an empty named portion will
still return `true` if the positional arguments match, even if `f` only has methods that
accept keyword arguments.
This ensures agreement with the behavior of `hasmethod` on `Tuple`s.

More generally, the names in the `FrankenTuple` must be a subset of the keyword argument
names in the matching method, _except_ when the method accepts a variable number of
keyword arguments (e.g. `kwargs...`).
In that case, the names in the method must be a subset of the `FrankenTuple`'s names.

# Examples
```jldoctest
julia> f(x::Int; y=3, z=4) = x + y + z;

julia> hasmethod(f, FrankenTuple{Tuple{Int},(:y,)})
true

julia> hasmethod(f, FrankenTuple{Tuple{Int},(:a,)}) # no keyword `a`
false

julia> g(; a, b, kwargs...) = +(a, b, kwargs...);

julia> hasmethod(g, FrankenTuple{Tuple{},(:a,:b,:c,:d)}) # g accepts arbitrarily many kwargs
true
```
"""
Base.hasmethod(::Any, ::Type{<:FrankenTuple})

"""
    ftuple(args...; kwargs...)

Construct a [`FrankenTuple`](@ref) from the given positional and keyword arguments.

# Examples
```jldoctest
julia> ftuple(1, 2)
FrankenTuple((1, 2), NamedTuple())

julia> ftuple(1, 2, a=3, b=4)
FrankenTuple((1, 2), (a = 3, b = 4))
```
"""
function ftuple(args...; kwargs...)
    @static if VERSION < v"1.7.0-DEV.1017"
        nt = kwargs.data
    else
        # NOTE: We don't use this unconditionally because the `NamedTuple` constructor
        # method that accepts arbitrary key-value iterators requires 1.6.0-DEV.877
        nt = NamedTuple(kwargs)
    end
    FrankenTuple(args, nt)
end

"""
    @ftuple (x...; y...)
    @ftuple (a, x=t, b, y=u)

Construct a [`FrankenTuple`](@ref) from the given tuple expression, which can contain
both positional and named elements. The tuple can be "sectioned" in the same manner as
a function signature, with positional elements separated from the named elements by a
semicolon, or positional and named elements can be intermixed, occurring in any order.

# Examples
```jldoctest
julia> @ftuple (1, 2; a=3, b=4)
FrankenTuple((1, 2), (a = 3, b = 4))

julia> @ftuple (1, a=3, 2, b=4)
FrankenTuple((1, 2), (a = 3, b = 4))
```
"""
macro ftuple(ex::Expr)
    ex.head === :tuple || throw(ArgumentError("@ftuple: expected tuple expression"))
    # ()
    if isempty(ex.args)
        t = Expr(:tuple)
        nt = Expr(:call, :NamedTuple)
    # (a, b; x=t, y=u)
    elseif ex.args[1] isa Expr && ex.args[1].head === :parameters
        t = Expr(:tuple)
        length(ex.args) > 1 && append!(t.args, ex.args[2:end])
        nt = Expr(:tuple)
        for kw in ex.args[1].args
            @assert kw isa Expr && kw.head === :kw
            push!(nt.args, Expr(:(=), esc(kw.args[1]), kw.args[2]))
        end
    # (a, x=t, b, y=u)
    else
        t = Expr(:tuple)
        nt = Expr(:tuple)
        for arg in ex.args
            if arg isa Expr && arg.head === :(=)
                push!(nt.args, esc(arg))
            else
                push!(t.args, arg)
            end
        end
        if nt == Expr(:tuple)
            # No named elements found
            nt = Expr(:call, :NamedTuple)
        end
    end
    :(FrankenTuple($t, $nt))
end

# XXX: I'm not convinced that I like this or think it's useful in any way, but it's cute
"""
    ftcall(f::Function, ft::FrankenTuple)

Call the function `f` using the unnamed portion of `ft` as its positional arguments and
the named portion of `ft` as its keyword arguments.

# Examples
```jldoctest
julia> ftcall(mapreduce, ftuple(abs2, -, 1:4; init=0))
-30
```
"""
ftcall(f, ft::FrankenTuple) = f(Tuple(ft)...; NamedTuple(ft)...)
ftcall(f, ft::FrankenTuple{Tuple{},(),Tuple{}}) = f()

end # module


text:
module JWebImageDemo

using Images
using Random

export permute_image

center_begin_end(actual::Int, required::Int) =
    round(Int, (actual - required) / 2) .+ [1, required]

function resize(img, maxSize::Int)
    sz = round.(Int, size(img) .* (maxSize / max(size(img)...)))
    σ = map((o,n) -> 0.05 * o / n, size(img), sz)
    kern = KernelFactors.gaussian(σ)   # from ImageFiltering
    return imresize(imfilter(img, kern, NA()), sz)
end

zone_to_coords(ind::Vector{Int}, step::Int) =
    ind .* step |>
    list -> map(x -> range(x + 1, stop = x + step), list)

function permute_image(img::Matrix, num::Int)::Matrix
    l = min(size(img)...)
    x = center_begin_end(size(img, 2), l)
    y = center_begin_end(size(img, 1), l)

    data_c = img[first(y):last(y), first(x):last(x)]
    img_round = resize(data_c, round(Int, ceil(l / num) * num))

    img_out = similar(img_round)
    step = round(Int, size(img_round, 1) / num)

    indices = [[i, j] for i = 0:(num-1) for j = 0:(num-1)]
    for (ind, ind_base) in zip(Random.shuffle(indices), indices)
        r_x, r_y = zone_to_coords(ind, step)
        rb_x, rb_y = zone_to_coords(ind_base, step)
        img_out[r_x, r_y] = img_round[rb_x, rb_y]
    end
    return img_out
end

end # module

ext: jl


text:
using Test, ExtraTrees, SharedArrays, DelimitedFiles, Statistics

Xconst = ones(Float32, 1000, 10)
Yconst = rand(size(Xconst, 1)) .> 0.5

model = Model(Xconst, Yconst)

@test model.trees[1].is_leaf == true
@test model.trees[1].probability .== Float32(mean(Yconst))
@test ExtraTrees.tree_depth(model.trees[1]) == 1
@test ExtraTrees.tree_nodes(model.trees[1]) == 1
@test ExtraTrees.tree_leaves(model.trees[1]) == 1

################

Xrand = rand(Float32, 1000, 3)
Yrand = rand(size(Xrand, 1)) .> 0.5

model = Model(Xrand, Yrand, max_depth=2)

@test model.trees[1].is_leaf == false
@test ExtraTrees.tree_depth(model.trees[1]) == 2
@test ExtraTrees.tree_leaves(model.trees[1]) == 2
@test ExtraTrees.tree_nodes(model.trees[1]) == 1+2

################

model = Model(Xrand, Yrand, min_samples_leaf=length(Yrand))

@test model.trees[1].is_leaf == true
@test model.trees[1].probability .== Float32(mean(Yrand))
@test ExtraTrees.tree_depth(model.trees[1]) == 1

################

model = Model(Xrand, Yrand, min_samples_split=length(Yrand))

@test model.trees[1].is_leaf == false
@test ExtraTrees.tree_depth(model.trees[1]) == 2

################

Xoptim = rand(Float32, 1000, 10)
Yoptim = rand(size(Xoptim, 1)) .> 0.5
Xoptim = hcat(Xoptim, Yoptim)

model = Model(Xoptim, Yoptim, n_subfeat=size(Xoptim, 2))

@test model.trees[1].is_leaf == false
@test ExtraTrees.tree_depth(model.trees[1]) == 2
@test mean(Yoptim .== model(Xoptim)) == 1.0

################

D = readdlm("./digits-01.csv", ',', Float32)

Ytrain = D[1:300, 1] .== 1.0
Xtrain = D[1:300, 2:end]

Ytest = D[301:end, 1] .== 1.0
Xtest = D[301:end, 2:end]

tree = Model(Xtrain, Ytrain)
forest = Model(Xtrain, Ytrain, n_trees=100)
flatten_forest = FlattenModel(Xtrain, Ytrain, n_trees=100)

acc_tree = mean(Ytest .== (tree(Xtest) .> 0.5))
acc_forest = mean(Ytest .== (forest(Xtest) .> 0.5))
acc_flatten_forest = mean(Ytest .== (flatten_forest(Xtest) .> 0.5))

@test acc_tree > 0.5
@test acc_forest > 0.98
@test acc_flatten_forest > 0.98
@test acc_forest >= acc_tree

#################

flatten_forest = FlattenModel(forest)

@test flatten_forest(Xtest[1, :]) == forest(Xtest[1, :])
@test all(flatten_forest(Xtest) .== forest(Xtest))
@test all(flatten_forest(SharedArray(Xtest)) .== forest(Xtest))


text:
"""
deflation criterion.
"""
function deflation_criterion(sub::T, da::T, db::T) where {T<:Number}
  deflation_criterion1(abs(sub), (abs(da) + abs(db) ))
end

@inline deflation_criterion1(sub::T, da::T) where {T<:Number} = sub <= eps(da) / 4

max_stars_repo_name: KlausC/LinAlgBigFloat

max_issues_repo_name: KlausC/LinAlgBigFloat

max_forks_repo_name: KlausC/LinAlgBigFloat


text:
"""
NearestNeighborsPBC is a module for using NearestNeighbors kdtree
implementation to find neighbors of points in a periodic square box.
This is done by simply mirroring particles across boundaries, using
the NearestNeighbors's kdtree on that larger set of points, then
relabeling indices to correspond to the original particles. The only
method currently implemented on the new periodic kdtrees is inrange(),
which can be used analogously to the inrange() defined in
NearestNeighbors. From a user perspective the only differences are
than KDTreePBC must be passed a box size, `L`, and may optionally be
passed a maximum distance, `r_max`, over which neighbors will be
sought. Also, `X` must be passed as a ndim x N matrix. Setting `r_max`
will not change the results but will avoid mirroring distant points if
    `r_max < L` increasing the speed of the calculation.

Currently NearestNeighborsPBC.KDTreePBC is only coded to use Euclidean
norms, however, it can accept sets of points in any dimension and `X`,
`L`, and `r_max` arguments of any AbstractFloat type, however, all
arguments must be of the same type.
"""
module NearestNeighborsPBC
using LinearAlgebra: norm_sqr
using NearestNeighbors
using StaticArrays

export KDTreePBC, inrange, inrangePBC_naive, smallest_separation_PBC

"""
    smallest_separation_PBC(X1, X2, L::Real)

Return the result of `X2 .- X1` subject to periodic boundary conditions in a box from `0` to `L`
along every axis. Since the elements of these points are only defined modulo `L`, the result
of the subtraction is multivalued. This function returns the separation with the smallest norm.
The elements of the result will all be bounded between `-L/2` and `L/2`.
"""
function smallest_separation_PBC(X1, X2, L::Real)
    return @. mod(X1 - X2 + L / 2, L) - L / 2
end

"""
    mirror_axes!(list::AbstractVector, mirror_pids::Vector{Int}, pid::Int, point::AbstractVector, mirror_along::Vector{Int}, L::Real, pad_distance::Real)

mirror `point` across `0` or `L` according to whether elements in `mirror_along` are -1 (mirror element across 0),
0 (don't mirror element), 1 (mirror element across `L`), pushing mirrored point to list and the pid of the point
it mirrors to `mirror_pids`. `L` sets the size of the container and d the distance from the boundary to be mirrored. Note that points in the interior are trivally mirrored when `mirror_along` is [0, 0, ...] an pushed once to list which is a flattened list of points.
"""
function mirror_axes!(list::AbstractVector, mirror_pids::Vector{Int}, pid::Int, point::AbstractVector, mirror_along::Vector{Int}, L::Real, pad_distance::Real)
    ndim = length(point)
    d = pad_distance
    mirror_point = copy(point)
    write = falses(ndim)

    for (ax, low_no_highQ) in enumerate(mirror_along)
        el = point[ax]
        if low_no_highQ == -1 && el < d
            mirror_point[ax] += L
            write[ax] = true
        elseif low_no_highQ == 0
            write[ax] = true
        elseif low_no_highQ == 1 && el > L - d
            mirror_point[ax] -= L
            write[ax] = true
        end
    end
    if all(write)
        push!(mirror_pids, pid)
        append!(list, mirror_point)
    end
end

"""
    function mirror_points_periodically(X::AbstractMatrix, L::Real, r_max::Real)

Apply periodic bounary conditions to points, `X=[[x0, y0,...] [x1, y1, ...] ...]`, which are forced by a modulo operation to live in a periodic square with `0<=xi<L`, by mirroring points within `r_max` of a boundary along an of the axes. The returned points will then have coordinates `-r_max < xi < L + r_max`.

Returns `mirror_pids`, `points_incl_mirrors` where `mirror_pids` are indices of the mirrored points in the original list of points and `points_incl_mirrors` is the set of points expanded to include mirrored points.
"""
function mirror_points_periodically(X::AbstractMatrix, L::Real, r_max::Real)
    ndim, n = size(X)

    X0 = mod.(X, L)
    mirror_points_list = collect(reshape(X0, :))
    mirror_pids = Vector{Int}(1:n)

    outer_prod = filter((x)->norm_sqr(x) > 0.5, map(collect, Iterators.product([[-1, 0, 1] for i in 1:ndim]...)))

    for pid in 1:n
        for mirror_along in outer_prod
            mirror_axes!(mirror_points_list, mirror_pids,
                         pid, X0[:, pid], mirror_along, L, r_max)
        end
    end

    mirror_points_list = reshape(mirror_points_list, (ndim, :))
    return mirror_pids, mirror_points_list
end

"""
    KDTreePBC(X::AbstractMatrix, L::Real; r_max::Real=-1.0)

Create a KD tree for points, `X=[[x0, y0,...] [x1, y1, ...] ...]`, subject to periodic bounary conditions within a periodic square with `0<=xi<L`. This is accomplished by mirroring points within `r_max` of a boundary along an of the axes. `r_max` should be specified if neighbors will be sought within distance significantly less than `L` as this minimizes the number of points that need to be mirrored.

To return a vector of vectors of neighbors of `X`, use `inrange(kdt::KDTreePBC, X, r_max)`.
"""
struct KDTreePBC{T, ndim}
    kdtree::KDTree{SArray{Tuple{ndim},T,1,ndim},Euclidean, T}
    indices_of_mirror_points::Vector{Int}
    r_max::T
    function KDTreePBC(X::Matrix{T}, L::T; r_max::T=T(-1)) where {T<:AbstractFloat}
        if r_max < 0
            r_max = L
        end
        if r_max > L
            throw(ErrorException("r_max must not be larger than L"))
        elseif L <=0
            throw(ErrorException("L must be > 0"))
        end
        mirror_pids, X_mirror = mirror_points_periodically(X, L, r_max)
        kdtree = KDTree(X_mirror)
        new{T, size(X, 1)}(kdtree, mirror_pids, r_max)
    end
end

"""
    inrange(kdt::KDTreePBC{T, ndim}, X::Matrix{T}, r_max::T) where {T, ndim}

Return a vector of vectors of neighbors of `X` within distance, `r_max`, anologously to NearestNeighbors' inrange() method only subject to periodic boundary conditions within a box of size `L` which is passed to KDTreePCB(). Passing `r_max` when building the KDTreePBC will not effect the result but will speed up the calculation when `r_max << L`.

For the naive N^2 version of the calculation, see
inrangePBC_naive(X::AbstractMatrix, L::Real, r_max::Real).
"""
function inrange(kdt::KDTreePBC{T, ndim}, X::Matrix{T}, r_max::T) where {T, ndim}
    mirror_pids = kdt.indices_of_mirror_points
    idxs = NearestNeighbors.inrange(kdt.kdtree, X, r_max)

    if r_max > kdt.r_max
        throw(ErrorException("r_max cannot be larger than kdt.r_max."))
    end

    contact_list = Vector{Int}[]
    for (pid, neigh_list) in enumerate(idxs)
        _list = Int[]
        for mi in neigh_list
            i = mirror_pids[mi]
            if i != pid
                push!(_list, i)
            end
        end
        push!(contact_list, _list)
    end
    return contact_list
end

"""
    inrangePBC_naive(X::AbstractMatrix, L::Real, r_max::Real)

calculate the neighbors subject to periodic boundary conditions
equivalent to `inrange(KDTreePBC(X, L, r_max=r_max), X, L, r_max)`
by the naive `N^2` algorithm of the distance between all pairs of points
(where `N=size(X, 2)`).

This algorithm should be faster for very small `N`, however,
it is primarily intended to be used to test `KDTreePBC` and
`inrange(kdt::KDTreePBC, L, rmax)`.
"""
function inrangePBC_naive(X::AbstractMatrix, L::Real, r_max::Real)
    r_max_sq = r_max^2

    neighbor_list = Vector{Int}[]
    for i in 1:size(X, 2)
        _list = Int[]
        for j in 1:size(X, 2)
            if i != j && norm_sqr(
                    smallest_separation_PBC(X[:, i], X[:, j], L)
                    ) < r_max_sq
                push!(_list, j)
            end
        end
        push!(neighbor_list, _list)
    end
    return neighbor_list
end

end#module

max_stars_repo_name: gregrf/Geo

max_issues_repo_name: gregrf/Geo

max_forks_repo_name: gregrf/Geo


text:
function write_message(communicator, result)
   start_result_message(communicator)
   write_bin(communicator, RESULT_INDICATOR)
   write_element(communicator, result)
   end_result_message(communicator)
end

function write_message(communicator, fail::Fail)
   start_fail_message(communicator)
   write_bin(communicator, FAIL_INDICATOR)
   write_string(communicator, fail.message)
   end_fail_message(communicator)
end

function write_callback_message(communicator, callbackid::String, args::ElementList)
   start_callback_message(communicator)
   write_bin(communicator, CALL_INDICATOR)
   write_string(communicator, callbackid)
   write_list(communicator, args)
   end_callback_message(communicator)
end

# attributes must be an iterable collection of keys (strings) and values.
function write_attributes(communicator, attributes)
   write_nattributes(communicator, length(attributes))
   for (key, value) in attributes
      write_string(communicator, key)
      write_element(communicator, value)
   end
end

function write_int32(communicator, i::Int)
   write_bin(communicator, Int32(i))
end

function write_int32s(communicator, intarr::Array{Int})
   for i in intarr
      write_int32(communicator, i)
   end
end

function write_nattributes(communicator, nattributes::Int)
   write_bin(communicator, convert(UInt8, nattributes))
end

function write_string(communicator, str::String)
   write_int32(communicator, ncodeunits(str))
   write_bin(communicator, str)
end

function write_dimensions(communicator, arr::AbstractArray)
   dims = size(arr)
   write_int32(communicator, length(dims))
   write_int32s(communicator, collect(dims))
end

function write_element(communicator, arr::Array{String}, attributes = ())
   write_bin(communicator, TYPE_ID_STRING)
   write_dimensions(communicator, arr)
   for i in eachindex(arr)
      if isassigned(arr, i)
         write_string(communicator, arr[i])
      else
         write_string(communicator, "")
      end
   end
   write_attributes(communicator, attributes)
end

function write_element(communicator, arr::Array{Union{String, Missing}})
   missings = findall(ismissing, vec(arr))
   arr2 = copy(arr)
   arr2[missings] .= ""
   strarr = convert(Array{String}, arr2)
   attributes = (("NA", missings), )
   write_element(communicator, strarr, attributes)
end

function write_element(communicator, arr::Array{F}) where {F <: SEND_AS_DOUBLE}
   attributes = (("JLTYPE", string(F)), )
   write_element(communicator, convert(Array{Float64}, arr), attributes)
end

function write_element(communicator, arr::Array{Union{Missing,F}}
      ) where {F <: Union{SEND_AS_DOUBLE}}
   attributes = (("JLTYPE", "Union{Missing," * string(F)* "}"), )
   arr2 = replace(arr, missing => R_NA_REAL)
   write_element(communicator, arr2, attributes)
end

function write_element(communicator, arr::Array{Union{Missing,Float64}})
   arr2 = replace(arr, missing => R_NA_REAL)
   write_element(communicator, arr2)
end

function write_element(communicator, arr::Array{Float64}, attributes = ())
   write_bin(communicator, TYPE_ID_FLOAT64)
   write_dimensions(communicator, arr)
   for f in arr
      write_bin(communicator, f)
   end
   write_attributes(communicator, attributes)
end

function write_uint8array_element(communicator, arr::AbstractArray{UInt8},
      attributes = ())

   if length(arr) == 1
      attributes = tuple(attributes..., ("JLDIM", 1))
   end

   write_bin(communicator, TYPE_ID_RAW)
   write_dimensions(communicator, arr)
   for d in arr
      write_bin(communicator, d)
   end
   write_attributes(communicator, attributes)
end

function write_element(communicator, arr::Array{UInt8}, attributes = ())
   write_uint8array_element(communicator, arr, attributes)
end

function write_element(communicator, arr::Array{T},
      ) where {T <: SEND_AS_INT32}

   attributes = (("JLTYPE", string(T)), )
   write_element(communicator, convert(Array{Int32}, arr), attributes, true)
end

function write_element(communicator, arr::Array{Union{Missing, T}},
      ) where {T <: Union{Int32, SEND_AS_INT32}}

   attributes = (("JLTYPE", "Union{Missing," * string(T)* "}"), )
   arr2 = convert(Array{Int32}, replace(arr, missing => R_NA_INTEGER))
   write_element(communicator, arr2, attributes, true)
end

function write_element(communicator, arr::Array{Int32},
      attributes = (), isothertype::Bool = false)

   if !isothertype
      attributes = (("JLTYPE", "Int32"), )
   end

   if length(arr) == 1
      attributes = tuple(attributes..., ("JLDIM", 1))
   end

   write_bin(communicator, TYPE_ID_INT)
   write_dimensions(communicator, arr)
   for i in arr
      write_bin(communicator, i)
   end
   write_attributes(communicator, attributes)
end

function write_element(communicator, arr::Array{Int64}; missings = Int[])

   arrmin, arrmax = extrema(arr)
   if arrmin >= typemin(Int32) && arrmax <= typemax(Int32)
      @static if Int == Int64
         attributes = ()
      else
         attributes = (("JLTYPE", "Int64"), )
      end
      if !isempty(missings)
         @static if Int != Int64
            attributes = (("JLTYPE", "Union{Missing,Int64}"), )
         end
         arr[missings] .= R_NA_INTEGER
      end
      write_element(communicator, convert(Array{Int32}, arr),
            attributes, true)
   else
      attributes = (("JLTYPE", "Int64"), )
      arr2 = convert(Array{Float64}, arr)
      if !isempty(missings)
         attributes = (("JLTYPE", "Union{Missing,Int64}"), )
         arr2[missings] .= R_NA_REAL
      end
      write_element(communicator, arr2, attributes)
      #TODO possible inexactness?
   end
end

function write_element(communicator, arr::Array{Union{Missing, Int64}})
   write_element(communicator, replace(arr, missing => 0);
         missings = findall(ismissing, arr))
end

function write_element(communicator, arr::Array{T},
      ) where {T <: SEND_AS_RAW_TYPES}

   if length(arr) == 1
      attributes = (("JLTYPE", string(T)), ("JLDIM", collect(size(arr))))
   else
      attributes = (("JLTYPE", string(T)), )
   end
   write_uint8array_element(communicator, reinterpret(UInt8, arr), attributes)
end

function write_element(communicator, arr::Array{Complex{Float64}},
      attributes = ())

   if length(arr) == 1
      attributes = tuple(attributes..., ("JLDIM", 1))
   end

   write_bin(communicator, TYPE_ID_COMPLEX)
   write_dimensions(communicator, arr)
   for d in arr
      write_bin(communicator, real(d))
      write_bin(communicator, imag(d))
   end
   write_attributes(communicator, attributes)
end

function write_element(communicator,
      arr::Array{Union{Missing,Complex{Float64}}})

   arr2 = replace(arr, missing => R_NA_COMPLEX)
   write_element(communicator, arr2)
end

function write_element(communicator, arr::Array{Union{Missing,C}},
      ) where {C <: SEND_AS_COMPLEX}

   arr2 = replace(arr, missing => R_NA_COMPLEX)
   attributes = (("JLTYPE", "Union{Missing," * string(C) * "}"), )
   write_element(communicator, arr2, attributes)
end

function write_element(communicator, arr::Array{C},
      ) where {C <: SEND_AS_COMPLEX}

   attributes = (("JLTYPE", string(C)), )
   write_element(communicator, convert.(Complex{Float64}, arr), attributes)
end

function write_element(communicator, arr::Array{Bool})
   write_bin(communicator, TYPE_ID_BOOL)
   write_dimensions(communicator, arr)
   for d in arr
      write_bin(communicator, d)
   end
end

function write_element(communicator, arr::Array{Union{Missing,Bool}})
   arr2 = replace(arr, true => Int32(1), false => Int32(0),
         missing => R_NA_INTEGER)
   attributes = (("R_LOGICAL", true), )
   write_element(communicator, arr2, attributes, true)
end

function write_array_element(communicator, arr::AbstractArray)

   attributes = Dict{String, Any}("JLTYPE" => string(typeof(arr)))

   arrsize = size(arr)

   arr2 = arr[:] # R lists support only 1 dimension
   if length(arrsize) > 1
      # we have a multidimensional array:
      # add the dimensions as attribute to be able to reconstruct it
      attributes["JLDIM"] = collect(arrsize)
   end

   ellist = ElementList(undefs_replaced(arr2), Vector{Symbol}(),
         Dict{Symbol, Any}(), attributes)
   write_element(communicator, ellist)
end

function write_element(communicator, arr::Array)
   if communicator.full_translation
      write_array_element(communicator, arr)
   else
      write_object_reference(communicator, arr, OBJECT_CLASS_ID_ARRAY)
   end
end

function write_element(communicator, arr::AbstractArray)
   if communicator.full_translation
      # translate as struct because you can never be sure of the implementation
      write_struct_element(communicator, arr)
   else
      # allow indexing as an array
      write_object_reference(communicator, arr, OBJECT_CLASS_ID_ARRAY)
   end
end

function write_element(communicator, dict::AbstractDict)

   if communicator.full_translation
      # Assure that keys and values are always in a list,
      # to allow a straightforward reconstruction later on.
      # Strings and numbers would otherwise be converted to vectors in R.
      alwaysaslist(x) = x
      alwaysaslist(x::Array{N}) where {N <: Union{Number, String}} =
            ElementList(Vector{Any}(x))

      ellist = ElementList(Vector{Any}(),
            [:keys, :values],
            Dict{Symbol, Any}(:keys => alwaysaslist(collect(keys(dict))),
                  :values => alwaysaslist(collect(values(dict)))),
            Dict{String, Any}("JLTYPE" => string(typeof(dict))))
      write_element(communicator, ellist)
   else
      write_object_reference(communicator, dict, OBJECT_CLASS_ID_STRUCT)
   end
end

function write_element(communicator, set::AbstractSet)
   if communicator.full_translation
      ellist = ElementList(Vector{Any}(collect(set)),
            Vector{Symbol}(), Dict{Symbol, Any}(),
            Dict{String, Any}("JLTYPE" => string(typeof(set))))
      write_element(communicator, ellist)
   else
      write_object_reference(communicator, set, OBJECT_CLASS_ID_STRUCT)
   end
end

function write_element(communicator, f::Function)

   callbackid = get(communicator.registered_callbacks, f, "")
   if callbackid == "" # it's not a callback function
      f_as_string = string(f)
      if endswith(f_as_string, "()") || startswith(f_as_string, '#')
         # Detected an anonymous function:
         # An anonymous function will have a string representation like
         # "getfield(Main, Symbol(\"##5#6\"))()" in Julia 1.0.
         # In Julia 1.3 it is something like "#3".
         ref = share_immutable_object!(communicator, f)
         write_bin(communicator, TYPE_ID_OBJECT_REFERENCE)
         write_bin(communicator, OBJECT_CLASS_ID_ANONYMOUS_FUNCTION)
         write_bin(communicator, ref)
         return

      else # Write an element for a named function.

         # In 1.0 f_as_string contains the full path, in 1.3 only the name.
         # We need the full path.
         # For operators, this does not work, because parsing fails then.
         if !('.' in f_as_string) && !Base.isoperator(Symbol(f_as_string))
            firstmodulepath = string(iterate(methods(f))[1].module)
            f_as_string = firstmodulepath * "." * f_as_string
         end

         write_bin(communicator, TYPE_ID_NAMED_FUNCTION)
         write_string(communicator, f_as_string)
         return
      end
   end

   write_bin(communicator, TYPE_ID_CALLBACK)
   write_string(communicator, callbackid)
end

function write_element(communicator, t::Tuple)
   if communicator.full_translation
      write_bin(communicator, TYPE_ID_LIST)
      attributes = Dict{String, Any}("JLTYPE" => string(typeof(t)))
      ellist = ElementList(
            Vector{Any}(collect(t)),
            Vector{Symbol}(), Dict{Symbol, Any}(),
            attributes)
      write_list(communicator, ellist)
   else
      write_object_reference(communicator, t, OBJECT_CLASS_ID_ARRAY)
   end
end

function write_element(communicator, ellist::ElementList)
   write_bin(communicator, TYPE_ID_LIST)
   write_list(communicator, ellist)
end

function write_element(communicator, obj::Symbol)
   write_bin(communicator, TYPE_ID_SYMBOL)
   write_string(communicator, string(obj))
end

function write_element(communicator, obj::Module)
   modulename = string(obj)
   if startswith(modulename, "Main.")
      modulename = String(modulename[6:end])
   end
   ellist = ElementList(
      Vector{Any}(),
      [:name],
      Dict{Symbol, Any}(:name => modulename),
      Dict{String, Any}("JLTYPE" => "Module"))
   write_element(communicator, ellist)
end

function write_element(communicator, obj::NamedTuple)
   if communicator.full_translation
      write_struct_element(communicator, obj)
   else
      write_object_reference(communicator, obj, OBJECT_CLASS_ID_STRUCT)
   end
end

function write_element(communicator, d::T) where {T2, T <: Type{T2}}
   write_expression(communicator, string(d))
end

function write_struct_element(communicator, obj::T,
      attributes::Dict{String, Any} = Dict{String, Any}()) where T

   names = fieldnames(T)
   if isempty(names) # type without members
      write_expression(communicator, string(T) * "()")
   else
      write_bin(communicator, TYPE_ID_LIST)
      fieldvalues = map(name -> getfield(obj, name), names)
      attributes["JLTYPE"] = string(T)
      ellist = ElementList(
            Vector(), collect(names),
            Dict{Symbol, Any}(zip(names, fieldvalues)),
            attributes)
      write_list(communicator, ellist)
   end
end

function write_mutable_struct_element(communicator, obj::T) where T

   ref = share_mutable_object!(communicator, obj)
   refexists = ref in communicator.objrefs
   push!(communicator.objrefs, ref)

   if refexists # recursion detected
      decrefcount!(communicator, ref) # not shared, actually
      write_struct_element(communicator, CircularReference(ref),
            Dict{String, Any}())
      return
   else
      # write object, giving reference to original object
      write_struct_element(communicator, obj,
            Dict{String, Any}("JLREF" => ref))
      return
   end
end

function write_element(communicator, obj::T) where T

   if isstructtype(T)
      if communicator.full_translation
         if !isimmutable(obj)
            write_mutable_struct_element(communicator, obj)
         else
            write_struct_element(communicator, obj)
         end
      else
         write_object_reference(communicator, obj, OBJECT_CLASS_ID_STRUCT)
      end
   else
      write_struct_element(communicator,
            Fail("Lost in translation: $(string(obj))"))
   end
end

function write_element(communicator, f::F) where {F <: SEND_AS_DOUBLE}
   attributes = (("JLTYPE", string(F)), )
   write_element(communicator, Float64(f), attributes)
end

function write_element(communicator, m::Missing)
   write_element(communicator, R_NA_REAL)
end

function write_element(communicator, m::Array{Missing})
   write_element(communicator, fill(R_NA_REAL, size(m)))
end

function write_element(communicator, f::Float64, attributes = ())
   write_bin(communicator, TYPE_ID_FLOAT64)
   write_int32(communicator, 0)
   write_bin(communicator, f)
   write_attributes(communicator, attributes)
end

function write_element(communicator, i::T) where {T <: SEND_AS_INT32}
   attributes = (("JLTYPE", string(T)), )
   write_element(communicator, Int32(i), attributes, true)
end

function write_element(communicator, i::Int32,
      attributes = (), isothertype::Bool = false)

   if !isothertype
      attributes = (("JLTYPE", "Int32"), )
   end

   write_bin(communicator, TYPE_ID_INT)
   write_int32(communicator, 0)
   write_bin(communicator, i)
   write_attributes(communicator, attributes)
end

function write_element(communicator, i::Int64)
   if typemin(Int32) <= i <= typemax(Int32)
      @static if Int == Int64
         attributes = ()
      else
         attributes = (("JLTYPE", "Int64"), )
      end
      write_element(communicator, Int32(i), attributes, true)
   else
      attributes = (("JLTYPE", "Int64"), )
      # TODO inexactness?
      write_element(communicator, Float64(i), attributes)
   end
end

function write_element(communicator, i::T) where {T <: SEND_AS_RAW_TYPES}
   attributes = (("JLTYPE", string(T)), )
   write_uint8array_element(communicator, reinterpret(UInt8, [i]), attributes)
end

function write_element(communicator, b::Bool)
   write_bin(communicator, TYPE_ID_BOOL)
   write_int32(communicator, 0)
   write_bin(communicator, b)
end

function write_element(communicator, str::String)
   write_bin(communicator, TYPE_ID_STRING)
   write_int32(communicator, 0)
   write_string(communicator, str)
   write_bin(communicator, 0x00) # no attributes
end

function write_element(communicator, c::Complex{Float64}, attributes = ())
   write_bin(communicator, TYPE_ID_COMPLEX)
   write_int32(communicator, 0)
   write_bin(communicator, Float64(real(c)))
   write_bin(communicator, Float64(imag(c)))
   write_attributes(communicator, attributes)
end

function write_element(communicator, c::C) where {C <: SEND_AS_COMPLEX}
   attributes = (("JLTYPE", string(C)), )
   write_element(communicator, convert(Complex{Float64}, c), attributes)
end

function write_element(communicator, u::UInt8)
   write_bin(communicator, TYPE_ID_RAW)
   write_int32(communicator, 0)
   write_bin(communicator, u)
   write_bin(communicator, 0x00) # no attributes
end

function write_element(communicator, n::Nothing)
   write_bin(communicator, TYPE_ID_NOTHING)
end

function write_element(communicator, p::EnforcedProxy{T}) where {T <: Array}
   write_object_reference(communicator, p.obj, OBJECT_CLASS_ID_SIMPLE_ARRAY)
end

function write_element(communicator, p::EnforcedProxy)
   write_object_reference(communicator, p.obj, OBJECT_CLASS_ID_NO_INFO)
end

function write_element(communicator, df::RDataFrame)
   if communicator.full_translation
      write_bin(communicator, TYPE_ID_LIST)
      write_list(communicator, getfield(df, :ellist))
   else
      write_object_reference(communicator, df, OBJECT_CLASS_ID_STRUCT)
   end
end

# This function and the struct serve only to communicate the CommunicatoR
# object of this client to R.
struct GetCommunicatoR
end
function write_element(communicator, _::GetCommunicatoR)
   ref = object_reference(communicator)
   write_object_reference(communicator, communicator, OBJECT_CLASS_ID_STRUCT)
end

function write_expression(communicator, str::AbstractString)
   write_bin(communicator, TYPE_ID_EXPRESSION)
   write_string(communicator, str)
end

function write_list(communicator, arr::AbstractArray)
   write_list(communicator, ElementList(vec(arr)))
end

function write_list(communicator, ellist::ElementList)
   write_int32(communicator, length(ellist.positionalelements))
   for el in ellist.positionalelements
      write_element(communicator, el)
   end

   write_int32(communicator, length(ellist.namedelements))
   for name in ellist.names # use the order of the names here
      write_string(communicator, string(name))
      write_element(communicator, ellist.namedelements[name])
   end

   write_attributes(communicator, ellist.attributes)
end

function write_object_reference(communicator, obj, object_class_id::UInt8)
   ref = shareobject!(communicator, obj)
   write_bin(communicator, TYPE_ID_OBJECT_REFERENCE)
   write_bin(communicator, object_class_id)
   write_bin(communicator, ref)
end


text:
const actions = Dict( "nop" => (i,x,p) -> (i+1, x),
                      "acc" => (i,x,p) -> (i+1, x+p),
                      "jmp" => (i,x,p) -> (i+p, x) )

function bootcode( filename )
    open(filename) do f
        return split.(readlines(f))
    end
end

function get_value!(instructions, executed, index, value)
    a, p = instructions[index]
    executed[index] = true
    index, value = actions[a]( index, value, parse(Int,p))
    if index > length(instructions)
        return nothing, value
    elseif executed[index]
        return index, nothing
    else
        return index, value
    end
end

function run_code( instructions )

    index = 1
    value = 0
    executed = [false for i in eachindex(instructions)]
    while true
        last_value = value
        index, value = get_value!(instructions, executed, index, value)
        index == nothing && ( return value, true )
        value == nothing && ( return last_value, false)
    end

end

function switch( instructions, position )

    action = instructions[position][1]
    new_instructions = [ copy(x) for x in instructions]

    if action == "nop"
        new_instructions[position][1] = "jmp"
    elseif action == "jmp"
        new_instructions[position][1] = "nop"
    end
    return new_instructions

end

instructions = bootcode("input08.txt")

println(" Part 1 : ", first(run_code( instructions)))
positions = findall( a -> ( a[1] in ["jmp","nop"]), instructions)

for position in positions
    new_instructions = switch(instructions, position)
    value, check = run_code(new_instructions)
    if check
       println(" Part 2 : ", value)
       break
    end
end

max_stars_repo_name: pnavaro/AdventOfCode

max_issues_repo_name: pnavaro/AdventOfCode

max_forks_repo_name: pnavaro/AdventOfCode


text:
export idx
export local_basis

idx(σ::Int) = (σ == -1) ? 1 : σ + 1
local_basis(d::Int) = union(-1, 1:d-1)

function all_states(rank::Union{Vector, NTuple})
    basis = [local_basis(r) for r ∈ rank]
    Iterators.product(basis...)
end


text:
module LangevinDynamics

using Random
using LinearAlgebra

using ..Utilities: get_index
using ..Models: AbstractModel, HolsteinModel, SSHModel, update_model!, mulMᵀ!, muldMdx!
using ..PhononAction: calc_dSbdx!
using ..FourierAcceleration: FourierAccelerator, fourier_accelerate!

using ..KPMPreconditioners: setup!

export evolve!, Dynamics, EulerDynamics, RungeKuttaDynamics, HeunsDynamics
export calc_dSdx!, calc_dSfdx!

"""
Abstract type for representing different way of evolving Langevin Dynamics.
"""
abstract type Dynamics end

##################
## EULER UPDATE ##
##################

struct EulerDynamics{T<:AbstractFloat} <: Dynamics

    """
    Number of degrees of freedom (DOF) to update.
    """
    Ndof::Int

    """
    Dimension of M matrix.
    """
    Ndim::Int

    """
    Timestep.
    """
    Δt::T

    """
    Derivative of action with respect to each DOF.
    """
    dSdx::Vector{T}

    """
    Noise vector in Langevin dynamics.
    """
    η::Vector{T}

    """
    Change in phonon fields.
    """
    Δx::Vector{T}

    """
    Random noise vector.
    """
    R::Vector{T}

    """
    M⁻¹⋅R
    """
    M⁻¹R::Vector{T}

    function EulerDynamics(model::AbstractModel, Δt::T) where {T<:AbstractFloat}

        Ndof     = model.Ndof
        Ndim     = model.Ndim
        dSdx     = zeros(T,Ndof)
        η        = zeros(T,Ndof)
        Δx       = zeros(T,Ndof)
        R        = zeros(T,Ndim)
        M⁻¹R     = zeros(T,Ndim)

        return new{T}(Ndof,Ndim,Δt,dSdx,η,Δx,R,M⁻¹R)
    end
end

function evolve!(model::AbstractModel{T1,T2,T3}, dyn::EulerDynamics{T1}, fa::FourierAccelerator{T1}, preconditioner=I)::Int  where {T1,T2,T3}

    Δt       = dyn.Δt
    dSdx     = dyn.dSdx
    QdSdx    = dyn.dSdx
    η        = dyn.η
    sqrtQη   = dyn.η
    Δx       = dyn.Δx
    R        = dyn.R
    M⁻¹R     = dyn.M⁻¹R

    # update the exponentiated interaction matrix so that it reflects the current
    # phonon field configuration.
    update_model!(model)

    # itialize η as vector of gaussian random number
    randn!(η,model)

    # calculate dSdx = [∂S/∂x₁(1),...,∂S/∂xₙ(1),...,∂S/∂x₁(τ),...,∂S/∂xₙ(τ),...,∂S/∂x₁(Lτ),...,∂S/∂xₙ(Lτ)]
    randn!(model.rng,R)
    iters = calc_dSdx!(dSdx, R, M⁻¹R, model, preconditioner)

    # dSdx ==> Q⋅dSdx
    fourier_accelerate!( QdSdx, fa , dSdx, 1.0 )

    # η = √Q⋅η
    fourier_accelerate!( sqrtQη, fa , η, 0.5 )

    # Δx = √(2⋅Δt⋅Q)⋅η - Δt⋅Q⋅dS/dx
    @. Δx = sqrt(2.0*Δt)*sqrtQη - Δt*QdSdx

    # update phonon fields
    @. model.x += Δx

    # update the exponentiated interaction matrix to reflect current phonon field configuration.
    update_model!(model)

    return iters
end

########################
## RUNGE-KUTTA UPDATE ##
########################

# Implemented Definition of Runge-Kutta update with Fourier Acceleration:
# 1. initialize η
# 2. calcualte dS/dx
# 3. x′ = x - Δt⋅dS/dx + √(2Δt)⋅η
# 4. calculate dS/dx′
# 5. x″ = x - Δt⋅[F⁻¹⋅Q⋅F]⋅(dS/dx′+dS/dx)/2 + √(2Δt)⋅[F⁻¹⋅√Q⋅F]⋅η
# Note: Fourier Acceleration Only Applied in Step 5
# Note: F is fourier transform from τ ⟶ ω
# Note: Q is the diagonal acceleration matrix

struct RungeKuttaDynamics{T<:AbstractFloat} <: Dynamics

    Ndof::Int
    Ndim::Int
    Δt::T
    dSdx::Vector{T}
    dSdx′::Vector{T}
    η::Vector{T}
    Δx::Vector{T}
    R::Vector{T}
    M⁻¹R::Vector{T}

    function RungeKuttaDynamics(model::AbstractModel, Δt::T) where {T<:AbstractFloat}

        Ndof     = model.Ndof
        Ndim     = model.Ndim
        dSdx     = zeros(T,Ndof)
        dSdx′    = zeros(T,Ndof)
        η        = zeros(T,Ndof)
        Δx       = zeros(T,Ndof)
        R        = zeros(T,Ndim)
        M⁻¹R     = zeros(T,Ndim)

        return new{T}(Ndof,Ndim,Δt,dSdx,dSdx′,η,Δx,R,M⁻¹R)
    end
end

function evolve!(model::AbstractModel{T1,T2,T3}, dyn::RungeKuttaDynamics{T1}, fa::FourierAccelerator{T1}, preconditioner=I)::Int  where {T1,T2,T3}

    Δt       = dyn.Δt
    dSdx     = dyn.dSdx
    QdSdx    = dyn.dSdx
    dSdx′    = dyn.dSdx′
    η        = dyn.η
    sqrtQη   = dyn.η
    Δx       = dyn.Δx
    R        = dyn.R
    M⁻¹R     = dyn.M⁻¹R
    x        = model.x
    x′       = model.x
    x″       = model.x

    # update the exponentiated interaction matrix to reflect current phonon field configuration.
    update_model!(model)

    # itialize η as vector of gaussian random number
    randn!(η,model)

    # calculate dSdx = [∂S/∂x₁(1),...,∂S/∂xₙ(1),...,∂S/∂x₁(τ),...,∂S/∂xₙ(τ),...,∂S/∂x₁(Lτ),...,∂S/∂xₙ(Lτ)]
    randn!(model.rng,R)
    iters = calc_dSdx!(dSdx, R, M⁻¹R, model, preconditioner)

    # get the update for the fields using euler method
    @. Δx = sqrt(2*Δt)*η - Δt*dSdx

    # update phonon fields
    @. x′ = x + Δx

    # update the exponentiated interaction matrix so that it reflects the current phonon field configuration.
    update_model!(model)

    # calculate dSdx = [∂S/∂x₁(1),...,∂S/∂xₙ(1),...,∂S/∂x₁(τ),...,∂S/∂xₙ(τ),...,∂S/∂x₁(Lτ),...,∂S/∂xₙ(Lτ)]
    randn!(model.rng,R)
    iters = calc_dSdx!(dSdx′, R, M⁻¹R, model, preconditioner)

    # revert back to original phonon fields
    @. x = x′ - Δx

    # update the exponentiated interaction matrix so that it reflects the current phonon field configuration.
    update_model!(model)

    # get the partial derivative for the second RK step
    @. dSdx = (dSdx′+dSdx)/2.0

    # dS/dx ==> Q⋅dS/dx
    fourier_accelerate!(QdSdx,fa,dSdx,1.0)

    # η ==> √Q⋅η
    fourier_accelerate!(sqrtQη,fa,η,0.5)

    # Δx = √(2⋅Δt⋅Q)⋅η - Δt⋅Q⋅dS/dx
    @. Δx = sqrt(2.0*Δt)*sqrtQη - Δt*QdSdx

    # update phonon fields
    @. x″ = x + Δx

    # update the exponentiated interaction matrix so that it reflects the current phonon field configuration.
    update_model!(model)

    return iters
end

#####################
## HEUN'S DYNAMICS ##
#####################

# Implemented Definition of Heun's update with Fourier Acceleration:
# 1.  initialize η
# 2.  ξ      = [F⁻¹⋅√Q⋅F]⋅η
# 3.  calcualte dS/dx
# 4.  dΓ/dx  = [F⁻¹⋅Q⋅F]⋅dS/dx
# 5.  Δx     = √(2Δt)⋅ξ - Δt⋅dΓ/dx
# 6.  x′     = x + Δx
# 7.  calculate dS/dx′
# 8.  dΓ/dx′ = [F⁻¹⋅Q⋅F]⋅dS/dx′
# 9.  x      = x′- Δx
# 10. x″     = x + √(2Δt)⋅ξ - Δt⋅(dΓ/dx+dΓ/dx′)/2
# Note: F is fourier transform from τ ⟶ ω
# Note: Q is the diagonal acceleration matrix

struct HeunsDynamics{T<:AbstractFloat} <: Dynamics

    Ndof::Int
    Ndim::Int
    Δt::T
    η::Vector{T}
    dSdx::Vector{T}
    dSdx′::Vector{T}
    Δx::Vector{T}
    R::Vector{T}
    M⁻¹R::Vector{T}

    function HeunsDynamics(model::AbstractModel,Δt::T) where {T<:AbstractFloat}

        Ndof  = model.Ndof
        Ndim  = model.Ndim
        η     = zeros(T,Ndof)
        dSdx  = zeros(T,Ndof)
        dSdx′ = zeros(T,Ndof)
        Δx    = zeros(T,Ndof)
        R     = zeros(T,Ndim)
        M⁻¹R  = zeros(T,Ndim)

        return new{T}(Ndof,Ndim,Δt,η,dSdx,dSdx′,Δx,R,M⁻¹R)
    end
end

function evolve!(model::AbstractModel{T1,T2,T3}, dyn::HeunsDynamics{T1}, fa::FourierAccelerator{T1}, preconditioner=I)::Int  where {T1,T2,T3}

    x        = model.x
    x′       = model.x
    x″       = model.x
    Δx       = dyn.Δx
    Δt       = dyn.Δt
    η        = dyn.η
    dSdx     = dyn.dSdx
    dΓdx     = dyn.dSdx
    dSdx′    = dyn.dSdx′
    dΓdx′    = dyn.dSdx′
    R        = dyn.R
    M⁻¹R     = dyn.M⁻¹R

    # 1. intialize η
    randn!(η,model)

    # 2. ξ = [F⁻¹⋅√Q⋅F]⋅η
    ξ = η
    fourier_accelerate!(ξ,fa,η,0.5)

    # 3. calcualte dS/dx
    update_model!(model)
    randn!(model.rng,R)
    iters1 = calc_dSdx!(dSdx, R, M⁻¹R, model, preconditioner)

    # 4. dΓ/dx  = [F⁻¹⋅Q⋅F]⋅dS/dx
    fourier_accelerate!(dΓdx,fa,dSdx,1.0)

    # 5. Δx = √(2Δt)⋅ξ - Δt⋅dΓ/dx
    @. Δx = sqrt(2*Δt)*ξ - Δt*dΓdx

    # 6. x′ = x + Δx
    @. x′ = x + Δx
    update_model!(model)

    # 7. calculate dS/dx′
    randn!(model.rng,R)
    iters2 = calc_dSdx!(dSdx′, R, M⁻¹R, model, preconditioner)

    # 8. dΓ/dx′ = [F⁻¹⋅Q⋅F]⋅dS/dx′
    fourier_accelerate!(dΓdx′,fa,dSdx′,1.0)

    # 9. x = x′- Δx
    @. x = x′ - Δx

    # 10. x″ = x + √(2Δt)⋅ξ - Δt⋅(dΓ/dx+dΓ/dx′)/2
    @. x″ = x + sqrt(2*Δt)*ξ - Δt*(dΓdx+dΓdx′)/2
    update_model!(model)

    return div(iters1+iters2,2)
end

#################################################
## FUNCTION FOR CALCULATING FORCE IN DYNAMICS ##
#################################################

"""
Calculates all of the partial derivatives ∂S/∂xᵢ(τ) and stores each partial derivative in a vector dSdx.
The expression we are evaluating is `∂S/∂xᵢ(τ) = ∂Sbose/∂xᵢ(τ) - 2gᵀ[∂M/∂xᵢ(τ)]⋅M⁻¹g`.
"""
function calc_dSdx!(dSdx::AbstractVector{T2},g::AbstractVector{T2},M⁻¹g::AbstractVector{T2},
                    model::AbstractModel{T1,T2,T3}, preconditioner=I)::Int where {T1,T2,T3}

    # ∂S/∂xᵢ(τ) = -2gᵀ⋅[∂M/∂xᵢ(τ)]⋅M⁻¹g
    iters = calc_dSfdx!(dSdx, g, M⁻¹g, model, preconditioner)

    # ∂S/∂xᵢ(τ) = ∂Sbose/∂xᵢ(τ) - 2gᵀ⋅[∂M/∂xᵢ(τ)]⋅M⁻¹g ==> All Done!
    calc_dSbdx!( dSdx , model , true )

    # returning number of iterations in solving for M⁻¹g
    return iters
end

"""
Calculate just the force associated with the fermionic part of the action.
"""
function calc_dSfdx!(dSfdx::AbstractVector{T2},g::AbstractVector{T2},M⁻¹g::AbstractVector{T2},
                     model::AbstractModel{T1,T2,T3}, preconditioner=I)::Int where {T1,T2,T3}

    # NOTE: The method I use below for calculating all of the partial derivatives {∂S/∂xᵢ(τ)} works only
    # because I am assuming a strictly local form of the electron-phonon interaction given by ∑(λᵢ⋅xᵢ⋅nᵢ).
    # If a longer-range electron-phonon interaction were added of the form ∑(λᵢⱼ⋅xᵢ⋅nⱼ), then this way of doing
    # things would no longer work, and a loop over each individual phonon field xᵢ(τ) would need to be added.

    # intialize random vector g.
    randn!(model.rng,g)

    # solve linear system to get M⁻¹⋅g
    iters = 0
    setup!(preconditioner) # setup block preconditioner
    fill!(M⁻¹g,0.0)
    if model.mul_by_M
        # solve M⋅x=g ==> x=M⁻¹⋅g
        model.transposed = false
        iters, err, flag = ldiv!(M⁻¹g,model,g,preconditioner)
    else
        model.transposed=false
        # solve MᵀM⋅x=Mᵀg ==> x=[MᵀM]⁻¹⋅Mᵀg=M⁻¹⋅g
        mulMᵀ!(model.v″,model,g)
        iters, err, flag = ldiv!(M⁻¹g,model,model.v″,preconditioner)
    end

    # ⟨∂M/∂xᵢ(τ)⟩=gᵀ⋅∂M/∂xᵢ(τ)⋅M⁻¹g
    muldMdx!( dSfdx , g, model , M⁻¹g )

    # ∂Sf/∂xᵢ(τ) = -2⟨∂M/∂xᵢ(τ)⟩ = -2gᵀ⋅∂M/∂xᵢ(τ)⋅M⁻¹g
    @. dSfdx = -2.0 * dSfdx

    return iters
end

end

max_stars_repo_name: cohensbw/Langevin

max_issues_repo_name: cohensbw/ElPhDynamics

max_forks_repo_name: cohensbw/ElPhDynamics


text:
import Test: @test, @testset, @test_throws

@testset "Radial" begin
    # mode average and radial integral for single mode and only Kerr should be identical
    using Luna
    import LinearAlgebra: norm
    a = 13e-6
    gas = :Ar
    pres = 5
    τ = 30e-15
    λ0 = 800e-9
    grid = Grid.RealGrid(5e-2, 800e-9, (160e-9, 3000e-9), 1e-12)
    m = Capillary.MarcatiliMode(a, gas, pres, loss=false)
    aeff(z) = Modes.Aeff(m, z=z)
    energyfun, energyfunω = Fields.energyfuncs(grid)

    dens0 = PhysData.density(gas, pres)
    densityfun(z) = dens0
    responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
    linop, βfun!, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0)
    inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6)
    Eω, transform, FT = Luna.setup(
        grid, densityfun, responses, inputs, βfun!, aeff)
    statsfun = Stats.collect_stats(grid, Eω,
                               Stats.ω0(grid),
                               Stats.energy(grid, energyfunω))
    output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
    Luna.run(Eω, grid, linop, transform, FT, output, status_period=5)

    modes = (
         Capillary.MarcatiliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false),
    )
    energyfun, energyfunω = Fields.energyfuncs(grid)
    inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6)
    Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs,
                                modes, :y; full=false)
    outputr = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
    linop = LinearOps.make_const_linop(grid, modes, λ0)
    Luna.run(Eω, grid, linop, transform, FT, outputr, status_period=10)

    Iω = abs2.(output.data["Eω"])
    Iωr = abs2.(dropdims(outputr.data["Eω"], dims=2))
    @test norm(Iω - Iωr)/norm(Iω) < 0.003
end

@testset "Full" begin
    # mode average and full integral for single mode and only Kerr should be identical
    using Luna
    import LinearAlgebra: norm
    a = 13e-6
    gas = :Ar
    pres = 5
    τ = 30e-15
    λ0 = 800e-9
    grid = Grid.RealGrid(5e-2, 800e-9, (160e-9, 3000e-9), 1e-12)
    m = Capillary.MarcatiliMode(a, gas, pres, loss=false)
    aeff(z) = Modes.Aeff(m, z=z)
    energyfun, energyfunω = Fields.energyfuncs(grid)

    dens0 = PhysData.density(gas, pres)
    densityfun = let dens0=PhysData.density(gas, pres)
        z -> dens0
    end
    responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
    linop, βfun!, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0)
    inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6)
    Eω, transform, FT = Luna.setup(
        grid, densityfun, responses, inputs, βfun!, aeff)
    statsfun = Stats.collect_stats(grid, Eω,
                               Stats.ω0(grid),
                               Stats.energy(grid, energyfunω))
    output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
    Luna.run(Eω, grid, linop, transform, FT, output, status_period=5)

    modes = (
         Capillary.MarcatiliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false),
    )
    energyfun, energyfunω = Fields.energyfuncs(grid)
    inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6)
    Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs,
                                modes, :y; full=true)
    outputf = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
    linop = LinearOps.make_const_linop(grid, modes, λ0)
    Luna.run(Eω, grid, linop, transform, FT, outputf, status_period=10)

    Iω = abs2.(output.data["Eω"])
    Iωf = abs2.(dropdims(outputf.data["Eω"], dims=2))
    @test norm(Iω - Iωf)/norm(Iω) < 0.003
end

@testset "Full, LP11 vs TM01" begin
    # TM01 and LP11 (which is TM01+HE21) have the same dispersion and a ratio of effective
    # areas of 2/3 -- mode averaged propagation in TM01 with 2/3x Aeff should be identical
    # to full modal propagation in LP11, as long as only Kerr is considered.
    using Luna
    import LinearAlgebra: norm
    a = 13e-6
    gas = :Ar
    pres = 5
    τ = 30e-15
    λ0 = 800e-9
    energy = 1e-6
    grid = Grid.RealGrid(5e-2, 800e-9, (160e-9, 3000e-9), 0.4e-12)
    m = Capillary.MarcatiliMode(a, gas, pres; kind=:TM, n=0, m=1, loss=false)
    aeff(z) = 2/3*Modes.Aeff(m, z=z) # Aeff of LP11 is 2/3x smaller than that of TM01
    energyfun, energyfunω = Fields.energyfuncs(grid)

    dens0 = PhysData.density(gas, pres)
    densityfun = let dens0=PhysData.density(gas, pres)
        z -> dens0
    end
    responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
    linop, βfun!, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0)
    inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=energy)
    Eω, transform, FT = Luna.setup(
        grid, densityfun, responses, inputs, βfun!, aeff)
    statsfun = Stats.collect_stats(grid, Eω,
                               Stats.ω0(grid),
                               Stats.energy(grid, energyfunω))
    output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
    Luna.run(Eω, grid, linop, transform, FT, output, status_period=5)

    # vertical LP11 double-lobe
    modes = (
        Capillary.MarcatiliMode(a, gas, pres, n=0, m=1, kind=:TM, ϕ=0.0, loss=false),
        Capillary.MarcatiliMode(a, gas, pres, n=2, m=1, kind=:HE, ϕ=-π/4, loss=false),
    )

    field = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=energy/2)
    inputs = ((mode=1, fields=(field,)), (mode=2, fields=(field,)))
    # inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=energy)
    Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs,
                                   modes, :xy; full=true)
    outputf = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
    linop = LinearOps.make_const_linop(grid, modes, λ0)
    Luna.run(Eω, grid, linop, transform, FT, outputf, status_period=10)

    Iω = abs2.(output["Eω"])
    Iωf = dropdims(sum(abs2.(outputf["Eω"]); dims=2), dims=2)
    @test norm(Iω - Iωf)/norm(Iω) < 0.0006
end

@testset "FieldInputs" begin
    using Luna
    import LinearAlgebra: norm
    a = 13e-6
    gas = :Ar
    pres = 5
    τ = 30e-15
    λ0 = 800e-9
    grid = Grid.RealGrid(5e-2, 800e-9, (160e-9, 3000e-9), 1e-12)
    modes = (
         Capillary.MarcatiliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false),
         Capillary.MarcatiliMode(a, gas, pres, n=1, m=2, kind=:HE, ϕ=0.0, loss=false)
    )
    aeff(z) = Modes.Aeff(m, z=z)
    energyfun, energyfunω = Fields.energyfuncs(grid)
    dens0 = PhysData.density(gas, pres)
    densityfun(z) = dens0
    responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
    inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6)
    Eω_single, transform, FT = Luna.setup(grid, densityfun, responses, inputs,
                                modes, :y; full=false)
    inputs = (Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6),)
    Eω_tuple, transform, FT = Luna.setup(grid, densityfun, responses, inputs,
                                modes, :y; full=false)
    inputs = ((mode=1, fields=(Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6),)),)
    Eω_tuple_of_namedtuples, transform, FT = Luna.setup(grid, densityfun, responses, inputs,
                                modes, :y; full=false)
    @test Eω_single ≈ Eω_tuple
    @test Eω_single ≈ Eω_tuple_of_namedtuples
end

@testset "Nonlinear coupling" begin
    using Luna
    import LinearAlgebra: norm
    a = 125e-6
    gas = :Ar
    pres = 0.167
    flength = 3

    τfwhm = 10e-15
    λ0 = 800e-9
    energy = 150e-6

    # HE11 and LP11 (consisting of HE21 + TM01)
    modes = (
        Capillary.MarcatiliMode(a, gas, pres, n=0, m=1, kind=:TM, ϕ=0.0),
        Capillary.MarcatiliMode(a, gas, pres, n=2, m=1, kind=:HE, ϕ=float(-π/4)),
        Capillary.MarcatiliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=float(0.0)),
    )

    grid = Grid.RealGrid(flength, λ0, (160e-9, 3000e-9), 0.4e-12)

    densityfun = let dens0=PhysData.density(gas, pres)
        z -> dens0
    end
    responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)

    field = Fields.GaussField(;λ0, τfwhm, energy=energy/2)
    inputs = ((mode=1, fields=(field,)), (mode=2, fields=(field,)))

    Eω, transform, FT = Luna.setup(
        grid, densityfun, responses, inputs, modes, :xy; full=true)
    nl = similar(Eω)

    transform(nl, Eω, 0.0);
    # nonlinear polarisation in LP11
    Inl12 = dropdims(sum(abs2.(nl[:, 1:2]); dims=2); dims=2)
    # nonlinear polarisation in HE11
    Inl3 = abs2.(nl[:, 3])
    # LP11 should not couple nonlinearly to HE11
    @test norm(Inl3)/norm(Inl12) < 1e-32
end


text:
immutable LedgerReport <: Report
    heading::Compat.UTF8String
    transactions::Vector{LedgerEntry}
end


text:
module DelaySSAdocs

using DocStringExtensions
# Write your package code here.
include("doc_funcs.jl")
include("aggregators.jl")

end


text:
# A transformation for testing the local affine mapping
struct SquareMe <: Transformation; end
(::SquareMe)(x) = x.^2
CoordinateTransformations.transform_deriv(::SquareMe, x0) = Matrix(Diagonal(2*x0))

struct Ray37
    origin::Vector
    direction::Vector
end
function (m::LinearMap)(r::Ray37)
    Ray37(m(r.origin),
          m(r.direction))
end
function (t::Translation)(r::Ray37)
    Ray37(t(r.origin),
          r.direction)
end

@testset "Common Transformations" begin
    @testset "AffineMap" begin
        @testset "Simple" begin
            M = [1 2; 3 4]
            v = [-1, 1]
            x = [1,0]
            y = [0.0,1.0]
            A = AffineMap(M,v)
            @test A(x) == M*x + v
        end

        @testset "composition " begin
            M1 = [1 2; 3 4]
            v1 = [-1, 1]
            A1 = AffineMap(M1,v1)
            M2 = [0 1; 1 0]
            v2 = [-2, 0]
            A2 = AffineMap(M2,v2)
            x = [1,0]
            y = [0,1]
            @test A1(x) == M1*x + v1
            @test A1(y) == M1*y + v1
            @test (A2∘A1)(x) == M2*(M1*x + v1) + v2
            @test (A2∘A1)(y) == M2*(M1*y + v1) + v2
        end

        @testset "inverse" begin
            M = [1.0 2.0; 3.0 4.0]
            v = [-1.0, 1.0]
            x = [1,0]
            y = [0.0,1.0]
            A = AffineMap(M,v)
            @test inv(A)(A(x)) ≈ x
            @test inv(A)(A(y)) ≈ y
        end

        @testset "Affine approximation" begin
            S = SquareMe()
            x0 = [1,2,3]
            dx = 0.1*[1,-1,1]
            A = AffineMap(S, x0)
            @test isapprox(S(x0 + dx), A(x0 + dx), atol=maximum(2*dx.^2))
        end
    end

    @testset "LinearMap" begin
        M = [1 2; 3 4]
        x = [1,0]
        y = [0,1]
        L = LinearMap(M)
        @test L(x) == M*x
        @test inv(L)(x) == inv(M)*x
        @test inv(L)(y) == inv(M)*y
        @test inv(L)(L(x)) ≈ x
        @test inv(L)(L(y)) ≈ y
        @test (L∘L)(x) == (M*M)*x
    end

    @testset "Translation" begin
        x = SVector(1.0, 2.0)
        trans = Translation(2.0, -1.0)

        # Inverse
        @test inv(trans) == Translation(-2.0, 1.0)

        # Composition
        @test trans ∘ trans == Translation(4.0, -2.0)

        # Transform
        @test trans(x) === SVector(3.0, 1.0)
        @test trans(collect(x)) == [3.0, 1.0]

        # Transform derivative
        m1 = transform_deriv(trans, x)
        m2 = [m1[i,j] for i=1:2, j=1:2] # m1 might be a UniformScaling
        @test m2 == Matrix(I, 2, 2)  # in v0.7 and above, this can just be `== I`

        # Transform parameter derivative
        m1 = transform_deriv_params(trans, x)
        m2 = [m1[i,j] for i=1:2, j=1:2] # m1 might be a UniformScaling
        @test m2 == Matrix(I, 2, 2)  # In v0.7 and above, this can just be `== I`
    end

    @testset "Recenter" begin
        M = [1 2; 3 4]
        for origin in ([5,-3], @SVector([5,-3]))
            c = recenter(M, origin)
            @test c(origin) == origin
            @test c(zero(origin)) == [6,-6]
        end
    end

    @testset "application of AffineMap in terms of LinearMap and Translation" begin
        origin = randn(3)
        direction = randn(3)
        linear = randn(3,3)
        trans = randn(3)
        r = Ray37(origin, direction)
        expected = Ray37(linear*origin+trans, linear*direction)
        result   = AffineMap(linear, trans)(r)
        @test expected.origin ≈ result.origin
        @test expected.direction ≈ result.direction
    end
end


text:
module CodecXz

export
    XzCompressor,
    XzCompressorStream,
    XzDecompressor,
    XzDecompressorStream

import TranscodingStreams:
    TranscodingStreams,
    TranscodingStream,
    Memory,
    Error,
    initialize,
    finalize,
    splitkwargs
using Libdl

const liblzmapath = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(liblzmapath)
    error("CodecXz.jl is not installed properly, run Pkg.build(\"CodecXz\") and restart Julia.")
end
include(liblzmapath)
check_deps()

include("liblzma.jl")
include("compression.jl")
include("decompression.jl")

end # module


text:
# ------------------------------------------------------------
# Exact Fourier transforms of Haar scaling function

@doc """
	FourHaarScaling(xi[, J, k])

The Fourier transform of the Haar scaling function on scale `J` and translation `k` evaluated at `xi`.
If not supplied, `J = 0` and `k = 0`.
"""->
function FourHaarScaling(xi::Real)
	if xi == zero(xi)
		return ComplexOne
	else
		return (1.0 - cis(-twoπ*xi)) / (twoπ*im*xi)
	end
end

# ------------------------------------------------------------
# Fourier transform of Daubechies scaling functions

@doc """
	feval(xi::Float64, C::Vector, offset::Int) -> Complex[{Float64}

*F*ilter *eval*uation at `xi` of the filter `C`, i.e., compute

	sum( C[n]*exp(-2*pi*(offset+n)*xi) )

The sum starts at `offset+1`.
"""->
function feval(xi::Float64, C::Vector{Float64}, offset::Integer)
	y = ComplexZero
	for n in 1:length(C)
		@inbounds y += C[n]*cis(-twoπ*(offset+n)*xi)
	end
	return y
end

@doc """
	FourDaubScaling(xi, C [, p]; ...)

The Fourier transform of the Daubechies scaling function with filter `C` evaluated at `xi`.
`p` is the lower bound of `C`'s support minus 1, which by default is 0.

The function is computed as an 'infinite' product;
to control this there are optional arguments:

- `prec`: Include factors that are numerically smaller than 1-prec.
- `maxcount`: The maximum number of factors.
"""->
function FourDaubScaling( xi::Real, C::Vector{Float64}, p::Integer=0; prec=SMALL_EPS, maxcount=100)
	isapprox(sum(C), 1.0) || throw(AssertionError("Filter must sum to one"))
	prec >= SMALL_EPS || throw(DomainError())
	maxcount >= 1 || throw(DomainError())

	almost1 = 1.0 - prec
	Y = ComplexOne
	count = 1
	while count <= maxcount
		xi *= 0.5
		y = feval(xi, C, p)
		Y *= y

		# Convergence check: |y(xi) - 1| is small for small xi. But y is
		# exactly 1 in all even integers, so prevent premature exit
		if abs(xi) <= 0.5
			abs(y) >= almost1 && break
			count += 1
		end
	end

	return Y
end

@doc """
	FourScalingFunc( xi, wavename, J=0; ... )
	FourScalingFunc( xi, wavename, side, J=0; ... )

Compute the Fourier transform at `xi` of the (boundary) scaling function `wavename` with scale `J` (at boundary `side`, which is either `'L'` or `'R'`).
Optional arguments are passed to the Fourier transform of `wavename` (if relevant).

For the boundary functions an input vector `xi` of length `M` returns a matrix of size `M`-by-`p`, where `p` is the number of vanishing moments for `wavename`.

**Note**:
This function is intended to return the Fourier transform used in generalized sampling with reconstruction on [-1/2, 1/2].
This has the following consequences which is **not** shared with the lower lever functions for Fourier transforms.

- When `side` is `'L'` the function is translated with -1/2 (in the time domain) and when `side` is `'R'` the function is translated with 1/2.
- When `side` is `'L'` the *first* column of the output is related to the function closest to the left boundary, but when `side` is 'R'` the *last* column is related to the function closest to the right boundary.

The internal scaling functions are not phase shifted since this is handled by the NFFT package.
"""->
function FourScalingFunc( xi, wavename::AbstractString, J::Integer=0; args... )
	J >= 0 || throw(AssertionError("Scale must be a non-negative integer"))

	if ishaar(wavename)
		return FourHaarScaling(xi, J)
	elseif isdaubechies(wavename)
		return FourDaubScaling(xi, van_moment(wavename), J; args...)
	else
		error(string("Fourier transform for ", wavename, " is not implemented"))
	end
end

function FourScalingFunc( xi, wavename::AbstractString, side::Char, J::Integer=0; args... )
	J >= 0 || throw(AssertionError("Scale must be a non-negative integer"))

	if !isdaubechies(wavename)
		error(string("Fourier transform for boundary ", wavename, " is not implemented"))
	end

	Y = FourDaubScaling(xi, van_moment(wavename), side, J; args...).'

	if side == 'L'
		phase_shift = cis( pi*xi )
		broadcast!(*, Y, Y, phase_shift)
	elseif side == 'R'
		phase_shift = cis( -pi*xi )
		broadcast!(*, Y, Y, phase_shift)
		Y = flipdim(Y,2)
	end

	return Y
end

# ------------------------------------------------------------------------
# Fourier transform of Daubechies boundary wavelets

@doc """
	UVmat(F::BoundaryFilter) -> Matrix, Matrix

Return the matrices `U` and `V` used in `FourDaubScaling` for boundary scaling functions.
"""->
function UVmat(B::BoundaryFilter)
	vm = van_moment(B)

	UV = zeros(Complex{Float64}, vm, 3*vm-1)

	for i in 1:vm
		@inbounds UV[i, 1:vm+2*i-1] = bfilter(B, i-1) / sqrt2
	end

	return UV[:,1:vm], UV[:,vm+1:end]
end

@doc """
	FourDaubScaling( xi, p [, J], side::Char; args... ) -> Matrix

The Fourier transform at `xi` of the Daubechies boundary scaling function with `p` vanishing moments at scale `J` at the left or right `side` ('L' or 'R').

The output is an `p`-by-`length(xi)` matrix where column `j` holds the Fourier transform of the `p` boundary functions at `xi[j]`.
"""->
function FourDaubScaling{T<:Real}( xi::AbstractVector{T}, p::Integer, side::Char; maxcount=40 )
	1 <= maxcount <= 100 || throw(DomainError())
	#=
	- F is a vector with the Fourier transforms of the boundary scaling
	functions at the current xi
	- F0 is F at the zero vector
	- phihat is the Fourier transform of the internal scaling function
	at the current xi (with different phaseshifts)

	For large j,
	F(xi) \approx U^j*F0 + sum_{l=0}^{j-1} U^l*V*phihat( xi/2^{l+1} )
	where U & V are from UVmat
	=#

	U, V = UVmat(bfilter(p, side))
	C = coef(ifilter(p, true))
	scale!(C, 1/sum(C))

	F0 = complex( (eye(U) - U) \ vec(sum(V,2)) )

	# ----------------------------------------
	# Precompute U^l and U^l*V

	Upower = Dict{Int, Matrix{Complex{Float64}}}()
	sizehint!(Upower, maxcount)
	Upower[0] = eye(U)

	UpowerV = Dict{Int, Matrix{Complex{Float64}}}()
	sizehint!(UpowerV, maxcount)
	UpowerV[0] = V

	for l in 1:maxcount
		Upower[l] = U*Upower[l-1]
		UpowerV[l] = U*UpowerV[l-1]
	end

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

	Nxi = length(xi)
	Y = Array{Complex{Float64}}(p, Nxi)

	# Arrays with temporary results
	F = Array{Complex{Float64}}(p)
	phi_phase = collect(( side == 'L' ? -(p:1.0:3*p-2) : 1+(p:1.0:3*p-2) ))
	scale!(phi_phase, twoπ)
	Nphihat = length(phi_phase)
	phihat = Array{Complex{Float64}}(Nphihat)

	for nxi in 1:Nxi
		fill!(F, ComplexZero)

		current_xi = xi[nxi]
		for l in 0:maxcount-1
			current_xi *= 0.5
			internal = FourDaubScaling(current_xi, C, -p)
			for m in 1:Nphihat
				@inbounds phihat[m] = internal*cis( phi_phase[m]*current_xi )
			end

			# F = F + U^l*V*phihat
			BLAS.gemv!('N', ComplexOne, UpowerV[l], phihat, ComplexOne, F)

			# TODO: Test for convergence after each pass through the loop?
			# http://julialang.org/blog/2013/09/fast-numeric for |F-Fold|^2
		end
		# F = F + U^maxcount*F0
		BLAS.gemv!('N', ComplexOne, Upower[maxcount], F0, ComplexOne, F)

		Y[:,nxi] = F
	end

	return Y
end

# ------------------------------------------------------------
# Dilation and translation.

@vectorize_1arg Real FourHaarScaling

function FourHaarScaling(xi::Real, J::Integer)
	sqrt(2.0^(-J)) * FourHaarScaling(2.0^(-J)*xi)
end

function FourHaarScaling(xi::Real, J::Integer, k::Integer)
	cis( -twoπ*2.0^(-J)*k*xi ) * FourHaarScaling(xi, J)
end

function FourHaarScaling{T<:Real}( xi::AbstractArray{T}, J::Integer)
	dilation = 2.0^(-J)
	scale_factor = sqrt(dilation)
	Y = Array{Complex{Float64}}(size(xi))
	for idx in eachindex(xi)
		@inbounds Y[idx] = scale_factor * FourHaarScaling( dilation*xi[idx] )
	end

	return Y
end

function FourHaarScaling{T<:Real}( xi::AbstractArray{T}, J::Integer, k::Integer)
	dilation = 2.0^(-J)
	scale_factor = sqrt(dilation)
	Y = Array{Complex{Float64}}(size(xi))
	for idx in eachindex(xi)
		@inbounds Y[idx] = cis(-twoπ*dilation*k*xi[idx]) * scale_factor * FourHaarScaling( dilation*xi[idx] )
	end

	return Y
end

function FourDaubScaling( xi, p::Integer, J::Integer=0; offset::Integer=-p, args... )
	C = coef( ifilter(p) )
	scale!(C, 1/sum(C))
	FourDaubScaling( xi, C, J; offset=offset, args... )
end

function FourDaubScaling( xi, p::Integer, J::Integer=0, k::Integer=0; offset::Integer=-p, args... )
	C = coef( ifilter(p) )
	scale!(C, 1/sum(C))
	FourDaubScaling( xi, C, J, k; offset=offset, args... )
end

function FourDaubScaling( xi::AbstractArray{Float64}, C::Vector{Float64}; offset::Integer=-div(length(C),2), args... )
	Y = Array{Complex{Float64}}(size(xi))
	for idx in eachindex(xi)
		@inbounds Y[idx] = FourDaubScaling( xi[idx], C, offset; args... )
	end

	return Y
end

function FourDaubScaling{T<:Real}( xi::AbstractArray{T}, C::Vector{Float64}, J::Integer; offset::Integer=-div(length(C),2), args... )
	dilation = 2.0^(-J)
	scale_factor = sqrt(dilation)
	Y = Array{Complex{Float64}}(size(xi))
	for idx in eachindex(xi)
		@inbounds Y[idx] = scale_factor*FourDaubScaling( dilation*xi[idx], C, offset; args... )
	end

	return Y
end

function FourDaubScaling{T<:Real}( xi::AbstractArray{T}, C::Vector{Float64}, J::Integer, k::Integer; offset::Integer=-div(length(C),2), args... )
	dilation = 2.0^(-J)
	scale_factor = sqrt(dilation)
	Y = Array{Complex{Float64}}(size(xi))
	for idx in eachindex(xi)
		@inbounds Y[idx] = cis(-twoπ*dilation*k*xi[idx]) * scale_factor * FourDaubScaling( dilation*xi[idx], C, offset; args... )
	end

	return Y
end

function FourDaubScaling( xi, p::Integer, side::Char, J::Integer; args... )
	J >= 0 || throw(AssertionError("Scale must be a non-negative integer"))

	xi *= 2.0^(-J)
	y = FourDaubScaling(xi, p, side; args...)
	scale!( y, sqrt(2.0^(-J)) )

	return y
end


text:
export install_geodynamics
import ADCME: get_gfortran, require_file, GFORTRAN, http_file, change_directory
function install_geodynamics()
    PWD = pwd()
    get_gfortran()
    GEODYNAMICS = joinpath(ADCME.PREFIXDIR, "geodynamics")
    change_directory(GEODYNAMICS)
    require_file("seismic_CPML_2D_pressure_second_order.f90") do
        http_file("https://raw.githubusercontent.com/kailaix/seismic_cpml/master/seismic_CPML_2D_pressure_second_order.f90", "seismic_CPML_2D_pressure_second_order.f90")
    end
    require_file("seismic_CPML_2D_isotropic_fourth_order.f90") do
        http_file("https://raw.githubusercontent.com/kailaix/seismic_cpml/master/seismic_CPML_2D_isotropic_fourth_order.f90", "seismic_CPML_2D_isotropic_fourth_order.f90")
    end
    # compiling
    run(`$GFORTRAN -std=gnu -fimplicit-none -frange-check -O3 -fmax-errors=10 -pedantic -pedantic-errors -Waliasing -Wampersand -Wcharacter-truncation -Wline-truncation -Wsurprising -Wno-tabs -Wunderflow
            -mcmodel=medium seismic_CPML_2D_pressure_second_order.f90  -o acoustic`)
    run(`$GFORTRAN -std=gnu -fimplicit-none -frange-check -O3 -fmax-errors=10 -pedantic -pedantic-errors -Waliasing -Wampersand -Wcharacter-truncation -Wline-truncation -Wsurprising -Wno-tabs -Wunderflow
            -mcmodel=medium seismic_CPML_2D_isotropic_fourth_order.f90 -o elastic`)
    exec1 = joinpath(pwd(), "acoustic"), joinpath(pwd(), "elastic")
    cd(PWD)
    exec1
end


text:
"""
    pprint_struct([io::IO=stdout], [mimetype], x; kw...)

Pretty print object `x` that is a struct type (`isstructtype` returns `true`).

# Keyword Arguments

See [`pprint`](@ref), they share the same keyword arguments.
"""
function pprint_struct end

pprint_struct(@nospecialize(x); kw...) = pprint_struct(stdout, x; kw...)

function pprint_struct(io::IO, @nospecialize(x);kw...)
    pprint_struct(GarishIO(io; kw...), x)
end

function pprint_struct(io::GarishIO, @nospecialize(x))
    pprint_struct(io, MIME"text/plain"(), x)
end

"""
    pprint_struct(io, ::MIME"text/plain", @nospecialize(x))

Print `x` as a struct type with mime type `MIME"text/plain"`.
"""
function pprint_struct(io::GarishIO, mime::MIME"text/plain", @nospecialize(x))
    t = typeof(x)
    isstructtype(t) || throw(ArgumentError("expect a struct type, got $t"))
    pprint_struct_type(io, mime, t); print(io.bland_io, "(")

    nf = nfields(x)::Int
    nf == 0 && return print(io.bland_io, ")")

    io.compact || println(io.bland_io)

    # findout fields to print
    fields_to_print = Int[]
    for i in 1:nf
        f = fieldname(t, i)
        value = getfield(x, i)
        if !io.include_defaults && is_option(x) && value == field_default(t, f)
        else
            push!(fields_to_print, i)
        end
    end

    within_nextlevel(io) do
        for i in fields_to_print
            f = fieldname(t, i)
            value = getfield(x, i)
            print_indent(io)
            print_token(io, :fieldname, f)
            if io.compact
                print(io.bland_io, "=")
            else
                print(io.bland_io, " = ")
                io.state.offset = 3 + length(string(f))
            end

            if !isdefined(x, f) # print undef as comment color
                pprint(io, undef)
            else
                new_io = GarishIO(io; limit=true)
                pprint_field(new_io, mime, value)
            end

            if !io.compact || i != last(fields_to_print)
                print(io.bland_io, ", ")
            end

            if i != last(fields_to_print)
                io.compact || println(io.bland_io)
            end
        end
    end
    io.compact || println(io.bland_io)
    print_indent(io)
    print(io.bland_io, ")")
end

function pprint_struct_type(io::GarishIO, mime::MIME"text/plain", @nospecialize(x))
    if x isa DataType
        # max_indent_reached = io.indent * io.state.level + io.state.offset + length(string(t)) + 2 > io.displaysize[2]
        if max_indent_reached(io, length(string(x)) + 2)
            if isempty(x.parameters)
                # nothing we can do here, let it break line
                print_token(io, :type, x)
            else
                print_token(io, :type, nameof(x)); print_token(io, :type, "{")
                println(io.bland_io)
                within_nextlevel(io) do
                    len_params = sum(length(string(p)) for p in x.parameters) + 2 * length(x.parameters)
                    if max_indent_reached(io, len_params)
                        for p in x.parameters
                            print_indent(io)
                            pprint_struct_type(io, mime, p)
                            print(io.bland_io, ", ")
                            println(io.bland_io)
                            # if p !== last(x.parameters)
                            # end
                        end
                    else
                        # try to print parameters in one line
                        print_indent(io)
                        print_token(io, :type, join(x.parameters, ", "))
                        println(io.bland_io)
                    end
                end
                print_indent(io)
                print_token(io, :type, "}")
            end
        else
            print_token(io, :type, x)
        end
    else
        print_token(io, :type, x)
    end
end

pprint_field(io::GarishIO, x) = pprint_field(io, MIME"text/plain"(), x)

function pprint_field(io::GarishIO, mime::MIME"text/plain", x)
    upperlevel_type = io.state.type
    upperlevel_noindent_in_first_line = io.state.noindent_in_first_line
    io.state.type = StructField
    io.state.noindent_in_first_line = true
    pprint(io, mime, x)
    io.state.noindent_in_first_line = upperlevel_noindent_in_first_line
    io.state.type = upperlevel_type
end


text:
## Benders Decomposition
## Edward J. Xu <edxu96@outlook.com>
## April 25th, 2019

module BendersDecomp

include("./milp.jl")
include("./l-shaped.jl")
include("./slave.jl")

"""
Make sure "result_q == obj_sub" in the final iteration,
while ((bu - bl > epsilon) && (ite <= timesIterationMax))  !!!
"""
function check_whe_continue(bu, bl, epsilon, result_q, obj_sub,
    ite, timesIterationMax)
  whe_continue = true
  if (bu - bl <= epsilon) & (bu - bl >= - epsilon)
    if (result_q - obj_sub <= epsilon) & (result_q - obj_sub >= - epsilon)
      whe_continue = false
    end
  end

  return (ite <= timesIterationMax) && whe_continue
end

"JuMP Model of Master Problem"
mutable struct ModMas
  expr::Model
  vec_y

  function ModMas(n_y, vec_min_y, vec_max_y)
    expr = Model(with_optimizer(GLPK.Optimizer))
    @variable(expr, q)
    @variable(expr, vec_y[1: n_y], Int)
    @objective(expr, Min, (transpose(vec_f) * vec_y)[1] + q)
    @constraint(expr, vec_y[1: n_y] .<= vec_max_y)
    @constraint(expr, vec_y[1: n_y] .>= vec_min_y)
    new(expr, vec_y)
  end
end

"Results of master problem."
struct LogMas
  obj
  vec_y_result
end

"Solve master problem."
function solve_mas!(mas::ModMas, q, vec_y, vec_b, mat_b)
  if bool_solution_sub
    @constraint(mas.mod, (transpose(vec_uBar) * (vec_b - mat_b * vec_y))[1] <= q)
  else  # Add feasible cut Constraints
    @constraint(mas.mod, (transpose(vec_uBar) * (vec_b - mat_b * vec_y))[1] <= q)
    @constraint(mas.mod, (transpose(vec_uBar) * (vec_b - mat_b * vec_y))[1] <= 0)
  end

  optimize!(mas.mod)
  log = LogMas(obj = objective_value(mas), vec_y = value_vec(vec_y))

  return log
end

"Linked lists for optimization results."
mutable struct Log
  ite::Int64
  log_slave::Union{LogRay,LogSub,Missing}
  log_mas::Union{LogMas,Missing}
  bu::Float64
  bl::Float64
  next::Log

  function Log(ite, bu=missing, bl=missing)
    new(missing, missing, bu, bl, missing)
  end
end

"Head of Linked List for Log"
struct LogHead
  next::Log

  function LogHead(n_cons, vec_y_bar, vec_b, mat_b, mat_a, vec_c)
    vec_y_bar = zeros(mod.n_y, 1)
    vec_result_x = zeros(mod.n_x, 1)
    new(
      next=Log(ite=0, bu=Inf, bl=- Inf)
      )
  end
end

"Append a new log after this log."
function append_log!(log::Log)
  log.next = Log(log.ite + 1)
  log.next.bu = log.bu
  log.next.bl = log.bl
end

"Find latest log."
function find_log_lastest(logs::LogHead)
  crt = LogHead
  while !(isequal(missing, crt.next))
    crt = crt.next
  end

  return crt
end

export ModMilpBenders, solveModMilpBenders!

end


text:
module RocketStorageActorTest

using Test
using Rocket

@testset "StorageActor" begin

    println("Testing: actor StorageActor")

    @testset begin
        source = from([ 1, 2, 3 ])
        actor  = StorageActor{Int}()

        @test actor.value === nothing

        subscribe!(source, actor)

        @test actor.value == 3
        @test getvalues(actor) == 3
    end

    @testset begin
        source = faulted(Int, "Error")
        actor  = StorageActor{Int}()

        @test actor.value === nothing
        @test_throws ErrorException subscribe!(source, actor)
        @test actor.value === nothing
        @test getvalues(actor) === nothing
    end

    @testset begin
        @test storage(Int) isa StorageActor{Int}
        @test storage(Real) isa StorageActor{Real}
    end
end

end


text:
module EvolutionaryStrategies

using Cambrian
using Random
using LinearAlgebra

abstract type ESState end

include("individual.jl")
include("exponential_nes.jl")
include("separable_nes.jl")
include("cmaes.jl")
include("random_search.jl")

end


text:
module TestMOIWrapper

using MadNLP
using Test

const MOI = MadNLP.MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

function test_MOI_Test()
    model = MOI.Bridges.full_bridge_optimizer(
        MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            MadNLP.Optimizer(),
        ),
        Float64,
    )
    MOI.set(model, MOI.Silent(), true)
    MOI.Test.runtests(
        model,
        MOI.Test.Config(
            atol = 1e-4,
            rtol = 1e-4,
            optimal_status = MOI.LOCALLY_SOLVED,
            exclude = Any[
                MOI.delete,
                MOI.ConstraintDual,
                MOI.ConstraintBasisStatus,
                MOI.DualObjectiveValue,
                MOI.ObjectiveBound,
            ]
        );
        exclude = String[
            "test_modification",
            # - Need to implement TimeLimitSec
            "test_attribute_TimeLimitSec",
            # - Wrong return type
            "test_model_UpperBoundAlreadySet",
            # - Final objective value is not equal to 0.0
            "test_objective_FEASIBILITY_SENSE_clears_objective",

            # TODO: Need to investigate why these tests are breaking
            #   get(model, MOI.ConstraintPrimal(), c) returns the
            #   opposite value: if 1.0 is expected, -1.0 is returned
            "test_constraint_ScalarAffineFunction_EqualTo",
            "test_quadratic_nonconvex_constraint_basic",
            "test_linear_integration",

            # TODO: there might be an issue with VectorAffineFunction/VectorOfVariables
            "test_conic_NormOneCone_VectorOfVariables",
            "test_conic_NormOneCone_VectorAffineFunction",
            "test_conic_NormInfinityCone_VectorOfVariables",
            "test_conic_NormInfinityCone_VectorAffineFunction",
            "test_conic_linear_VectorAffineFunction",
            "test_conic_linear_VectorOfVariables",

            # Tests excluded on purpose
            # - Excluded as MadNLP returns LOCALLY_INFEASIBLE instead of INFEASIBLE
            "INFEASIBLE",
            "test_solve_DualStatus_INFEASIBILITY_CERTIFICATE_",
            # - Excluded because Hessian information is needed
            "test_nonlinear_hs071_hessian_vector_product",
            # - Excluded because Hessian information is needed
            "test_nonlinear_hs071_no_hessian",
            # - Excluded because Hessian information is needed
            "test_nonlinear_invalid",

            #  - Excluded because this test is optional
            "test_model_ScalarFunctionConstantNotZero",
            #  - Excluded because MadNLP returns INVALID_MODEL instead of LOCALLY_SOLVED
            "test_linear_VectorAffineFunction_empty_row",
        ]
    )
    return
end

end

TestMOIWrapper.runtests()


text:
include("../src/linearMeshIO.jl")

gmsh_file = joinpath(pwd(), "tests/meshes/gmsh/t1.msh")
# gmsh_file = joinpath(pwd(), "tests/meshes/gmsh/t2.msh")
# gmsh_file = joinpath(pwd(), "tests/meshes/gmsh/t3.msh")

gmsh = linearMeshIO.GMSHreader.readGMSH(gmsh_file);

max_stars_repo_name: c-schubert/linearMeshIO

max_issues_repo_name: c-schubert/linearMeshIO

max_forks_repo_name: c-schubert/linearMeshIO


text:
using Dolo
import Dolo: Euler
model = Model("examples/models/rbc.yaml")

F = Euler(model);

function L_p(F, x0, x1, dp; ε=1e-8)

    p0 = F.p

    dp = SVector(dp...)

    u = maximum(abs, dp)

    dp_ε = dp/u*1e-8

    f0 = F(x0, x1)

    F.p = p0 + dp_ε

    f0_ϵ = F(x0, x1)

    F.p = p0

    f = (f0_ε-f0)/ε*u

    return f

end


text:
using Test
using RootSolvers
using StaticArrays

struct RootSolvingProblem{S,F,F′,FT,FTA,N}
  "Name"
  name::S
  "Root equation"
  f::F
  "Derivative of f"
  f′::F′
  "Exact solution"
  x̃::FT
  "Initial guess, for iterative methods"
  x_init::FTA
  "Lower bound for bracketed methods"
  x_lower::FTA
  "Upper bound for bracketed methods"
  x_upper::FTA
end
function RootSolvingProblem(args...)
  S = typeof(args[1])
  F = typeof(args[2])
  F′ = typeof(args[3])
  FT = eltype(args[end])
  FTA = typeof(args[end])
  N = length(args[end])
  RootSolvingProblem{S,F,F′,FT,FTA,N}(args...)
end

problem_size() = 5
problem_size(::RootSolvingProblem{S,F,F′,FT,FTA,N}) where {S,F,F′,FT,FTA,N} = N
float_types() = [Float32, Float64]
get_tolerances(FT) = [ResidualTolerance{FT}(1e-6), SolutionTolerance{FT}(1e-3), nothing]

test_verbose!(::CompactSolution, sol, problem, tol, converged) = nothing
function test_verbose!(::VerboseSolution, sol, problem, tol, converged)
  # Not sure what else to test here
  sol isa AbstractArray || @test length(sol.root_history) == length(sol.err_history)

  if converged
    if tol isa ResidualTolerance
      if sol isa AbstractArray
        @test all(map(x -> problem.f(x.root), sol) .< tol.tol)
      else
        @test problem.f(sol.root) < tol.tol
      end
    end
  end
end

"""
    expand_data_inputs!(problem_list, N)

Expands data inputs in `problem_list` from scalars to SArray's
"""
function expand_data_inputs!(problem_list, ε, N)
  for problem in deepcopy(problem_list)
    FT = typeof(problem.x̃)
    push!(problem_list,
      RootSolvingProblem(
        problem.name,
        problem.f,
        problem.f′,
        problem.x̃,
        SArray{Tuple{N, N}, FT}(problem.x_init  .+ FT(ε)*rand(FT, N, N)),
        SArray{Tuple{N, N}, FT}(problem.x_lower .+ FT(ε)*rand(FT, N, N)),
        SArray{Tuple{N, N}, FT}(problem.x_upper .+ FT(ε)*rand(FT, N, N))
        ))
  end
end

#####
##### Construct problem list
#####

"""
    get_methods(x_init, x_lower, x_upper, f′)

Numerical methods to test, given arguments from `RootSolvingProblem`.
"""
function get_methods(x_init, x_lower, x_upper, f′)
    return (
    SecantMethod(x_lower, x_upper),
    RegulaFalsiMethod(x_lower, x_upper),
    NewtonsMethodAD(x_init),
    NewtonsMethod(x_init, f′)
    )
end

# Convenience types for dispatching
# (since instances of SecantMethod
# are not `isbits` with `CuArray`s).
struct SecantMethodType end
struct RegulaFalsiMethodType end
struct NewtonsMethodADType end
struct NewtonsMethodType end

# Convenience methods for unifying interfaces
# in test suite:
get_method(::SecantMethodType, x_init, x_lower, x_upper, f′) = SecantMethod(x_lower, x_upper)
get_method(::RegulaFalsiMethodType, x_init, x_lower, x_upper, f′) = RegulaFalsiMethod(x_lower, x_upper)
get_method(::NewtonsMethodADType, x_init, x_lower, x_upper, f′) = NewtonsMethodAD(x_init)
get_method(::NewtonsMethodType, x_init, x_lower, x_upper, f′) = NewtonsMethod(x_init, f′)

#####
##### Construct problem list
#####

problem_list = RootSolvingProblem[]

for FT in float_types()
  push!(problem_list, RootSolvingProblem(
    "simple quadratic",
    x -> x^2 - 100^2,
    x -> 2x,
    FT(100),
    FT(1),
    FT(0),
    FT(1000),
    ))
end

expand_data_inputs!(problem_list, 0.1, problem_size())