text:
using Polynomials
using Plots
include("funkcje.jl")
return

prec = 0.01
idealne = [(cos(t), sin(t)) for t in -pi:0.001:pi]
return

function makeSpline(coordinates)
    function CalculateCs()
        μ = [NaN64 for k in 1:n]
        α = vcat([0.], [3. * (A[k+1] - A[k])/h[k] - 3. * (A[k] - A[k-1])/h[k-1] for k in 2:n])

        C = [0. for k in 1:nPlus1]
        l = [0. for k in 1:nPlus1]
        z = [0. for k in 1:nPlus1]

        l[1] = 1.
        μ[1] = 0.
        z[1] = 0.

        for k in 2:n
            l[k] = 2 * (xs[k+1] - xs[k-1]) - h[k-1] * μ[k-1]
            μ[k] = h[k] / l[k]
            z[k] = (α[k] - h[k-1]*z[k-1]) / l[k]
        end

        l[nPlus1] = 1.
        z[nPlus1] = 0.
        C[nPlus1] = 0.

        for k in n:(-1):1
            C[k] = z[k] - μ[k]*C[k+1]
        end
        return C
    end
    nPlus1 = length(coordinates)
    n = nPlus1 - 1
    # x1 x2 ... xk+1
    xs = map(pair->pair[1], coordinates)
    # y1 y2 ... yk+1
    # ys = map(pair->pair[2], coordinates)
    A = map(pair->pair[2], coordinates)

    h = [xs[k+1] - xs[k] for k in 1:n]
    C = CalculateCs()

    B = [(A[k+1] - A[k]) / h[k] - h[k]*(C[k+1] + 2. * C[k]) / 3. for k in 1:n]
    D = [(C[k+1] - C[k]) / (3. * h[k]) for k in 1:n]

    outputSet = [Poly(NaN) for k in 1:n]
    intervals = myzip(xs,xs[2:end])

    for k in 1:n
        xk = xs[k]
        polynom = A[k] + B[k]*poly([xk]) + C[k]*(poly([xk])^2) + D[k]*(poly([xk]) ^ 3)
        outputSet[k] = polynom
    end
    return myzip(outputSet, intervals)
end

function calculateSpline(splajn)
#     prec = 0.01
    coors = []
    angles = []

    for i in 1:length(splajn)
        xk = splajn[i][2][1]
        xkPlus1 = splajn[i][2][2]
        wielomian = splajn[i][1]
        przedzial = [x for x in xk:prec:(xkPlus1-prec)]

        coors = vcat(coors, map(wielomian, przedzial))
        angles = vcat(angles, przedzial)
    end

    return angles, coors
end

function circle(n, r = 1)
    ps = [t for t in -pi:2*pi/n:pi]
    xs = r * cos.(ps)
    ys = r * sin.(ps)

    xspline = makeSpline(myzip(ps, xs))
    yspline = makeSpline(myzip(ps, ys))

    angles, xvalues = calculateSpline(xspline)
    angles, yvalues = calculateSpline(yspline) #jak zrobi się na odwrót kolejności to wychodzą ładne rysunki błędów

    xerrors = abs.(cos.(angles) - xvalues)
    yerrors = abs.(sin.(angles) - yvalues)

    errors = sqrt.([x^2 for x in xerrors] + [y^2 for y in yerrors])

    return [(angles, errors), (xvalues, yvalues), (xs, ys)]
end

function plotCircle(c, n)
    plot!(c[2], lab = "n = " * string(n), linestyle = :dash)
#     scatter!(c[3], lab = "Zadane punkty")
end

function mainCircle(ns :: Array{Int})

    Plots.plot(size = (600*1.1, 600))

    title!("Splajn")

    xlabel!("Os x")
    ylabel!("Os y")

    plot!(idealne, lab = "idealne kolo", color = :black)
    cns = [(circle(n), n) for n in ns]
    for cn in cns
        plotCircle(cn[1], cn[2])
    end

    savefig("splajn")
end
function mainCircle(n :: Int)

    Plots.plot(size = (600*1.1, 600))

    title!("Splajn")

    xlabel!("Os x")
    ylabel!("Os y")

    plot!(idealne, lab = "idealne kolo", color = :black)
    c = circle(n)
    plotCircle(c, n)

    savefig("splajn")
end

function plotError(c, n)
    plot!(c[1], lab = "n = " * string(n))
end
function mainError(ns :: Array{Int})
    Plots.plot(size = (600*1.1, 600))

    title!("Splajn Bledy")

    xlabel!("Os x")
    ylabel!("Os y")

    cns = [(circle(n), n) for n in ns]
    for cn in cns
        plotError(cn[1], cn[2])
    end

    savefig("splajnbledy")
end
function mainError(n :: Int)

    Plots.plot(size = (600*1.1, 600))

    title!("Splajn Blad")

    xlabel!("Os x")
    ylabel!("Os y")

    c = circle(n)
    plotError(c, n)

    savefig("splajnblad")
end

function helix(n, p, r = 1)
    ps = [t for t in -pi:pi/n:pi]
    xs = r * cos.(ps)
    ys = r * sin.(ps)

    xspline = makeSpline(myzip(ps, xs))
    yspline = makeSpline(myzip(ps, ys))

    angles, xvalues = calculateSpline(xspline)
    angles, yvalues = calculateSpline(yspline)

    pts = [p*t for t in angles]

    Plots.plot(size = (600*1.1, 600), leg = false)
    plot!(pts, xvalues, yvalues)
    plot!(xvalues, pts, yvalues)
    plot!(xvalues, yvalues, pts)
    title!("Helisa Splajn")
    xlabel!("Os x")
    ylabel!("Os y")
end

function B(t)
    c = (4/3)*(sqrt(2) - 1)

    P0 = [0,1]
    P1 = [c,1]
    P2 = [1,c]
    P3 = [1,0]

    return (1-t)^3*P0 + 3*(1-t)^2*t*P1 + 3*(1-t)*t^2*P2 + t^3*P3
end

function circle()
    n = 1000
    prec = 1/n
    ts = [t for t in 0:prec:1]

    xvalues, yvalues = myunzip([B(t) for t in ts])

    angles = [t for t in 0:pi/(2*n+1):pi/2]

    xerrors = abs.(sin.(angles) - xvalues)
    yerrors = abs.(cos.(angles) - yvalues)

    errors = sqrt.([x^2 for x in xerrors] + [y^2 for y in yerrors])

    xvs = vcat(xvalues, yvalues, -xvalues, -yvalues) #cosinus
    yvs = vcat(-yvalues, xvalues, yvalues, -xvalues) #sinus

    xvalues = xvs
    yvalues = yvs

    return [(angles, errors), (xvalues, yvalues)]
end

function silnia(k)::BigInt
    if k == 0
        return 1
    else
        return k * silnia(k-1)
    end
end

function ilorazysilni(k)::Tuple{Array{BigFloat,1}, Array{BigFloat,1}}
    xis = [BigFloat(1.) for i in 0:k]
    yis = [BigFloat(1.) for i in 0:k]
    xis[1] = yis[1] = BigFloat(1.)
    for j in 2:k
        i = j - 1
        xis[j] = BigFloat(xis[i] * BigFloat(k+i-1) * BigFloat(k - i) / BigFloat(i))
        yis[j] = xis[j] * (k+i) / (k-i)
    end
    yis[k+1] = yis[k] * 2
    return xis, yis
end

function Hk(k, α, facts, e, f)
    suma = .0
    for i in 0:e
        hi = f(pi*(k+i)/2 + α) / (2*α)^(k+i+1)
        fact = facts[i+1]
        suma += (hi*fact)
    end
    return suma / silnia(k)
end

function Xk(k, α)
    return Hk(k, α, ilorazysilni(k)[1], k-1, cos)
end
function Yk(k, α)
    return Hk(k, α, ilorazysilni(k)[2], k, sin)
end

function Hn(t, n, α, Hks, c1, c2)
    suma = .0
    c = t^2 - α^2
    for k in 1:n
        suma += Hks[k]*c^k
    end
    suma *= c1
    suma += c2
    return suma
end

function Xn(t, n, α, Xks)
    return Hn(t, n, α, Xks, 2*α, cos(α))
end
function Yn(t, n, α, Xks)
    return Hn(t, n, α, Xks, 2*t, t*sin(α)/α)
end

function arc(n, α, r)
    Xks = [Xk(k, α) for k in 1:n]
    Yks = [Yk(k, α) for k in 1:n]

    angles = [t for t in -α:prec:α]

    xvalues = r * [Xn(t, n, α, Xks) for t in angles]
    yvalues = r * [Yn(t, n, α, Yks) for t in angles]

    xerrors = abs.(r * cos.(angles) - xvalues)
    yerrors = abs.(r * sin.(angles) - yvalues)

    angles = convert.(Float64, round.(angles; digits = 10))
    xvalues = convert.(Float64, round.(xvalues; digits = 10))
    yvalues = convert.(Float64, round.(yvalues; digits = 10))
    xerrors = convert.(Float64, round.(xerrors; digits = 10))
    yerrors = convert.(Float64, round.(yerrors; digits = 10))

    errors = sqrt.([x^2 for x in xerrors] + [y^2 for y in yerrors])

    return [(angles,errors), (xvalues, yvalues)]
end
function circle(n, α = pi, r = 1)
    a = arc(n, α, r)
    angles, errors = a[1]
    xvalues, yvalues = a[2]

    if α == pi
        xpoints = [cos(pi), cos(-pi)]
        ypoints = [sin(pi), sin(-pi)]
        return [(angles, errors), (xvalues, yvalues), (xpoints, ypoints)]
    end
    if α == pi/2
        xvalues = vcat(xvalues, reverse(-xvalues))
        yvalues = vcat(yvalues, reverse(yvalues))

        xpoints = [cos(pi/2), cos(-pi/2)]
        ypoints = [sin(pi/2), sin(-pi/2)]
        return [(angles, errors), (xvalues, yvalues), (xpoints, ypoints)]
    end
    if α == pi/4
        xvs = vcat(reverse(xvalues), reverse(yvalues), -xvalues, yvalues)
        yvs = vcat(reverse(yvalues), reverse(-xvalues), yvalues, xvalues)
        xvalues = xvs
        yvalues = yvs

        xpoints = [cos(pi/4), cos(3*pi/4), cos(5*pi/4), cos(7*pi/4)]
        ypoints = [sin(pi/4), sin(3*pi/4), sin(5*pi/4), sin(7*pi/4)]

        return [(angles, errors), (xvalues, yvalues), (xpoints, ypoints)]
    end
end

function helix(n, p, α = pi, r = 1)
    Xks = [Xk(k, α) for k in 1:n]
    Yks = [Yk(k, α) for k in 1:n]

    angles = [t for t in -α:prec:α]

    xvalues = r * [Xn(t, n, α, Xks) for t in angles]
    yvalues = r * [Yn(t, n, α, Yks) for t in angles]
    zvalues = [p * t for t in angles]

    return [xvalues, yvalues, zvalues]
end

max_stars_repo_path: II rok/I semestr/ANM/Pracownia 3/prog/program.jl

max_stars_repo_name: antonitomaszewski/Studia

max_issues_repo_name: antonitomaszewski/Studia

max_forks_repo_name: antonitomaszewski/Studia


text:
module fTst

export f

function f(x::AbstractArray)::Real
    f = 2.0*x[1]^2 + 3.0*x[2]
    return f
end

end

ext: jl

max_stars_repo_name: westrl/tst

max_issues_repo_name: westrl/tst

max_forks_repo_name: westrl/tst


text:
"""

Update the data with the data accumulated in a frame

"""
function update_counters_frame!(
    R::Result,
    rdf_count_random_frame::Vector{Float64},
    md_count_random_frame::Vector{Float64},
    volume_frame::Volume,
    solute::Selection,
    solvent::Selection,
    n_solvent_in_bulk::Float64,
)

    # Volume of each bin shell and of the solute domain
    @. volume_frame.shell =
        volume_frame.total *
        (rdf_count_random_frame / (R.options.n_random_samples * solvent.nmols))
    volume_frame.domain = sum(volume_frame.shell)

    # Volume and density of the bulk region in this frame
    if !R.options.usecutoff
        volume_frame.bulk = volume_frame.total - volume_frame.domain
    else
        ibulk = setbin(R.dbulk + 0.5 * R.options.binstep, R.options.binstep)
        for i = ibulk:R.nbins
            volume_frame.bulk = volume_frame.bulk + volume_frame.shell[i]
        end
    end

    bulk_density_at_frame = n_solvent_in_bulk / volume_frame.bulk

    # The number of "random" molecules generated must be that to fill the complete
    # volume with the density of the bulk. Now that we know the density of the bulk
    # and the volume of the frame, we can adjust that random count
    fillup_factor = bulk_density_at_frame * volume_frame.total / solvent.nmols
    @. R.rdf_count_random = R.rdf_count_random + rdf_count_random_frame * fillup_factor
    @. R.md_count_random = R.md_count_random + md_count_random_frame * fillup_factor

    # Update final counters, which have to be normalized by the number of frames
    # at the end
    @. R.volume.shell = R.volume.shell + volume_frame.shell
    R.volume.bulk = R.volume.bulk + volume_frame.bulk
    R.volume.domain = R.volume.domain + volume_frame.domain
    R.density.solvent_bulk = R.density.solvent_bulk + bulk_density_at_frame

    return nothing
end


text:
Dims{N} = NTuple{N,Integer}

include("impl/pool.jl")
include("impl/conv.jl")

# Convolutions

function cdims(x::NTuple{N}, w::NTuple{N}, pad, stride) where N
  ntuple(Val(N)) do i
    if i < N-1
      1 + div(x[i] - w[i] + 2*pad[i], stride[i])
    elseif i == N-1
      w[N]
    else # i == N
      x[N]
    end
  end
end

# Conv Transpose dims

function ctdims(x::NTuple{N}, w::NTuple{N}, pad, stride, dilation) where N
  ntuple(Val(N)) do i
    if i < N-1
      (x[i] - 1) * stride[i] + dilation[i] * (w[i] - 1) - 2*pad[i] + 1
    elseif i == N-1
      w[N-1]
    else # i == N
      x[N]
    end
  end
end

# Kernel dims

function wdims(x::NTuple{N}, y::NTuple{N}, pad, stride, dilation) where N
  ntuple(Val(N)) do i
    if i < N-1
      1 + div((1 - y[i]) * stride[i] + x[i] + 2pad[i] - 1, dilation[i])
    elseif i == N-1
      x[i]
    else # i == N
      y[i-1]
    end
  end
end

# Interface

head(x) = reverse(Base.tail(reverse(x)))
padtuple(x::Tuple,p::Integer) = map(_->p, head(head(x)))
padtuple(x::Tuple,p::Tuple) = p
padtuple(x::AbstractArray,p) = padtuple(size(x),p)

function conv(x::AbstractArray, w::AbstractArray; size=nothing, pad = 0, stride = 1, dilation = 1)
  pad_, stride_ = padtuple(x, pad), padtuple(x, stride)
  if size === nothing
    size = cdims(Base.size(x), dilation_dims(w, dilation), pad_, stride_)
  end
  conv!(similar(x, size), x, w, pad = pad_, stride = stride_, dilation = dilation)
end

function crosscor(x::A, w::A; size=nothing, pad = 0, stride = 1, dilation = 1) where A<:AbstractArray
  pad_, stride_ = padtuple(x, pad), padtuple(x, stride)
  if size === nothing
    size = cdims(Base.size(x), dilation_dims(w, dilation), pad_, stride_)
  end
  crosscor!(similar(x, size), x, w, pad = pad_, stride = stride_, dilation = dilation)
end

function ∇conv_data(dy::AbstractArray, w::AbstractArray; size=nothing, pad = 0, stride = 1, dilation = 1, flipkernel = 0)
  pad_, stride_, dilation_ = padtuple(dy, pad), padtuple(dy, stride), padtuple(dy, dilation)
  if size === nothing
    size = ctdims(Base.size(dy), Base.size(w), pad_, stride_, dilation_)
  end
  ∇conv_data!(similar(dy, size), dy, w, pad = pad_, stride = stride_, dilation = dilation_, flipkernel=flipkernel)
end

function ∇conv_filter(dy::AbstractArray, x::AbstractArray; size = nothing, pad = 0, stride = 1, dilation = 1, flipkernel=0)
  pad_, stride_, dilation_ = padtuple(dy, pad), padtuple(dy, stride), padtuple(dy, dilation)
  if size === nothing
    size = wdims(Base.size(x), Base.size(dy), pad_, stride_, dilation_)
  end
  ∇conv_filter!(zero(similar(dy, size)), dy, x; pad = pad, stride = stride, dilation = dilation, flipkernel=flipkernel)
end

# N-D dispatch

function conv!(y::AbstractArray{T,3}, x::AbstractArray{T,3}, w::AbstractArray{T,3};
               pad = 0, stride = 1, dilation = 1, flipkernel =0) where T
    args = map(x -> reshape(x, size(x,1),1,size(x,2),size(x,3)), (y, x, w))
    conv!(args..., pad = (pad...,0), stride = (stride...,1), dilation = (dilation...,1), flipkernel=flipkernel)
    return y
end

function crosscor!(y::AbstractArray, x::AbstractArray, w::AbstractArray;
               pad = 0, stride = 1, dilation = 1)
    conv!(y, x, w, pad=pad, stride=stride, dilation=dilation, flipkernel=1)
end

function ∇conv_filter!(dw::AbstractArray{T,3}, dy::AbstractArray{T,3}, x::AbstractArray{T,3};
                       pad = 0, stride = 1, dilation = 1, flipkernel=0) where T
    args = map(x -> reshape(x, size(x,1),1,size(x,2),size(x,3)), (dw, dy, x))
    ∇conv_filter!(args..., pad = (pad...,0), stride = (stride...,1), dilation = (dilation...,1), flipkernel=flipkernel)
    return dw
end

function ∇conv_data!(dx::AbstractArray{T,3}, dy::AbstractArray{T,3}, w::AbstractArray{T,3};
		     pad = 0, stride = 1, dilation = 1, flipkernel = 0) where T
    args = map(x -> reshape(x, size(x,1),1,size(x,2),size(x,3)), (dx, dy, w))
    ∇conv_data!(args..., pad = (pad...,0), stride = (stride...,1), dilation = (dilation..., 1), flipkernel = flipkernel)
    return dx
end

conv!(y::AbstractArray{T,4}, x::AbstractArray{T,4}, w::AbstractArray{T,4};
      pad = 0, stride = 1, dilation = 1, flipkernel=0) where T =
  conv2d!(y, x, w, padding = pad, stride = stride, dilation = dilation, mode=flipkernel)

∇conv_filter!(dw::AbstractArray{T,4}, dy::AbstractArray{T,4}, x::AbstractArray{T,4};
              pad = 0, stride = 1, dilation = 1, flipkernel=0) where T =
  conv2d_grad_w!(dw, x, dy, padding = pad, stride = stride, dilation = dilation, mode=flipkernel)

∇conv_data!(dx::AbstractArray{T,4}, dy::AbstractArray{T,4}, w::AbstractArray{T,4};
            pad = 0, stride = 1, dilation = 1, flipkernel=0) where T =
  conv2d_grad_x!(dx, w, dy, padding = pad, stride = stride, dilation = dilation, mode=flipkernel)

conv!(y::AbstractArray{T,5}, x::AbstractArray{T,5}, w::AbstractArray{T,5};
      pad = 0, stride = 1, dilation = 1, flipkernel=0) where T =
  conv3d!(y, x, w, padding = pad, stride = stride, dilation = dilation, mode=flipkernel)

∇conv_filter!(dw::AbstractArray{T,5}, dy::AbstractArray{T,5}, x::AbstractArray{T,5};
              pad = 0, stride = 1, dilation = 1, flipkernel=0) where T =
  conv3d_grad_w!(dw, x, dy, padding = pad, stride = stride, dilation = dilation, mode=flipkernel)

∇conv_data!(dx::AbstractArray{T,5}, dy::AbstractArray{T,5}, w::AbstractArray{T,5};
            pad = 0, stride = 1, dilation = 1, flipkernel=0) where T =
  conv3d_grad_x!(dx, w, dy, padding = pad, stride = stride, dilation = dilation, mode=flipkernel)

  # Depthwise Conv

function dcdims(x::NTuple{4,Int}, w::NTuple{4,Int}, pad, stride)
  ((x[1] + 2 * pad[1] - w[1])÷stride[1] + 1,(x[2] + 2 * pad[2] - w[2])÷stride[2] + 1,w[3]*w[4],x[4])
end

function depthwiseconv(x::AbstractArray, w::AbstractArray; pad = 0, stride = 1)
  pad_, stride_ = padtuple(x, pad), padtuple(x, stride)
  depthwiseconv!(similar(x, dcdims(size(x), size(w), pad_, stride_)), x, w, pad = pad_, stride = stride_)
end

function depthwisecrosscor(x::A, w::A; pad = 0, stride = 1) where A<:AbstractArray
  pad_, stride_ = padtuple(x, pad), padtuple(x, stride)
  depthwisecrosscor!(similar(x, dcdims(size(x), size(w), pad_, stride_)), x, w, pad = pad_, stride = stride_)
end

depthwiseconv!(y::AbstractArray{T,4}, x::AbstractArray{T,4}, w::AbstractArray{T,4};
      pad = 0, stride = 1, flipkernel=0) where T =
  depthwiseconv2d!(y, x, w, padding = pad, stride = stride, mode= flipkernel)

depthwisecrosscor!(y::AbstractArray{T,4}, x::AbstractArray{T,4}, w::AbstractArray{T,4};
      pad = 0, stride = 1) where T =
  depthwiseconv!(y, x, w, pad = pad, stride = stride, flipkernel=1)

∇depthwiseconv_data(dy::AbstractArray, x::AbstractArray, w::AbstractArray; pad = 0, stride = 1, flipkernel=0) =
  ∇depthwiseconv_data!(zero(x), dy, x, w; pad = pad, stride = stride, flipkernel=flipkernel)

∇depthwiseconv_filter(dy::AbstractArray, x::AbstractArray, w::AbstractArray; pad = 0, stride = 1, flipkernel=0) =
  ∇depthwiseconv_filter!(zero(w), dy, x, w; pad = pad, stride = stride, flipkernel=flipkernel)

∇depthwiseconv_filter!(dw::AbstractArray{T,4}, dy::AbstractArray{T,4}, x::AbstractArray{T,4}, w::AbstractArray{T,4};
              pad = 0, stride = 1, flipkernel=0) where T =
  depthwiseconv2d_grad_w!(dw, x, w, dy, padding = pad, stride = stride, mode=flipkernel)

∇depthwiseconv_data!(dx::AbstractArray{T,4}, dy::AbstractArray{T,4}, x::AbstractArray{T,4}, w::AbstractArray{T,4};
            pad = 0, stride = 1, flipkernel=0) where T =
  depthwiseconv2d_grad_x!(dx, x, w, dy, padding = pad, stride = stride, mode=flipkernel)

# Pooling

function pdims(dims::Dims{N}, window, padding, stride) where N
  ntuple(Val(N)) do i
    if i < N-1
      1 + (dims[i] + 2*padding[i] - window[i])÷stride[i]
    else
      dims[i]
    end
  end
end

expand(::Type{Val{N}}, i::Integer) where N = ntuple(_ -> i, Val(N))
expand(::Type{Val{N}}, i::NTuple{N, Integer}) where N = i

# Interface

maxpool(x::AbstractArray, k; pad = map(_->0,k), stride = k) =
  maxpool!(similar(x, pdims(size(x), k, expand(Val{length(k)}, pad),
            expand(Val{length(k)}, stride))), x, k, pad = expand(Val{length(k)}, pad),
            stride = expand(Val{length(k)}, stride))

maxpool!(y::A, x::A, k; kw...) where A<:AbstractArray =
  maxpool_cpu!(y, x, k; kw...)

∇maxpool(dy::A, y::A, x::A, k; pad = map(_->0,k), stride = k) where A<:AbstractArray =
  ∇maxpool!(similar(x), dy, y, x, k, pad = expand(Val{length(k)}, pad),
            stride = expand(Val{length(k)}, stride))

∇maxpool!(dx::A, dy::A, y::A, x::A, k; kw...) where A<:AbstractArray =
  ∇maxpool_cpu!(dx, dy, y, x, k; kw...)

meanpool(x::AbstractArray, k; pad = map(_->0,k), stride = k) =
  meanpool!(similar(x, pdims(size(x), k, expand(Val{length(k)}, pad),
            expand(Val{length(k)}, stride))), x, k, pad = expand(Val{length(k)}, pad),
            stride = expand(Val{length(k)}, stride))

meanpool!(y::A, x::A, k; kw...) where A<:AbstractArray =
  meanpool_cpu!(y, x, k; kw...)

∇meanpool(dy::A, y::A, x::A, k; pad = map(_->0,k), stride = k) where A<:AbstractArray =
  ∇meanpool!(similar(x), dy, y, x, k, pad = expand(Val{length(k)}, pad),
            stride = expand(Val{length(k)}, stride))

∇meanpool!(dx::A, dy::A, y::A, x::A, k; kw...) where A<:AbstractArray =
  ∇meanpool_cpu!(dx, dy, y, x, k; kw...)

# N-D dispatch
# We use a separate function to avoid ambiguity issues
# (more specific array types vs. more specific dimensions)

maxpool_cpu!(y::A, x::A, k::Dims{2}; pad = (0,0), stride = k) where A<:AbstractArray{<:Real,4} =
  maxpool2d!(y, x, window = k, padding = pad, stride = stride)

∇maxpool_cpu!(dx::AbstractArray{<:Real,4}, dy::AbstractArray{<:Real,4}, y::AbstractArray{<:Real,4}, x::AbstractArray{<:Real,4},
              k::Dims{2}; pad = (0,0), stride = k) =
  maxpool2d_grad!(dx, dy, y, x,
                  window = k, padding = pad, stride = stride)

maxpool_cpu!(y::AbstractArray{<:Real,5}, x::AbstractArray{<:Real,5}, k::Dims{3}; pad = (0,0), stride = k) =
  maxpool3d!(y, x, window = k, padding = pad, stride = stride)

∇maxpool_cpu!(dx::AbstractArray{<:Real,5}, dy::AbstractArray{<:Real,5}, y::AbstractArray{<:Real,5}, x::AbstractArray{<:Real,5},
              k::Dims{3}; pad = (0,0), stride = k) =
  maxpool3d_grad!(dx, dy, y, x,
                  window = k, padding = pad, stride = stride)

meanpool_cpu!(y::AbstractArray{<:Real,4}, x::AbstractArray{<:Real,4}, k::Dims{2}; pad = (0,0), stride = k) =
  meanpool2d!(y, x, window = k, padding = pad, stride = stride)

∇meanpool_cpu!(dx::AbstractArray{<:Real,4}, dy::AbstractArray{<:Real,4}, y::AbstractArray{<:Real,4}, x::AbstractArray{<:Real,4},
              k::Dims{2}; pad = (0,0), stride = k) =
  meanpool2d_grad!(dx, dy, y, x,
                   window = k, padding = pad, stride = stride)

meanpool_cpu!(y::AbstractArray{<:Real,5}, x::AbstractArray{<:Real,5}, k::Dims{3}; pad = (0,0), stride = k) =
  meanpool3d!(y, x, window = k, padding = pad, stride = stride)

∇meanpool_cpu!(dx::AbstractArray{<:Real,5}, dy::AbstractArray{<:Real,5}, y::AbstractArray{<:Real,5}, x::AbstractArray{<:Real,5},
              k::Dims{3}; pad = (0,0), stride = k) =
  meanpool3d_grad!(dx, dy, y, x,
                   window = k, padding = pad, stride = stride)

# Deprecated

# 0.4.2
@deprecate ∇conv_data(dy::A, x::A, w::A; kw...) where A<:AbstractArray ∇conv_data(dy, w; size=size(x), kw...)
@deprecate ∇conv_filter(dy::A, x::A, w::A; kw...) where A<:AbstractArray ∇conv_filter(dy, x; size=size(w), kw...)


text:
## =============================================================================
# Hidden markov model with zero inflated likelihood

struct ZIPoissonProblem{Y,N,I,P} <: SurveyProblem
    y::Y
    n::N
    indices::I
    prior::P
end

function ZIPoissonProblem(y,n,indices; c=nothing)
    prior = ZIPoissonPrior(n; c)
    ZIPoissonProblem(y,n,indices, prior)
end

struct ZIPoissonPrior{P,L,N}
    prior_p::P
    prior_μ::L
    n::N
end

function ZIPoissonPrior(n; c=nothing)
    if isnothing(c)
        prior_p  = filldist(Beta(), n)
        prior_μ = filldist(GammaMeanCv(15,.8), n)
        return ZIPoissonPrior( prior_p, prior_μ, n)
    elseif c == :family
        prior_p  = filldist(Beta(), n)
        # prior_p  = truncated(Normal(0.51,0.05), 0, 1)
        prior_μ = filldist(truncated(Normal(3.5,1), 0, Inf), n)
        return ZIPoissonPrior( prior_p, prior_μ, n)
    elseif c == :colleagues
        prior_p  = filldist(Beta(), n)
        # prior_p  = truncated(Normal(0.33,0.05), 0, 1)
        prior_μ = filldist(truncated(Normal(5,1.5), 0, Inf), n)
        return ZIPoissonPrior( prior_p, prior_μ, n)
    elseif c == :friends
        prior_p  = filldist(Beta(), n)
        # prior_p  = truncated(Normal(0.46,0.05), 0, 1)
        prior_μ = filldist(truncated(Normal(4.5,2), 0, Inf), n)
        return ZIPoissonPrior( prior_p, prior_μ, n)
    elseif c == :strangers
        prior_p  = filldist(Beta(), n)
        # prior_p  = truncated(Normal(0.32,0.05), 0, 1)
        prior_μ = filldist(truncated(Normal(20,8), 0, Inf), n)
        return ZIPoissonPrior( prior_p, prior_μ, n)
    else
        error( "wrong choice" )
    end

end

function Base.rand(priors::ZIPoissonPrior)
    @unpack prior_p, prior_μ, n = priors
    p  = rand(prior_p)
    μ  = rand(prior_μ)
    return (; p, μ)
end

function Distributions.logpdf(prior::ZIPoissonPrior, θ)
    @unpack p, μ = θ
    @unpack prior_p, prior_μ = prior
    ℓ  = logpdf(prior_p, p)
    ℓ += logpdf(prior_μ, μ)
    return ℓ
end

# likelihood_dist(p, μ, ϕ) =
#     ZeroInflatedDistribution.(
#         Ref(Bernoulli(p)),
#         NegativeBinomial2.(μ, Ref(ϕ)
#     )
# )
likelihood_dist(p, μ) =
    ZeroInflatedDistribution.(
        Bernoulli.(p),
        Poisson.(μ)
)

function (problem::ZIPoissonProblem)(θ)
    @unpack p, μ  = θ               # extract the parameters
    @unpack y, indices, prior = problem
    # !(ϕ > 0) && (@error "not ϕ>0"; return -Inf)
    # !(0 <= p <= 1) && (@error "not 0 <= p <= 1"; return -Inf)

    ℓ  = logpdf(prior, θ)
    # !all(0 .<= μ ) && (@error "not all(0 .<= μ )"; return -Inf)

    d = likelihood_dist(p, μ)
    ℓ += logpdf(arraydist(d[indices]),y)
    return ℓ
end

function problem_transformation(problem::ZIPoissonProblem)
    as((
        p = as(Vector, as𝕀, problem.n),
        μ = as(Vector, asℝ₊, problem.n)
    ))
end

function ZIPoissonProblem(data, col)
    y = data[:,col]
    indices = date2index( data.date )
    n = last(indices)
    return ZIPoissonProblem(y, n, indices; c=col)
end

function StatsBase.predict(problem::ZIPoissonProblem, θ)
    @unpack p, μ = θ               # extract the parameters
    @unpack y, indices = problem

    d = likelihood_dist(p, μ)
    y = rand(arraydist(d[indices]))
    return y
end

generate_quantities(problem::ZIPoissonProblem, θ) = θ.μ

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
# ================================
# Abstract types
# ================================

"""
    AbstractSolution

Abstract supertype of all solution types of a rechability problem.
"""
abstract type AbstractSolution end

# ================================
# Property checking problem
# ================================

"""
    CheckSolution{T} <: AbstractSolution

Type that wraps the solution of a verification problem.

### Fields

TODO

### Notes

This type contains the answer if the property is satisfied, and if not, it
contains the index at which the property might be violated for the first time.
"""
struct CheckSolution{PT, ST<:AbstractPost} <: AbstractSolution
    property::PT
    satisfied::Bool
    vidx::Int
    vtspan::TimeInterval
    alg::ST
end

property(sol::CheckSolution) = sol.property
satisfied(sol::CheckSolution) = sol.satisfied
violation_index(sol::CheckSolution) = sol.vidx
violation_tspan(sol::CheckSolution) = sol.vtspan
alg(sol::CheckSolution) = sol.alg

# ================================
# Reachability problem
# ================================

"""
    ReachSolution{FT<:AbstractFlowpipe, ST<:AbstractPost} <: AbstractSolution

Type that wraps the solution of a reachability problem as a sequence of lazy
sets, and a dictionary of options.

### Fields

- `Xk`       -- the list of [`AbstractReachSet`](@ref)s
- `alg`      -- algorihm used
- `options`  -- the dictionary of options
"""
struct ReachSolution{FT<:AbstractFlowpipe, ST<:AbstractPost} <: AbstractSolution
    F::FT
    alg::ST
    ext::Dict{Symbol, Any} # dictionary used by extensions
end

# constructor from empty extension dictionary
function ReachSolution(F::FT, alg::ST) where {FT<:AbstractFlowpipe, ST<:AbstractPost}
    return ReachSolution(F, alg, Dict{Symbol, Any}())
end

# getter functions
flowpipe(sol::ReachSolution) = flowpipe(sol.F)
tstart(sol::ReachSolution) = tstart(sol.F)
tend(sol::ReachSolution) = tend(sol.F)
tspan(sol::ReachSolution) = tspan(sol.F)
tstart(sol::ReachSolution, arr::AbstractVector) = tstart(sol.F, arr)
tend(sol::ReachSolution, arr::AbstractVector) = tend(sol.F, arr)
tspan(sol::ReachSolution, arr::AbstractVector) = tspan(sol.F, arr)
ensemble(sol::ReachSolution) = get(sol.ext, :ensemble, nothing)

# NOTE: using sol.alg for setrep and rsetrep may be inaccurate as solutions
# do not change their sol.alg (eg. overapproximate(sol, Zonotope)) where sol is
# a Taylor model flowpipe keeps sol.alg being TMJets => we have to check
# for the rep's of the flowpipe
setrep(sol::ReachSolution{FT, ST}) where {FT, ST} = setrep(FT)
rsetrep(sol::ReachSolution{FT, ST}) where {FT, ST} = rsetrep(FT)
vars(sol::ReachSolution) = vars(sol.F)
numrsets(sol::ReachSolution) = numrsets(sol.F)

# iteration and indexing iterator interface
array(sol::ReachSolution) = array(sol.F)
Base.iterate(sol::ReachSolution) = iterate(sol.F)
Base.iterate(sol::ReachSolution, state) = iterate(sol.F, state)
Base.length(sol::ReachSolution) = length(sol.F)
Base.first(sol::ReachSolution) = first(sol.F)
Base.last(sol::ReachSolution) = last(sol.F)
Base.firstindex(sol::ReachSolution) = 1
Base.lastindex(sol::ReachSolution) = length(sol.F)
Base.getindex(sol::ReachSolution, i::Int) = getindex(sol.F, i)
Base.getindex(sol::ReachSolution, i::Number) = getindex(sol.F, i)
Base.getindex(sol::ReachSolution, I::AbstractVector) = getindex(sol.F, I)

# indexing: fp[j, i] returning the j-th reach-set of the i-th flowpipe
Base.getindex(sol::ReachSolution, I::Int...) = getindex(sol.F, I...)
Base.eachindex(sol::ReachSolution) = eachindex(sol.F)
Base.keys(sol::ReachSolution) = keys(sol.F)

# evaluation interface
Base.getindex(sol::ReachSolution, t::Float64) = getindex(sol.F, t)
(sol::ReachSolution)(t::Float64) = sol.F(t)
(sol::ReachSolution)(t::Number) = sol.F(t)
(sol::ReachSolution)(dt::IA.Interval{Float64}) = sol.F(dt)

function overapproximate(sol::ReachSolution{FT}, args...) where {FT<:AbstractFlowpipe}
    return ReachSolution(overapproximate(sol.F, args...), sol.alg, sol.ext)
end

function box_approximation(sol::ReachSolution{FT}) where {FT<:AbstractFlowpipe}
    return ReachSolution(overapproximate(sol.F, Hyperrectangle), sol.alg, sol.ext)
end

function project(sol::ReachSolution{FT}, args...) where {FT<:AbstractFlowpipe}
    return project(sol.F, args...)
end

# convenience alias to match the usage in the plot recipe
function project(sol::ReachSolution{FT}; vars) where {FT<:AbstractFlowpipe}
    return project(sol.F, Tuple(vars))
end

function shift(sol::ReachSolution{<:AbstractFlowpipe}, t0::Number)
    return ReachSolution(shift(sol.F, t0), sol.alg, sol.ext)
end

# LazySets interface falls back to the associated flowpipe
LazySets.dim(sol::ReachSolution) = dim(sol.F)
LazySets.ρ(d, sol::ReachSolution{FT}) where {FT<:AbstractFlowpipe} = ρ(d, sol.F)
LazySets.σ(d, sol::ReachSolution{FT}) where {FT<:AbstractFlowpipe} = σ(d, sol.F)

# further setops functions acting on solutions' flowpipes
∈(x::AbstractVector, sol::ReachSolution) = ∈(x, sol.F)

# ------------------------------
# Methods to check disjointness
# ------------------------------

# interface
is_intersection_empty(sol::ReachSolution, Y::SetOrReachSet, method::AbstractDisjointnessMethod=FallbackDisjointness()) = is_intersection_empty(sol.F, Y, method)
is_intersection_empty(Y::SetOrReachSet, sol::ReachSolution, method::AbstractDisjointnessMethod=FallbackDisjointness()) = is_intersection_empty(sol.F, Y, method)
Base.isdisjoint(sol::ReachSolution, Y::SetOrReachSet, method::AbstractDisjointnessMethod=FallbackDisjointness()) = is_intersection_empty(sol.F, Y, method)
Base.isdisjoint(Y::SetOrReachSet, sol::ReachSolution, method::AbstractDisjointnessMethod=FallbackDisjointness()) = is_intersection_empty(sol.F, Y, method)

is_intersection_empty(sol1::ReachSolution, sol2::ReachSolution, method::AbstractDisjointnessMethod=FallbackDisjointness()) = is_intersection_empty(sol1.F, sol2.F, method)
Base.isdisjoint(sol1::ReachSolution, sol2::ReachSolution, method::AbstractDisjointnessMethod=FallbackDisjointness()) = is_intersection_empty(sol1.F, sol2.F, method)

# ------------------------------
# Methods to check inclusion
# ------------------------------

# inclusion checks
Base.:⊆(sol::ReachabilityAnalysis.ReachSolution, X::LazySet) = ⊆(sol.F, X)
Base.:⊆(sol::ReachabilityAnalysis.ReachSolution, Y::AbstractLazyReachSet) = ⊆(sol.F, set(Y))

# -------------------------------------
# Methods for concrete set operations
# -------------------------------------

linear_map(M::AbstractMatrix, sol::ReachSolution) = linear_map(M, sol.F)

# -------------------------------------------
# Specialized methods for template solutions
# -------------------------------------------

function support_function_matrix(sol::ReachSolution{<:Flowpipe{N, <:TemplateReachSet}}) where {N}
    return support_function_matrix(flowpipe(sol))
end

function flatten(sol::ReachSolution{<:Flowpipe{N, <:TemplateReachSet}}, rows=(1, 2)) where {N}
    return flatten(flowpipe(sol), rows)
end


text:
### INITIALIZE

@cdef :Py_Initialize Cvoid ()
@cdef :Py_InitializeEx Cvoid (Cint,)
@cdef :Py_Finalize Cvoid ()
@cdef :Py_FinalizeEx Cint ()
@cdef :Py_AtExit Cint (Ptr{Cvoid},)
@cdef :Py_IsInitialized Cint ()

@cdef :Py_SetPythonHome Cvoid (Cwstring,)
@cdef :Py_SetProgramName Cvoid (Cwstring,)
@cdef :Py_GetVersion Cstring ()

### REFCOUNT

@cdef :Py_IncRef Cvoid (PyPtr,)
@cdef :Py_DecRef Cvoid (PyPtr,)
Py_RefCnt(o) = GC.@preserve o UnsafePtr(Base.unsafe_convert(PyPtr, o)).refcnt[]

Py_DecRef(f::Function, o::Ptr, dflt = PYERR()) =
    isnull(o) ? dflt : (r = f(o); Py_DecRef(o); r)

Py_Is(o1, o2) = Base.unsafe_convert(PyPtr, o1) == Base.unsafe_convert(PyPtr, o2)

### EVAL

@cdef :PyEval_EvalCode PyPtr (PyPtr, PyPtr, PyPtr)
@cdef :Py_CompileString PyPtr (Cstring, Cstring, Cint)
@cdef :PyEval_GetBuiltins PyPtr ()
@cdef :PyRun_InteractiveOne Cint (Ptr{Cvoid}, Cstring) # (FILE* file, const char* filename)

### GIL & THREADS

@cdef :PyEval_SaveThread Ptr{Cvoid} ()
@cdef :PyEval_RestoreThread Cvoid (Ptr{Cvoid},)
@cdef :PyGILState_Ensure PyGILState_STATE ()
@cdef :PyGILState_Release Cvoid (PyGILState_STATE,)

### IMPORT

@cdef :PyImport_ImportModule PyPtr (Cstring,)
@cdef :PyImport_Import PyPtr (PyPtr,)
@cdef :PyImport_GetModuleDict PyPtr ()

PyImport_GetModule(name) = begin
    ms = PyImport_GetModuleDict()
    ok = PyMapping_HasKeyString(ms, name)
    ism1(ok) && return PyNULL
    ok != 0 ? PyMapping_GetItemString(ms, name) : PyNULL
end

### MODULES

@cdef :PyModule_GetDict PyPtr (PyPtr,)

### ERRORS

@cdef :PyErr_Occurred PyPtr ()
@cdef :PyErr_GivenExceptionMatches Cint (PyPtr, PyPtr)
@cdef :PyErr_Clear Cvoid ()
@cdef :PyErr_SetNone Cvoid (PyPtr,)
@cdef :PyErr_SetString Cvoid (PyPtr, Cstring)
@cdef :PyErr_SetObject Cvoid (PyPtr, PyPtr)
@cdef :PyErr_Fetch Cvoid (Ptr{PyPtr}, Ptr{PyPtr}, Ptr{PyPtr})
@cdef :PyErr_NormalizeException Cvoid (Ptr{PyPtr}, Ptr{PyPtr}, Ptr{PyPtr})
@cdef :PyErr_Restore Cvoid (PyPtr, PyPtr, PyPtr)

PyErr_IsSet() = !isnull(PyErr_Occurred())

PyErr_IsSet(t) =
    (o = PyErr_Occurred(); !isnull(o) && PyErr_GivenExceptionMatches(o, t) != 0)

function PyErr_FetchTuple(normalize::Bool = false)
    t = Ref{PyPtr}()
    v = Ref{PyPtr}()
    b = Ref{PyPtr}()
    PyErr_Fetch(t, v, b)
    normalize && PyErr_NormalizeException(t, v, b)
    (t[], v[], b[])
end

### EXCEPTIONS

for x in [
    :BaseException,
    :Exception,
    :StopIteration,
    :GeneratorExit,
    :ArithmeticError,
    :LookupError,
    :AssertionError,
    :AttributeError,
    :BufferError,
    :EOFError,
    :FloatingPointError,
    :OSError,
    :ImportError,
    :IndexError,
    :KeyError,
    :KeyboardInterrupt,
    :MemoryError,
    :NameError,
    :OverflowError,
    :RuntimeError,
    :RecursionError,
    :NotImplementedError,
    :SyntaxError,
    :IndentationError,
    :TabError,
    :ReferenceError,
    :SystemError,
    :SystemExit,
    :TypeError,
    :UnboundLocalError,
    :UnicodeError,
    :UnicodeEncodeError,
    :UnicodeDecodeError,
    :UnicodeTranslateError,
    :ValueError,
    :ZeroDivisionError,
    :BlockingIOError,
    :BrokenPipeError,
    :ChildProcessError,
    :ConnectionError,
    :ConnectionAbortedError,
    :ConnectionRefusedError,
    :FileExistsError,
    :FileNotFoundError,
    :InterruptedError,
    :IsADirectoryError,
    :NotADirectoryError,
    :PermissionError,
    :ProcessLookupError,
    :TimeoutError,
    :EnvironmentError,
    :IOError,
    :WindowsError,
    :Warning,
    :UserWarning,
    :DeprecationWarning,
    :PendingDeprecationWarning,
    :SyntaxWarning,
    :RuntimeWarning,
    :FutureWarning,
    :ImportWarning,
    :UnicodeWarning,
    :BytesWarning,
    :ResourceWarning,
]
    f = Symbol(:PyExc_, x)
    r = Symbol(f, :__ref)
    @eval const $r = Ref(PyNULL)
    @eval $f() = pyloadglobal($r, $(QuoteNode(f)))
end

### INPUT HOOK

function PyOS_RunInputHook()
    hook = unsafe_load(Ptr{Ptr{Cvoid}}(pyglobal(:PyOS_InputHook)))
    isnull(hook) || ccall(hook, Cint, ())
    nothing
end


text:
convert(::Type{Array}, f::ArrayFun) = reshape(vec(f), size(space(f))...)
Array(f::ArrayFun) = convert(Array, f)
Vector(f::VectorFun) = Array(f)
Matrix(f::MatrixFun) = Array(f)

map(f,A::ArrayFun) = Base.collect_similar(A, Base.Generator(f,A))

similar(a::VectorFun, S::Type) = Array{S,1}(undef, size(a,1))
similar(a::MatrixFun, S::Type) = Array{S,2}(undef, size(a,1), size(a,2))

getindex(f::MatrixFun,
               k::Union{Integer,AbstractRange,Colon},
               j::Union{Integer,AbstractRange,Colon}) =
    Fun(Array(f)[k,j])

const FunTypes = Union{Fun,Number}
const ScalarFunTypes = Union{ScalarFun,Number}
function vcat(vin::FunTypes...)
    #  remove tuple spaces
    v=Vector{Fun}(undef,0)
    for f in vin
        if rangetype(space(f)) <: AbstractVector
            append!(v,vec(f))
        else
            push!(v,f)
        end
    end

    S = Space(space.(v))
    Fun(S,interlace(v,S))
end

function hcat(v::ScalarFunTypes...)
    ff = vcat(v...)  # A vectorized version
    transpose(ff)
end

hvcat(rows::Tuple{Vararg{Int}},v::FunTypes...) = Fun(hvnocat(rows,v...))

function hcat(v::VectorFun...)
    N = length(v[1])
    M = length(v)

    V = Array{Fun}(undef, N,M)
    for J=1:M
        V[:,J] = vec(v[J])
    end
    Fun(V)
end

function Fun(v::AbstractVector{F}) where F<:Fun
    S = Space(space.(v))
    Fun(S,interlace(v,S))
end

#TODO: rewrite
function Fun(v::AbstractArray{<:Fun})
    ff=Fun(vec(v))  # A vectorized version
    Fun(Space(map(space,v)),coefficients(ff))
end

Fun(v::AbstractArray{NN}) where {NN<:Number} =
    Fun(v,Space(fill(ConstantSpace(NN),size(v))))
Fun(v::AbstractArray) = Fun(Fun.(Array(v)) :: AbstractArray{<:Fun})

Fun(f::ArrayFun, d::Space{D,R}) where {D,R<:AbstractArray} = space(f)==d ? f : Fun(d,coefficients(f,d))
Fun(f::ArrayFun, d::Space) = Fun(f,Space(fill(d,size(space(f)))))

Fun(M::AbstractMatrix{<:Number},sp::Space) = Fun([Fun(M[:,k],sp) for k=1:size(M,2)])

for OP in (:(transpose),)
    @eval begin
        $OP(f::ArrayFun) = Fun($OP(Array(f)))
        $OP(sp::Space{D,R}) where {D,R<:AbstractArray} = Space($OP(Array(sp)))
    end
end

## calculus

for op in (:differentiate,:integrate,:(cumsum),:(real),:(imag),:(conj))
    @eval $op(f::ArrayFun) = Fun(map($op,f))
end

# TODO: use QR
function det(f::MatrixFun)
    @assert size(space(f))==(2,2)
    m=Array(f)
    m[1,1]*m[2,2]-m[1,2]*m[2,1]
end

function inv(V::MatrixFun)
    n,m = size(space(V))
    if n ≠ m
        throw(DimensionMismatch("space $(space(V)) is not square"))
    end

    # TODO: This assumes other columns have same spaces
    M=Multiplication(V,Space(space(V).spaces[:,1]))
    # convert I to the rangespace of M
    M\Fun(Matrix(I,m,m), repeat(rangespace(M),1,m))
end

## Algebra

for OP in (:*,:+,:-)
    @eval begin
        $OP(A::AbstractArray{<:Number}, f::ArrayFun) = Fun($OP(A,Array(f)))
        $OP(f::ArrayFun,                A::AbstractArray{<:Number}) = Fun($OP(Array(f),A))
        $OP(A::AbstractArray{<:Fun},    f::ArrayFun) = Fun($OP(A,Array(f)))
        $OP(f::ArrayFun,                A::AbstractArray{<:Fun}) = Fun($OP(Array(f),A))
        $OP(A::UniformScaling,          f::ArrayFun) = Fun($OP(A,Array(f)))
        $OP(f::ArrayFun,                A::UniformScaling) = Fun($OP(Array(f),A))
        $OP(A::Number,                  f::ArrayFun) = Fun($OP(A,Array(f)))
        $OP(f::ArrayFun,                A::Number) = Fun($OP(Array(f),A))

        $OP(f::ScalarFun,               A::AbstractArray) = Fun(broadcast($OP,f,A))
        $OP(A::AbstractArray,           f::ScalarFun) = Fun(broadcast($OP,A,f))

        $OP(f::ScalarFun,               A::ArrayFun) = $OP(f,Array(A))
        $OP(A::ArrayFun,                f::ScalarFun) = $OP(Array(A),f)
    end
end

# use standard +, -
*(A::ArrayFun,f::ArrayFun) = Fun(Array(A)*Array(f))

norm(A::VectorFun, p::Real) = norm(norm.(Array(A)),p)

## Vector of fun routines

function coefficientmatrix(::Type{N},f::AbstractVector{F},o...) where {N,F}
    if isempty(f)
        return Matrix{N}(0,0)
    end

    n=mapreduce(ncoefficients,max,f)
    m=length(f)
    R=zeros(N,n,m)
    for k=1:m
        R[1:ncoefficients(f[k]),k]=coefficients(f[k],o...)
    end
    R
end

scalarorfuntype(::Fun{S,T}) where {S,T<:Number} = T
scalarorfuntype(::T) where {T<:Number} = T
scalarorfuntype(b::AbstractVector{T}) where {T<:Number} = T
scalarorfuntype(b::AbstractVector{Any}) = promote_type(map(scalarorfuntype,b)...)
scalarorfuntype(b::AbstractVector{F}) where {F<:Fun} = promote_type(map(scalarorfuntype,b)...)

coefficientmatrix(Q::AbstractVector{F},o...) where {F<:Fun}=coefficientmatrix(scalarorfuntype(Q),Q,o...)
coefficientmatrix(Q::AbstractVector{Any})=(@assert isempty(Q); zeros(0,0))

function values(f::AbstractVector{Fun{D,N,VN}}) where {D,N,VN}
    n=mapreduce(ncoefficients,max,f)
    m=length(f)
    R=zeros(N,n,m)
    for k=1:m
        R[:,k] = values(pad(f[k],n))
    end
    R
end

function values(p::AbstractMatrix{Fun{D,T,VT}}) where {D,T,VT}
    @assert size(p)[1] == 1

   values(vec(p))
end

## evaluation

#TODO: fix for complex
evaluate(A::AbstractArray{T},x::Number) where {T<:Fun} =
    typeof(first(A)(x))[Akj(x) for Akj in A]


text:
using Graphs

"""
    SimpleGraph(R::AbstractRecurrenceMatrix)

Create a recurrent network as a `SimpleGraph`, from the symmetric recurrence matrix `R`.
That matrix is typically the result of calculating a [`RecurrenceMatrix`](@ref)
or a [`JointRecurrenceMatrix`](@ref).

The recurrence structure of `R` is interpreted as the adjacency matrix of an
undirected complex network, where two different vertices are connected if they are
neighbors in the embedded phase space, i.e.

```math
A_{i,j} = R_{i,j} - \\delta_{i,j}
```

Following this definition, diagonal points of `R` are ommited, i.e.
the graph does not contain self-connected nodes.

See the package [Graphs.jl](https://juliagraphs.org/Graphs.jl/stable/)
for further options to work with `SimpleGraph` objects, besides the functions
for Recurrence Network Analysis provided in this package.

# References

[1] : R.V. Donner *et al.* "Recurrence networks — a novel paradigm for nonlinear time series analysis",
*New Journal of Physics* 12, 033025 (2010).
[DOI:10.1088/1367-2630/12/3/033025](https://doi.org/10.1088/1367-2630/12/3/033025)

[2] : R.V. Donner *et al.* "Complex Network Analysis of Recurrences", in:
Webber, C.L. & Marwan N. (eds.) *Recurrence Quantification Analysis.
Theory and Best Practices*, Springer, pp. 101-165 (2015).
"""
function Graphs.SimpleGraphs.SimpleGraph(R::AbstractRecurrenceMatrix)
    graph = SimpleGraph(R.data)
    delta = SimpleGraphFromIterator(Edge(v,v) for v in 1:size(graph, 1))
    graph = difference(graph, delta)
    return graph
end


text:
"
`FitnessFunction` represents a fitness function that is going to be evaluated
over the individuals of a population.

# Fields
- `_fitnessFunction::Function`: the actual function that is going to be run.
- `_varArgs::Array`: arguments of the fitness function (if any).
- `_weight::Float32`: the weight of the fitness value. If this value is positive,
    the fitness function will be maximized, if it is negative, the fitness
    function will be minimized.

Several considerations must be taken in order to create a fitness function.

Every fitness function must **always** recieve the individual's genotype as its
first argument. If the function needs other arguments besides the genotype, they
are specified in the `_varArgs` field of the fitness function (see Fields). The
user does **not** need to specify the genotype as one of the fitness function
arguments, as the system already knows it is one of them. Here is an example of
how to create a `FitnessFunction` instance:

```jldoctest
julia> function countBinaryOnes(genotype::BinaryGenotype)
    representation = getRepresentation(genotype)
    tam = length(representation)
    fitness = 0
    for i=1:tam
        if representation[i] == 1
            fitness += 1
        end
    end
    return fitness
end
countBinaryOnes (generic function with 1 method)

julia> FitnessFunction(countBinaryOnes)
FitnessFunction(countBinaryOnes, Any[], 1.0)
```

A `FitnessFunction` instance has been created with `countBinaryOnes` as the
actual function, which receives two arguments (2 and 0.5), and a weight of -1.0,
so it is to be minimized. Notice that, in order to get the actual representation
of the genotype (which is BinaryGenotype), [`getRepresentation`](@ref) function
is called, which returns in this case an array of Bool.

If it is a Genetic Programming problem, the fitness functions must receive
another mandatory argument as their second argument, even if they do not use it,
which is the information about the GP experiment (see [`GPExperimentInfo`](@ref)).
Again, the user does **not** need to specify de GP experiment information as one
of the fitness function arguments. Here is an example of how to create a fitness
function for a GP problem:

```jldoctest
julia> function myFitnessFunction(genotype::STGPGenotype, gpExperimentInfo::STGPInfo, arg1, arg2, arg3)
    #···Code···
    return fitness
end
myFitnessFunction (generic function with 1 method)

julia> FitnessFunction(myFitnessFunction, [1,0,0,1], [5,2,6,1], [6,7,2,1], weight=2)
FitnessFunction(myFitnessFunction, Any[[1,0,0,1], [5,2,6,1], [6,7,2,1]], 2.0)
```

Here, a `FitnessFunction` instance has been created with `myFitnessFunction` as the
actual function, which receives three arrays as arguments, and a weight of 2.0,
so it is to be maximized and it weights twice as much.

For the GP case, the obtention of the genotype's representation can be most of
the time meaningless, so the user might want to obtain the phenotype instead.
The user can obtain the phenotype right from the genotype by calling the
[`getPhenotype`](@ref) function (click to see documentation), and once they
get the phenotype, it can be evaluated by calling the [`evaluate`](@ref) function
(click to see documentation), which returns the result of the phenotype for a
specific set of values. In this example, a fitness function for a STGP problem
that creates decision trees for binary classification is shown:

```jldoctest
julia> function CCR(genotype::STGPGenotype, gpExperimentInfo::STGPInfo, classes::Array{Bool}, xValues::Array{Int64}, yValues::Array{Int64})

    rightGuesses = 0
    phenotype = getPhenotype(genotype); #Obtains the phenotype
    nValues = length(xValues)

    for i=1:nValues
        # Evaluates the phenotype for a set of the values of the variables x and y
        estimated_class = evaluate(phenotype, getVariables(gpExperimentInfo), xValues[i], yValues[i])
        if typeof(estimated_class) <: Bool
            if classes[i] == estimated_class
                rightGuesses += 1
            end
        end
    end

    # Calculates CCR, which is the fitness value of the individual
    fitness = (rightGuesses/nValues)*100
    return fitness
end
CCR (generic function with 1 method)

julia> classes = [true, true, false, true]
4-element Array{Bool,1}:
[...]

julia> xValues = [2, 1, 8, 5]
4-element Array{Bool,1}:
[...]

julia> yValues = [4, 6, 2, 5]
4-element Array{Bool,1}:
[...]

julia> FitnessFunction(CCR, classes, xValues, yValues, weight=1)
FitnessFunction(CCR, Any[Bool[1, 1, 0, 1], [2, 1, 8, 5], [4, 6, 2, 5]], 1.0)
```

As the reader can already guess, the fitness of the individuals of the population
will be their CCR, which is the number of patterns correctly classified. In order
to do so, we have created a fitness function that receives the mandatory arguments,
which are the individual's genotype and the GP experiment information, and, in
addition, it also receives the true classes and the values for the two variable
nodes of the problem, which are \"x\" and \"y\".

To evaluate the code generated from a GP individual, first, the user shall obtain
its phenotype, which is the Julia code translated from a GP tree. Then, the user
can evaluate the phenotype using the actual values as it is shown in the fitness
function by calling `evaluate`.

Once the fitness function and the arguments it is going to receive are created,
the user can create its instance of `FitnessFunction` as it is shown.

See also: [`Evaluator`](@ref)
"
struct FitnessFunction
    _fitnessFunction::Function
    _varArgs::Array
    _weight::Float32

    function FitnessFunction(fitnessFunction::Function, varArgs...;
                             weight::Real = 1.0)
        varArgs = [i for i in varArgs]
        new(fitnessFunction, varArgs, weight)
    end
end # struct

"""
    getFunction(fitnessfunction::FitnessFunction)::Function

Obtains the function of a fitness function.
"""
getFunction(fitnessfunction::FitnessFunction)::Function = fitnessfunction._fitnessFunction
# function

"""
    getFunctionArgs(fitnessfunction::FitnessFunction)::Array

Obtains the arguments of a fitness function.
"""
getFunctionArgs(fitnessfunction::FitnessFunction)::Array = fitnessfunction._varArgs
# function

"""
    getWeight(fitnessFunction::FitnessFunction)::Float32

Obtains the weight of a fitness function.
"""
getWeight(fitnessFunction::FitnessFunction)::Float32 = fitnessFunction._weight
# function

"""
    getWeight(fitnessFunctions::Array{FitnessFunction}, index::Integer=1)::Float32

Obtains the weight of a fitness function within an array of fitness functions.
"""
getWeight(fitnessFunctions::Array{FitnessFunction}, index::Integer=1)::Float32 = getWeight(fitnessFunctions[index])
# function

"""
    getWeights(fitnessFunctions::Array{FitnessFunction})::Array{Float32}

Obtains the weights of all the fitness function of an array.
"""
function getWeights(fitnessFunctions::Array{FitnessFunction})::Array{Float32}
    len = length(fitnessFunctions)
    weights = Array{Float32}(undef, len)
    for i=eachindex(weights)
        @inbounds weights[i] = getWeight(fitnessFunctions, i)
    end
    return weights
end # function

"""
    fitnessEval(fitnessFunc::FitnessFunction, genotype::Genotype)::Float64

Returns the result of a fitness function for a determined genotype.
"""
function fitnessEval(fitnessFunc::FitnessFunction, genotype::Genotype)::Float64
    ff = getFunction(fitnessFunc)
    args = getFunctionArgs(fitnessFunc)
    ff(genotype, args...)
end # function

"""
    fitnessEval(fitnessFunc::FitnessFunction,
                gpExperimentInfo::GPExperimentInfo, genotype::Genotype)::Float64

Returns the result of a fitness function for a determined genotype.
"""
function fitnessEval(fitnessFunc::FitnessFunction,
                     gpExperimentInfo::GPExperimentInfo, genotype::Genotype)::Float64
    ff = getFunction(fitnessFunc)
    args = getFunctionArgs(fitnessFunc)
    ff(genotype, gpExperimentInfo, args...)
end # function

"
This struct evaluates each individual according to a fitness function and an
assigned weight to each value of fitness.

# Fields
- `_fitnessFunctions::Array{FitnessFunction}`: set of fitness functions that are
    going to be evaluated for every individual of the population.
- `_globalFitnessFunction::FitnessFunction`: second-phase fitness function. If set,
    this function will be evaluated once all the individuals of the population
    have been evaluated with their raw fitness values.
- `_compareFunction::Function`: method used for comparison between individuals.
- `_compareFunctionArgs::String`: tells the evaluator what is going to be
    compared. \"raw\" is for fitness comparison, \"rawW\" is for weighted fitness comparison,
    \"global\" is for second-phase fitness comparison, \"globalW\" is for weighted
    second-phase fitness comparison, \"both\" is for both fitness and second-phase
    fitness comparison and \"bothW\" is for both weighted fitness and weighted
    second-phase fitness comparison.

If the user wants to set a second-phase fitness function, there are some
considerations to be taken.

The global fitness function must **always** receive the population of individuals
as its first argument and the fitness weights as its second argument, as they are
used to determine which fitness values weight more and which ones weight less.
The user can specify other arguments needed when creating the `FitnessFunction`
instance, but they do **not** need to specify the population or the weights as
one of them, as the system already know they are. Also, the function must return
an array with the second-phase fitness values of all the individuals in the
population. Here is an example of how a second-phase fitness function should be
created:

```jldoctest
julia> function myGlobalFitnessFunction(population::Array{Individual}, weights::Array{Real})
    #···Code···
    return globalFitnessValues
end
myGlobalFitnessFunction (generic function with 1 method)

julia> FitnessFunction(myGlobalFitnessFunction)
FitnessFunction(myGlobalFitnessFunction, Any[], 1.0)
```

Here, a second-phase fitness function has been created using `myGlobalFitnessFunction`
as the actual function, which takes no arguments, and assigning 1.0 as its weight.
The way the system know this is the second-phase fitness function is
by specifying so when calling `setEvaluator` with the keyword argument
`globalFitnessFunction`.

See also: [`setEvaluator`](@ref), [`FitnessFunction`](@ref)
"
struct Evaluator
    _fitnessFunctions::Array{FitnessFunction}
    _globalFitnessFunction::FitnessFunction
    _compareFunction::Function
    _compareFunctionArgs::String
end # struct

"""
    getNumberFitness(evaluator::Evaluator)

Obtains the number of fitness functions of the problem
"""
getNumberFitness(evaluator::Evaluator) = length(evaluator._fitnessFunctions)
# function

"""
    getFitnessFunction(evaluator::Evaluator, index::Integer)::FitnessFunction

Obtains a fitness function of the evaluator.
"""
getFitnessFunction(evaluator::Evaluator, index::Integer)::FitnessFunction = evaluator._fitnessFunctions[index]
# function

"""
    hasGlobalFitnessFunction(evaluator::Evaluator)::Bool

Checks if a global fitness function has been assigned.
"""
function hasGlobalFitnessFunction(evaluator::Evaluator)::Bool
    return getFunction(evaluator._globalFitnessFunction) != noFunc
end
# function

"""
    getGlobalFitnessFunction(evaluator::Evaluator)::FitnessFunction

Obtains the global fitness function of the evaluator.
"""
getGlobalFitnessFunction(evaluator::Evaluator)::FitnessFunction = evaluator._globalFitnessFunction
# function

"""
    getWeight(evaluator::Evaluator, index::Integer=1)::Float32

Obtains the weight of a fitness function of the evaluator.
"""
getWeight(evaluator::Evaluator, index::Integer=1)::Float32 = getWeight(evaluator._fitnessFunctions, index)
# function

"""
    getWeights(evaluator::Evaluator)::Array{Float32}

Obtains the weights of all the fitness function of the evaluator.
"""
getWeights(evaluator::Evaluator)::Array{Float32} = getWeights(evaluator._fitnessFunctions)
# function

"""
    getCompareFunction(evaluator::Evaluator)::Function

Obtains the method used for comparison of the evaluator.
"""
getCompareFunction(evaluator::Evaluator)::Function = evaluator._compareFunction
# function

"""
    getParametrizedCompareFunction(compareFunction::Function, args::String,
                                   weights::Array{Real})::Function

Obtains the comparison function based on the weights and the number of fitness
functions of the evaluator.
"""
function getParametrizedCompareFunction(compareFunction::Function, args::String,
                                        weights::Array{Float32})::Function

    local compare
    #compare = () -> nothing
    if compareFunction == noFunc
        if length(weights) == 1
            if sign(weights[1]) == 1
                compare = (ind1::Individual, ind2::Individual) -> compareFitnessMaximize(getFitness(ind1, all=true), getFitness(ind2, all=true))
            else
                compare = (ind1::Individual, ind2::Individual) -> compareFitnessMinimize(getFitness(ind1, all=true), getFitness(ind2, all=true))
            end
        else
            compare = (ind1::Individual, ind2::Individual) -> compareFitnessLinearCombination(getFitness(ind1, all=true), getFitness(ind2, all=true), weights)
        end
    else
        weights = length(weights) == 1 ? weights[1] : weights
        if args == "raw"
            compare = (ind1::Individual, ind2::Individual) -> compareFunction(getFitness(ind1, all=true), getFitness(ind2, all=true))
        elseif args == "rawW"
            compare = (ind1::Individual, ind2::Individual) -> compareFunction(getFitness(ind1, all=true), getFitness(ind2, all=true), weights)
        elseif args == "global"
            compare = (ind1::Individual, ind2::Individual) -> compareFunction(getGlobalFitness(ind1), getGlobalFitness(ind2))
        elseif args == "globalW"
            compare = (ind1::Individual, ind2::Individual) -> compareFunction(getGlobalFitness(ind1), getGlobalFitness(ind2), weights)
        elseif args == "both"
            compare = (ind1::Individual, ind2::Individual) -> compareFunction(getFitness(ind1, all=true), getGlobalFitness(ind1), getFitness(ind2, all=true), getGlobalFitness(ind2))
        else
            compare = (ind1::Individual, ind2::Individual) -> compareFunction(getFitness(ind1, all=true), getGlobalFitness(ind1), getFitness(ind2, all=true), getGlobalFitness(ind2), weights)
        end
    end
    return compare
end # function

"""
    computeFitness(evaluator::Evaluator, experimentInfo::ExperimentInfo,
                   genotype::Genotype, stopCondition::StopCondition)::Array{Float64}

Computes the fitness of a single individual.
"""
function computeFitness(evaluator::Evaluator, experimentInfo::ExperimentInfo,
                        genotype::Genotype, stopCondition::StopCondition)::Array{Float64}

    len = length(evaluator._fitnessFunctions)
    fitnesses = Array{Float64}(undef,len)

    if experimentInfo._individualType <: GAGenotype
        for i=eachindex(fitnesses)
            @inbounds fitnesses[i] = fitnessEval(getFitnessFunction(evaluator, i), genotype)
        end
    else
        for i=eachindex(fitnesses)
            @inbounds fitnesses[i] = fitnessEval(getFitnessFunction(evaluator, i), experimentInfo._GPExperimentInfo, genotype)
        end
    end

    if stopCondition._maxEvaluations != -1
        notifyEvaluation(stopCondition)
    end

    return fitnesses
end # function

"""
    globalFitnessEval(fitnessFunc::FitnessFunction,
                      population::Array{Individual},
                      weights::Array{Real})::Float64

Returns the result of a global fitness function for a population.
"""
function globalFitnessEval(fitnessFunc::FitnessFunction,
                           population::Array{Individual},
                           weights::Array{Float32})::Array{Float64}

    ff = getFunction(fitnessFunc)
    args = getFunctionArgs(fitnessFunc)
    ff(population, weights, args...)
end # function

"""
    evaluate!_(evaluator::Evaluator, experimentInfo::ExperimentInfo,
               population::Array{Individual}, stopCondition::StopCondition)

Evaluates the fitness of the individuals of a population and sets their values.
"""
function evaluate!_(evaluator::Evaluator, experimentInfo::ExperimentInfo,
                    population::Array{Individual}, stopCondition::StopCondition)
    for ind in population
        if hasNullFitness(ind)
            ind._fitness = computeFitness(evaluator, experimentInfo, ind._genotype, stopCondition)
        end
    end

    globalFitnessFunc = getGlobalFitnessFunction(evaluator)
    if getFunction(globalFitnessFunc) != noFunc
        fitnesses = getFitness(population)
        weights = getWeights(evaluator)
        globalFitnesses = globalFitnessEval(globalFitnessFunc, population, weights)
        setGlobalFitness!(population, globalFitnesses)
    end
    nothing
end # function

"""
    compareFitnessMaximize(fitness1::Float64, fitness2::Float64)

Compares two values of fitness and returns true if the first one is greater than
the second one, returns false otherwise.
"""
function compareFitnessMaximize(fitness1::Float64, fitness2::Float64)
    return fitness1 > fitness2 ? true : false
end # function

"""
    compareFitnessMinimize(fitness1::Float64, fitness2::Float64)

Compares two values of fitness and returns true if the first one is less than
the second one, returns false otherwise.
"""
function compareFitnessMinimize(fitness1::Float64, fitness2::Float64)
    return fitness1 < fitness2 ? true : false
end # function

"""
    compareFitnessLinearCombination(fitness1::Array{Float64},
                                    fitness2::Array{Float64},
                                    weights::Array{Float32})

Compares two sets of fitness values calculating the dot product between each one
of them and the weights assigned. Returns true if the first set of fitness values
is greater than the second one, returns false otherwise.
"""
function compareFitnessLinearCombination(fitness1::Array{Float64},
                                         fitness2::Array{Float64},
                                         weights::Array{Float32})

    function dot(array1::Array, array2::Array)
        sum(array1 .* array2)
    end

    conjunction1 = dot(fitness1, weights)
    conjunction2 = dot(fitness2, weights)

    return conjunction1 > conjunction2 ? true : false
end # function

"""
    getBestIndividual(population::Array{Individual}, compare::Function,
                      n::Integer = 1)::Array{Int}

Obtains the `n` best individuals of a population according to a comparing function.
"""
function getBestIndividual(population::Array{Individual}, compare::Function,
                           n::Integer = 1)::Array{Int}

    bestIndividual = collect(1:n)
    len = length(population)

    @inbounds for i=1:n
        for j=i+1:n
            if compare(population[bestIndividual[i]], population[bestIndividual[j]])
                bestIndividual[j], bestIndividual[i] = bestIndividual[i],  bestIndividual[j]
            end
        end
    end

    @inbounds for i=n+1:len
        if compare(population[i], population[bestIndividual[1]])
            bestIndividual[1] = i
            for j=1:n-1
                if compare(population[bestIndividual[j]], population[bestIndividual[j+1]]) || break
                    bestIndividual[j], bestIndividual[j+1] = bestIndividual[j+1],  bestIndividual[j]
                end
            end
        end
    end

    return bestIndividual
end # function

"""
    getWorstIndividual(population::Array{Individual}, compare::Function,
                       n::Integer=1)::Array{Int}

Obtains the `n` worst individuals of a population according to a comparing function.
"""
function getWorstIndividual(population::Array{Individual}, compare::Function,
                            n::Integer=1)::Array{Int}

    worstIndividual = collect(1:n)
    len = length(population)

    @inbounds for i=1:n
        for j=i+1:n
            if !compare(population[worstIndividual[i]], population[worstIndividual[j]])
                worstIndividual[j], worstIndividual[i] = worstIndividual[i],  worstIndividual[j]
            end
        end
    end

    @inbounds for i=n+1:len
        if !compare(population[i], population[worstIndividual[1]])
            worstIndividual[1] = i
            for j=1:n-1
                if !compare(population[worstIndividual[j]], population[worstIndividual[j+1]]) || break
                    worstIndividual[j], worstIndividual[j+1] = worstIndividual[j+1],  worstIndividual[j]
                end
            end
        end
    end

    if n == 1
        worstIndividual = worstIndividual[1]
    end

    return worstIndividual
end # function

"""
    pareto(population::Array{Individual}, weights::Array{Float32})::Array{Float64}

Assigns every individual of a population its second-phase fitness which is the
level they are in the pareto front.

# Arguments
- `population::Array{Individual}`: population of individuals that are going to be
    second-phase evaluated.
- `weights::Array{Real}`: array of weights for every single "raw fitness"
    function.

# Returns
An array with the levels of the individuals in the pareto front.
"""
function pareto(population::Array{Individual}, weights::Array{Float32})::Array{Float64}

    weightsSign = sign.(weights)

    function dominates(fitness1, fitness2)
        nFitness = length(fitness1)
        nDominations1 = 0
        nDominations2 = 0

        fitness1 = fitness1 .* weightsSign
        fitness2 = fitness2 .* weightsSign

        for i=1:nFitness
            if fitness1[i] > fitness2[i]
                nDominations1 += 1
            elseif fitness1[i] < fitness2[i]
                nDominations2 += 1
            end
        end
        if nDominations1 == nFitness
            return 1
        elseif nDominations2 == nFitness
            return -1
        end

        return 0
    end # function

    popsize = length(population)
    fitnesses = getFitness(population)
    dominateMatrix = zeros(Integer, popsize, popsize+1)
    dominateMatrix[:,end] = collect(1:popsize)

    @inbounds for i=1:popsize
        for j=i+1:popsize
            dominateMatrix[i,j] = dominates(fitnesses[i,:], fitnesses[j,:])
            dominateMatrix[j,i] = dominateMatrix[i,j] * (-1)
        end
    end

    paretoLevel = zeros(Integer, popsize)

    level = 0

    while (true)
        level+=1
        len = size(dominateMatrix)[1]
        if len == 0
            break
        end
        rows = ones(Bool, len)
        cols = ones(Bool, len+1)
        @inbounds for i=1:len
            n = count(i->(i==-1), view(dominateMatrix, i, :))
            if n == 0
                paretoLevel[dominateMatrix[i,end]] = level
                rows[i] = false
                cols[i] = false
            end
        end
        dominateMatrix = dominateMatrix[rows, cols]
    end

    return paretoLevel
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
# development script for EconDatasets

using DynAssMgmt
using EconDatasets
using TimeSeries

cd("/home/chris/scalable/julia/EconDatasets.jl/")

# download data
getDataset("IndustryPfs")

# load data
xxRets = dataset("IndustryPfs")

# store with information regarding the return type
retType = ReturnType(true, false, Dates.Day(1), false)
rets = Returns(xxRets, retType)

# derive associated prices
synthPrices = rets2prices(rets, 1.0, true)
logSynthPrices = getLogPrices(synthPrices)

function getLogPrices(prices::TimeSeries.TimeArray)
    return TimeSeries.TimeArray(prices.timestamp, log.(prices.values), prices.colnames)
end

# plot prices over time
DynAssMgmt.tsPlot(logSynthPrices)

# visualize universe
ewmaEstimator = EWMA(1, 1)
thisUniv = apply(ewmaEstimator, rets)
Plots.plot(thisUniv, rets.data.colnames)

## define efficient frontier / diversfication frontier strategies
DynAssMgmt.getUnivExtrema(thisUniv)
sigTargets = [linspace(1., 1.2, 15)...]

# get efficient frontier
effFrontStrats = EffFront(10)
effFrontWgts = apply(effFrontStrats, thisUniv)

# get as strategy types
diversTarget = 0.8
divFrontStrats = DivFront(diversTarget, sigTargets)
divFrontWgts = apply(divFrontStrats, thisUniv)

## mu/sigma results for full series of portfolios
DynAssMgmt.vizPfSpectrum(thisUniv, effFrontWgts[:])
DynAssMgmt.vizPfSpectrum!(thisUniv, divFrontWgts[:])

##
DynAssMgmt.wgtsOverStrategies(divFrontWgts)


text:
function mark_labyrinth!(r)

    if isborder(r,side)

    else
        putmarker!(r)
        for side in (Nord, West, Sud, Ost)
            try_move!(r, side)
            mark_labyrinth!(r)
            move!(r, inverse(side))
        end
    end
end

function try_move!(r,side)
    if isborder(r,side)
        return
    else
        move!(r,side)
    end
end

max_stars_repo_name: NikitaDestroyer/Basics

max_issues_repo_name: NikitaDestroyer/Basics

max_forks_repo_name: NikitaDestroyer/Basics


text:
using TiledIteration: TileIterator
using FFTW: fft, dct

function upsample(x::AbstractArray{T,D}, factor::NTuple{D}, offset::NTuple{D} = (fill(0,D)...,)) where {T,D}
    @assert all(0 .<= offset .< factor) "offset is out of range"
    szout = size(x) .* factor
    setindex!(zeros(T, szout), x, StepRange.(offset .+ 1, factor, szout)...)
end

function downsample(x::AbstractArray{T,D}, factor::NTuple{D}, offset::NTuple{D}=(fill(0,D)...,)) where {T,D}
    @assert all(0 .<= offset .< factor) "offset is out of range"
    x[StepRange.(offset .+ 1, factor, size(x))...]
end

# function downsampledview(x::AbstractArray{T,D}, factor::NTuple{D}, offset::NTuple{D}=(fill(0,D)...,)) where {T,D}
#     @assert all(0 .<= offset .< factor) "offset is out of range"
#     @view x[StepRange.(offset .+ 1, factor, size(x))...]
# end

representationmatrix(f, sz::NTuple) = representationmatrix(f, sz...)
function representationmatrix(f::Function, sz::Integer...)
    hcat([ setindex!(zeros(sz), 1, idx) |> f |> vec for idx in 1:prod(sz) ]...)
end

Base.@pure function haarbasis(d::Integer)
    w = walsh(d)
    @views [ reshape(w[p,:], fill(2,d)...) |> Array for p in 1:size(w, 1)]
end

Base.@pure function walsh(n::Integer)
    # ifelse(n >= 0, sub_walsh(Val(n)), error("n must to be a positive")) # this code is not work correctly.
    if n >= 0; sub_walsh(Val(n)) else error("n must to be a positive") end
end

function sub_walsh(::Val{N}) where {N}
    w = sub_walsh(Val(N-1))
    return [ w w ; w -w ]
end
sub_walsh(::Val{0}) = 1

# matrix-formed CDFT operator for D-dimensional signal
cdftmtx(sz::NTuple) = cdftmtx(sz...)
cdftmtx(sz::Integer...) = cdftmtx(Float64, sz...)
cdftmtx(T::Type, sz::NTuple) = cdftmtx(T, sz...)
cdftmtx(::Type, sz::Integer...) = cdftmtx(Float64, sz...)

cdftmtx(::Type{Complex{T}}, sz...) where {T} = cdftmtx(T, sz...)

Base.@pure function cdftmtx(::Type{T}, sz::Integer...) where T<:AbstractFloat
    len = prod(sz)

    mtx = representationmatrix(x->fft(T.(x)), sz)
    rm = Diagonal(Complex{T}[ cis(-angle(mtx[n,end])/2) for n in 1:len ])

    complex(T).(rm * mtx / sqrt(len))
end

permdctmtx(sz::NTuple) = permdctmtx(sz...)
permdctmtx(sz::Integer...) = permdctmtx(Float64, sz...)
permdctmtx(T::Type, sz::NTuple) = permdctmtx(T, sz...)

permdctmtx(::Type, sz::Integer...) = permdctmtx(Float64, sz...)

Base.@pure function permdctmtx(::Type{T}, sz::Integer...) where T<:AbstractFloat
    mtx = representationmatrix(x->dct(T.(x)), sz)

    isevenids = map(ci->iseven(sum(ci.I .- 1)), CartesianIndices(sz)) |> vec
    permids = sortperm(isevenids; rev=true, alg=Base.DEFAULT_STABLE)

    @views vcat([ transpose(mtx[pi,:]) for pi in permids ]...)
end

function getMatrixB(P::Integer, angs::AbstractVector{T}) where T
    @assert (length(angs) == fld(P,4)) "mismatch number of channels"
    hP = fld(P,2)
    psangs = (2 .* angs .+ pi) ./ 4
    ss, cs = sin.(psangs), cos.(psangs)

    LC = map(ss, cs) do s, c
        [ (-1im*c) (-1im*s); c (-s) ]
    end
    LS = map(ss, cs) do s, c
        [ s c; (1im*s) (-1im*c) ]
    end

    pbm = ones(fill(hP % 2,2)...)
    C = cat(LC..., pbm; dims=[1,2])
    S = cat(LS..., 1im*pbm; dims=[1,2])

    [ C conj(C); S conj(S) ] / sqrt(convert(T,2))
end

function analysisbank(nsolt::AbstractNsolt)
    M = prod(decimations(nsolt))
    ord = orders(nsolt)

    # create inpulse signal matrix
    mtx0 = reverse(Matrix(I, M, M .* prod(ord .+ 1) ), dims=1)
    krncenter = initialStep(nsolt, mtx0 )

    nStrides = (cumprod([ M, (ord[1:end-1] .+ 1)... ])...,)
    rotdimsfcns = (fill(identity, ndims(nsolt))...,)
    krnsym = extendAtoms(nsolt, krncenter, nStrides, rotdimsfcns, border=:circular_traditional)

    return shiftFilterSymmetry(nsolt, krnsym)
end

# compatible mode for SaivDr
# function analysisbank_compatible(nsolt::AbstractNsolt)
# #function analysisbank(nsolt::AbstractNsolt)
#     M = prod(decimations(nsolt))
#     ord = orders(nsolt)
#
#     # create inpulse signal matrix
#     mtx0 = reverse(Matrix(I, M, M .* prod(ord .+ 1) ), dims=1)
#     # mtx0 = circshift(mtx0, (0, -M))
#     krncenter = initialStep(nsolt, mtx0 )
#
#     nStrides = (cumprod([ M, (ord[1:end-1] .+ 1)... ])...,)
#     rotdimsfcns = (fill(identity, ndims(nsolt))...,)
#     krnsym = extendAtoms(nsolt, krncenter, nStrides, rotdimsfcns, border=:circular_traditional)
#
#     return shiftFilterSymmetry(nsolt, krnsym)
# end

kernels(pfb::PolyphaseFB) = (analysiskernels(pfb), synthesiskernels(pfb))

function analysiskernels(pfb::PolyphaseFB)
    df = decimations(pfb)
    afb = analysisbank(pfb)

    @views map([ reshape(afb[p,:], prod(df), :) for p in 1:size(afb, 1)]) do vf
        out = similar(vf, kernelsize(pfb)...)
        for (idx, tile) in enumerate(TileIterator(axes(out), df))
            out[tile...] = reshape(vf[:, idx], df...)
        end
        out
    end
end

function synthesiskernels(cc::AbstractNsolt)
    map(analysiskernels(cc)) do af
        reshape(af .|> conj |> vec |> reverse, size(af))
    end
end

function mdarray2polyphase(x::AbstractArray{TX,D}, szBlock::NTuple{D,TS}) where {TX,D,TS<:Integer}
    nBlocks = fld.(size(x), szBlock)

    @assert all(size(x) .% szBlock .== 0) "size error. input data: $(size(x)), block size: $(szBlock)."

    # outdata = hcat([ vec(@view x[tile...]) for tile in TileIterator(axes(x), szBlock)]...)
    outdata = similar(x, prod(szBlock), prod(nBlocks))
    @views for (idx, tile) in enumerate(TileIterator(axes(x), szBlock))
        outdata[:,idx] = vec(x[tile...])
    end
    PolyphaseVector(outdata, nBlocks)
end

function polyphase2mdarray(x::PolyphaseVector{TX,D}, szBlock::NTuple{D,TS}) where {TX,D,TS<:Integer}
    @assert (size(x.data, 1) == prod(szBlock)) "size mismatch! 'prod(szBlock)' must be equal to $(size(x.data,1))."

    out = similar(x.data, (x.nBlocks .* szBlock)...)
    @views for (idx, tile) in enumerate(TileIterator(axes(out), szBlock))
        out[tile...] = reshape(x.data[:,idx], szBlock...)
    end
    out
end

function rotatedimspv(x::PolyphaseVector{T,D}) where {T,D}
    data = rotatedimspv(x.data, x.nBlocks[1])
    nBlocks = (x.nBlocks[2:end]..., x.nBlocks[1])
    return PolyphaseVector(data, nBlocks)
end

function irotatedimspv(x::PolyphaseVector{T,D}) where {T,D}
    data = irotatedimspv(x.data, x.nBlocks[1])
    nBlocks = (x.nBlocks[end], x.nBlocks[2:end]...)
    return PolyphaseVector(data, nBlocks)
end

function rotatedimspv(x::AbstractMatrix, nBlocks::Integer)
    @views hcat([ x[:, (1:nBlocks:end) .+ idx] for idx = 0:nBlocks-1 ]...)
end

irotatedimspv(x::AbstractMatrix, nBlocks::Integer) = rotatedimspv(x, fld(size(x, 2), nBlocks))

@inline function unnormalized_butterfly!(xu::T, xl::T) where {T<:AbstractMatrix}
    tu, tl = (xu + xl, xu - xl)
    xu .= tu
    xl .= tl
    nothing
end

@inline function half_butterfly!(xu::T, xl::T) where {T<:AbstractMatrix}
    tu, tl = (xu + xl, xu - xl) ./ 2
    xu .= tu
    xl .= tl
    nothing
end

function shiftcoefs!(V::Val{:circular}, k::Integer, mtxup::AbstractMatrix, mtxlw::AbstractMatrix, nShift::Integer)
    if isodd(k)
        shiftcoefs_odd!(V, mtxlw, nShift)
    else
        shiftcoefs_even!(V, mtxup, nShift)
    end
    nothing
end

function shiftcoefs_odd!(::Val{:circular}, mtx::AbstractMatrix, nShift::Integer)
    mtx .= circshift(mtx, (0, nShift))
end

shiftcoefs_even!(V::Val{:circular}, mtx, nShift) = shiftcoefs_odd!(V, mtx, -nShift)

adjshiftcoefs!(v::Val{:circular}, k, mtxup, mtxlw, nShift::Integer) = shiftcoefs!(v, k, mtxup, mtxlw, -nShift)

# function shiftcoefs!(::Val{:zero}, k::Integer, mtxup::AbstractMatrix, mtxlw::AbstractMatrix, nShift::Integer)
#     if isodd(k)
#         mtxlw[:, 1+nShift:end] .= @view mtxlw[:, 1:end-nShift]
#         mtxlw[:, 1:nShift] .= 0
#     else
#         mtxup[:, 1:end-nShift] .= @view mtxup[:, 1+nShift:end]
#         mtxup[:, end-nShift+1:end] .= 0
#     end
#     nothing
# end
#
# function adjshiftcoefs!(::Val{:zero}, k::Integer, mtxup::AbstractMatrix, mtxlw::AbstractMatrix, nShift::Integer)
#     if isodd(k)
#         mtxlw[:, 1:end-nShift] .= @view mtxlw[:, 1+nShift:end]
#         mtxlw[:, end-nShift+1:end] .= 0
#     else
#         mtxup[:, 1+nShift:end] .= @view mtxup[:, 1:end-nShift]
#         mtxup[:, 1:nShift] .= 0
#     end
#     nothing
# end

function shiftcoefs!(::Val{:circular_traditional}, ::Integer, ::Any, mtxlw::AbstractMatrix, nShift::Integer)
    mtxlw .= circshift(mtxlw, (0, nShift))
end

adjshiftcoefs!(v::Val{:circular_traditional}, k, mtxup, mtxlw, nShift::Integer) = shiftcoefs!(v, k, mtxup, mtxlw, -nShift)

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
# I.Y. Kim and O.L. de Weck.  Adaptive weighted-sum method for bi-objective optimization
# Pareto front generation. Structural and Multidisciplinary Optimization.
# Vol. 29, Num. 2, pp. 149–158, 2005.
# https://link.springer.com/article/10.1007/s00158-004-0465-1

using NSGAIII, PyPlot

function plot_pop(P)
    clf()
    p = plot(map(x -> x.y[1], P), map(x -> x.y[2], P), "bo", markersize=1)
    !isinteractive() && show()
    sleep(0.2)
end

const d = RealCoding(8, [-3, -3], [3, 3])
z1(x1, x2) = -(3(1-x1)^2 * exp(-x1^2 - (x2+1)^2) - 10(x1/5 - x1^3 - x2^5) * exp(-x1^2-x2^2) -3exp(-(x1+2)^2 - x2^2) + 0.5(2x1 + x2))
z2(x1, x2) = -(3(1+x2)^2 * exp(-x2^2 - (1-x1)^2) - 10(-x2/5 + x2^3 + x1^5) * exp(-x1^2-x2^2) - 3exp(-(2-x2)^2 - x1^2))
z(x) = begin
    x1, x2 = decode(x, d)
    z1(x1, x2), z2(x1, x2)
end
nsga(300, 20, ()->bitrand(d.nbbitstotal), z, 10, fplot = plot_pop)


text:
g = EvolvingGraph()

add_edge!(g, 1, 2, 1)
add_edge!(g, 1, 3, 2)
add_edge!(g, 2, 3, 3)

ns = breadth_first_impl(g, 2, 1)
nns = nodes(g)

g = IntAdjacencyList(4,3)

add_edge!(g, 1, 2, 1)
add_edge!(g, 1, 3, 2)
add_edge!(g, 2, 3, 3)
add_edge!(g, 4, 1, 1)
add_edge!(g, 4, 3, 2)
ns = breadth_first_impl(g, 1, 1)
@test (1,1) in keys(ns)
@test (3,2) in keys(ns)
@test (2,3) in keys(ns)
@test (3,3) in keys(ns)
@test !((4,2) in keys(ns))


text:
"""
    get_nbResPairs()

Returns the number of resonant pairs (n,np).
Here, the resonant pairs (n,np) must satisfy the 3 conditions:
- `l <= n <= lmax`.
- `-lmax <= np <= lmax`.
- `n` and `np` have the same parity.
"""
function get_nbResPairs()
    s = 0 # Initialising the counter
    for n=1:lmax # Loop over the index n !! ATTENTION, we have have n>0
        for np=-lmax:lmax # Loop over the index np
            if (np != 0 && mod(n-np,2) == 0) # For np to be allowed, we must have np != 0, and n,np of the same parity
                s += 1 # Updating the counter
            end
        end
    end
    return s # Output
end

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

"""
    nbResPairs

Total number of resonant pairs (n,np).
Use `get_nbResPairs()` to determine this quantity.
"""
const nbResPairs = get_nbResPairs()

"""
    tabResPairs

Static container for the resonance pairs (n,np)
Use `tabResPairs!()` to fill it in.
"""
const tabResPairs = zeros(Int64,2,nbResPairs)

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

"""
    tabResPairs!()

Fills in the table tabResPairs, of size 2*nbResPairs, of the resonance pairs (n,np)
"""
function tabResPairs!()
    iPair = 1 # Counter for the resonant pair
    for n=1:lmax # Loop over the index n !! ATTENTION, we have have n>0
        for np=-lmax:lmax # Loop over the index np
            if (np != 0 && mod(n-np,2) == 0) # For np to be allowed, we must have np != 0, and n,np of the same parity
                tabResPairs[1,iPair], tabResPairs[2,iPair] = n, np # Filling in the array with (n,np)
                iPair += 1 # Updating the pair
            end
        end
    end
end

##################################################
tabResPairs!() # Preparing the table of the resonance pairs to consider
##################################################

"""
    bisection(fun, xl, xu, <keyword arguments>)

Algorithm which looks for zeroes of the functions `fun` given in argument.
It is a simple bisection algorithm, but it makes no allocations and is sufficiently fast.
It allows us not to have to use the Roots library that makes a lot of allocations.
The optional tolerances are set to the same as the ones found in Roots.jl, and are most likely overkill.

The zero is searched within the interval `[xl, xu]`, whose bounds are given as arguments.
"""
function bisection(fun, xl::Float64, xu::Float64, tolx::Float64=4.0*eps(Float64), tolf::Float64=4.0*eps(Float64), iterMAX::Int64=50)
    if (xl > xu)
        xl, xu = xu, xl # Ordering the input arguments
    end
    #####
    fl, fu = fun(xl), fun(xu) # Evaluating the function on the bracket
    #####
    if (abs(fl) <= tolf) # We have already found a solution on the left bracket
        return xl # Returning the left bracket
    end
    #####
    if (abs(fu) <= tolf) # We have already found a solution on the right bracket
        return xu # Returning the right bracket
    end
    #####
    @assert fl*fu < 0.0 "bisection: NOT A BRACKET"
    #####
    iter = 0 # Counter for the iterations
    #####
    while true # Bisection loop
        #####
        xm = (xl+xu)*0.5 # Middle value
        #####
        if ((abs(xu-xl) <= tolx) || (iter > iterMAX)) # The considered bracket is smaller than the tolerance, or we have made too many iterations
            return xm # Returning the middle value
        end
        #####
        fm = fun(xm) # Value at the midpoint
        #####
        iter += 1 # Updating the counter of iterations
        #####
        if (abs(fm) <= tolf) # The middle value is below the threshold
            return xm # Returning the middle value
        end
        #####
        # Otherwise, we iterate the bisection
        if (fm*fl < 0.0) # One root can be found between the left point and the middle point
            xu, fu = xm, fm # The upper point becomes the midpoint
        else
            xl, fl = xm, fm # The lower point becomes the midpoint
        end
    end
end

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

"""
    tabSMARes_serial

Static container for the resonant lines, contains the edge of the domain [amin,amax] of the resonant line.
Use `tabSMARes!()` to fill it in.
"""
const tabSMARes_serial = zeros(Float64,2,2)
"""
    amin

Minimum semi-major axis for which the resonance condition is solved, corresponds to jlc(amin)=0.97 !!
ATTENTION, we make sure that jlc(amin)<1.0 to always have meaningful brackets in SMAResQ().

Set as `17*rg` where `rg` is the Schwartzschild radius.
"""
const amin = 17.0*rg # Minimum SMA for which the resonance condition is solved, corresponds to jlc(amin)=0.97 !! ATTENTION, we make sure that jlc(amin)<1.0 to always have meaningful brackets in SMAResQ()

"""
    amax

Maximum semi-major axis for which the resonance condition is solved
Set as `rh` where `rh` is given as an argument of the initial command line as `--rh`.
"""
const amax = rh

"""
    nbSMAResSEARCH

Number of semi-major axis used for the search of the range of SMA solving the resonance condition.
"""
const nbSMAResSEARCH = 1000

"""
    tabSMAResSEARCH

Table of semi-major axis, uniform in log-space, where we determine if the resonance condition can be solved in j-space.
ATTENTION, the table is sorted!
"""
const tabSMAResSEARCH = exp.(range(log(amin),length=nbSMAResSEARCH,log(amax)))

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

"""
    SMAResQ(np,ap,omega)

Checks if a resonance is possible for the semi-major axis ap with resonance number np.
Returns `true` if a given semi-major axis `ap` is such that there exists `jp` with `jlc(a)<jp<1.0` such that `np*nup(ap,jp) = omega`.
This function is very simple because jp->nup(ap,jp) is always decreasing.
"""
function SMAResQ(np::Int64,ap::Float64,omega::Float64)
    omegaleft = np*nup(ap,jlc(ap)) # Value of the resonant frequency on the very left
    omegaright = np*nup(ap,1.0) # Value of the resonant frequency on the very right
    return ( (omegaleft-omega)*(omegaright-omega) < 0.0) # If this product is negative, there is indeed a solution in the range jlc(ap)<j<1.0
end

"""
    getjRes(np,ap,omega)

Returns the value of `j` such that `np*nup(ap,jp) = omega`.
This function is very simple because jp->nup(ap,jp) is always decreasing.
This function will fail if the bisection condition is not satisfied: check with SMAResQ(np,ap,omega).
"""
function getjRes(np::Int64,ap::Float64,omega::Float64)
    return bisection(jp -> np*nup(ap,jp)-omega,jlc(ap),1.0) # We search for the root, by bisection, within the range jlc(ap)<j<1.0
end

"""
    tabSMARes!(n,np,a,j,tabSMARes=tabSMARes_serial)

Computes, if it exists, a given resonance line as an array of (ap,jp) locations in orbital space.
The values are written the the array tabSMARes, which is set to be tabSMARes_serial by default.

This function returns true if a line has been constructed and false otherwise.

# Remarks
- Part of the code could be refactored.
- The scan in SMAResSearch sometimes makes a lot of useless calculations.
"""
function tabSMARes!(n::Int64,np::Int64,a::Float64,j::Float64,tabSMARes::Array{Float64,2}=tabSMARes_serial)
    #####
    omega = n*nup(a,j) # Value of the resonant frequency to match
    #####
    # First, we scroll over tabSMAResSEARCH, to determine in that array
    # the smallest and largest SMA for which the resonance condition is solvable
    iSMAmin = [nbSMAResSEARCH+1, nbSMAResSEARCH+1]
    iSMAmax = [0, 0]

    #####
    take = -1 # Index to check which resonance line we are on
    nbLines = 0 # Number of resonance lines
    #####
    for iSMA = 1:nbSMAResSEARCH # Loop over the table of SMA used for the search of the SMA where the condition can be solved
        ap = tabSMAResSEARCH[iSMA] # Current value of ap for which we try to solve the resonance condition
        #####
        if(SMAResQ(np,ap,omega)) # The resonance condition can be solved for this value of ap
            #####
            if (take == -1) # Enter first resonance line
                take = 0
                nbLines = 1
            elseif (take == 1) # Enter second resonance line
                nbLines = 2
            end
            if (take == 0) # First resonance line
                if (iSMA < iSMAmin[1]) # We found a smaller value of ap for which the resonance condition can be solved
                    iSMAmin[1] = iSMA # Updating iSMAmin
                end
            elseif (take == 1) # Second resonance line
                if (iSMA < iSMAmin[2]) # We found a smaller value of ap for which the resonance condition can be solved
                    iSMAmin[2] = iSMA # Updating iSMAmin
                end
            end
            #####
            if (take == 0) # First resonance line
                if (iSMA > iSMAmax[1]) # We found a larger value of ap for which the resonance condition can be solved
                    iSMAmax[1] = iSMA # Updating iSMAmax
                end
            elseif (take == 1) # Second resonance line
                if (iSMA > iSMAmax[2]) # We found a smaller value of ap for which the resonance condition can be solved
                    iSMAmax[2] = iSMA # Updating iSMAmin
                end
            end
        else
            if (take == 0) # End of the first (main) resonance line
                take = 1
            end
        end
    end

    #####
    # We have not been able to find SMA for which the resonance condition can be solved
    if (nbLines == 0)
        return 0
    end
    #####
    # In all the other cases, we have been able to find a SMA-region

    for iLine=1:2 # Loop over the resonance lines (there are 1 or 2)

        if !(iLine == 2 && nbLines != 2) # There are two resonance lines

            if (iSMAmin[iLine] == 1) # The minimum SMA is exactly at the edge of the region
                tabSMARes[iLine,1] = amin # Lower edge of the considered region
            #####
            else # The lower edge is not at the edge of the SMA domain
                # We must now find on which side of the (ap,jp)-triangle the resonant line ends, in the range given by tabSMAResSEARCH
                SMAlower, SMAupper = tabSMAResSEARCH[iSMAmin[iLine]-1], tabSMAResSEARCH[iSMAmin[iLine]] # Range in SMA over which we must search for a solution
                #####
                # First, we try to see if we can find a solution on the right, i.e. along jp=1.0
                omegalower, omegaupper = np*nup(SMAlower,1.0)-omega, np*nup(SMAupper,1.0)-omega # Value of the resonance frequency along jp=1.0, in the SMA range [SMAlower,SMAupper]
                #####
                if (omegalower*omegaupper < 0.0) # We found a bracketing interval along jp=1.0, i.e. the resonant location is on the right
                    tabSMARes[iLine,1] = bisection(ap -> np*nup(ap,1.0)-omega,SMAlower,SMAupper) # Finding the solution along jp=1.0, in the SMA range [SMAlower,SMAupper]
                else # Else, we must search for the resonant location along jp=jlc(ap), i.e. the resonant location is on the left
                    tabSMARes[iLine,1] = bisection(ap -> np*nup(ap,jlc(ap))-omega,SMAlower,SMAupper) # Finding the solution along jp=jlc(ap), in the SMA range [SMAlower,SMAupper]
                end
            end
            #####
            if (iSMAmax[iLine] == nbSMAResSEARCH) # The maximum SMA is exactly at the edge of the region
                tabSMARes[iLine,2] = amax # Upper edge of the considered region
            #####
            else # The upper edge is not at the edge of the SMA domain
                # We must now find on which side of the (ap,jp)-triangle the resonant line ends, in the range given by tabSMA ResSEARCH
                SMAlower, SMAupper = tabSMAResSEARCH[iSMAmax[iLine]], tabSMAResSEARCH[iSMAmax[iLine]+1] # Range in SMA over which we must search for a solution
                #####
                # First, we try to see if we can find a solution on the right, i.e. along jp=1.0
                omegalower, omegaupper = np*nup(SMAlower,1.0)-omega, np*nup(SMAupper,1.0)-omega # Value of the resonance frequency along jp=1.0, in the SMA range [SMAlower,SMAupper]
                #####
                if (omegalower*omegaupper < 0.0) # We found a bracketing interval along jp=1.0, i.e. the resonant location is on the right
                    tabSMARes[iLine,2] = bisection(ap -> np*nup(ap,1.0)-omega,SMAlower,SMAupper) # Finding the solution along jp=1.0, in the SMA range [SMAlower,SMAupper]
                else # Else, we must search for the resonant location along jp=jlc(ap), i.e. the resonant location is on the left
                    tabSMARes[iLine,2] = bisection(ap -> np*nup(ap,jlc(ap))-omega,SMAlower,SMAupper) # Finding the solution along jp=jlc(ap), in the SMA range [SMAlower,SMAupper]
                end
            end

        else # There is only one resonance line
            tabSMARes[2,1] = 0.0
            tabSMARes[2,2] = 0.0
        end
    end

end