text:
#############################################################################
#Copyright (c) 2016-2019 Charles Le Losq
#
#The MIT License (MIT)
#
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the #Software without restriction, including without limitation the rights to use, copy, #modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, #and to permit persons to whom the Software is furnished to do so, subject to the #following conditions:
#
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#
# functions.jl contains several mathematic functions
#
#
#############################################################################

"""
	bootsample(x::Array{Float64}, y::Array{Float64}; boottype::String = "np", ese::Array{Float64} = [0.0])

Inputs
------

	x: Array{Float64}, the x axis. It can have multiple columns.

	y: Array{Float64}, the y axis. It can have multiple columns.

Options
-------

	boottype: ASCIIString, either "np" or "p", this is the type of bootstrapping performed. "np" performes a non-parametric resampling fo the dataset with replacement. "p" performs a parametric resampling. The data are resample from a gaussian distribution centered on the y values with errors that should be provided in the ese variable.

	ese: Array{Float64}, containing the errors affecting the y values that are used during parametric bootstrapping.

Outputs
-------

	b_x_f: Array{Float64}, the bootstrapped x values

	b_y_f: Array{Float64}, the bootstrapped y values

The bootstrap function can be embedded in a for loop, and will each time produce a different dataset. Performing K times the bootstrapping and fitting each time the model will allow to estimate the error distribution on the peak parameters. This technic has the advantage of making no prior assumption on the probability distribution functions of parameters errors. However, it is  much more time consuming that using the covariance matrix.

"""
function bootsample(x::Array{Float64}, y::Array{Float64}; boottype::String = "np", ese::Array{Float64} = [0.0])
    testx = size(x,1)
    testy = size(y,1)
    teste = size(ese,1)
    if boottype == "np" && testx == testy
        vect = collect(1:size(x)[1]) # for real bootstrapping
        idx = sample(vect,size(vect)[1],replace=true) #resampling data with replacement...
        b_x_f = x[idx,:] # we pick the right x
        b_y_f = y[idx,:] # we pick the right y
    elseif boottype == "p" && testx == testy && testx == teste
	    b_x_f = x
	    b_y_f = ones(size(y))
	    b_y_f = randn!(b_y_f[:,:]) .* ese[:,:] + y[:,:]
    else
	    error("Something is wrong. Check size of x, y and ese as well as the boottype (either p or np). Providing an ese array is mandatory for parametric bootstrapping")
    end

    return b_x_f, b_y_f
end

"""
	bootperf(params_boot::Array{Float64}; plotting::String = "True", parameter::Int64 = 1, feature::Int64 = 1, histogram_step::Int64 = 100, savefigures::String = "False", save_bootrecord::String = "Boot_record.pdf", save_histogram::String = "Boot_histogram.pdf")

Inputs
------

	params_boot[i,j] or [i,j,k]
		array with i the bootstrrap experiment, j the parameter and k the feature being calculated (e.g., a peak during peak fitting)
	plotting: String
		switch to plotting mode ("True") or not ("False"). If true, parameter and feature must be provided, otherwise an error message is returned.

	histogram_step
		integer value to control the histogram X axis division.
	savefigures: "True" or "False"
		explicit, save in the current working directory.
	save_bootrecord:
		Name for the graphic showing the bootstrap mean and std evolutions, with extension.
	save_histogram:
		Name for the graphic showing the histogram for the parameter and feature of interest, with extension.

Outputs
-------

	std_record, mean_record, the arrays recording how the standard deviation and mean of the parameters as a function of the bootstrap advance.

"""
function bootperf(params_boot::Array{Float64}; plotting::String = "True", parameter::Int64 = 1, feature::Int64 = 1, histogram_step::Int64 = 100, savefigures::String = "False", save_bootrecord::String = "Boot_record.pdf", save_histogram::String = "Boot_histogram.pdf")

    # test
    if plotting != "True" && plotting != "False"
        error("Error: plotting should be set to True or False")
    end

    nb_boot::Int64 = size(params_boot)[1]
    nb_param::Int64 = size(params_boot)[2]

   if ndims(params_boot) == 2 #### TWO DIMENSIONS PARAMS_BOOT

       mean_record::Array{Float64} = zeros(size(params_boot))
       std_record::Array{Float64} = zeros(size(params_boot))
       for k = 1:nb_boot
           for j = 1:size(params_boot)[2]
       	        if k == 1
       	            mean_record[k,j] = params_boot[k,j]
       		    std_record[k,j] = 0
       	        else
       	            mean_record[k,j] = sum(mean(params_boot[1:k,j],1))
       		    std_record[k,j] = sum(std(params_boot[1:k,j],1))
       	        end
       	    end
       	end

       if plotting == "True" # Graphics
           # tests:
           if savefigures != "True" && savefigures != "False"
               error("savefigure should be set to True or False.")
           end
           if isa(save_bootrecord ,String) != true
               error("Error: save_bootrecord is not a valid savename. Check it is a string of characters.")
           end
           if isa(save_histogram ,String) != true
               error("Error: save_histogram is not a valid savename. Check it is a string of characters.")
           end

           if parameter == 0 || histogram_step == 0
               error("parameter or histogram_step appear to have been set to 0...")
           elseif parameter > size(params_boot)[2]
               error("Value entered for parameter is too high.")
           end

           fig = figure()
           #title("Mean and Std of parameter $(parameter), feature $(feature)",fontsize = 10, fontweight = "bold", fontname="Arial")
           subplot(121)
           plot(1:nb_boot, mean_record[:,parameter],label="Mean")
           xlabel("Number of iterations during bootstrap",fontsize = 10, fontweight = "bold", fontname="Arial")
           ylabel("Mean value",fontsize = 10, fontweight = "bold", fontname="Arial")

           subplot(122)
           plot(1:nb_boot, std_record[:,parameter],label="Std")
           xlabel("Number of iterations during bootstrap",fontsize = 10, fontweight = "bold", fontname="Arial")
           ylabel("Standard deviation",fontsize = 10, fontweight = "bold", fontname="Arial")

           if savefigures == "True"
               savefig(save_bootrecord)
           end

           fig2= figure()
           h = plt[:hist](params_boot[:,parameter],100)
           #title("Histogram of parameter $(parameter), feature $(feature)",fontsize = 10, fontweight = "bold", fontname="Arial")
           xlabel("Parameter value",fontsize = 10, fontweight = "bold", fontname="Arial")
           ylabel("Number of bootstrap values",fontsize = 10, fontweight = "bold", fontname="Arial")

           if savefigures == "True"
               savefig(save_histogram)
       	   end

       end # end if plotting

       return mean_record, std_record

   elseif ndims(params_boot) == 3 #### THREE DIMENSIONS PARAMS_BOOT
       mean_record = zeros(size(params_boot))
       std_record = zeros(size(params_boot))
       for k = 1:nb_boot
           for j = 1:size(params_boot)[2]
   	    for i = 1:size(params_boot)[3]
   	        if k == 1
   	            mean_record[k,j,i] = params_boot[k,j,i]
   		    std_record[k,j,i] = 0
   	        else
   	            mean_record[k,j,i] = sum(mean(params_boot[1:k,j,i],1))
   		    std_record[k,j,i] = sum(std(params_boot[1:k,j,i],1))
   	        end
   	    end
   	end
       end

       if plotting == "True" # Graphics
           # tests:
           if savefigures != "True" && savefigures != "False"
               error("savefigure should be set to True or False.")
           end
           if isa(save_bootrecord ,String) != true
               error("Error: save_bootrecord is not a valid savename. Check it is a string of characters.")
           end
           if isa(save_bootrecord ,String) != true
               error("Error: save_histogram is not a valid savename. Check it is a string of characters.")
           end

           if parameter == 0 || feature == 0 || histogram_step == 0
               error("parameter, feature or histogram_step appear to have been set to 0...")
           elseif parameter > size(params_boot)[2]
               error("Value entered for parameter is too high.")
           elseif feature > size(params_boot)[3]
   	    error("Value entered for feature is too high.")
           end

           fig = figure()
           #title("Mean and Std of parameter $(parameter), feature $(feature)",fontsize = 10, fontweight = "bold", fontname="Arial")
           subplot(121)
           plot(1:nb_boot, mean_record[:,parameter,feature],label="Mean")
           xlabel("Number of iterations during bootstrap",fontsize = 10, fontweight = "bold", fontname="Arial")
           ylabel("Mean value",fontsize = 10, fontweight = "bold", fontname="Arial")

           subplot(122)
           plot(1:nb_boot, std_record[:,parameter,feature],label="Std")
           xlabel("Number of iterations during bootstrap",fontsize = 10, fontweight = "bold", fontname="Arial")
           ylabel("Standard deviation",fontsize = 10, fontweight = "bold", fontname="Arial")

           if savefigures == "True"
               savefig(save_bootrecord)
           end

           fig2= figure()
           h = plt[:hist](params_boot[:,parameter,feature],100)
           #title("Histogram of parameter $(parameter), feature $(feature)",fontsize = 10, fontweight = "bold", fontname="Arial")
           xlabel("Parameter value",fontsize = 10, fontweight = "bold", fontname="Arial")
           ylabel("Number of bootstrap values",fontsize = 10, fontweight = "bold", fontname="Arial")

           if savefigures == "True"
               savefig(save_histogram)
   	end

       end # end if plotting

       return mean_record, std_record

    elseif ndims(params_boot) > 3
	error("Not implemented, params_boot should have 3 dimensions, maximum.")
    end

end # end function


text:
"""
Launch IPython in Julia

## Usage

Run `using IPython` and then type `.` in empty `julia>` prompt or run
`IPython.start_ipython()`.  If you are using IPython 7.0 or above, you
can switch back to Julia REPL by `backspace` or `ctrl-h` key (like
other REPL modes).  For older versions of IPython, exiting IPython as
usual (e.g., `ctrl-d`) brings you back to the Julia REPL.  Re-entering
IPython keeps the previous state.  Use pre-defined `Main` object to
access Julia namespace from IPython.  Use `py"..."` string macro to
access Python namespace from Julia.

**Note:**
First launch of IPython may be slow.

# Examples
```
julia> xs = [1, 2, 3];  # some Julia variable

julia> using IPython

julia> IPython.start_ipython()  # or type . (a dot)

In [1]: Main.xs  # accessing Julia variable
array([1, 2, 3], dtype=int64)

In [2]: Main.xs += 10

In [3]: ε = "varepsilon"

In [4]: exit()  # or type backspace or ctrl-d

julia> xs
3-element Array{Int64,1}:
 11
 12
 13

julia> using PyCall

julia> py"ε"  # you can access Python variables with PyCall.@py_str macro
"varepsilon"

julia> py"ϵ"  # but be careful when using Unicode!
"varepsilon"
```
"""
module IPython

using Compat
using Compat: @warn, @info
include("julia_api.jl")
include("core.jl")
include("convenience.jl")
include("julia_repl.jl")

end # module

ext: jl


text:
"""
    get_observable_timers() :: Nothing

Test performance of current observable implementation by Fourier transforming test correlations
on a few 2D and 3D lattices with 50 x 50 x 0 (2D) or 50 x 50 x 50 (3D) momenta.
"""
function get_observable_timers() :: Nothing

    # init test dummys
    l1   = get_lattice("square",      6, verbose = false); r1 = get_reduced_lattice("heisenberg", [[0.0]], l1, verbose = false)
    l2   = get_lattice("cubic",       6, verbose = false); r2 = get_reduced_lattice("heisenberg", [[0.0]], l2, verbose = false)
    l3   = get_lattice("kagome",      6, verbose = false); r3 = get_reduced_lattice("heisenberg", [[0.0]], l3, verbose = false)
    l4   = get_lattice("hyperkagome", 6, verbose = false); r4 = get_reduced_lattice("heisenberg", [[0.0]], l4, verbose = false)
    χ1   = Float64[exp(-norm(r1.sites[i].vec)) for i in eachindex(r1.sites)]
    χ2   = Float64[exp(-norm(r2.sites[i].vec)) for i in eachindex(r2.sites)]
    χ3   = Float64[exp(-norm(r3.sites[i].vec)) for i in eachindex(r3.sites)]
    χ4   = Float64[exp(-norm(r4.sites[i].vec)) for i in eachindex(r4.sites)]
    k_2D = get_momenta((0.0, 1.0 * pi), (0.0, 1.0 * pi), (0.0, 0.0), (50, 50, 0))
    k_3D = get_momenta((0.0, 1.0 * pi), (0.0, 1.0 * pi), (0.0, 1.0 * pi), (50, 50, 50))

    # init timer
    to = TimerOutput()

    @timeit to "=> Fourier transform" begin
        for rep in 1 : 10
            @timeit to "-> square"      compute_structure_factor(χ1, k_2D, l1, r1)
            @timeit to "-> cubic"       compute_structure_factor(χ2, k_3D, l2, r2)
            @timeit to "-> kagome"      compute_structure_factor(χ3, k_2D, l3, r3)
            @timeit to "-> hyperkagome" compute_structure_factor(χ4, k_3D, l4, r4)
        end
    end

    show(to)

    return nothing
end


text:
function assess(α::T, β::T, x::T, tol::T) where T <: AbstractFloat
    Nopt = ceil(Integer, x^(1/α))
    term = NaN
    Nmin = -1
    for n = 1:Nopt
        term = exp( loggamma(n*α-β+1) - n*log(x) )
        if term < tol
            Nmin = n
            break
        end
    end
    return Nmin, term
end


text:
"""
Provides Genie with encryption and decryption capabilities.
"""
module Encryption

import Genie, Nettle

const ENCRYPTION_METHOD = "AES256"

"""
    encrypt{T}(s::T) :: String

Encrypts `s`.
"""
function encrypt(s::T)::String where T
  (key32, iv16) = encryption_sauce()
  encryptor = Nettle.Encryptor(ENCRYPTION_METHOD, key32)

  Nettle.encrypt(encryptor, :CBC, iv16, Nettle.add_padding_PKCS5(Vector{UInt8}(s), 16)) |> bytes2hex
end

"""
    decrypt(s::String) :: String

Decrypts `s` (a `string` previously encrypted by Genie).
"""
function decrypt(s::String) :: String
  (key32, iv16) = encryption_sauce()
  decryptor = Nettle.Decryptor(ENCRYPTION_METHOD, key32)
  deciphertext = Nettle.decrypt(decryptor, :CBC, iv16, s |> hex2bytes)

  String(Nettle.trim_padding_PKCS5(deciphertext))
end

"""
    encryption_sauce() :: Tuple{Vector{UInt8},Vector{UInt8}}

Generates a pair of key32 and iv16 with salt for encryption/decryption
"""
function encryption_sauce() :: Tuple{Vector{UInt8},Vector{UInt8}}
  if length(Genie.secret_token()) < 64
    if !Genie.Configuration.isprod()
      @error "Invalid Genie.secret_token() with less than 64 characters; using a temporary token"
      Genie.secret_token!()
    else
      error("Can't encrypt - make sure that secret_token!(token) is called in config/secrets.jl")
    end
  end
  token = Genie.secret_token()
  passwd = token[1:32]
  salt = hex2bytes(token[33:64])
  Nettle.gen_key32_iv16(Vector{UInt8}(passwd), salt)
end

end


text:
#################### DistributionStruct ####################

#################### Base Methods ####################

dims(d::DistributionStruct) = size(d)

function dims(D::Array{MultivariateDistribution})
  size(D)..., mapreduce(length, max, D)
end

#################### List Fallbacks ####################

unlist(d::Distribution, x::AbstractArray) = vec(x)

unlist_sub(d::Distribution, x::AbstractArray) = unlist(d, x)

unlist_sub(d::UnivariateDistribution, X::AbstractArray) = vec(X)

unlist_sub(D::Array{UnivariateDistribution}, X::AbstractArray) = vec(X)

function unlist_sub(D::Array{MultivariateDistribution}, X::AbstractArray)
  y = similar(X, length(X))
  offset = 0
  for sub in CartesianRange(size(D))
    n = length(D[sub])
    inds = 1:n
    y[offset + inds] = X[sub, inds]
    offset += n
  end
  resize!(y, offset)
end

function relistlength(d::UnivariateDistribution, x::AbstractArray)
  (x[1], 1)
end

function relistlength(d::MultivariateDistribution, x::AbstractArray)
  n = length(d)
  value = x[1:n]
  (value, n)
end

function relistlength(d::MatrixDistribution, x::AbstractArray)
  n = length(d)
  value = reshape(x[1:n], size(d))
  (value, n)
end

relistlength_sub(d::Distribution, s::AbstractStochastic, x::AbstractArray) =
  relistlength(d, x)

function relistlength_sub(d::UnivariateDistribution, s::ArrayStochastic,
                          X::AbstractArray)
  n = length(s)
  value = reshape(X[1:n], size(s))
  (value, n)
end

function relistlength_sub(D::Array{UnivariateDistribution}, s::ArrayStochastic,
                          X::AbstractArray)
  n = length(s)
  value = reshape(X[1:n], size(s))
  (value, n)
end

function relistlength_sub(D::Array{MultivariateDistribution},
                          s::ArrayStochastic, X::AbstractArray)
  Y = similar(X, size(s))
  offset = 0
  for sub in CartesianRange(size(D))
    n = length(D[sub])
    inds = 1:n
    Y[sub, inds] = X[offset + inds]
    offset += n
  end
  (Y, offset)
end

#################### Link Fallbacks ####################

link(d::Distribution, x) = x

link_sub(d::Distribution, x) = link(d, x)

function link_sub(d::UnivariateDistribution, X::AbstractArray)
  Y = similar(X, Float64)
  map!(x -> link_sub(d, x), Y, X)
end

function link_sub(D::Array{UnivariateDistribution}, X::AbstractArray)
  Y = similar(X, Float64)
  map!(i -> link_sub(D[i], X[i]), Y, 1:length(D))
end

function link_sub(D::Array{MultivariateDistribution}, X::AbstractArray)
  Y = similar(X, Float64)
  for sub in CartesianRange(size(D))
    d = D[sub]
    inds = 1:length(d)
    Y[sub, inds] = link_sub(d, X[sub, inds])
  end
  Y
end

invlink(d::Distribution, x) = x

invlink_sub(d::Distribution, x) = invlink(d, x)

function invlink_sub(d::UnivariateDistribution, X::AbstractArray)
  Y = similar(X, Float64)
  map!(x -> invlink_sub(d, x), Y, X)
end

function invlink_sub(D::Array{UnivariateDistribution}, X::AbstractArray)
  Y = similar(X, Float64)
  map!(i -> invlink_sub(D[i], X[i]), Y, 1:length(D))
end

function invlink_sub(D::Array{MultivariateDistribution}, X::AbstractArray)
  Y = similar(X, Float64)
  for sub in CartesianRange(size(D))
    d = D[sub]
    inds = 1:length(d)
    Y[sub, inds] = invlink_sub(d, X[sub, inds])
  end
  Y
end

#################### Logpdf Fallbacks ####################

logpdf(d::Distribution, x, transform::Bool) = logpdf(d, x)

function logpdf_sub(d::Distribution, x, transform::Bool)
  insupport(d, x) ? logpdf(d, x, transform) : -Inf
end

function logpdf_sub(d::UnivariateDistribution, X::AbstractArray,
                    transform::Bool)
  lp = 0.0
  for x in X
    lp += logpdf_sub(d, x, transform)
  end
  lp
end

function logpdf_sub(D::Array{UnivariateDistribution}, X::AbstractArray,
                    transform::Bool)
  lp = 0.0
  for i in 1:length(D)
    lp += logpdf_sub(D[i], X[i], transform)
  end
  lp
end

function logpdf_sub(D::Array{MultivariateDistribution}, X::AbstractArray,
                    transform::Bool)
  lp = 0.0
  for sub in CartesianRange(size(D))
    d = D[sub]
    lp += logpdf_sub(d, vec(X[sub, 1:length(d)]), transform)
  end
  lp
end

#################### Rand Fallbacks ####################

rand_sub(d::Distribution, x) = rand(d)

rand_sub(d::UnivariateDistribution, X::AbstractArray) = rand(d, size(X))

rand_sub(D::Array{UnivariateDistribution}, X::AbstractArray) = map(rand, D)

function rand_sub(D::Array{MultivariateDistribution}, X::AbstractArray)
  Y = fill(NaN, size(X))
  for sub in CartesianRange(size(D))
    d = D[sub]
    Y[sub, 1:length(d)] = rand(d)
  end
  Y
end


text:
using Distributions: MvNormal

struct GaussianDataset{D} <: AbstractDataset{D}
    X
    Xt
end

n_display(::GaussianDataset) = 1_000

function GaussianDataset(
    n_data::Int,
    mean::AbstractVector{T},
    cov::AbstractVecOrMat{T};
    seed::Int=1,
    test_ratio=1/6,
    n_test::Int=ratio2num(n_data, test_ratio),
    rng=MersenneTwister(seed),
) where {T<:AbstractFloat}
    X = rand(rng, MvNormal(mean, cov), n_data)
    Xt = rand(rng, MvNormal(mean, cov), n_test)
    D = length(mean)
    @info "Oh you just get the $(D)D Gaussian dataset" n_data=n_data n_test=n_test mean=mean cov=cov
    return GaussianDataset{D}(X, Xt)
end

GaussianDataset(n_data; kwargs...) = GaussianDataset(n_data, 2f0 * ones(Float32, 2), [1 81f-2; 81f-2 1])


text:
using Flux, Test, Mill
using Mill: reflectinmodel, length2bags
using Combinatorics

BAGS = [
        length2bags([1 for _ in 1:10]),
        length2bags([2 for _ in 1:5]),
        length2bags([5, 5]),
        length2bags([10]),
length2bags([3, 4, 3]),
AlignedBags([1:3, 0:-1, 0:-1, 4:7, 0:-1, 8:10]),
AlignedBags([0:-1, 1:5, 0:-1, 0:-1, 0:-1, 6:10]),
ScatteredBags([collect(1:3), collect(7:10), collect(4:6)]),
ScatteredBags([collect(7:10), [], collect(1:3), [], collect(4:6), []]),
ScatteredBags([[], collect(1:10), []]),
]

BAGS2 = [
         (AlignedBags([1:2, 3:4, 0:-1]), ScatteredBags([[2,3,4], [1], []]), AlignedBags([1:4, 0:-1, 5:8, 0:-1])),
         (AlignedBags([0:-1, 1:2, 3:4]), ScatteredBags([[1], [2], [3, 4]]), AlignedBags([0:-1, 1:7, 0:-1, 8:8])),
         (AlignedBags([0:-1, 0:-1, 1:2, 3:4]), ScatteredBags([[2,4], [], [3, 1], []]), AlignedBags([1:1, 2:2, 0:-1, 3:8])),
         (AlignedBags([0:-1, 1:2, 3:4, 0:-1]), ScatteredBags([[], [1,3], [2,4], []]), AlignedBags([0:-1, 1:2, 3:6, 7:8]))
        ]

@testset "aggregation grad check w.r.t. input" begin
    for bags in BAGS
        d = rand(1:20)
        x = randn(d, 10)
        w = abs.(randn(size(x, 2))) .+ 0.1
        w_mat = abs.(randn(size(x))) .+ 0.1

        # generate all combinations of aggregations
        anames = ["Mean", "Max", "PNorm", "LSE"]
        for idxs in powerset(collect(1:length(anames)))
            !isempty(idxs) || continue
            # not a thorough testing of all functions, but fast enough
            # for idxs in permutations(idxs)

            s = Symbol("Segmented", anames[idxs]...)
            a = @eval $s($d)
            @test mgradtest(x) do x
                a(x, bags)
            end
            @test mgradtest(x) do x
                a(x, bags, w)
            end
            @test mgradtest(x) do x
                a(x, bags, w_mat)
            end
        end
    end
end

@testset "aggregation grad check w.r.t. agg params" begin
    for bags in BAGS
        d = rand(1:20)
        x = randn(d, 10)
        w = abs.(randn(size(x, 2))) .+ 0.1
        w_mat = abs.(randn(size(x))) .+ 0.1

        fs = [:SegmentedMax, :SegmentedMean, :SegmentedPNorm, :SegmentedLSE]
        params = (:C1,), (:C2,), (:ρ, :c, :C3), (:p, :C4)

        for idxs in powerset(collect(1:length(fs)))
            !isempty(idxs) || continue;
            rs = []; as = []; cs = [];
            for (f, ps) in zip(fs[idxs], params[idxs])
                # p too small causes numerical issues
                push!(rs, (p == :p ? :(clamp.(randn($d), 0.1, Inf)) : :(randn($d)) for p in ps)...)
                push!(as, ps...)
                push!(cs, Expr(:call, f, ps...))
            end
            @eval begin
                @test mgradtest($(map(eval, rs)...)) do $(as...)
                    a = Aggregation($(cs...))
                    a(missing, $bags)
                end
                @test mgradtest($(map(eval, rs)...)) do $(as...)
                    a = Aggregation($(cs...))
                    a($x, $bags)
                end
                @test mgradtest($(map(eval, rs)...)) do $(as...)
                    a = Aggregation($(cs...))
                    a(missing, $bags, nothing)
                end
                @test mgradtest($(map(eval, rs)...)) do $(as...)
                    a = Aggregation($(cs...))
                    a($x, $bags, $w)
                end
                @test mgradtest($(map(eval, rs)...)) do $(as...)
                    a = Aggregation($(cs...))
                    a($x, $bags, $w_mat)
                end
            end
        end
    end
end

@testset "derivative w.r.t weights in aggregations" begin
    for bags in BAGS
        d = rand(1:20)
        x = randn(d, 10)
        w = abs.(randn(size(x, 2))) .+ 0.1
        w_mat = abs.(randn(size(x))) .+ 0.1

        a1 = f64(SegmentedMean(d))
        a2 = f64(SegmentedMax(d))
        a3 = f64(SegmentedPNorm(d))
        a4 = f64(SegmentedLSE(d))
        for g in [
                  w -> a1(x, bags, w),
                  w -> a2(x, bags, w),
                  w -> a4(x, bags, w)
                 ]
            @test mgradtest(g, w)
        end
        for g in [
                  w_mat -> a1(x, bags, w_mat),
                  w_mat -> a2(x, bags, w_mat),
                  w_mat -> a4(x, bags, w_mat)
                 ]
            @test mgradtest(g, w_mat)
        end
        # for g in [
        #           w -> a3(x, bags, w_mat)
        #          ]
        #     # NOT IMPLEMENTED YET
        #     @test_throws Exception mgradtest(g, w_mat)
        # end
        # for g in [
        #           w -> a3(x, bags, w_mat)
        #          ]
        #     # NOT IMPLEMENTED YET
        #     @test_throws Exception mgradtest(g, w_mat)
        # end
    end
end

@testset "model aggregation grad check w.r.t. inputs" begin
    for (bags1, bags2, bags3) in BAGS2
        layerbuilder(k) = Dense(k, 2, NNlib.relu)
        x = randn(4, 4)
        y = randn(3, 4)
        z = randn(2, 8)

        n = ArrayNode(x)
        m = f64(reflectinmodel(n, layerbuilder))
        @test mgradtest(x) do x
            n = ArrayNode(x)
            m(n).data
        end

        bn = BagNode(ArrayNode(x), bags1)
        abuilder = d -> SegmentedPNormLSE(d)
        m = f64(reflectinmodel(bn, layerbuilder))
        @test mgradtest(x) do x
            bn = BagNode(ArrayNode(x), bags1)
            m(bn).data
        end

        tn = TreeNode((ArrayNode(x), ArrayNode(y)))
        m = f64(reflectinmodel(tn, layerbuilder))
        @test mgradtest(x, y) do x, y
            tn = TreeNode((ArrayNode(x), ArrayNode(y)))
            m(tn).data
        end

        tn = TreeNode((BagNode(ArrayNode(y), bags1), BagNode(ArrayNode(x), bags2)))
        abuilder = d -> SegmentedMeanMax(d)
        m = f64(reflectinmodel(tn, layerbuilder, abuilder))
        @test mgradtest(x, y) do x, y
            tn = TreeNode((BagNode(ArrayNode(y), bags1), BagNode(ArrayNode(x), bags2)))
            m(tn).data
        end

        bn = BagNode(ArrayNode(z), bags3)
        bnn = BagNode(bn, bags1)
        abuilder = d -> SegmentedPNormLSEMeanMax(d)
        m = f64(reflectinmodel(bnn, layerbuilder, abuilder))
        @test mgradtest(z) do z
            bn = BagNode(ArrayNode(z), bags3)
            bnn = BagNode(bn, bags1)
            m(bnn).data
        end
    end
end

@testset "model aggregation grad check w.r.t. inputs weighted" begin
    for (bags1, bags2, bags3) in BAGS2
        layerbuilder(k) = Dense(k, 2, NNlib.relu)
        x = randn(4, 4)
        y = randn(3, 4)
        z = randn(2, 8)
        w = abs.(randn(4)) .+ 0.1
        w2 = abs.(randn(4)) .+ 0.1
        w3 = abs.(randn(8)) .+ 0.1

        bn = BagNode(ArrayNode(x), bags1, w)
        abuilder = d -> SegmentedPNormLSE(d)
        m = f64(reflectinmodel(bn, layerbuilder))
        @test mgradtest(x) do x
            bn = BagNode(ArrayNode(x), bags1, w)
            m(bn).data
        end

        tn = TreeNode((BagNode(ArrayNode(y), bags1, w), BagNode(ArrayNode(x), bags2, w2)))
        abuilder = d -> SegmentedMeanMax(d)
        m = f64(reflectinmodel(tn, layerbuilder, abuilder))
        @test mgradtest(x, y) do x, y
            tn = TreeNode((BagNode(ArrayNode(y), bags1, w), BagNode(ArrayNode(x), bags2, w2)))
            m(tn).data
        end

        bn = BagNode(ArrayNode(z), bags3, w3)
        bnn = BagNode(bn, bags1)
        abuilder = d -> SegmentedPNormLSEMeanMax(d)
        m = f64(reflectinmodel(bnn, layerbuilder, abuilder))
        @test mgradtest(z) do z
            bn = BagNode(ArrayNode(z), bags3, w3)
            bnn = BagNode(bn, bags1, w)
            m(bnn).data
        end
    end
end

@testset "model aggregation grad check w.r.t. params" begin
    for (bags1, bags2, bags3) in BAGS2
        layerbuilder(k) = Dense(k, 2, NNlib.relu)
        x = randn(4, 4)
        y = randn(3, 4)
        z = randn(2, 8)

        n = ArrayNode(x)
        m = f64(reflectinmodel(n, layerbuilder))
        @test mgradtest(params(m)...) do W, b
            m = ArrayModel(Dense(W, b, relu))
            m(n).data
        end

        bn = BagNode(ArrayNode(x), bags1)
        abuilder = d -> SegmentedPNormLSE(d)
        m = f64(reflectinmodel(bn, layerbuilder, abuilder))
        @test mgradtest(params(m)...) do W1, b1, ρ, c, C1, p, C2, W2, b2
            m = BagModel(Dense(W1, b1, relu),
                         Aggregation(
                                     SegmentedPNorm(ρ, c, C1),
                                     SegmentedLSE(p, C2)
                                    ),
                         Dense(W2, b2, σ))
            m(bn).data
        end

        tn = TreeNode((ArrayNode(x), ArrayNode(y)))
        m = f64(reflectinmodel(tn, layerbuilder))
        @test mgradtest(params(m)...) do W1, b1, W2, b2, W3, b3
            m = ProductModel(ArrayModel.((
                                          Dense(W1, b1, σ),
                                          Dense(W2, b2, relu)
                                         )), Dense(W3, b3, σ))
            m(tn).data
        end

        tn = TreeNode((BagNode(ArrayNode(y), bags1), BagNode(ArrayNode(x), bags2)))
        abuilder = d -> SegmentedPNormLSEMeanMax(d)
        m = f64(reflectinmodel(tn, layerbuilder, abuilder))
        @test mgradtest(params(m)...) do W1, b1, ρ1, c1, C11, p1, C12, C13, C14,
            W2, b2, W3, b3, ρ2, c2, C21, p2, C22, C23, C24, W4, b4, W5, b5
            m = ProductModel((
                              BagModel(
                                       Dense(W1, b1, σ),
                                       Aggregation(
                                                   SegmentedPNorm(ρ1, c1, C11),
                                                   SegmentedLSE(p1, C12),
                                                   SegmentedMean(C13),
                                                   SegmentedMax(C14)
                                                  ),
                                       Dense(W2, b2, relu)
                                      ),
                              BagModel(
                                       Dense(W3, b3, relu),
                                       Aggregation(
                                                   SegmentedPNorm(ρ2, c2, C21),
                                                   SegmentedLSE(p2, C22),
                                                   SegmentedMean(C23),
                                                   SegmentedMax(C24)
                                                  ),
                                       Dense(W4, b4, σ)
                                      ),
                             ), Dense(W5, b5, relu))
            m(tn).data
        end

        bn = BagNode(ArrayNode(z), bags3)
        bnn = BagNode(bn, bags1)
        abuilder = d -> SegmentedMeanMax(d)
        m = f64(reflectinmodel(bnn, layerbuilder, abuilder))
        @test mgradtest(params(m)...) do W1, b1, C11, C12, W2, b2, C21, C22, W3, b3
            m = BagModel(
                         BagModel(
                                  Dense(W1, b1),
                                  Aggregation(
                                              SegmentedMean(C11),
                                              SegmentedMax(C12)
                                             ),
                                  Dense(W2, b2)
                                 ),
                         Aggregation(
                                     SegmentedMean(C21),
                                     SegmentedMax(C22)
            ),
                         Dense(W3, b3)
                        )
            m(bnn).data
        end
    end
end

@testset "model aggregation grad check w.r.t. params weighted" begin
    for (bags1, bags2, bags3) in BAGS2
        layerbuilder(k) = Dense(k, 2, NNlib.relu)
        x = randn(4, 4)
        y = randn(3, 4)
        z = randn(2, 8)
        w = abs.(randn(4)) .+ 0.1
        w2 = abs.(randn(4)) .+ 0.1
        w3 = abs.(randn(8)) .+ 0.1

        bn = BagNode(ArrayNode(x), bags1, w)
        abuilder = d -> SegmentedPNormLSE(d)
        m = f64(reflectinmodel(bn, layerbuilder, abuilder))
        @test mgradtest(params(m)...) do W1, b1, ρ, c, C1, p, C2, W2, b2
            m = BagModel(Dense(W1, b1, relu),
                         Aggregation(
                                     SegmentedPNorm(ρ, c, C1),
                                     SegmentedLSE(p, C2)
                                    ),
                         Dense(W2, b2, σ))
            m(bn).data
        end

        tn = TreeNode((BagNode(ArrayNode(y), bags1, w), BagNode(ArrayNode(x), bags2, w2)))
        abuilder = d -> SegmentedPNormLSEMeanMax(d)
        m = f64(reflectinmodel(tn, layerbuilder, abuilder))
        @test mgradtest(params(m)...) do W1, b1, ρ1, c1, C11, p1, C12, C13, C14,
            W2, b2, W3, b3, ρ2, c2, C21, p2, C22, C23, C24, W4, b4, W5, b5
            m = ProductModel((
                              BagModel(
                                       Dense(W1, b1, σ),
                                       Aggregation(
                                                   SegmentedPNorm(ρ1, c1, C11),
                                                   SegmentedLSE(p1, C12),
                                                   SegmentedMean(C13),
                                                   SegmentedMax(C14)
                                                  ),
                                       Dense(W2, b2, relu)
                                      ),
                              BagModel(
                                       Dense(W3, b3, relu),
                                       Aggregation(
                                                   SegmentedPNorm(ρ2, c2, C21),
                                                   SegmentedLSE(p2, C22),
                                                   SegmentedMean(C23),
                                                   SegmentedMax(C24)
                                                  ),
                                       Dense(W4, b4, σ)
                                      ),
                             ), Dense(W5, b5, relu))
            m(tn).data
        end

        bn = BagNode(ArrayNode(z), bags3, w3)
        bnn = BagNode(bn, bags1, w)
        abuilder = d -> SegmentedMeanMax(d)
        m = f64(reflectinmodel(bnn, layerbuilder, abuilder))
        @test mgradtest(params(m)...) do W1, b1, C11, C12, W2, b2, C21, C22, W3, b3
            m = BagModel(
                         BagModel(
                                  Dense(W1, b1),
                                  Aggregation(
                                              SegmentedMean(C11),
                                              SegmentedMax(C12)
                                             ),
                                  Dense(W2, b2)
                                 ),
                         Aggregation(
                                     SegmentedMean(C21),
                                     SegmentedMax(C22)
            ),
                         Dense(W3, b3)
                        )
            m(bnn).data
        end
    end

  @testset "A gradient of ProductNode with a NamedTuple " begin
    a = TreeNode((a = ArrayNode(randn(2,4)), b = ArrayNode(randn(3,4))))
    m = f64(ProductModel((a = ArrayModel(Dense(2,2)), b = ArrayModel(Dense(3,1))), Dense(3,2)))
    ps = params(m)
    gradient(() -> sum(m(a).data), ps)
    @test mgradtest(params(m)...) do W1, b1, W2, b2, W3, b3
      m = ProductModel((a = ArrayModel(Dense(W1, b1)), b = ArrayModel(Dense(W2, b2))), Dense(W3, b3))
      sum(m(a).data)
    end
  end
end


text:
# Scenario 1, peak of the storm
θ₁ = [0.7, 0.8, 3.3, 5, π / 2, 4, 2.7, 0.55, 0.26]
# Scenario 2, higher frequency constant spread
θ₂ = [0.7, 1.1, 3.3, 5, π / 2, 4, 2.7, 0.55, 0]
# Scenario 3, PM style
θ₃ = [0.7, 1.0, 1, 5, π / 2, 4, 2.7, 0.55, 0.26]
#
Θ = [θ₁, θ₂, θ₃]

cutoff = [0.6, 0.75, 0.7]

##
include("src/simulation/PaperSimulation.jl")
import Random
Random.seed!(1234)

## Check on small example
@time sim_small = [
    PaperSimulation.SimulationStudy(1, θᵢ; ωₗ = c, name = "Scenario $i ") for
    ((i, θᵢ), c) ∈ zip(enumerate(Θ), cutoff)
]

## run simulation study
@time sim = [
    PaperSimulation.SimulationStudy(1000, θᵢ; ωₗ = c, name = "Scenario $i ") for
    ((i, θᵢ), c) ∈ zip(enumerate(Θ), cutoff)
]
using JLD2
@save "dat/simulation.jld" sim


text:
using ReinforcementLearningEnvironmentViZDoom, Test
import ReinforcementLearningEnvironmentViZDoom: interact!, reset!, getstate

function test1()
    env = ViZDoomEnvironment("", "map01")
    s, r, done = interact!(1, env)
    @test typeof(s) == Array{UInt8, 1}
    s = reset!(env)
    @test typeof(s) == Array{UInt8, 1}
    s, done = getstate(env)
    @test typeof(s) == Array{UInt8, 1}
end
test1()

max_stars_repo_name: jbrea/RLEnvVizDoom

max_issues_repo_name: jbrea/RLEnvVizDoom

max_forks_repo_name: jbrea/RLEnvVizDoom


text:
# # Poisson with Neumann conditions

#md # ```@meta
#md # CurrentModule = ImmersedLayers
#md # ```

#=
Here, we'll demonstrate a solution of Laplace's equation,
with Neumann boundary conditions on a surface. Similar to the
Dirichlet problem in [A Dirichlet Poisson problem](@ref), we will solve the
problem with one Neumann condition external to the surface,
and another Neumann value internal to the surface.

Our underlying problem is still

$$\nabla^2\varphi^+ = 0,\qquad \nabla^2\varphi^- = 0$$

where $+$ denotes the exterior and $-$ the interior of the surface. (We will
consider a circle of radius 1.) The boundary conditions on this surface are

$$\mathbf{n}\cdot\nabla\varphi^+ = v^+_n, \qquad \mathbf{n}\cdot\nabla\varphi^- = v^-_n$$

In other words, we seek to set the value on the exterior normal derivative to $v_n$
of the local normal vector on the surface, while the interior should have zero normal
derivative.

Discretizing this problem by the usual techniques, we seek to solve

$$\begin{bmatrix} L & D_s \\ G_s & R_n^T R_n \end{bmatrix} \begin{pmatrix} f \\ -[\phi] \end{pmatrix} = \begin{pmatrix} R [v_n] \\ \overline{v}_n \end{pmatrix}$$

where $\overline{v}_n = (v^+_n + v^-_n)/2$ and $[v_n] = v^+_n - v^-_n$. The resulting
$[\phi]$ is $f^+-f^-$.

As with the Dirichlet problem, this saddle-point problem can be solved by block-LU decomposition. First solve

$$L f^{*} = R [v_n]$$

for $f^*$. Then solve

$$-S [\phi] = \overline{v}_n - G_s f^{*}$$

for $[\phi]$, where $S = R_n^T R_n - G_s L^{-1} D_s = -C_s L^{-1}C_s^T$, and finally, compute

$$f = f^{*} + L^{-1}D_s [\phi]$$

It should be remembered that, for any scalar potential field, there is
a corresponding streamfunction $\psi$ that generates the same flow. We
can get that field, as well, with only a little bit more effort:

$$S [\psi] = C_s f^{*}$$

and then solve

$$L s = C_s^T [\phi] - \hat{C}_s^T [\psi]$$.

for the streamfunction $s$.
=#

using ImmersedLayers
using Plots
using LinearAlgebra
using UnPack

#=
## Set up the extra cache and solve function
=#
#=
The problem type takes the usual basic form
=#
@ilmproblem NeumannPoisson scalar

#=
The extra cache holds additional intermediate data, as well as
the Schur complement. We don't bother creating a filtering matrix here.
=#
struct NeumannPoissonCache{SMT,DVT,VNT,FT,ST} <: AbstractExtraILMCache
   S :: SMT
   dvn :: DVT
   vn :: VNT
   fstar :: FT
   sstar :: ST
end
#=
The function `prob_cache`, as before, constructs the operators and extra
cache data structures
=#
function ImmersedLayers.prob_cache(prob::NeumannPoissonProblem,base_cache::BasicILMCache)
    S = create_CLinvCT(base_cache)
    dvn = zeros_surface(base_cache)
    vn = zeros_surface(base_cache)
    fstar = zeros_grid(base_cache)
    sstar = zeros_gridcurl(base_cache)
    NeumannPoissonCache(S,dvn,vn,fstar,sstar)
end
nothing #hide

#=
And finally, here's the steps we outlined above, used to
extend the `solve` function
=#
function ImmersedLayers.solve(vnplus,vnminus,prob::NeumannPoissonProblem,sys::ILMSystem)
    @unpack extra_cache, base_cache = sys
    @unpack S, dvn, vn, fstar, sstar = extra_cache

    fill!(fstar,0.0)
    fill!(sstar,0.0)

    f = zeros_grid(base_cache)
    s = zeros_gridcurl(base_cache)
    df = zeros_surface(base_cache)
    ds = zeros_surface(base_cache)

    vn .= 0.5*(vnplus+vnminus)
    dvn .= vnplus - vnminus

    ## Find the potential
    regularize!(fstar,dvn,base_cache)
    inverse_laplacian!(fstar,base_cache)

    surface_grad!(df,fstar,base_cache)
    df .= vn - df
    df .= -(S\df);

    surface_divergence!(f,df,base_cache)
    inverse_laplacian!(f,base_cache)
    f .+= fstar

    ## Find the streamfunction
    surface_curl!(sstar,df,base_cache)

    surface_grad_cross!(ds,fstar,base_cache)
    ds .= S\ds

    surface_curl_cross!(s,ds,base_cache)
    s .-= sstar
    s .*= -1.0

    inverse_laplacian!(s,base_cache)

    return f, df, s, ds
end
nothing #hide

#=
## Solve the problem
Here, we will demonstrate the solution on a circular shape of radius 1,
with $v_n^+ = n_x$ and $v_n^- = 0$. This is actually the set of conditions
used to compute the unit scalar potential field (and, as we will see, the added mass) in potential flow.
=#

#=
Set up the grid
=#
Δx = 0.01
Lx = 4.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx)
Δs = 1.4*cellsize(g)
body = Circle(1.0,Δs);

#=
Create the system
=#
prob = NeumannPoissonProblem(g,body,scaling=GridScaling)
sys = ImmersedLayers.__init(prob)
nothing #hide

#=
Set the boundary values
=#
nrm = normals(sys)
vnplus = zeros_surface(sys)
vnminus = zeros_surface(sys)
vnplus .= nrm.u
nothing #hide

#=
Solve it
=#
solve(vnplus,vnminus,prob,sys) #hide
@time f, df, s, ds = solve(vnplus,vnminus,prob,sys);

#=
and plot the field
=#
plot(plot(f,sys,layers=true,levels=30,title="ϕ"),
plot(s,sys,layers=true,levels=30,title="ψ"))

#=
and the Lagrange multiplier field, $[\phi]$, on the surface
=#
plot(df)

#=
If, instead, we set the inner boundary condition to $n_x$ and the
outer to zero, then we get the flow *inside* of a translating circle
=#
vnplus = zeros_surface(sys)
vnminus = zeros_surface(sys)
vnminus .= nrm.u
f, df, s, ds = solve(vnplus,vnminus,prob,sys);
plot(plot(f,sys,layers=true,levels=30,title="ϕ"),
plot(s,sys,layers=true,levels=30,title="ψ"))

#=
## Multiple bodies
The cache and solve function we created above can be applied for
any body or set of bodies. Let's apply it here to a circle of radius 0.25 inside of a
square of half-side length 2, where our goal is to find the effect of the enclosing square
on the motion of the circle. As such, we will set the Neumann conditions
both internal to and external to the square to be 0, but for the exterior of the
circle, we set it to $n_x$.
=#
bl = BodyList();
push!(bl,Square(1.0,Δs))
push!(bl,Circle(0.25,Δs))
nothing #hide

#=
We don't actually have to transform these shapes, but it is
illustrative to show how we would move them.
=#
t1 = RigidTransform((0.0,0.0),0.0)
t2 = RigidTransform((0.0,0.0),0.0)
tl = RigidTransformList([t1,t2])
tl(bl)
nothing #hide

#=
Create the problem and system
=#
prob = NeumannPoissonProblem(g,bl,scaling=GridScaling)
sys = ImmersedLayers.__init(prob)
nothing #hide

#=
Set the boundary conditions. We set only the exterior Neumann value
of body 2 (the circle), using [`copyto!`](@ref)
=#
nrm = normals(sys)
vnplus = zeros_surface(sys)
vnminus = zeros_surface(sys)
copyto!(vnplus,nrm.u,sys,2)

#=
Solve it and plot
=#
f, df, s, ds  = solve(vnplus,vnminus,prob,sys)
plot(plot(f,sys,layers=true,levels=30,title="ϕ"),
plot(s,sys,layers=true,levels=30,title="ψ"))

#=
Now, let's compute the added mass components of the circle associated
with this motion. We are approximating

$$M = -\int_{C_2} f^+ \mathbf{n}\mathrm{d}s$$

where $C_2$ is shape 2 (the circle), and $f^+$ is simply $[\phi]$ on body 2.
=#
M = -integrate(df∘nrm,sys,2)

#=
As one would expect, the circle has added mass in the $x$ direction
associated with moving in that direction.
=#


text:
module Experimental
using RvSpectMLBase
using EchelleInstruments
using EchelleCCFs
using DataFrames, Query
using Dates
using Statistics

CCFs = EchelleCCFs

include("pipeline/pipeline.jl")

include("interp/interp.jl")
#=
include("interp/gp/temporalgps.jl")
export TemporalGPInterpolation
export construct_gp_posterior, gp_marginal, predict_gp
export predict_mean, predict_deriv, predict_deriv2, predict_mean_and_deriv, predict_mean_and_derivs
=#

include("line_finder/line_finder.jl")
export LineFinderPlan

include("alg/project_flux_common_wavelengths.jl")
include("alg/make_template_spectrum.jl")

#include("alg/dcpca.jl")
#include("alg/ppcca.jl")
#include("alg/rvs_from_gp_pairs.jl")
end


text:
using BenchmarkTools
using Augmentor, Images, TestImages

const SUITE = BenchmarkGroup()

pattern = testpattern()
affpattern = Augmentor.prepareaffine(pattern)
camera = testimage("cameraman")
cameras = similar(camera, size(camera)..., 2)
copy!(view(cameras,:,:,1), camera)
copy!(view(cameras,:,:,2), camera)

# compile the gauss filter (takes a couple seconds)
Augmentor.applyeager(ElasticDistortion(10), pattern)

shortname(op) = string(typeof(op).name.name)

# missing: CombineChannels, CacheImageInto
SUITE["applyeager"] = BenchmarkGroup(["eager", "operations"])
for op in [
        NoOp(),
        FlipX(), FlipY(),
        ShearX(5), ShearY(5),
        Either((Rotate90(), Rotate270()), (1,0)),
        Rotate90(), Rotate180(), Rotate270(),
        Rotate(45),
        Crop(1:100,1:100),
        CropNative(1:100,1:100),
        CropSize(100,100),
        CropRatio(1), RCropRatio(1),
        Scale(1.1), Zoom(1.1),
        Resize(100,100),
        ConvertEltype(Gray),
        CacheImage(),
        ElasticDistortion(10),
        SplitChannels(),
        Reshape(400,300),
        PermuteDims(2,1),
    ]
    SUITE["applyeager"][shortname(op)] = @benchmarkable Augmentor.applyeager($op, $pattern)
end

SUITE["applylazy"] = BenchmarkGroup(["lazy", "operations"])
for op in [
        NoOp(),
        FlipX(), FlipY(),
        ShearX(5), ShearY(5),
        Either((Rotate90(), Rotate270()), (1,0)),
        Rotate90(), Rotate180(), Rotate270(),
        Rotate(45),
        Crop(1:100,1:100),
        CropNative(1:100,1:100),
        CropSize(100,100),
        CropRatio(1), RCropRatio(1),
        Scale(1.1), Zoom(1.1),
        Resize(100,100),
        ConvertEltype(Gray),
        ElasticDistortion(10),
        SplitChannels(),
        Reshape(400,300),
        PermuteDims(2,1),
    ]
    SUITE["applylazy"][shortname(op)] = @benchmarkable Augmentor.plain_array(Augmentor.applylazy($op, $pattern))
end

SUITE["applyaffine"] = BenchmarkGroup(["affine", "lazy", "operations"])
for op in [
        NoOp(),
        FlipX(), FlipY(),
        ShearX(5), ShearY(5),
        Either((Rotate90(), Rotate270()), (1,0)),
        Rotate90(), Rotate180(), Rotate270(),
        Rotate(45),
        Scale(1.1),
    ]
    SUITE["applyaffine"][shortname(op)] = @benchmarkable Augmentor.plain_array(Augmentor.applyaffine($op, $affpattern))
end

SUITE["applyaffineview"] = BenchmarkGroup(["affine", "lazy", "operations"])
for op in [
        NoOp(),
        FlipX(), FlipY(),
        ShearX(5), ShearY(5),
        Either((Rotate90(), Rotate270()), (1,0)),
        Rotate90(), Rotate180(), Rotate270(),
        Rotate(45),
        Crop(1:100,1:100),
        CropNative(1:100,1:100),
        CropSize(100,100),
        CropRatio(1), RCropRatio(1),
        Scale(1.1), Zoom(1.1),
        Resize(100,100),
    ]
    SUITE["applyaffineview"][shortname(op)] = @benchmarkable Augmentor.plain_array(Augmentor.applyaffineview($op, $affpattern))
end

grp = SUITE["augment"] = BenchmarkGroup(["-"])
pl = ShearX(10)
grp["affine1"] = @benchmarkable augment($pattern, $pl)
pl = ShearX(10) |> ShearX(-10)
grp["affine2"] = @benchmarkable augment($pattern, $pl)
pl = ShearX(10) |> ShearX(-10) |> ShearY(10) |> ShearY(-10)
grp["affine3"] = @benchmarkable augment($pattern, $pl)
pl = Rotate180() |> Rotate90(1)
grp["lazyrotate"] = @benchmarkable augment($pattern, $pl)
pl = Rotate180() |> Crop(5:200,100:400) |> Rotate90(1) |> Crop(1:250,1:150)
grp["lazycrop"] = @benchmarkable augment($pattern, $pl)
pl = Resize(100,100) |> Resize(200,400) |> Resize(20,20) |> Resize(100,100)
grp["resize"] = @benchmarkable augment($pattern, $pl)

grp = SUITE["augment!"] = BenchmarkGroup(["-"])
pl = Rotate180() |> Rotate90()
out = similar(pattern, 400, 300)
grp["lazyrotate"] = @benchmarkable augment!($out, $pattern, $pl)
pl = Resize(100,100) |> Resize(200,400) |> Resize(20,20) |> Resize(100,100)
out = similar(pattern, 100, 100)
grp["resize"] = @benchmarkable augment!($out, $pattern, $pl)


text:
#=
    Absorb behavior into the medium (Float64).

    Project state-free value as an unmarked stateful quantity (the projection).
    Mark the stateful quantity, or do not mark it (the stateful projection).

    Unmark the stateful projection.
    Reflect the projection as state-free value (the reflection).
=#

#=
    givenFloat
         FloatExponent  -->  stationedExponent
                                      stationedFloat

                                      stationedFloat
                                               stationedExponent  -->  stationaryExponent
          FloatExponent                                                stationaryExponent
               givenFloat                  =                 stationaryFloat

=#
@inline function pushout{F<:AbstractFloat}(fp::F)
   stationedExponent = get_exponent(fp) - Bias(F)
   reinterpret(F,put_exponent(fp, stationedExponent))
end                                                    # pushout of/from a systemFloat, pushout a stationedFloat

@inline function pullback{F<:AbstractFloat}(fp::F)
   stationaryExponent = get_exponent(fp) + Bias(F)
   reinterpret(F,put_exponent(fp, stationaryExponent))
end                                                    # pullback from/of a stationedFloat, pullback a systemFloat

function project{F<:AbstractFloat}(fp::F)
    if Tiny(F) < fp < Huge(F)
      pushout(fp)
    elseif fp == 0.0
       fp
    elseif isfinite(fp)
       if fp == 0.0
           fp
       elseif fp < Huge(F)    # Tiny and Huge are formed in the enchaned state (collective values)
           # signbit(fp) ? -set_ebit(TinyProjected(F)) : set_ebit(TinyProjected(F))
           signbit(fp) ? -TinyProjected(F) : TinyProjected(F)
       else
           # signbit(fp) ? -set_ebit(HugeProjected(F)) : set_ebit(HugeProjected(F))
           signbit(fp) ? -HugeProjected(F) : HugeProjected(F)
       end
    else  # ±Inf or NaN
       isinf(fp) ? set_ebit(fp) : fp
    end
end

function reflect{F<:AbstractFloat}(fp::F)
    if iszero(fp) | isnotfinite(fp)
       signbit(fp) ? -fp : fp
    else
       cafp = abs(clr_ebit(fp))
       if TinyProjected(F) < cafp < HugeProjected(F)
           flipsign( pullback(cafp), fp)
       elseif cafp <= TinyProjected(F)
           flipsign( Tiny(F), fp )
       else
           flipsign( Huge(F), fp )
       end
    end
end


text:
module SwitchGate

using ..Ahorn, Maple

function gateFinalizer(entity)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    entity.data["nodes"] = [(x + width, y)]
end

const textures = ["block", "mirror", "temple", "stars"]

const placements = Ahorn.PlacementDict(
    "Switch Gate (Stone)" => Ahorn.EntityPlacement(
        Maple.SwitchGate,
        "rectangle",
        Dict{String, Any}(
            "sprite" => "block"
        ),
        gateFinalizer
    ),
    "Switch Gate (Mirror)" => Ahorn.EntityPlacement(
        Maple.SwitchGate,
        "rectangle",
        Dict{String, Any}(
            "sprite" => "mirror"
        ),
        gateFinalizer
    ),
    "Switch Gate (Temple)" => Ahorn.EntityPlacement(
        Maple.SwitchGate,
        "rectangle",
        Dict{String, Any}(
            "sprite" => "temple"
        ),
        gateFinalizer
    ),
    "Switch Gate (Moon)" => Ahorn.EntityPlacement(
        Maple.SwitchGate,
        "rectangle",
        Dict{String, Any}(
            "sprite" => "stars"
        ),
        gateFinalizer
    ),

    "Touch Switch" => Ahorn.EntityPlacement(
        Maple.TouchSwitch
    )
)

Ahorn.editingOptions(entity::Maple.SwitchGate) = Dict{String, Any}(
    "sprite" => textures
)

Ahorn.nodeLimits(entity::Maple.SwitchGate) = 1, 1

Ahorn.minimumSize(entity::Maple.SwitchGate) = 16, 16
Ahorn.resizable(entity::Maple.SwitchGate) = true, true

function Ahorn.selection(entity::Maple.SwitchGate)
    x, y = Ahorn.position(entity)
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return [Ahorn.Rectangle(x, y, width, height), Ahorn.Rectangle(stopX, stopY, width, height)]
end

function Ahorn.selection(entity::Maple.TouchSwitch)
    x, y = Ahorn.position(entity)

    return  Ahorn.Rectangle(x - 7, y - 7, 14, 14)
end

iconResource = "objects/switchgate/icon00"

function renderGateSwitch(ctx::Ahorn.Cairo.CairoContext, x::Int, y::Int, width::Int, height::Int, sprite::String)
    iconSprite = Ahorn.getSprite(iconResource, "Gameplay")

    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    frame = "objects/switchgate/$sprite"

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y, 8, 0, 8, 8)
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y + height - 8, 8, 16, 8, 8)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, x, y + (i - 1) * 8, 0, 8, 8, 8)
        Ahorn.drawImage(ctx, frame, x + width - 8, y + (i - 1) * 8, 16, 8, 8, 8)
    end

    for i in 2:tilesWidth - 1, j in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y + (j - 1) * 8, 8, 8, 8, 8)
    end

    Ahorn.drawImage(ctx, frame, x, y, 0, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, x + width - 8, y, 16, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, x, y + height - 8, 0, 16, 8, 8)
    Ahorn.drawImage(ctx, frame, x + width - 8, y + height - 8, 16, 16, 8, 8)

    Ahorn.drawImage(ctx, iconSprite, x + div(width - iconSprite.width, 2), y + div(height - iconSprite.height, 2))
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.SwitchGate, room::Maple.Room)
    sprite = get(entity.data, "sprite", "block")
    startX, startY = Int(entity.data["x"]), Int(entity.data["y"])
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    renderGateSwitch(ctx, stopX, stopY, width, height, sprite)
    Ahorn.drawArrow(ctx, startX + width / 2, startY + height / 2, stopX + width / 2, stopY + height / 2, Ahorn.colors.selection_selected_fc, headLength=6)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.TouchSwitch, room::Maple.Room)
    Ahorn.drawSprite(ctx, "objects/touchswitch/container.png", 0, 0)
    Ahorn.drawSprite(ctx, "objects/touchswitch/icon00.png", 0, 0)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.SwitchGate, room::Maple.Room)
    sprite = get(entity.data, "sprite", "block")

    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    renderGateSwitch(ctx, x, y, width, height, sprite)
end

end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
# Maximum values of model parameters to ensure that the same set of data is used across the case studies

# max number of covariates we'll ever consider
const maxNumCovariates = 200

# max number of samples we'll ever consider
const maxNumDataSamples = 12000

# max number of MC replicates
const maxNumMCReplicates = 1000

# max number of data replicates
const maxNumDataReplicates = 1000

# max number of model replicates
const maxNumModelReplicates = 1000

max_stars_repo_name: rohitkannan/DD-SAA

max_issues_repo_name: rohitkannan/DD-SAA

max_forks_repo_name: rohitkannan/DD-SAA


text:
###############################################################################
#
# Calculate zenith angle
#
###############################################################################
"""
    zenith_angle(
                latd::FT,
                decd::FT,
                lhad::FT
    ) where {FT<:AbstractFloat}
    zenith_angle(
                latd::FT,
                day::FT,
                hour::FT,
                minute::FT
    ) where {FT<:AbstractFloat}

Calculate the zenith angle, given
- `latd` Latitude in degree
- `decd` Declination of the Sun in degree
- `lhad` Local hour angle in degree
- `day` Day of year
- `hour` Hour of day
- `minute` Minute of hour
"""
function zenith_angle(
            latd::FT,
            decd::FT,
            lhad::FT
) where {FT<:AbstractFloat}
    cosz = sind(latd) * sind(decd) + cosd(latd) * cosd(decd) * cosd(lhad);

    return acosd(cosz)
end

function zenith_angle(
            latd::FT,
            day::FT,
            hour::FT,
            minute::FT
) where {FT<:AbstractFloat}
    _deg::FT = 360 / YEAR_D(FT) * (day + (hour+minute/60) / 24 + 10);
    decd::FT = -FT(23.44) * cosd(_deg);
    lhad::FT = (hour-12) * 15;

    return zenith_angle(latd, decd, lhad)
end

function zenith_angle(
            latd::FT,
            toy::FT
) where {FT<:AbstractFloat}
    _deg::FT = 360 / YEAR_D(FT) * (toy + 10);
    decd::FT = -FT(23.44) * cosd(_deg);
    lhad::FT = ((toy%1) - FT(0.5)) * 360;

    return zenith_angle(latd, decd, lhad)
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
# Helper functions to get node numbers and points
@inline getleft(i::Int) = 2i
@inline getright(i::Int) = 2i + 1
@inline getparent(i::Int) = div(i, 2)
@inline isleaf(n_internal_nodes::Int, idx::Int) = idx > n_internal_nodes

function show(io::IO, tree::NNTree{V}) where {V}
    println(io, typeof(tree))
    println(io, "  Number of points: ", length(tree.data))
    println(io, "  Dimensions: ", length(V))
    println(io, "  Metric: ", tree.metric)
    print(io,   "  Reordered: ", tree.reordered)
end

# We split the tree such that one of the sub trees has exactly 2^p points
# and such that the left sub tree always has more points.
# This means that we can deterministally (with just some comparisons)
# find if we are at a leaf node and how many
function find_split(low, leafsize, n_p)

    # The number of leafs node left in the tree,
    # use `ceil` to count a partially filled node as 1.
    n_leafs = ceil(Int, n_p / leafsize)

    # Number of leftover nodes needed
    k = floor(Integer, log2(n_leafs))
    rest = n_leafs - 2^k

    # The conditionals here fulfill the desired splitting procedure but
    # can probably be written in a nicer way

    # Can fill less than two nodes -> leafsize to left node.
    if n_p <= 2 * leafsize
        mid_idx = leafsize

    # The last leaf node will be in the right sub tree -> fill the left
    # sub tree with
    elseif rest > 2^(k - 1) # Last node over the "half line" in the row
        mid_idx = 2^k * leafsize

    # Perfectly filling both sub trees -> half to left and right sub tree
    elseif rest == 0
        mid_idx = 2^(k - 1) * leafsize

    # Else we fill the right sub tree -> send the rest to the left sub tree
    else
        mid_idx = n_p - 2^(k - 1) * leafsize
    end
    return mid_idx + low
end

# Gets number of points in a leaf node, this is equal to leafsize for every node
# except the last node.
@inline function n_ps(idx::Int, td::TreeData)
    if idx != td.last_full_node
        return td.leafsize
    else
        return td.last_node_size
    end
end

# Returns the index for the first point for a given leaf node.
@inline function point_index(idx::Int, td::TreeData)
    if idx >= td.cross_node
        return td.offset_cross + idx * td.leafsize
    else
        return td.offset + idx * td.leafsize
    end
end

# Returns a range over the points in a leaf node with a given index
@inline function get_leaf_range(td::TreeData, index)
    p_index = point_index(index, td)
    n_p =  n_ps(index, td)
    return p_index:p_index + n_p - 1
end

# Store all the points in a leaf node continuously in memory in data_reordered to improve cache locality.
# Also stores the mapping to get the index into the original data from the reordered data.
function reorder_data!(data_reordered::Vector{V}, data::Vector{V}, index::Int,
                         indices::Vector{Int}, indices_reordered::Vector{Int}, tree_data::TreeData) where {V}

    for i in get_leaf_range(tree_data, index)
        idx = indices[i]
        data_reordered[i] = data[idx]
        # Saves the inverse n
        indices_reordered[i] = idx
    end
end

# Checks the distance function and add those points that are among the k best.
# Uses a heap for fast insertion.
@inline function add_points_knn!(best_dists::Vector, best_idxs::Vector{Int},
                                 tree::NNTree, index::Int, point::AbstractVector,
                                 do_end::Bool, skip::F) where {F}
    for z in get_leaf_range(tree.tree_data, index)
        @POINT 1
        idx = tree.reordered ? z : tree.indices[z]
        dist_d = evaluate(tree.metric, tree.data[idx], point, do_end)
        if dist_d <= best_dists[1]
            if skip(tree.indices[z])
                continue
            end

            best_dists[1] = dist_d
            best_idxs[1] = idx
            percolate_down!(best_dists, best_idxs, dist_d, idx)
        end
    end
end

# Add those points in the leaf node that are within range.
# TODO: If we have a distance function that is incrementally increased
# as we sum over the dimensions (like the Minkowski norms) then we could
# stop computing the distance function as soon as we reach the desired radius.
# This will probably prevent SIMD and other optimizations so some care is needed
# to evaluate if it is worth it.
@inline function add_points_inrange!(idx_in_ball::Vector{Int}, tree::NNTree,
                                     index::Int, point::AbstractVector, r::Number, do_end::Bool)
    for z in get_leaf_range(tree.tree_data, index)
        @POINT 1
        idx = tree.reordered ? z : tree.indices[z]
        dist_d = evaluate(tree.metric, tree.data[idx], point, do_end)
        if dist_d <= r
            push!(idx_in_ball, idx)
        end
    end
end

# Add all points in this subtree since we have determined
# they are all within the desired range
function addall(tree::NNTree, index::Int, idx_in_ball::Vector{Int})
    tree_data = tree.tree_data
    @NODE 1
    if isleaf(tree.tree_data.n_internal_nodes, index)
        for z in get_leaf_range(tree.tree_data, index)
            @POINT_UNCHECKED 1
            idx = tree.reordered ? z : tree.indices[z]
            push!(idx_in_ball, idx)
        end
        return
    else
        addall(tree, getleft(index), idx_in_ball)
        addall(tree, getright(index), idx_in_ball)
    end
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

# compute graph related operators

################################################################################
# SymbolicNode attribute getter and setter
################################################################################

macro _list_symbol_info(self, func_name)
  quote
    ref_sz    = Ref{MX_uint}(0)
    ref_names = Ref{char_pp}(0)
    @mxcall($func_name, (MX_handle, Ref{MX_uint}, Ref{char_pp}),
            $(esc(self)), ref_sz, ref_names)
    narg = ref_sz[]
    names = unsafe_wrap(Array, ref_names[], narg)
    names = [Symbol(unsafe_string(x)) for x in names]
    return names
  end
end

"""
    list_arguments(s::SymbolicNode)

List all the arguments of this node. The argument for a node contains both
the inputs and parameters. For example, a `FullyConnected` node will
have both data and weights in its arguments. A composed node (e.g. a MLP) will
list all the arguments for intermediate nodes.

Returns a list of symbols indicating the names of the arguments.
"""
list_arguments(s::SymbolicNode) = @_list_symbol_info(s, :MXSymbolListArguments)

"""
    list_outputs(s::SymbolicNode)

List all the outputs of this node.

Returns a list of symbols indicating the names of the outputs.
"""
list_outputs(s::SymbolicNode) = @_list_symbol_info(s, :MXSymbolListOutputs)

"""
    list_auxiliary_states(s::SymbolicNode)

List all auxiliary states in the symbool.

Auxiliary states are special states of symbols that do not corresponds to an argument,
and do not have gradient. But still be useful for the specific operations.
A common example of auxiliary state is the moving_mean and moving_variance in BatchNorm.
Most operators do not have Auxiliary states.

Returns a list of symbols indicating the names of the auxiliary states.
"""
list_auxiliary_states(s::SymbolicNode) =
  @_list_symbol_info(s, :MXSymbolListAuxiliaryStates)

"""
    get_internals(s::SymbolicNode)

Get a new grouped `SymbolicNode` whose output contains all the internal outputs of
this `SymbolicNode`.
"""
function get_internals(s::SymbolicNode)
  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXSymbolGetInternals, (MX_handle, Ref{MX_handle}), s, ref_hdr)
  return SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end

"""
    get_children(x::SymbolicNode)

Gets a new grouped `SymbolicNode` whose output contains inputs to output
nodes of the original symbol.

```julia
julia> x, y = @mx.var x y
(SymbolicNode x, SymbolicNode y)

julia> z = x + y
SymbolicNode _plus0

julia> z |> mx.get_children |> mx.list_outputs
2-element Array{Symbol,1}:
 :x
 :y
```
"""
function get_children(x::SymbolicNode)
  hdl = Ref{MX_handle}(C_NULL)
  @mxcall(:MXSymbolGetChildren, (MX_handle, Ref{MX_handle}), x, hdl)
  sym = hdl[] |> MX_SymbolHandle |> SymbolicNode
  isempty(list_outputs(sym)) ? nothing : sym
end

"""
    get_attr(s::SymbolicNode, key::Symbol)

Get attribute attached to this `SymbolicNode` belonging to key.

Returns the value belonging to key as a `String`.
If not available, returns `missing`.
"""
function get_attr(s::SymbolicNode, key::Symbol)
  key_s = string(key)
  ref_out = Ref{Cstring}()
  ref_success = Ref{Cint}(-1)
  @mxcall(:MXSymbolGetAttr, (MX_handle, Cstring, Ref{Cstring}, Ref{Cint}),
          s, key_s, ref_out, ref_success)
  if ref_success[] == 1
    unsafe_string(ref_out[])
  else
    missing
  end
end

"""
    list_attr(s::SymbolicNode)

Get all attributes from a symbol.

Returns a dictionary of attributes.
"""
function list_attr(s::SymbolicNode)
  ref_sz    = Ref{MX_uint}(0)
  ref_strings = Ref{char_pp}(0)
  @mxcall(:MXSymbolListAttrShallow, (MX_handle, Ref{MX_uint}, Ref{char_pp}),
            s, ref_sz, ref_strings)
  narg = 2*ref_sz[]
  strings = unsafe_wrap(Array, ref_strings[], narg)
  out = Dict{Symbol, String}()
  for i in 1:2:narg
    key = Symbol(unsafe_string(strings[i]))
    value = unsafe_string(strings[i+1]) # Creates a copy of string
    out[key] = value
  end
  return out
end

"""
    list_all_attr(s::SymbolicNode)

Get all attributes from the symbol graph.

Returns a dictionary of attributes.
"""
function list_all_attr(s::SymbolicNode)
  ref_sz    = Ref{MX_uint}(0)
  ref_strings = Ref{char_pp}(0)
  @mxcall(:MXSymbolListAttr, (MX_handle, Ref{MX_uint}, Ref{char_pp}),
            s, ref_sz, ref_strings)
  narg = 2*ref_sz[]
  strings = unsafe_wrap(Array, ref_strings[], narg)
  out = Dict{Symbol, String}()
  for i in 1:2:narg
    key = Symbol(unsafe_string(strings[i]))
    value = unsafe_string(strings[i+1])
    out[key] = value
  end
  return out
end

"""
    set_attr(s::SymbolicNode, key::Symbol, value::AbstractString)

Set the attribute key to value for this `SymbolicNode`.

!!! note
    It is encouraged not to call this function directly, unless you know exactly what you are doing. The
    recommended way of setting attributes is when creating the `SymbolicNode`. Changing
    the attributes of a `SymbolicNode` that is already been used somewhere else might
    cause unexpected behavior and inconsistency.
"""
function set_attr(s::SymbolicNode, key::Symbol, value::AbstractString)
  key_s = string(key)
  value_s = String(value)

  @mxcall(:MXSymbolSetAttr, (MX_handle, Cstring, Cstring), s, key_s, value_s)
end

"""
    get_name(s::SymbolicNode)

Get the name of the symbol.

    julia> x = mx.Variable(:data)
    julia> mx.get_name(x)
    :data

    julia> y = mx.FullyConnected(x, num_hidden = 128)
    julia> mx.get_name(y)
    :fullyconnected0
"""
function get_name(s::mx.SymbolicNode)
    name = Ref{mx.char_p}(C_NULL)
    success = Ref(0)
    @mxcall(:MXSymbolGetName, (MX_handle, Ref{char_p}, Ref{Int}), s.handle.value, name, success)
    @assert success[] != -1

    str = name[]
    if str == C_NULL  # e.g. the symbol returned via get_internals
        string(s.handle.value)
    else
        Symbol(unsafe_string(str))
    end
end

################################################################################
# Atomic SymbolicNode functions dynamically imported from libmxnet
################################################################################

@inline function _create_atomic_symbol(creator::MX_handle, keys::Vector{String},
                                       vals::Vector{String})
  ref_sym_hdr = Ref{MX_handle}(C_NULL)
  @mxcall(:MXSymbolCreateAtomicSymbol,
          (MX_handle, MX_uint, Ptr{char_p}, Ptr{char_p}, Ref{MX_handle}),
          creator, length(keys), keys, vals, ref_sym_hdr)
  SymbolicNode(MX_SymbolHandle(ref_sym_hdr[]))
end

@inline function _create_atomic_symbol(creator::MX_handle, keys::Vector{String},
                                       vals::Vector{String},
                                       attrs::Dict{Symbol, String})
  node = _create_atomic_symbol(creator, keys, vals)
  # set attrs
  for (k, v) in attrs
    set_attr(node, k, v)
  end
  node
end

function _define_atomic_symbol_creator(name::String)
  handle = _get_libmx_op_handle(name)
  f_desc, key_narg = _get_libmx_op_description(name, handle)

  f_desc *= "* `name::Symbol`: The name of the `SymbolicNode`. (e.g. `:my_symbol`), optional.\n"
  f_desc *= "* `attrs::Dict{Symbol,String}`: The attributes associated with this `SymbolicNode`.\n\n"

  func_name = Symbol(name)
  import_expr = _import_expr(func_name)

  func_def = quote
  function $func_name(::Type{SymbolicNode}, args::SymbolicNode...; name = "", kwargs...)

    # NOTE: hacky way of solving the problem that the arguments of `dot` should be swapped
    # See https://github.com/dmlc/MXNet.jl/issues/55
    if $name == "dot"
      args = reverse(args)
    end

    # NOTE: hacky way of solving the semantic difference of the axes parameter in Julia
    # and in libmxnet.
    # See https://github.com/dmlc/MXNet.jl/pull/123
    if $name == "transpose"
      kwargs = Any[key != :axes ? (key, arg) : (key, reverse(map(i->length(arg)-i, arg))) for (key, arg) in kwargs]
    end

    param_keys = String[]
    param_vals = String[]
    symbol_kws = Dict{Symbol,SymbolicNode}()
    attrs = Dict{Symbol,String}()

    $(if key_narg != ""
      quote
        if !in($key_narg, param_keys)
          push!(param_keys, $key_narg)
          push!(param_vals, string(length(args)))
        end
      end
    end)

    for (k,v) in kwargs
      if k == :name; continue; end
      if isa(v, SymbolicNode)
        symbol_kws[k] = v
      elseif k == :attrs
        if isa(v, Dict)
          attrs = convert(Dict{Symbol, String}, v)
        else
          throw(ArgumentError("attrs needs to be a Dictionary"))
        end
      else
        push!(param_keys, string(k))
        push!(param_vals, dump_mx_param(v))
      end
    end

    if length(args) > 1 && length(symbol_kws) != 0
      @assert(false, $name * " only accepts SymbolicNode either as positional or keyword arguments with optional positional `data` argument, not both.")
    end
    $(if key_narg != ""
      quote
        if length(symbol_kws) > 0
          @assert(false, $name * " takes variable number of SymbolicNode arguments, " *
                         "please pass input Symbols via positional arguments, instead of keyword arguments.")
        end
      end
    end)

    local op = _get_cached_libmx_op_handle($name)
    node = _create_atomic_symbol(op.value, param_keys, param_vals, attrs)

    # generate a new name for the new symbol if user not provided in kwargs
    hint = lowercase($name)
    name = get!(DEFAULT_NAME_MANAGER, name, hint)

    if length(symbol_kws) == 0
      _compose!(node, name, args...)
    elseif length(args) == 1
      _compose!(node; name=name, data=args[1], symbol_kws...)
    else
      _compose!(node; name=name, symbol_kws...)
    end

    return node
  end # function
  end # quote

  func_def2 = quote
  @doc $f_desc
  function $func_name(args::SymbolicNode...; kwargs...)
    $func_name(SymbolicNode, args...; kwargs...)
  end # function
  end # quote

  return quote
    $import_expr
    $func_def
    $func_def2
  end
end

macro _import_atomic_symbol_creators()
  # NOTE: those are operators defined for NDArray, we exclude them here
  # because the calling convention for the type signature is not strong
  # enough to disambiguate the method for NDArray and SymbolicNode
  ignored_ops = ("_set_value", "reshape")  # in lowercase

  op_names = _get_libmx_op_names()
  func_exprs = map(op_names) do name
    if lowercase(name) ∉ ignored_ops
      expr = _define_atomic_symbol_creator(name)
    end
  end

  esc(quote
    $(func_exprs...)
  end)
end

@_import_atomic_symbol_creators

################################################################################
# Utility macros to chain up symbols
################################################################################

macro chain(layers)
    exprs = []
    last_layer = nothing

    function _chain_layer(layer, last_layer)
        if last_layer ≡ nothing
            return esc(layer)
        else
            if @capture(layer, f_(x__))
                x′ = esc.(x)
                return :($f($last_layer, $(x′...)))
            else
                throw(AssertionError("$layer is not a valid function call and cannot be chained."))
            end
        end
    end

    while true
        if @capture(layers, l1_=>l2_)
            new_layer = gensym()
            push!(exprs, :($new_layer = $(_chain_layer(l1, last_layer))))
            last_layer = new_layer
            layers = l2
        else
            push!(exprs, _chain_layer(layers, last_layer))
            break
        end
    end
    Expr(:block, exprs...)
end

################################################################################
# compose
################################################################################

function _compose!(node::SymbolicNode; kwargs...)
  name     = char_p(C_NULL)
  arg_keys = AbstractString[]  # FIXME: can it be String[] ?
  arg_vals = MX_handle[]

  for (k, v) in kwargs
    if k == :name
      name = string(v)
    else
      @assert(isa(v, SymbolicNode), "Compose expect `SymbolicNode` as arguments")
      push!(arg_keys, string(k))
      push!(arg_vals, v)
    end
  end

  @mxcall(:MXSymbolCompose,
          (MX_handle, char_p, MX_uint, Ptr{char_p}, Ptr{MX_handle}),
          node, name, length(arg_keys), arg_keys, arg_vals)
  node
end

_compose!(node::SymbolicNode, args::SymbolicNode...) =
  _compose!(node, char_p(0), args...)

function _compose!(node::SymbolicNode, name::Union{Symbol, char_p}, args::SymbolicNode...)
  if name isa Symbol
    name = string(name)
  end
  arg_keys = Ptr{char_p}(C_NULL)
  arg_vals = MX_handle[args...]

  @mxcall(:MXSymbolCompose,
          (MX_handle, char_p, MX_uint, Ptr{char_p}, Ptr{MX_handle}),
          node, name, length(arg_vals), arg_keys, arg_vals)
  node
end