text:
# ------------------------------------------------------------------------------------------
# # Learning with a single neuron using Flux.jl
#
# In this notebook, we'll use `Flux` to create a single neuron and teach it to learn, as we
# did by hand in notebook 10!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Read in data and process it
#
# Let's start by reading in our data
# ------------------------------------------------------------------------------------------

using CSV
using TextParse
using DataFrames

applecols, applecolnames = TextParse.csvread("data/Apple_Golden_1.dat", '\t')
bananacols, bananacolnames = TextParse.csvread("data/bananas.dat", '\t')

apples = DataFrame(Dict(strip(name)=>col for (name, col) in zip(applecolnames, applecols)))
bananas = DataFrame(Dict(strip(name)=>col for (name, col) in zip(bananacolnames, bananacols)));

# ------------------------------------------------------------------------------------------
# and processing it to extract information about the red and green coloring in our images:
# ------------------------------------------------------------------------------------------

col1 = :red
col2 = :green

x_apples  = [ [apples[i, col1], apples[i, col2]] for i in 1:size(apples)[1] ]
x_bananas = [ [bananas[i, col1], bananas[i, col2]] for i in 1:size(bananas)[1] ]

xs = vcat(x_apples, x_bananas)

ys = vcat( zeros(size(x_apples)[1]), ones(size(x_bananas)[1]) );

# ------------------------------------------------------------------------------------------
# The input data is in `xs` and the labels (true classifications as bananas or apples) in
# `ys`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Using `Flux.jl`
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Now we can load `Flux` to really get going!
# ------------------------------------------------------------------------------------------

using Flux

# ------------------------------------------------------------------------------------------
# We saw in the last notebook that σ is a built-in function in `Flux`.
#
# Another function that is used a lot in neural networks is called `ReLU`; in Julia, the
# function is called `relu`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 1
#
# Use the docs to discover what `ReLU` is all about.
#
# `relu.([-3, 3])` returns
#
# A) [-3, 3]

# B) [0, 3]

# C) [0, 0]

# D) [3, 3]

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

# ------------------------------------------------------------------------------------------
# ### Making a single neuron in Flux
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Let's use `Flux` to build our neuron with 2 inputs and 1 output:
#
#
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We previously put the two weights in a vector, $\mathbf{w}$. Flux instead puts weights in
# a $1 \times 2$ matrix (i.e. a matrix with 1 *row* and 2 *columns*).
#
# Previously, to compute the dot product of $\mathbf{w}$ and $\mathbf{x}$ we had to use
# either the `dot` function, or we had to transpose the vector $\mathbf{w}$:
#
# ```julia
# # transpose w
# b = w' * x
# # or use dot!
# b = dot(w, x)
# ```
# If the weights are instead stored in a $1 \times 2$ matrix, `W`, then we can simply
# multiply `W` and `x` together instead!
#
# We start off with random values for our parameters and data:
# ------------------------------------------------------------------------------------------

W = rand(1, 2)

x = rand(2)

# ------------------------------------------------------------------------------------------
# Note that the product of `W` and `x` will now be an array (vector) with a single element,
# rather than a single number:
# ------------------------------------------------------------------------------------------

W * x

# ------------------------------------------------------------------------------------------
# This means that our bias `b` is treated as an array when we're using `Flux`:
# ------------------------------------------------------------------------------------------

b = rand(1)

# ------------------------------------------------------------------------------------------
# #### Exercise 2
#
# Write a function `mypredict` that will take a single input, array `x` and use `W`, `b`,
# and built-in `σ` to generate an output prediction (stored as an array). This function
# defines our neural network!
#
# Hint: This function will look very similar to $f_{\mathbf{w},\mathbf{b}}$ from the last
# notebook but has changed since our data structures to store our parameters have changed!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 3
#
# Define a loss function called `loss`.
#
# `loss` should take two inputs: a vector storing data, `x`, and a vector storing the
# correct "labels" for that data. `loss` should return the sum of the squares of differences
# between the predictions and the correct labels.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Calculating gradients using Flux: backpropagation
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# For learning, we know that what we need is a way to calculate derivatives of the `loss`
# function with respect to the parameters `W` and `b`. So far, we have been doing that using
# finite differences.
#
# `Flux.jl` instead implements a numerical method called **backpropagation** that calculates
# gradients (essentially) exactly, in an automatic way, by indirectly applying the rules of
# calculus.
# To do so, it provides a new type of object called "tracked" arrays. These are arrays that
# store not only their current value, but also information about gradients, which is used by
# the backpropagation method.
#
# [If you want to understand the maths behind backpropagation, we recommend e.g. [this
# lecture](https://www.youtube.com/watch?v=i94OvYb6noo).]
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# To do so, `Flux` provides a function `param` to define such objects that will contain the
# information for a *param*eter.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Let's start, as usual, by setting up some random initial values for the parameters:
# ------------------------------------------------------------------------------------------

W_data = rand(1, 2)
b_data = rand(1)

W_data, b_data

# ------------------------------------------------------------------------------------------
# We now set up `Flux.jl` objects that will contain these values *and* their derivatives,
# and allow to propagate
# this information around:
# ------------------------------------------------------------------------------------------

W = param(W_data)
b = param(b_data)

# ------------------------------------------------------------------------------------------
# Here, `param` is a function that `Flux` provides to create an object that represents a
# parameter of a machine learning model, i.e. an object which has both a value and
# derivative information, and such that other objects know how to *keep track* of when it is
# used in an expression.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 4
#
# What type does `W` have?
#
# A) Array (1D)

# B) Array (2D)

# C) TrackedArray (1D)

# D) TrackedArray (2D)

# E) Parameter (1D)

# F) Parameter (2D)

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

# ------------------------------------------------------------------------------------------
# #### Exercise 5
#
# `W` stores not only its current value, but also has space to store gradient information.
# You can access the values and gradient of the weights as follows:
#
# ```julia
# W.data
# W.grad
# ```
#
# At this point, are the values of the weights or the gradient of the weights larger?
#
# A) the values of the weights

# B) the gradient of the weights
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 6
#
# For data `x` and `y` where
#
# ```julia
# x, y = [0.413759, 0.692204], [0.845677]
# ```
# apply the loss function to `x` and `y` to give a new variable `l`. What is the type of
# `l`? (How many dimensions does it have?)
#
# A) Array (0D)

# B) Array (1D)

# C) TrackedArray (0D)

# D) TrackedArray (1D)

# E) Float64

# F) Int64

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

# ------------------------------------------------------------------------------------------
# ### Stochastic gradient descent
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can now use these features to reimplement stochastic gradient descent, following the
# method we used in the previous notebook, but now using backpropagation!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 7
#
# Modify the code from the previous notebook for stochastic gradient descent to use Flux
# instead.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Investigating stochastic gradient descent
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Let's look at the values stored in `b` before we run stochastic gradient descent:
# ------------------------------------------------------------------------------------------

b

# ------------------------------------------------------------------------------------------
# After running `stochastic_gradient_descent`, we find the following:
# ------------------------------------------------------------------------------------------

W_final, b_final = stochastic_gradient_descent(loss, W, b, xs, ys, 100000)

# ------------------------------------------------------------------------------------------
# we can look at the values of `W_final` and `b_final`, which our machine learned to
# generate our desired classification.
# ------------------------------------------------------------------------------------------

W_final

b_final

# ------------------------------------------------------------------------------------------
# #### Exercise 8
#
# Plot the data and the learned function.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 9
#
# Do this plot every so often as the learning process is proceeding in order to have an
# animation of the process.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Automation with Flux.jl
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We will need to repeat the above process for a lot of different systems.
# Fortunately, `Flux.jl` provides us with tools to automate this!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Flux allows to create a neuron in a simple way:
# ------------------------------------------------------------------------------------------

using Flux

model = Dense(2, 1, σ)

# ------------------------------------------------------------------------------------------
# The `2` and `1` refer to the number of inputs and outputs, and the neuron is defined using
# the $\sigma$ function.
# ------------------------------------------------------------------------------------------

typeof(model)

# ------------------------------------------------------------------------------------------
# We have made an object of type `Dense`, defined by `Flux`, with the name `model`. This
# represents a "dense neural network layer" (see later for more on neural network layers).
# The parameters that will be modified during the learning process live *inside* the `model`
# object.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 10
#
# Investigate which variables live inside the `model` object and what type they are. How
# does that compare to the call to create the `Dense` object that we started with?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Model object as a function
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can apply the `model` object to data just as if it were a standard function:
# ------------------------------------------------------------------------------------------

model(rand(2))

# ------------------------------------------------------------------------------------------
# #### Exercise 11
#
# Prove to yourself that you understand what is going on when we call `model`. Create two
# arrays `W` and `b` with the same elements as `model.W` and `model.b`. Use `W` and `b` to
# generate the same answer that you get when we call `model([.5, .5])`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Using Flux
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We now need to provide Flux with three pieces of information:
#
# 1. A loss function
# 2. Some training data
# 3. An optimization method
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Loss functions
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Flux has various loss functions built in, for example the mean-squared error (`mse`) that
# we have been using:
# ------------------------------------------------------------------------------------------

loss(x, y) = Flux.mse(model(x), y)

# ------------------------------------------------------------------------------------------
# Another common one is the cross entropy, `Flux.crossentropy`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Data
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# The data can take a couple of different forms.
# One form is a single **iterator**, consisting of pairs $(x, y)$ of data and labels.
# We can achieve this with `zip`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 12
#
# Use `zip` to "zip together" `xs` and `ys`. Then use the `collect` function to check what
# `zip` actually does.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Optimization routine
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Now we need to tell Flux what kind of optimization routine to use. It has several built
# in; the standard stochastic gradient descent algorithm that we have been using is called
# `SGD`. We must pass it two things: a list of parameter objects which will be modified by
# the optimization routine, and a step size:
# ------------------------------------------------------------------------------------------

opt = SGD([model.W, model.b], 0.01)
# give a list of the parameters that will be modified

# ------------------------------------------------------------------------------------------
# The gradient calculations and parameter updates will be carried out by this optimizer
# function; we do not see those details, but if you are curious, you can, of course, look at
# the `Flux.jl` source code!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Training
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We now have all the pieces in place to actually **train** our model (a single neuron) on
# the data.
# "Training" refers to using pre-labeled data to learn the function that relates the input
# data to the desired output data given by the labels.
#
# `Flux` provides the function `train!`, which performs a single pass through the data and
# does a single step of optimization using the partial cost function for each data point:
# ------------------------------------------------------------------------------------------

Flux.train!(loss, data, opt)

# ------------------------------------------------------------------------------------------
# We can then just repeat this several times to train the network more and coax it towards
# the minimum of the cost function:
# ------------------------------------------------------------------------------------------

for i in 1:100
    Flux.train!(loss, data, opt)
end

# ------------------------------------------------------------------------------------------
# Now let's look at the parameters after training:
# ------------------------------------------------------------------------------------------

model.W

model.b

# ------------------------------------------------------------------------------------------
# Instead of writing out a list of parameters to modify, `Flux` provides the function
# `params`, which extracts all available parameters from a model:
# ------------------------------------------------------------------------------------------

opt = SGD(params(model), 0.01)

params(model)

# ------------------------------------------------------------------------------------------
# ## Adding more features
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 13
#
# So far we have just used two features, red and green.
#
# (i) Add a third feature, blue. Plot the new data.
#
# (ii) Train a neuron with 3 inputs and 1 output on the data.
#
# (iii) Can you find a good way to visualize the result?
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/14. Tools - Learning with a single neuron using Flux.jl.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
Nemo.fit!(::fmpq_rel_series, Int) = nothing
Nemo.fit!(::fmpq_abs_series, Int) = nothing

@doc Markdown.doc"""
    integral(f::RelSeriesElem{T}) -> RelSeriesElem

Return the integral of the power series $f$.
"""
function Nemo.integral(f::RelSeriesElem{T}) where T
  g = parent(f)()
  fit!(g, precision(f)+1)
  set_precision!(g, precision(f)+1)
  v = valuation(f)
  Nemo.set_valuation!(g, v+1)
  for i=0:Nemo.pol_length(f)
    c = Nemo.polcoeff(f, i)
    if !iszero(c)
      setcoeff!(g, i, divexact(c, i+v+1))
    end
  end
  Nemo.renormalize!(g)
  return g
end

function *(f::PolyElem{<:SeriesElem{qadic}}, g::PolyElem{<:SeriesElem{qadic}})
  if degree(f) > 2 &&  degree(g) > 2
    fg = mymul_ks(f, g)
#    @hassert :AbsFact 2 fg == Nemo.mul_classical(f, g)
# This cannot be asserted unfortunately: mymul_ks works in the
# capped abs. prec world, while mul_classical does not:
# f = t*((p+O(p^2) + O(s))), then in the classical world
# f^2 = t^2*p^2+O(p^3), while in the capped world we get 0
# (assuming the cap is at 2)
    return fg
  else
    return Nemo.mul_classical(f, g)
  end
end

#=
function *(f::RelSeriesElem{qadic}, g::RelSeriesElem{qadic})
  return mymul_ks(f, g)
  if pol_length(f) > 2 &&  pol_length(g) > 2
    fg = mymul_ks(f, g)
#    @hassert :AbsFact 2 fg == Nemo.mul_classical(f, g)
# This cannot be asserted unfortunately: mymul_ks works in the
# capped abs. prec world, while mul_classical does not:
# f = t*((p+O(p^2) + O(s))), then in the classical world
# f^2 = t^2*p^2+O(p^3), while in the capped world we get 0
# (assuming the cap is at 2)
    return fg
  else
    return Nemo.mul_classical(f, g)
  end
end
=#

function Base.minimum(::typeof(precision), a::Vector{<:SeriesElem})
  return minimum(map(precision, a))
end

function Base.maximum(::typeof(precision), a::Vector{<:SeriesElem})
  return maximum(map(precision, a))
end
Base.length(a::qadic) = a.length

@inline function coeffraw(q::qadic, i::Int)
  @assert i < length(q)
  return reinterpret(Ptr{fmpz}, q.coeffs)+i*sizeof(Ptr{Int})
end

@inline function coeffraw(q::fmpz_poly, i::Int)
  @assert i < length(q)
  return reinterpret(Ptr{fmpz}, q.coeffs)+i*sizeof(Ptr{Int})
end

@inline function Hecke.setcoeff!(z::fmpz_poly, n::Int, x::Ptr{fmpz})
   ccall((:fmpz_poly_set_coeff_fmpz, Hecke.libflint), Nothing,
                    (Ref{fmpz_poly}, Int, Ptr{fmpz}), z, n, x)
   return z
end

@inline function Hecke.mul!(a::Ref{fmpz}, b::Ref{fmpz}, c::fmpz)
  ccall((:fmpz_mul, Hecke.libflint), Cvoid, (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}),a, b, c)
end
@inline function Hecke.iszero(a::Ref{fmpz})
  return unsafe_load(reinterpret(Ptr{Int}, a))==0
end

#=
function mul!(C::Generic.RelSeries{qadic}, f::Generic.RelSeries{qadic}, g::Generic.RelSeries{qadic})
  return f*g
end
=#

function mymul_ks(f::SeriesElem{qadic}, g::SeriesElem{qadic})
  rf = precision(f)
  rg = precision(g)
  S = parent(f)
  Qq = base_ring(S)
  h = degree(Qq)
  mp = typemax(Int)

  #assume no denominator!!!
  F = Hecke.Globals.Zx()
  fit!(F, 1+rf*(2*h-1))
  for j=0:rf-1
    d = coeff(f, j)
    mp = min(mp, precision(d))
    v = valuation(d)
    if v > 0
      sc = prime(Qq)^v
    end
    for k=0:length(d)-1
      Base.GC.@preserve d setcoeff!(F, (j)*(2*h-1)+k, coeffraw(d, k))
      if v > 0 && length(F)-1 >= (j)*(2*h-1)+k
        #problem: if the new coeff is zero, the length isn't increased
        Base.GC.@preserve F Hecke.mul!(coeffraw(F, (j)*(2*h-1)+k), coeffraw(F, (j)*(2*h-1)+k), sc)
      end
    end
  end

  G = Hecke.Globals.Zx()
  fit!(G, 1+rg*(2*h-1))
  for j=0:rg-1
    d = coeff(g, j)
    mp = min(mp, precision(d))
    v = valuation(d)
    if v > 0
      sc = prime(Qq)^v
    end
    for k=0:length(d)-1
      Base.GC.@preserve d setcoeff!(G, (j)*(2*h-1)+k, coeffraw(d, k))
      if v > 0 && length(G)-1 >= (j)*(2*h-1)+k
        #problem: if the new coeff is zero, the length isn't increased
        Base.GC.@preserve F Hecke.mul!(coeffraw(G, (j)*(2*h-1)+k), coeffraw(G, (j)*(2*h-1)+k), sc)
      else
      end
    end
  end

#  @show density(F), density(G), mp
  FG = mullow(F, G, min(rf, rg)*(2*h-1))

  c = qadic[]
  for j=0:min(rf,rg)-1
    H = Hecke.Globals.Zx()
    for x = 0:2*h-2
      if x + (j)*(2*h-1) < length(FG)
        Base.GC.@preserve FG, setcoeff!(H, x, coeffraw(FG, x + (j)*(2*h-1)))
      end
    end
    push!(c, Qq(H, mp))
    @assert valuation(c[end]) >= 0
  end
  while iszero(c[end]) && length(c) > 1
    pop!(c)
  end
  s = S(c, length(c), min(rf, rg), 0)
  Hecke.renormalize!(s)
  return s
end

function mymul_ks(f::PolyElem{<:SeriesElem{qadic}}, g::PolyElem{<:SeriesElem{qadic}})
  nf = degree(f)
  ng = degree(g)
  rf = minimum(precision, coefficients(f))
  rg = minimum(precision, coefficients(g))
  S = base_ring(f)
  Qq = base_ring(S)
  h = degree(Qq)
  mp = typemax(Int)

  #assume no denominator!!!
  nfg = nf+ng+1

  F = Hecke.Globals.Zx()
  fit!(F, 1+(rf*nfg+nf)*(2*h-1))
  for i=0:degree(f)
    c = coeff(f, i)
    for j=0:rf-1
      d = coeff(c, j)
      mp = min(mp, precision(d))
      v = valuation(d)
      if v > 0
        sc = prime(Qq)^v
      end
      for k=0:length(d)-1
        Base.GC.@preserve d setcoeff!(F, (j*nfg+i)*(2*h-1)+k, coeffraw(d, k))
        if v > 0 && length(F)-1 >= (j*nfg+i)*(2*h-1)+k
          #problem: if the new coeff is zero, the length isn't increased
          Base.GC.@preserve F Hecke.mul!(coeffraw(F, (j*nfg+i)*(2*h-1)+k), coeffraw(F, (j*nfg+i)*(2*h-1)+k), sc)
        end
      end
    end
  end
  G = Hecke.Globals.Zx()
  fit!(G, 1+(rg*nfg+ng)*(2*h-1))
  for i=0:degree(g)
    c = coeff(g, i)
    for j=0:rg-1
      d = coeff(c, j)
      mp = min(mp, precision(d))
      v = valuation(d)
      if v > 0
        sc = prime(Qq)^v
      end
      for k=0:length(d)-1
        Base.GC.@preserve d setcoeff!(G, (j*nfg+i)*(2*h-1)+k, coeffraw(d, k))
        if v > 0 && length(G)-1 >= (j*nfg+i)*(2*h-1)+k
          #problem: if the new coeff is zero, the length isn't increased
          Base.GC.@preserve F Hecke.mul!(coeffraw(G, (j*nfg+i)*(2*h-1)+k), coeffraw(d, k), sc)
        else
        end
      end
    end
  end

#  @show density(F), density(G), mp
  FG = mullow(F, G, min(rf, rg)*nfg*(2*h-1))

  fg = parent(f)()

  for i=0:degree(f)+degree(g)
    c = qadic[]
    for j=0:min(rf,rg)-1
      H = Hecke.Globals.Zx()
      for x = 0:2*h-2
        if x + (j*nfg+i)*(2*h-1) < length(FG)
          Base.GC.@preserve FG, setcoeff!(H, x, coeffraw(FG, x + (j*nfg+i)*(2*h-1)))
        end
      end
      push!(c, Qq(H, mp))
      @assert valuation(c[end]) >= 0
    end
    while iszero(c[end]) && length(c) > 1
      pop!(c)
    end
    s = S(c, length(c), min(rf, rg), 0)
    Hecke.renormalize!(s)
    if !iszero(s)
      setcoeff!(fg, i, s)
    end
  end
  return fg
end

function Nemo.canonical_unit(a::SeriesElem)
  iszero(a) && return one(parent(a))
  v = valuation(a)
  v == 0 && return a
  v > 0 && return shift_right(a, v)
  return shift_left(a, -v)
end

#TODO: this is for rings, not for fields, maybe different types?
function Base.gcd(a::T, b::T) where {T <: SeriesElem}
  iszero(a) && iszero(b) && return a
  iszero(a) && return gen(parent(a))^valuation(b)
  iszero(b) && return gen(parent(a))^valuation(a)
  return gen(parent(a))^min(valuation(a), valuation(b))
end

function Base.lcm(a::T, b::T) where {T <: SeriesElem}
  iszero(a) && iszero(b) && return a
  iszero(a) && return a
  iszero(b) && return b
  return gen(parent(a))^max(valuation(a), valuation(b))
end

function Hecke.ResidueField(S::SeriesRing{T}) where {T <: Nemo.RingElem} #darn nmod/gfp
  k = base_ring(S)
  return k, MapFromFunc(x -> coeff(x, 0), y -> set_precision(S(y), 1), S, k)
end

#TODO: in Nemo, rename to setprecision
#      fix/report series add for different length
function set_precision(a::SeriesElem, i::Int)
  b = deepcopy(a)
  set_precision!(b, i)
  return b
end

# should be Nemo/AA
# TODO: symbols vs strings
#       lift(PolyRing, Series)
#       lift(FracField, Series)
#       (to be in line with lift(ZZ, padic) and lift(QQ, padic)
#TODO: some of this would only work for Abs, not Rel, however, this should be fine here
function Hecke.map_coefficients(f, a::RelSeriesElem; parent::SeriesRing)
  c = typeof(f(coeff(a, 0)))[]
  for i=0:Nemo.pol_length(a)-1
    push!(c, f(Nemo.polcoeff(a, i)))
  end
  b = parent(c, length(c), precision(a), valuation(a))
  return b
end

#=
function Hecke.map_coefficients(f, a::RelSeriesElem)
  d = f(coeff(a, 0))
  T = parent(a)
  if parent(d) == base_ring(T)
    S = T
  else
    S = PowerSeriesRing(parent(d), max_precision(T), string(var(T)), cached = false)[1]
  end
  c = typeof(d)[d]
  for i=1:Nemo.pol_length(a)-1
    push!(c, f(Nemo.polcoeff(a, i)))
  end
  b = S(c, length(c), precision(a), valuation(a))
  return b
end
=#
function lift(R::PolyRing{S}, s::SeriesElem{S}) where {S}
  t = R()
  for x = 0:pol_length(s)
    setcoeff!(t, x, polcoeff(s, x))
  end
  return shift_left(t, valuation(s))
end

function rational_reconstruction(a::SeriesElem; parent::PolyRing = PolynomialRing(base_ring(a), cached = false)[1])
  C = base_ring(a)
  Ct = parent
  t = gen(Ct)
  b = Ct()
  v = valuation(a)
  for i=0:Nemo.pol_length(a)
    setcoeff!(b, i+v, Nemo.polcoeff(a, i))
  end
  return rational_reconstruction(b, t^precision(a))
end

function rational_reconstruction(a::padic)
  return rational_reconstruction(Hecke.lift(a), prime(parent(a), precision(a)))
end

Hecke.gcd_into!(a::PolyElem, b::PolyElem, c::PolyElem) = gcd(b, c)
Base.copy(a::PolyElem) = deepcopy(a)
Base.copy(a::SeriesElem) = deepcopy(a)

Base.inv(f::MapFromFunc) = MapFromFunc(x->preimage(f, x), codomain(f), domain(f))

function Hecke.squarefree_part(a::PolyElem)
  return divexact(a, gcd(a, derivative(a)))
end

ext: jl


text:
using ProtoStructs, Test

@proto TestMe
test_me = @test_nowarn TestMe(A=1, B="2", C=complex(1))

@testset "Construction" begin
    @test TestMe((A=1,)) isa TestMe
    @test TestMe(A=1) == TestMe((A=1,))
end # testset

@testset "Access" begin
    @test test_me.A == 1
    @test test_me.B == "2"
    @test test_me.C == complex(1)
end # testset

@testset "Properties" begin
    @test propertynames( test_me ) == (:A, :B, :C)
end # testset


text:
"""
    AbstractDomain
An abstract super type for any domain type. A domain type `D <: AbstractDomain` must implement the following methods to properly interface `AbstractDomain`.
- `Base.∈(val, ::D)`
- `Base.rand(::D)`
- `Base.length(::D)` that is the number of elements in a discrete domain, and the distance between bounds or similar for a continuous domain
Addtionally, if the domain is used in a dynamic context, it can extend
- `add!(::D, args)`
- `delete!(::D, args)`
where `args` depends on `D`'s structure
"""
abstract type AbstractDomain end

"""
    EmptyDomain
A struct to handle yet to be defined domains.
"""
struct EmptyDomain <: AbstractDomain end

"""
    domain()
Construct an [`EmptyDomain`](@ref).
"""
domain() = EmptyDomain()

"""
    Base.in(value, d <: AbstractDomain)
Fallback method for `value ∈ d` that returns `false`.
"""
Base.in(value, ::D) where D <: AbstractDomain = false

"""
    Base.rand(d <: AbstractDomain)
Fallback method for `length(d)` that return `0`.
"""
Base.length(::D) where D <: AbstractDomain = 0

"""
    Base.isempty(d <: AbstractDomain)
Fallback method for `isempty(d)` that return `length(d) == 0` which default to `0`.
"""
Base.isempty(d::D) where D <: AbstractDomain = length(d) == 0

"""
    domain_size(d <: AbstractDomain)
Fallback method for `domain_size(d)` that return `length(d)`.
"""
domain_size(d::D) where D <: AbstractDomain = length(d)

"""
    get_domain(::AbstractDomain)
Access the internal structure of any domain type.
"""
get_domain(d::D) where {D <: AbstractDomain} = d.domain


text:
using Test
using GeophysicalModelGenerator

# Create 1D dataset with lat/lon/depth
Lat         =   1.0:10.0;
Lon         =   11.0:20.0;
Depth       =   (-20:-11)*km;
Data        =   zeros(size(Lon));
Data_set    =   GeoData(Lat,Lon,Depth,(FakeData=Data,Data2=Data.+1.))     # create GeoData without attributes
@test Data_set.depth[2]==-19km

Depth1      =   (-20.:-11.)*m;            # depth has units of m
Data_set1   =   GeoData(Lat,Lon,Depth1,(FakeData=Data,Data2=Data.+1.))
@test Data_set1.depth[2]==-0.019km
@test Data_set.atts["note"]=="No attributes were given to this dataset" # check whether the default attribute assignment works

Depth2      =   -20.:-11.;              # no units
Data_set2   =   GeoData(Lat,Lon,Depth2,(FakeData=Data,Data2=Data.+1.))
@test Data_set2.depth[2]==-19km

# test that it works if we give a Data array, rather than a NamedTuple, as input (we add a default name)
Data_set3 = GeoData(Lat,Lon,Depth,Data)
@test keys(Data_set3.fields)[1] == :DataSet1

# test that it works if we give a Tuple, rather than a NamedTuple, as input (we add a default name)
Data_set4 = GeoData(Lat,Lon,Depth,(Data,))
@test keys(Data_set4.fields)[1] == :DataSet1

# Throw an error if we supply a Tuple with 2 fields (the user should really supply names in that case)
@test_throws ErrorException GeoData(Lat,Lon,Depth,(Data,Data))

# test assignment of attributes
att_dict = Dict("author"=>"Marcel", "year"=>2021)
Data_set4 = GeoData(Lat,Lon,Depth,(Data,),att_dict)
@test Data_set4.atts["author"]=="Marcel"
@test Data_set4.atts["year"]==2021

# check that an error is thrown if a different size input is given for depth
Depth2      =   (-100:10:1.0)               # no units
@test_throws ErrorException GeoData(Lat,Lon,Depth2,(FakeData=Data,Data2=Data.+1.))

# Convert 1D vector to cartesian structure
Data_cart = convert(ParaviewData,Data_set)

@test Data_cart.x[3] ≈ 6189.685604255086
@test Data_cart.y[3] ≈ 324.3876769792181
@test Data_cart.z[3] ≈ 1421.35608984477

# Create Lon/Lat/Depth and X/Y/Z grids from given numbers or 1D vectors
Lon,Lat,Depth =  LonLatDepthGrid(10:20,30:40,(-10:-1)km);
X,Y,Z         =  XYZGrid(10:20,30:40,(-10:-1)km);
@test size(Lon)==(11, 11, 10)
@test Lat[2,2,2]==31.0
@test size(X)==(11, 11, 10)
@test Y[2,2,2]==31.0

Lon,Lat,Depth   =  LonLatDepthGrid(10:20,30:40,-50km);
X,Y,Z           =  XYZGrid(10:20,30:40,-50km);
@test Lon[2,2] == 11.0
@test X[2,2]   == 11.0

Lon,Lat,Depth   =  LonLatDepthGrid(10,30,(-10:-1)km); # 1D line @ given lon/lat
X,Y,Z           =  XYZGrid(10,30,(-10:-1)km);
@test size(Lon)==(1,1,10)
@test Lat[2]==30.0
@test size(X)==(1,1,10)
@test Y[2]==30.0

# throw an error if a 2D array is passed as input
@test_throws ErrorException LonLatDepthGrid(10:20,30:40,[20 30; 40 50]);
@test_throws ErrorException XYZGrid(10:20,30:40,[20 30; 40 50]);

# Create 3D arrays & convert them
Lon,Lat,Depth   =  LonLatDepthGrid(10:20,30:40,(-10:-1)km); # 3D grid
Data            =   ustrip(Depth);

Data_set1       =   GeoData(Lon,Lat,Depth,(FakeData=Data,Data2=Data.+1.))
@test size(Data_set1.depth)==(11, 11, 10)
@test Data_set1.depth[1,2,3]==-8.0km

# double-check that giving 3D arrays in the wrong ordering triggers a warning message
@test_throws ErrorException("It appears that the lon array has a wrong ordering")    GeoData(Lat,Lon,Depth,(FakeData=Data,Data2=Data.+1.))

# Create 2D arrays & convert them
Lon,Lat,Depth   =   LonLatDepthGrid(10:20,30:40,-50km);
Data            =   ustrip(Depth);
Data_set2       =   GeoData(Lon,Lat,Depth,(FakeData=Data,Data2=Data.+1.))
@test Data_set2.depth[2,2]== -50.0km

# Convert the 2D and 3D arrays to their cartesian counterparts
Data_cart1      = convert(ParaviewData,Data_set1)
@test size(Data_cart1.z)==(11, 11, 10)
@test Data_cart1.z[2,2,2] ≈ 3261.2581739797533

Data_cart2      = convert(ParaviewData,Data_set2)
@test size(Data_cart2.z)==(11, 11, 1)
@test Data_cart2.z[2,2] ≈ 3240.141612908441

# Test projection points (used for map projections)
p1 = ProjectionPoint();
@test  p1.Lat==49.9929
@test  p1.Lon==8.2473
@test  p1.EW ≈ 446048.5158750616
@test  p1.NS ≈ 5.53811274482716e6

p2 = ProjectionPoint(p1.EW,p1.NS,p1.zone,p1.isnorth)
@test p1.EW-p2.EW + p1.NS-p2.NS + p1.zone-p2.zone == 0.0

# Create UTM Data structure
ew          =   422123.0:100:433623.0
ns          =   4.514137e6:100:4.523637e6
depth       =   -5400:250:600
EW,NS,Depth =   XYZGrid(ew, ns, depth);
Data        =   ustrip.(Depth);
Data_set    =   UTMData(EW,NS,Depth,33, true, (FakeData=Data,Data2=Data.+1.))

@test Data_set.EW[3,4,2]==422323.0
@test Data_set.NS[3,4,2]==4.514437e6
@test Data_set.depth[3,4,2]==-5150m
@test Data_set.northern[1] == true
@test Data_set.zone[1] == 33
@test Data_set.atts["note"]=="No attributes were given to this dataset" # check whether the default attribute assignment works

# convert from UTMData -> GeoData
Data_set1 = convert(GeoData, Data_set)
@test Data_set1.lon[20] == 14.099668158564413
@test Data_set1.lat[20] == 40.77470011887963
@test Data_set1.depth[20] == -5400m

# Convert from GeoData -> UTMData
Data_set2 = convert(UTMData, Data_set1)
@test sum(abs.(Data_set2.EW.val-Data_set.EW.val)) < 1e-5

# Test Projection point for negative values
proj = ProjectionPoint(Lat= -2.8, Lon=36)
@test proj.zone == 37
@test proj.isnorth == false

# Convert from GeoData -> UTMData, but for a fixed zone (used for map projection)
proj = ProjectionPoint(Lat= 40.77470011887963, Lon=14.099668158564413)
Data_set3 = Convert2UTMzone(Data_set1, proj)
@test Data_set3.EW.val[100] ≈  432022.99999999994

# Create CartData structure
x        =   0:2:10
y        =   -5:5
z        =   -10:2:2
X,Y,Z    =   XYZGrid(x, y, z);
Data     =   Z
Data_setC =   CartData(X,Y,Z, (FakeData=Data,Data2=Data.+1.))
@test sum(abs.(Data_setC.x.val)) ≈ 2310.0km

# Convert from CartData -> UTMData
Data_set4 = Convert2UTMzone(Data_setC, proj)

# Convert from  UTMData -> CartData with shifting (useful to create a LaMEM model, for example)
using Statistics
Data_set5 = Convert2CartData(Data_set, proj)
@test  Data_set5.x.val[22] == 0.2km
@test  Data_set5.y.val[22] ≈ -9.313225746154785e-13km
@test  Data_set5.z.val[22] == -5.4km

# Convert result back to UTM (convert LaMEM results back to UTM & afterwards to GeoData)
Data_set6 = Convert2UTMzone(Data_set5, proj)
@test sum(Data_set.EW.val-Data_set6.EW.val) ≈ 0.0
@test sum(Data_set.NS.val-Data_set6.NS.val) ≈ 0.0


text:
function write_doc(docformat::Pandoc2HTML, doc, rendered, out_path)
    _, weave_source = splitdir(abspath(doc.source))
    weave_version, weave_date = weave_info()

    # Header is inserted from displayed plots
    header_script = doc.header_script
    self_contained = (header_script ≠ "") ? [] : "--self-contained"

    if haskey(doc.header, "bibliography")
        filt = "--filter"
        citeproc = "pandoc-citeproc"
    else
        filt = []
        citeproc = []
    end

    out_path = get_out_path(doc, out_path, "html")
    cd_back = let d = pwd(); () -> cd(d); end
    cd(dirname(out_path))
    try
        out = basename(out_path)
        highlight_stylesheet = get_highlight_stylesheet(MIME("text/html"), docformat.highlight_theme)
        cmd = `pandoc -f markdown+raw_html -s --mathjax=""
        $filt $citeproc $(docformat.pandoc_options)
        --template $(docformat.template_path)
        -H $(docformat.stylesheet_path)
        $(self_contained)
        -V highlight_stylesheet=$(highlight_stylesheet)
        -V weave_version=$(weave_version)
        -V weave_date=$(weave_date)
        -V weave_source=$(weave_source)
        -V headerscript=$(header_script)
        -o $(out)`
        proc = open(cmd, "r+")
        println(proc.in, rendered)
        close(proc.in)
        proc_output = read(proc.out, String)
    catch
        rethrow() # TODO: just show error content instead of rethrow the err
    finally
        cd_back()
    end

    return out_path
end

function write_doc(docformat::Pandoc2PDF, doc, rendered, out_path)
    if haskey(doc.header, "bibliography")
        filt = "--filter"
        citeproc = "pandoc-citeproc"
    else
        filt = []
        citeproc = []
    end

    out_path = get_out_path(doc, out_path, "pdf")
    cd_back = let d = pwd(); () -> cd(d); end
    cd(dirname(out_path))
    try
        out = basename(out_path)
        cmd = `pandoc -f markdown+raw_tex -s  --pdf-engine=xelatex --highlight-style=tango
         $filt $citeproc $(docformat.pandoc_options)
         --include-in-header=$(docformat.header_template)
         -V fontsize=12pt -o $(out)`
        proc = open(cmd, "r+")
        println(proc.in, rendered)
        close(proc.in)
        proc_output = read(proc.out, String)
    catch
        rethrow()
    finally
        cd_back()
    end

    return out_path
end


text:
export Bag, enable_block!, disable_block!, setcontent!, isenabled
"""
    Bag{N}<:TagBlock{AbstractBlock, N}

A bag is a trivil container, but can
    * `setcontent!(bag, content)`
    * `disable_block!(bag)`
    * `enable_block!(bag)`
"""
mutable struct Bag{N}<:TagBlock{AbstractBlock, N}
    content::AbstractBlock{N}
    mask::Bool
end
Bag(b::AbstractBlock) = Bag(b,true)

YaoBlocks.content(bag) = bag.content
YaoBlocks.chcontent(bag::Bag, content) = Bag(content)
YaoBlocks.mat(bag::Bag{N}) where N = bag.mask ? mat(bag.content) : IMatrix{1<<N}()
YaoBlocks.apply!(reg::AbstractRegister, bag::Bag) = bag.mask ? apply!(reg, bag.content) : reg
YaoBlocks.ishermitian(bag::Bag) = bag.mask ? ishermitian(bag.content) : true
YaoBlocks.isreflexive(bag::Bag) = bag.mask ? isreflexive(bag.content) : true
YaoBlocks.isunitary(bag::Bag) = bag.mask ? isunitary(bag.content) : true
YaoBlocks.occupied_locs(bag::Bag) = bag.mask ? occupied_locs(bag.content) : ()
setcontent!(bag::Bag, content) = (bag.content = content; bag)
disable_block!(b::Bag) = (b.mask = false; b)
enable_block!(b::Bag) = (b.mask = true; b)
isenabled(b::Bag) = b.mask

function YaoBlocks.print_annotation(io::IO, bag::Bag)
    printstyled(io, isenabled(bag) ? "[⊙] " : "[⊗] "; bold=true, color=isenabled(bag) ? :green : :red)
end

function Base.show(io::IO, ::MIME"plain/text", blk::Bag)
    return print_tree(io, blk; title=false, compact=false)
end


text:
# entities immediately available as private members in IIF, but requires other packages for actual use

# only export once the convenience constructors are available along with conditional Interpolations dependency

struct FluxModelsDistribution{ID,OD,P,D<:AbstractArray}
  # shape of the input data
  inputDim::NTuple{ID,Int}
  # shape of the output data
  outputDim::NTuple{OD,Int}
  # actual Flux models
  models::Vector{P}
  # the data used for prediction, must be <: AbstractArray
  data::D
  # shuffle model predictions relative to particle index at each sampling
  shuffle::Base.RefValue{Bool}
  # false for default serialization with model info, set true for separate storage of models
  serializeHollow::Base.RefValue{Bool}
    # # TODO remove requirement and standardize sampler API
    # specialSampler::Function
end

"""
    $TYPEDEF

Generate a `<:SamplableBelief` from a heatmap, e.g. a digital elevation model.

Notes
- Give in heatmap and grid, and object becomes a density function that can also be sampled.
- Sampling can be more nuanced by injecting a hint, or location of interest:
  - Mostly aimed at limiting compute when faced with massive heatmaps, e.g. nav units are 10's but map is ~1e6.
- Density approximation is constructed on Guassian measurement assumption of level set and sigma variation.
- Assume data is on a regular grid on TranslationGroup(2)
  - Assume on early implementation `x_grid, y_grid = domain`
- Serialization currently does not store the hint callback.
- To save space, serialization does not store the internal density, but rather reconstructs at unpacking.

DevNotes:
- Generalize to scalar fields on any Manifold.
- Generalize to vector fields if interpolation is sensible.
- TODO standardize with AliasingScalarSampler see IIF #1341
- TODO store the hint function (at least any easy cases)
"""
struct HeatmapGridDensity{T <: Real, H <: Union{<:Function, Nothing}, B <: Union{ManifoldKernelDensity, BallTreeDensity}}
  """intensity data, on regular grid"""
  data::Matrix{T}
  """domain as grid or locations at which scalar intensity elements exist"""
  domain::Tuple{<:AbstractVector{T},<:AbstractVector{T}}
  """use location hint to focus sampling to specific area of data, requires additional info at `getSample`
      assumed the callback will return _____ NOT ACTIVE YET"""
  hint_callback::H
  """general rule for kernel bandwidths used in construction of density, e.g. 0.7 of domain grid spacing"""
  bw_factor::T
  """density function as samplable representation of the data over the domain"""
  densityFnc::B # TODO change to ::ManifoldKernelDensity{TranslationGroup(2),BallTreeDensity}
end

(hmd::HeatmapGridDensity)(w...;kw...) = hmd.densityFnc(w...;kw...)

function sampleTangent(M::AbstractManifold, hms::HeatmapGridDensity)
  sampleTangent(M, hms.densityFnc)
end

function Base.show(io::IO, x::HeatmapGridDensity{T,H,B}) where {T,H,B}
  printstyled(io, "HeatmapGridDensity{", bold=true, color=:blue)
  println(io)
  printstyled(io, "    T", color=:magenta, bold=true )
  println(io, "      = ", T)
  printstyled(io, "    H", color=:magenta, bold=true )
  println(io, "`int  = ", H )
  printstyled(io, "    B", color=:magenta, bold=true )
  println(io, "      = ", B )
  printstyled(io, " }", color=:blue, bold=true)
  println(io, "(")
  println(io, "  data:       ", size(x.data))
  println(io, "    min/max:    ", round(minimum(x.data),digits=5), " / ", round(maximum(x.data),digits=5))
  println(io, "  domain:     ", size(x.domain[1]), ", ", size(x.domain[2]))
  println(io, "    min/max:    ", round(minimum(x.domain[1]),digits=5), " / ", round(maximum(x.domain[1]),digits=5))
  println(io, "    min/max:    ", round(minimum(x.domain[2]),digits=5), " / ", round(maximum(x.domain[2]),digits=5))
  println(io, "  bw_factor:  ", x.bw_factor)
  print(io, "  ")
  show(io, x.densityFnc)
  nothing
end

Base.show(io::IO, ::MIME"text/plain", x::HeatmapGridDensity) = show(io, x)
Base.show(io::IO, ::MIME"application/prs.juno.inline", x::HeatmapGridDensity) = show(io,x)

##

"""
    $TYPEDEF

Generate a `<:SamplableBelief` by selecing normal (Gaussian) deviation from a Level Set of a heatmap, e.g. known altitude on a digital elevation model (DEM).

Notes
- Give in heatmap and grid, a level set gets generated, and the object becomes a density function for sampling.
- Sampling can be more nuanced by injecting a hint, or location of interest:
  - Mostly aimed at limiting compute when faced with massive heatmaps, e.g. nav units are 10's but map is ~1e6.
- Density approximation is constructed on Guassian measurement assumption of level set and sigma variation.
- Assume data is on a regular grid on TranslationGroup(2)

DevNotes:
- Generalize to scalar fields on any Manifold.
- Generalize to vector fields if interpolation is sensible.

See also: [`HeatmapGridDensity`](@ref), [`ManifoldKernelDensity`](@ref)
"""
struct LevelSetGridNormal{T <: Real, H <: HeatmapGridDensity}
  level::T
  """one sigma value associated with measurement noise of `level` against `data`"""
  sigma::T
  """make samplible region of interest from data be `sigma_scale` from `level`, e.g. 3*sigma."""
  sigma_scale::T
  """HeatmapDensityGrid is used to sample the LevelSet regions of interest"""
  heatmap::H
end

(lsg::LevelSetGridNormal)(w...;kw...) = lsg.densityFnc(w...;kw...)

function sampleTangent(M::AbstractManifold, lsg::LevelSetGridNormal)
  sampleTangent(M, lsg.heatmap.densityFnc)
end

function Base.show(io::IO, x::LevelSetGridNormal{T,H}) where {T,H}
  printstyled(io, "LevelSetGridNormal{", bold=true, color=:blue)
  println(io)
  printstyled(io, "    T", color=:magenta, bold=true )
  println(io, "      = ", T)
  printstyled(io, "    H", color=:magenta, bold=true )
  println(io, "`int  = ", H )
  printstyled(io, " }", color=:blue, bold=true)
  println(io, "(")
  println(io, "  level:      ", x.level)
  println(io, "  sigma:      ", x.sigma)
  println(io, "  sig.scale:  ", x.sigma_scale)
  println(io, "  heatmap:    ")
  show(io, x.heatmap)
  nothing
end

Base.show(io::IO, ::MIME"text/plain", x::LevelSetGridNormal) = show(io, x)
Base.show(io::IO, ::MIME"application/prs.juno.inline", x::LevelSetGridNormal) = show(io,x)

#


text:
module Filterbanks

"""
    Filterbank(length::Int, windows::Vector{Vector{Number}}, windowstarts::Vector{Int})

A representation of a filterbank to be applied to a signal with `length`
samples. The `windows` is a vector of vectors, each representing the nonzero
part of a filter and the `windowstarts` is a vector containing the start of
the nonzero portion of the filter.

The `*` operator is overloaded to apply the filterbank to a signal as if the
filterbank were a matrix and an multiplication was performed, but in a more
optimized way.

See also: [`getfilters`](@ref), [`applyfilters`](@ref)
"""
struct Filterbank
	length::Int
	windows::Vector{Vector{Number}}
	windowstarts::Vector{Int}
end

Base.:*(filterbank::Filterbank, signal::AbstractArray) = begin

	return map((window, winstart)->
	window'*signal[winstart:winstart + length(window) - 1]
	, filterbank.windows, filterbank.windowstarts)

end

"""
    getfilters(filterbank::Filterbank)

Return a vector containing vectors filled with the values of each filter of the
given filterbank, including the zero portion.

See also: [`applyfilters`](@ref)
"""
function getfilters(filterbank::Filterbank)

	return hcat(map((window, winstart)->
		[
			zeros(winstart)
			window
			zeros(filterbank.length - (winstart + length(window)))
		], filterbank.windows, filterbank.windowstarts)...)

end

function filtermesh(filterbank::Filterbank)

	filtermesh = zeros(filterbank.length)

	foreach((window, winstart)->
		filtermesh[winstart:winstart + length(window) - 1] += window
	, filterbank.windows, filterbank.windowstarts)

	return filtermesh

end

"""
    applyfilters(signal, filterbank::Filterbank)

Apply a filterbank to a signal and return the result.

See also: [`getfilters`](@ref)
"""
function applyfilters(signal, filterbank::Filterbank)

	return hcat(map((window, winstart)->
		[
			zeros(winstart)
			window.*signal[winstart:winstart + length(window) - 1]
			zeros(filterbank.length - (winstart + length(window)))
		], filterbank.windows, filterbank.windowstarts)...)

end

export Filterbank, getfilters, filtermesh, applyfilters
end  # module Filterbanks


text:
#=
This file is auto-generated. Do not edit.
=#

"""Data Structure for a proportional reserve product for system simulations."""
mutable struct ProportionalReserve <: Reserve
    name::String
    contributingdevices::Vector{Device}  # devices from which the product can be procured
    timeframe::Float64  # the relative saturation timeframe
    internal::InfrastructureSystemsInternal
end

function ProportionalReserve(name, contributingdevices, timeframe, )
    ProportionalReserve(name, contributingdevices, timeframe, InfrastructureSystemsInternal())
end

function ProportionalReserve(; name, contributingdevices, timeframe, )
    ProportionalReserve(name, contributingdevices, timeframe, )
end

# Constructor for demo purposes; non-functional.

function ProportionalReserve(::Nothing)
    ProportionalReserve(;
        name="init",
        contributingdevices=[ThermalStandard(nothing)],
        timeframe=0.0,
    )
end

"""Get ProportionalReserve name."""
get_name(value::ProportionalReserve) = value.name
"""Get ProportionalReserve contributingdevices."""
get_contributingdevices(value::ProportionalReserve) = value.contributingdevices
"""Get ProportionalReserve timeframe."""
get_timeframe(value::ProportionalReserve) = value.timeframe
"""Get ProportionalReserve internal."""
get_internal(value::ProportionalReserve) = value.internal


text:
function checkReturnedQueryIds(xDf::D, ids::T) where {D<:DataFrames.DataFrame, T<:Vector{String}}

   # Check for symbols that were not found in the DB
   uniqueIds = unique(xDf[:InputId]);

   diffIds = setdiff(ids, uniqueIds);
   if !isempty(diffIds)
      warn("No values found for:", prefix = "WARNING: ");
      foreach(el -> warn(el, prefix = ""), diffIds);
   end

   # Check if there are multiple records per input value
   idsOne2M = @> begin
      xDf
      @select(:InputId, :InputSourceDb)
      unique()
      @by(:InputId, nRows = length(:InputSourceDb))
      @where(:nRows .> 1)
   end

   if !isempty(idsOne2M)
      warn("Some input IDs have multiple mappings", prefix = "WARNING: ");
      # foreach(el -> warn(el, prefix = ""), idsOne2M[:InputId]);
   end

end

function postCheckMirnas(x::DataFrames.DataFrame, check::Symbol, input::Vector{String})
   if DataFrames.nrow(x) > 0
      # Check for multi matches
      # multiMatches = @> begin
      #    x[[:InputId, check]]
      #    unique()
      #    @by(:InputId, nRows = length(names()[2]))
      #    @where(:nRows .> 1)
      # end

      # if (DataFrames.nrow(multiMatches) != 0)
      #    warn("Some input identifiers match to more than one MIMAT accession!");
      #    show(multiMatches);
      # end

      dif = setdiff(input, x[:InputId]);
      if length(dif) != 0
         warn("Some input identifiers were not found in the database.");
         show(dif);
      end

      dif = setdiff(x[:InputId], input);
      if length(dif) != 0
         warn("Some input identifiers have multiple matches.");
         show(dif);
      end

   end
end

max_stars_repo_name: StefanHaunsberger/GeneNameGenieJ

max_issues_repo_name: StefanHaunsberger/GeneNameGenieJ

max_forks_repo_name: StefanHaunsberger/GeneNameGenieJ


text:
@testsuite "random" (AT, eltypes)->begin
    rng = if AT <: AbstractGPUArray
        GPUArrays.default_rng(AT)
    else
        Random.default_rng()
    end
    cpu_rng = Random.default_rng()

    @testset "rand" begin  # uniform
        for T in eltypes, d in (10, (10,10))
            A = AT{T}(undef, d)
            B = copy(A)
            rand!(rng, A)
            rand!(rng, B)
            @test Array(A) != Array(B)

            Random.seed!(rng)
            Random.seed!(rng, 1)
            rand!(rng, A)
            Random.seed!(rng, 1)
            rand!(rng, B)
            @test all(Array(A) .== Array(B))

            if rng != cpu_rng
                rand!(cpu_rng, A)
            end
        end

        A = AT{Bool}(undef, 1024)
        fill!(A, false)
        rand!(rng, A)
        @test true in Array(A)
        fill!(A, true)
        rand!(rng, A)
        @test false in Array(A)
    end

    @testset "randn" begin  # normally-distributed
        # XXX: randn calls sqrt, and Base's sqrt(::Complex) performs
        #      checked type conversions that throw boxed numbers.
        for T in filter(isrealfloattype, eltypes), d in (2, (2,2))
            A = AT{T}(undef, d)
            B = copy(A)
            randn!(rng, A)
            randn!(rng, B)
            @test !any(Array(A) .== Array(B))

            Random.seed!(rng)
            Random.seed!(rng, 1)
            randn!(rng, A)
            Random.seed!(rng, 1)
            randn!(rng, B)
            @test all(Array(A) .== Array(B))

            if rng != cpu_rng
                randn!(cpu_rng, A)
            end
        end
    end
end


text:
function is_unix(os::Symbol)
    if (os==:Windows) return false;
    elseif (os==:Linux) return true;
    elseif (os==:FreeBSD) return true;
    elseif (os==:Darwin) return true;
    else error("unknown operating system")
    end
end

function _os_test(qm,ex,test)
    @assert qm == :?
    @assert isa(ex,Expr)
    @assert ex.head == :(:)
    @assert length(ex.args) == 2
    if test
        return esc(ex.args[1])
    else
        return esc(ex.args[2])
    end
end

macro windows(qm,ex)
    _os_test(qm, ex, OS_NAME===:Windows)
end
macro unix(qm,ex)
    _os_test(qm, ex, is_unix(OS_NAME))
end
macro osx(qm,ex)
    _os_test(qm, ex, OS_NAME===:Darwin)
end
macro linux(qm,ex)
    _os_test(qm, ex, is_unix(OS_NAME) && OS_NAME!==:Darwin)
end

macro windows_only(ex)
    @windows? esc(ex) : nothing
end
macro unix_only(ex)
    @unix? esc(ex) : nothing
end
macro osx_only(ex)
    @osx? esc(ex) : nothing
end
macro linux_only(ex)
    @linux? esc(ex) : nothing
end

# Windows version macros

@windows_only function windows_version()
    verinfo = ccall(:GetVersion, Uint32, ())
    (verinfo & 0xFF, (verinfo >> 8) & 0xFF)
end
@unix_only windows_version() = (0,0)

WINDOWS_XP_VER = (5,1)

macro windowsxp(qm,ex)
    _os_test(qm, ex, OS_NAME===:Windows && windows_version() <= WINDOWS_XP_VER)
end

macro windowsxp_only(ex)
    @windowsxp? esc(ex) : nothing
end

macro non_windowsxp_only(ex)
    @windowsxp? nothing : esc(ex)
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using PseudoPointParametrisationExperiments
using Test

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


text:
"""
qqmap(obs::ClimGrid, ref::ClimGrid, fut::ClimGrid; method="Additive", detrend=true, window::Int=15, rankn::Int=50, thresnan::Float64=0.1, keep_original::Bool=false, interp::Function = Linear(), extrap::Function = Flat())

Quantile-Quantile mapping bias correction. For each julian day of the year (+/- **window** size), a transfer function is estimated through an empirical quantile-quantile mapping.

The quantile-quantile transfer function between **ref** and **obs** is etimated on a julian day (and grid-point) basis with a moving window around the julian day. Hence, for a given julian day, the transfer function is then applied to the **fut** dataset for a given julian day.

**Options**

**method::String = "Additive" (default) or "Multiplicative"**. Additive is used for most climate variables. Multiplicative is usually bounded variables such as precipitation and humidity.

**detrend::Bool = true (default)**. A 4th order polynomial is adjusted to the time series and the residuals are corrected with the quantile-quantile mapping.

**window::Int = 15 (default)**. The size of the window used to extract the statistical characteristics around a given julian day.

**rankn::Int = 50 (default)**. The number of bins used for the quantile estimations. The quantiles uses by default 50 bins between 0.01 and 0.99. The bahavior between the bins is controlled by the interp keyword argument. The behaviour of the quantile-quantile estimation outside the 0.01 and 0.99 range is controlled by the extrap keyword argument.

**thresnan::Float64 = 0.1 (default)**. The fraction is missing values authorized for the estimation of the quantile-quantile mapping for a given julian days. If there is more than **treshnan** missing values, the output for this given julian days returns NaNs.

**keep_original::Bool = false (default)**. If **keep_original** is set to true, the values are set to the original values in presence of too many NaNs.

**interp = Interpolations.Linear() (default)**. When the data to be corrected lies between 2 quantile bins, the value of the transfer function is linearly interpolated between the 2 closest quantile estimation. The argument is from Interpolations.jl package.

**extrap = Interpolations.Flat() (default)**. The bahavior of the quantile-quantile transfer function outside the 0.01-0.99 range. Setting it to Flat() ensures that there is no "inflation problem" with the bias correction. The argument is from Interpolation.jl package.
"""
function qqmap(obs::ClimGrid, ref::ClimGrid, fut::ClimGrid; method::String = "Additive", detrend::Bool = true, window::Int = 15, rankn::Int = 50, thresnan::Float64 = 0.1, keep_original::Bool = false, interp = Linear(), extrap = Flat())

    # Consistency checks # TODO add more checks for grid definition
    @argcheck size(obs[1], 1) == size(ref[1], 1) == size(fut[1], 1)
    @argcheck size(obs[1], 2) == size(ref[1], 2) == size(fut[1], 2)

    # Modify dates (e.g. 29th feb are dropped/lost by default)
    obs = ClimateTools.dropfeb29(obs)
    ref = ClimateTools.dropfeb29(ref)
    fut = ClimateTools.dropfeb29(fut)

    # Remove trend if specified
    if detrend == true
        # Obs
        obs_polynomials = ClimateTools.polyfit(obs)
        obs = obs - ClimateTools.polyval(obs, obs_polynomials)
        # Ref
        ref_polynomials = ClimateTools.polyfit(ref)
        ref = ref - ClimateTools.polyval(ref, ref_polynomials)
        # Fut
        fut_polynomials = ClimateTools.polyfit(fut)
        poly_values = ClimateTools.polyval(fut, fut_polynomials)
        fut = fut - poly_values
    end

    #Get date vectors
    datevec_obs = get_timevec(obs)
    datevec_ref = get_timevec(ref)
    datevec_fut = get_timevec(fut)

    # Julian days vectors
    obs_jul = Dates.dayofyear.(datevec_obs)
    ref_jul = Dates.dayofyear.(datevec_ref)
    fut_jul = Dates.dayofyear.(datevec_fut)

    # Prepare output array (replicating data type of fut ClimGrid)
    dataout = fill(convert(typeof(fut[1].data[1]), NaN), (size(fut[1], 1), size(fut[1], 2), size(fut[1], 3)))::Array{typeof(fut[1].data[1]),T} where T

    if minimum(ref_jul) == 1 && maximum(ref_jul) == 365
        days = 1:365
    else
        days = minimum(ref_jul) + window:maximum(ref_jul) - window
        start = Dates.monthday(minimum(datevec_obs))
        finish = Dates.monthday(maximum(datevec_obs))
    end

    # Reshape. Allows multi-threading on grid points.
    obsin = reshape(obs[1].data, (size(obs[1].data, 1) * size(obs[1].data, 2), size(obs[1].data, 3)))
    refin = reshape(ref[1].data, (size(ref[1].data, 1) * size(ref[1].data, 2), size(ref[1].data, 3)))
    futin = reshape(fut[1].data, (size(fut[1].data, 1) * size(fut[1].data, 2), size(fut[1].data, 3)))
    dataoutin = reshape(dataout, (size(dataout, 1) * size(dataout, 2), size(dataout, 3)))

    # Looping over grid points using multiple-dispatch calls to qqmap
    Threads.@threads for k = 1:size(obsin, 1)
        # for k = 1:size(obsin, 1)

        obsvec = obsin[k,:]
        refvec = refin[k,:]
        futvec = futin[k,:]

        dataoutin[k, :] = qqmap(obsvec, refvec, futvec, days, obs_jul, ref_jul, fut_jul, method = method, detrend = detrend, window = window, rankn = rankn, thresnan = thresnan, keep_original = keep_original, interp = interp, extrap = extrap)

    end

    # Apply mask
    dataout = applymask(dataout, obs.msk)

    # Rebuilding ClimGrid
    lonsymbol = Symbol(fut.dimension_dict["lon"])
    latsymbol = Symbol(fut.dimension_dict["lat"])

    dataout2 = AxisArray(dataout, Axis{lonsymbol}(fut[1][Axis{lonsymbol}].val), Axis{latsymbol}(fut[1][Axis{latsymbol}].val), Axis{:time}(datevec_fut))

    timeattrib = fut.timeattrib
    timeattrib["calendar"] = "365_day"

    globalattribs = fut.globalattribs
    if haskey(globalattribs, "history")
        globalattribs["history"] = string(globalattribs["history"], " - ", "Bias-corrected with ClimateTools.jl")
    else
        globalattribs["history"] = string("Bias-corrected with ClimateTools.jl")
    end

    C = ClimGrid(dataout2; longrid = fut.longrid, latgrid = fut.latgrid, msk = fut.msk, grid_mapping = fut.grid_mapping, dimension_dict = fut.dimension_dict, timeattrib = timeattrib, model = fut.model, frequency = fut.frequency, experiment = fut.experiment, run = fut.run, project = fut.project, institute = fut.institute, filename = fut.filename, dataunits = fut.dataunits, latunits = fut.latunits, lonunits = fut.lonunits, variable = fut.variable, typeofvar = fut.typeofvar, typeofcal = "365_day", varattribs = fut.varattribs, globalattribs = globalattribs)

    if detrend == true
        C = C + poly_values
    end

    return C

end

"""
qqmap(obsvec::Array{N,1}, refvec::Array{N,1}, futvec::Array{N,1}, days, obs_jul, ref_jul, fut_jul; method::String="Additive", detrend::Bool=true, window::Int64=15, rankn::Int64=50, thresnan::Float64=0.1, keep_original::Bool=false, interp=Linear(), extrap=Flat())

Quantile-Quantile mapping bias correction for single vector. This is a low level function used by qqmap(A::ClimGrid ..), but can work independently.

"""
function qqmap(obsvec::Array{N,1} where N, refvec::Array{N,1} where N, futvec::Array{N,1} where N, days, obs_jul, ref_jul, fut_jul; method::String = "Additive", detrend::Bool = true, window::Int64 = 15, rankn::Int64 = 50, thresnan::Float64 = 0.1, keep_original::Bool = false, interp = Linear(), extrap = Flat())

    # range over which quantiles are estimated
    P = range(0.01, stop = 0.99, length = rankn)
    obsP = similar(obsvec, length(P))
    refP = similar(refvec, length(P))
    sf_refP = similar(refvec, length(P))

    # Prepare output array
    dataout = similar(futvec, (size(futvec)))

    # LOOP OVER ALL DAYS OF THE YEAR
    for ijulian in days

        # idx for values we want to correct
        idxfut = (fut_jul .== ijulian)

        # Find all index of moving window around ijulian day of year
        idxobs = ClimateTools.find_julianday_idx(obs_jul, ijulian, window)
        idxref = ClimateTools.find_julianday_idx(ref_jul, ijulian, window)

        obsval = obsvec[idxobs]
        refval = refvec[idxref]
        futval = futvec[idxfut]

        # Random perturbation is added to avoid division by zero
        if lowercase(method) == "multiplicative"
            # Seed for reproducible results
            Random.seed!(42)

            init     = 0.000001
            interval = 0.00001
            stop     = 0.001
            R = init:interval:stop

            obsval = obsval .+ rand(R, length(obsvec[idxobs]))
            refval = refval .+ rand(R, length(refvec[idxref]))
            futval = futval .+ rand(R, length(futvec[idxfut]))
        end

        obsvalnan = isnan.(obsval)
        refvalnan = isnan.(refval)
        futvalnan = isnan.(futval)

        if (sum(obsvalnan) < (length(obsval) * thresnan)) & (sum(refvalnan) < (length(refval) * thresnan)) & (sum(futvalnan) < (length(futval) * thresnan))

            # Estimate quantiles for obs and ref for ijulian
            quantile!(obsP, obsval[.!obsvalnan], P, sorted = false)
            quantile!(refP, refval[.!refvalnan], P, sorted = false)

            if lowercase(method) == "additive" # used for temperature
                sf_refP .= obsP .- refP
                itp = interpolate((refP,), sf_refP, Gridded(interp))
                itp = extrapolate(itp, extrap) # add extrapolation
                # dataout[idxfut] .= itp(futval) .+ futval
                futval .= itp(futval) .+ futval

                futval[futvalnan] .= 0.0

            elseif lowercase(method) == "multiplicative" # used for precipitation
                sf_refP .= obsP ./ refP
                sf_refP[sf_refP .< 0] .= eps(1.0)
                sf_refP[isnan.(sf_refP)] .= 1.0#eps(1.0)#1.0
                itp = interpolate((refP,), sf_refP, Gridded(interp))
                itp = extrapolate(itp, extrap) # add extrapolation
                futval .= itp(futval) .* futval

                futval[futvalnan] .= 0.0

            else
                error("Wrong method")
            end
            # Replace values with new ones
            dataout[idxfut] .= futval
        else

            if keep_original
                # # Replace values with original ones (i.e. too may NaN values for robust quantile estimation)
                dataout[idxfut] .= futval
            else
                # DO NOTHING (e.g. if there is no reference, we want NaNs and not original values)
                dataout[idxfut] .= NaN
            end
        end
    end

    return dataout

end

"""
biascorrect_extremes(obs::ClimGrid, ref::ClimGrid, fut::ClimGrid; method::String="Multiplicative", detrend::Bool=true, window::Int=15, rankn::Int=50, thresnan::Float64=0.1, keep_original::Bool=false, interp=Linear(), extrap=Flat(), gev_params::DataFrame)

Correct the tail of the distribution with a paramatric method, using the parameters μ, σ and ξ contained in gevparams DataFrame.
"""
function biascorrect_extremes(obs::ClimGrid, ref::ClimGrid, fut::ClimGrid; method::String = "Multiplicative", P::Real = 0.95, detrend::Bool = true, window::Int = 15, rankn::Int = 50, thresnan::Float64 = 0.1, keep_original::Bool = false, interp = Linear(), extrap = Flat(), gevparams::DataFrame)

    # Consistency checks # TODO add more checks for grid definition
    @argcheck size(obs[1], 1) == size(ref[1], 1) == size(fut[1], 1)
    @argcheck size(obs[1], 2) == size(ref[1], 2) == size(fut[1], 2)

    qqmap_base = qqmap(obs, ref, fut, method=method, detrend=detrend, window=window, rankn=rankn, thresnan=thresnan, keep_original=keep_original, interp=interp, extrap=extrap)

    # Modify dates (e.g. 29th feb are dropped/lost by default)
    obs = ClimateTools.dropfeb29(obs)
    ref = ClimateTools.dropfeb29(ref)
    fut = ClimateTools.dropfeb29(fut)

    # Remove trend if specified
    if detrend == true
        # Obs
        obs_polynomials = ClimateTools.polyfit(obs)
        obs = obs - ClimateTools.polyval(obs, obs_polynomials)
        # Ref
        ref_polynomials = ClimateTools.polyfit(ref)
        ref = ref - ClimateTools.polyval(ref, ref_polynomials)
        # Fut
        fut_polynomials = ClimateTools.polyfit(fut)
        poly_values = ClimateTools.polyval(fut, fut_polynomials)
        fut = fut - poly_values
    end

    #Get date vectors
    datevec_obs = get_timevec(obs)
    datevec_ref = get_timevec(ref)
    datevec_fut = get_timevec(fut)

    # Condition for a moving window
    movingwindow = maximum(year.(datevec_fut)) - minimum(year.(datevec_fut)) > (maximum(year.(datevec_ref)) - minimum(year.(datevec_ref)))*1.5

    # Prepare output array (replicating data type of fut ClimGrid)
    dataout = fill(convert(typeof(fut[1].data[1]), NaN), (size(fut[1], 1), size(fut[1], 2), size(fut[1], 3)))::Array{typeof(fut[1].data[1]),T} where T

    # Reshape. Allows multi-threading on grid points.
    obsin = reshape(obs[1].data, (size(obs[1].data, 1) * size(obs[1].data, 2), size(obs[1].data, 3)))
    refin = reshape(ref[1].data, (size(ref[1].data, 1) * size(ref[1].data, 2), size(ref[1].data, 3)))
    futin = reshape(fut[1].data, (size(fut[1].data, 1) * size(fut[1].data, 2), size(fut[1].data, 3)))
    dataoutin = reshape(dataout, (size(dataout, 1) * size(dataout, 2), size(dataout, 3)))
    latgrid = reshape(obs.latgrid, (size(obs.latgrid, 1) * size(obs.latgrid, 2)))
    longrid = reshape(obs.longrid, (size(obs.longrid, 1) * size(obs.longrid, 2)))
    qqvecin = reshape(qqmap_base[1].data, (size(qqmap_base[1].data, 1) * size(qqmap_base[1].data, 2), size(qqmap_base[1].data, 3)))

    # Looping over grid points using multiple-dispatch calls to qqmap
    # Threads.@threads for k = 1:size(obsin, 1)
    for k = 1:size(obsin, 1)

        obsvec = obsin[k,:]
        refvec = refin[k,:]
        futvec = futin[k,:]

        obsvalnan = isnan.(obsvec)
        refvalnan = isnan.(refvec)
        futvalnan = isnan.(futvec)

        if (sum(obsvalnan) < (length(obsvec) * thresnan)) & (sum(refvalnan) < (length(refvec) * thresnan)) & (sum(futvalnan) < (length(futvec) * thresnan))

            # Estimate threshold
            thres = ClimateTools.get_threshold(obsvec, refvec, thres = P)

            GPD_obs = ClimateTools.estimate_gpd(obs.latgrid[k], obs.longrid[k], gevparams, thres)

            if movingwindow
                # window_length will be roughly the length of REF
                ref_l = (maximum(year.(datevec_ref)) - minimum(year.(datevec_ref)) + 1)*365
                # in case of multiple members, we approximate ref length
                fut_l = length(futvec)
                w_length = fut_l/(floor(fut_l/ref_l))
                # The fut series is separated is a few time windows.
                # We also make the windows overlap to replicate a moving window
                R = w_length/2:w_length/3:fut_l
                for c in R
                    idxi = round(Int, max(c-w_length/2+1,1))
                    idxf = round(Int, min(c+w_length/2,fut_l))
                    futvec_tmp = futvec[idxi:idxf]

                    # However, we want each moving window to correct a smaller portion of the time series
                    idxi_corr = round(Int, w_length/2-(w_length/2)/3+1)
                    if c == w_length/2
                        idxi_corr = 1
                    end
                    idxf_corr = round(Int, w_length/2+(w_length/2)/3)
                    if c == maximum(w_length/2:w_length/3:fut_l) || idxi+idxf_corr-1 > fut_l
                        idxf_corr = min(round(Int, w_length-w_length/3),fut_l-idxi+1)
                    end
                    # nbyears = round(Int,length(futvec_tmp)/365)
                    # N = nbyears*2 # By default, we correct roughly 2 values per year

                    # Get clusters
                    futclusters = getcluster(futvec_tmp, thres, 1.0)
                    # Estimate GPD parameters of fut clusters
                    GPD_fut = gpdfit(futclusters[:Max], threshold = thres)
                    # Calculate quantile values of maximum clusters values
                    fut_cdf = Extremes.cdf.(GPD_fut, futclusters[:Max])
                    # Estimation of fut values based on provided GPD parameters
                    newfut = quantile.(GPD_obs, fut_cdf)

                    # Linear transition over a quarter of the distance
                    exIDX = futclusters[:Position]

                    target = (maximum(futvec_tmp[exIDX]) - minimum(futvec_tmp[exIDX]))/4.0
                    transition = (futvec_tmp[exIDX] .- minimum(futvec_tmp[exIDX]))/target
                    transition[transition .> 1.0] .= 1.0
                    # Apply linear transition
                    newfut_trans = (newfut .* transition) .+ (qqvecin[k, exIDX] .* (1.0 .- transition))

                    # We put the new data in the bias-corrected vector
                    exIDX_corr = findall(futclusters[:Position] .< idxf_corr)
                    dataoutin[k, idxi+idxi_corr-1:idxi+idxf_corr-1] = qqvecin[k, idxi_corr:idxf_corr]
                    dataoutin[k, idxi .+ futclusters[exIDX_corr,:Position]] = newfut_trans[exIDX_corr]

                end
            else
                # refclusters = getcluster(refvec[k,:], thres, 1.0)
                futclusters = getcluster(futvec, thres, 1.0)
                GPD_fut = gpdfit(futclusters[!,:Max], threshold = thres)
                fut_cdf = Extremes.cdf.(GPD_fut, futclusters[!,:Max])
                newfut = quantile.(GPD_obs, fut_cdf)

                # Linear transition over a quarter of the distance
                exIDX = futclusters[!,:Position]

                target = (maximum(futvec[exIDX]) - minimum(futvec[exIDX]))/4.0
                transition = (futvec[exIDX] .- minimum(futvec[exIDX]))/target
                transition[transition .> 1.0] .= 1.0
                # Apply linear transition
                newfut_trans = (newfut .* transition) .+ (qqvecin[k, exIDX] .* (1.0 .- transition))

                # We put the new data in the bias-corrected vector
                dataoutin[k, :] = qqvecin[k, :]
                dataoutin[k, futclusters[!,:Position]] = newfut_trans

            end

        # else
        #
        #     if keep_original
        #         # # Replace values with original ones (i.e. too may NaN values for robust quantile estimation)
        #         dataoutin[k, idxfut] .= futvec
        #     else
        #         # DO NOTHING (e.g. if there is no reference, we want NaNs and not original values)
        #         dataoutin[k, idxfut] .= NaN
        #     end
        end
    end

    # Apply mask
    dataout = applymask(dataout, obs.msk)

    # Rebuilding ClimGrid
    lonsymbol = Symbol(fut.dimension_dict["lon"])
    latsymbol = Symbol(fut.dimension_dict["lat"])

    dataout2 = AxisArray(dataout, Axis{lonsymbol}(fut[1][Axis{lonsymbol}].val), Axis{latsymbol}(fut[1][Axis{latsymbol}].val), Axis{:time}(datevec_fut))

    timeattrib = fut.timeattrib
    timeattrib["calendar"] = "365_day"

    globalattribs = fut.globalattribs
    if haskey(globalattribs, "history")
        globalattribs["history"] = string(globalattribs["history"], " - ", "Bias-corrected with ClimateTools.jl")
    else
        globalattribs["history"] = string("Bias-corrected for extreme values with ClimateTools.jl")
    end

    C = ClimGrid(dataout2; longrid = fut.longrid, latgrid = fut.latgrid, msk = fut.msk, grid_mapping = fut.grid_mapping, dimension_dict = fut.dimension_dict, timeattrib = timeattrib, model = fut.model, frequency = fut.frequency, experiment = fut.experiment, run = fut.run, project = fut.project, institute = fut.institute, filename = fut.filename, dataunits = fut.dataunits, latunits = fut.latunits, lonunits = fut.lonunits, variable = fut.variable, typeofvar = fut.typeofvar, typeofcal = "365_day", varattribs = fut.varattribs, globalattribs = globalattribs)

    if detrend == true
        C = C + poly_values
    end

    return C

end

"""
gev2gpd(μ, σ, ξ, thres)

Transform GEV parameters to GPD parameters (Coles, 2001)
"""
function gev2gpd(μ, σ, ξ, thres)

    return σ₂ = σ + ξ * (thres - μ)
end

function estimate_gpd(lat, lon, gevparams, thres)

    # Find closest gev parameters
    l1 = (Float64(lat), Float64(lon))
    l2 = (Array(gevparams[!,:lat]), Array(gevparams[!,:lon]))
    dist, idx = findmindist(l1, l2)

    # Extract GEV parameters
    μ = gevparams[idx[1], :mu]
    σ = gevparams[idx[1], :sigma]
    ξ = gevparams[idx[1], :xi]

    # Transform to GPD parameters and build reference GP distribution.
    μ₂ = ClimateTools.gev2gpd(μ, σ, ξ, thres)
    GPD_obs = Extremes.GeneralizedPareto(μ₂, σ, ξ)

    return GPD_obs

end

# """
#     biascorrect_extremes(obs::AxisArray, ref::AxisArray, fut::AxisArray, μ, σ, ξ)
#
# Correct the tail of the distribution with a paramatric method, using the parameters μ, σ and ξ.
# """
# function biascorrect_extremes(obsvec::AxisArray, refvec::AxisArray, futvec::AxisArray, days, μ::Real, σ::Real, ξ::Real)
#
#     obs_jul = dayofyear.(obsvec[Axis{:time}][:])
#     ref_jul = dayofyear.(refvec[Axis{:time}][:])
#     fut_jul = dayofyear.(futvec[Axis{:time}][:])
#
#     # Prepare output array (replicating data type of fut ClimGrid)
#     dataout = fill(convert(typeof(futvec.data[1]), NaN), length(futvec))::Array{typeof(fut[1].data[1]),T} where T
#
#     if minimum(ref_jul) == 1 && maximum(ref_jul) == 365
#         days = 1:365
#     else
#         days = minimum(ref_jul) + window:maximum(ref_jul) - window
#         start = Dates.monthday(minimum(datevec_obs))
#         finish = Dates.monthday(maximum(datevec_obs))
#     end
#
#     threshold = mean([quantile(obsvec[obsvec .>= 1.0], P) quantile(refvec[refvec .>= 1.0], P)])
#
#     idx_obs = obsvec .> threshold
#     idx_ref = refvec .> threshold
#     idx_fut = futvec .> threshold
#
#     dataout .= qqmap(obsvec.data, refvec.data, futvec.data, days, obs_jul, ref_jul, fut_jul, method = "multiplicative", detrend = false)
#
#
# end


text:
export CayleyConfiguration, CircuitTable, MixedCell, circuitdot

const CellIndices = Vector{Tuple{Int,Int}} # Maybe also SVector{N, Tuple{T, T}}

"""
    CayleyConfiguration{I<:Integer}

The Cayley matrix of a set of matrices ``A_1, …, A_n`` where ``A_i ∈ 𝐑^{n×mᵢ}``.
"""
struct CayleyConfiguration{I<:Integer} <: AbstractMatrix{I}
    A::Matrix{I}
    offsets::Vector{Int} # Maybe also SVector{N, Int}
end

function CayleyConfiguration(Aᵢs::Matrix{I}...) where {I<:Integer}
    n = size(Aᵢs[1], 1)
    @assert all(Aᵢ -> size(Aᵢ, 1) == n, Aᵢs) "Matrices do not have the same number of rows."

    A = fill(zero(I), 2n, sum(Aᵢ -> size(Aᵢ, 2), Aᵢs))

    offsets = Int[]
    offset = 0
    for (i, Aᵢ) in enumerate(Aᵢs)
        push!(offsets, offset)
        mᵢ = size(Aᵢ, 2)

        for j=1:mᵢ
            jj = offset + j
            for k=1:n
                A[k, jj] = Aᵢ[k, j]
            end
            A[n+i, jj] = one(I)
        end

        offset += mᵢ
    end
    CayleyConfiguration(A, offsets)
end

"""
    linearindex(CC::CayleyConfiguration, i::Integer, col::Integer)

Get the linear index into the Caylex matrix corresponding to the configuration `i` and column `col`.
"""
function linearindex(CC::CayleyConfiguration, configuration::Integer, col::Integer)
    CC.offsets[configuration] + col
end

"""
    linearindices(CC::CayleyConfiguration, cellindices::CellIndices)

Get the linear indexes into the Caylex matrix corresponding to the cell indices.
"""
function linearindices(A::CayleyConfiguration, cellindices::CellIndices)
    linearindices = Int[]
    for i=1:length(cellindices)
        aᵢ, bᵢ = cellindices[i]
        push!(linearindices, linearindex(A, i, aᵢ), linearindex(A, i, bᵢ))
    end
    linearindices
end

Base.getindex(CC::CayleyConfiguration, i::Int) = getindex(CC.A, i)
Base.getindex(CC::CayleyConfiguration, I::Vararg{Int, 2}) = getindex(CC.A, I...)
Base.size(CC::CayleyConfiguration) = size(CC.A)
Base.size(CC::CayleyConfiguration, I) = size(CC.A, I)


text:
# This file is a part of ROOTFramework.jl, licensed under the MIT License (MIT).

using Cxx

cxxinclude("TObject.h")

import Base: show

export TObject, TObjectPtr, TObjectInst
const TObject = cxxt"TObject"
const TObjectPtr = pcpp"TObject"
const TObjectInst = Union{TObject, TObjectPtr}

tobject_name(x::TObjectInst) = unsafe_string((@cxx x->GetName()))

show(io::IO, x::TObject) =
    print(io, "TObject(\"$(escape_string(tobject_name(x)))\")")

show(io::IO, x::TObjectPtr) =
    print(io, "TObjectPtr(\"$(escape_string(tobject_name(x)))\")")

function dyncast(x::TObjectPtr)
    tstring = icxx"""dynamic_cast<TString*>($x);"""
    !isnullptr(tstring) && return tstring

    tfile = icxx"""dynamic_cast<TFile*>($x);"""
    !isnullptr(tfile) && return tfile

    tdirectory = icxx"""dynamic_cast<TDirectory*>($x);"""
    !isnullptr(tdirectory) && return tdirectory

    tchain = icxx"""dynamic_cast<TChain*>($x);"""
    !isnullptr(tchain) && return tchain

    ttree = icxx"""dynamic_cast<TTree*>($x);"""
    !isnullptr(ttree) && return ttree

    tcanvas = icxx"""dynamic_cast<TCanvas*>($x);"""
    !isnullptr(tcanvas) && return tcanvas

    th1d = icxx"""dynamic_cast<TH1D*>($x);"""
    !isnullptr(th1d) && return th1d
    th1f = icxx"""dynamic_cast<TH1F*>($x);"""
    !isnullptr(th1f) && return th1f
    th1i = icxx"""dynamic_cast<TH1I*>($x);"""
    !isnullptr(th1i) && return th1i

    th2d = icxx"""dynamic_cast<TH2D*>($x);"""
    !isnullptr(th2d) && return th2d
    th2f = icxx"""dynamic_cast<TH2F*>($x);"""
    !isnullptr(th2f) && return th2f
    th2i = icxx"""dynamic_cast<TH2I*>($x);"""
    !isnullptr(th2i) && return th2i

    th3d = icxx"""dynamic_cast<TH3D*>($x);"""
    !isnullptr(th3d) && return th3d
    th3f = icxx"""dynamic_cast<TH3F*>($x);"""
    !isnullptr(th3f) && return th3f
    th3i = icxx"""dynamic_cast<TH3I*>($x);"""
    !isnullptr(th3i) && return th3i
end


text:
using Plots
using Distributions
using Random

1 + 1
a = 2.0

d = Normal(0,1)
data = rand(d,5)
plot(data)


text:
using st16SseqJuliaTools
if abspath(PROGRAM_FILE) == @__FILE__
    map = read_jc_clusters(ARGS[1])
    sizes = Array{Tuple{String,Int64},1}()
    for k in keys(map)
        push!(sizes,(k,length(map[k])))
    end
    sort!(sizes,by=x->x[2],rev=true)
    for p in map[sizes[1][1]]
        println(p)
    end
end