text:
################################################################################
#  The simplified computation of Im parts. we have integrated out qs & cos(θ)  #
################################################################################
@doc raw"""
    dkVIm(p0, ps, q0, k, m, T, Npi, lam4pik)

compute $\tilde{\partial_k}\mathrm{Im}V_k(q_0)$.

`dkVIm` only contains $V(q_0)$, no $V(-q_0)$

`dkVIm` contains type-1 delta function

`dkVIm` doesn't contains type-2 delta function
"""
function dkVIm(p0, ps, q0, k, m, T, Npi, lam4pik)
    lam4pik^2 *
    (2 + Npi) *
    π *
    (
        3 * (
            dkF1All(p0 - q0, ps, k, m, T) +
            dkF1All(p0 + q0, ps, k, m, T) +
            dkF2All(p0 - q0, ps, k, m, T) +
            dkF2All(p0 + q0, ps, k, m, T)
        )
        # +
        # (2 + Npi) * (
        #     dkF1All(1e-8 - 1e-14, 1e-8, k, m, T) +
        #     dkF2All(1e-8 - 1e-14, 1e-8, k, m, T)
        # )
    )
end

@doc raw"""
    dkVImintqs(p0, ps, q0, qsmax, k, m, T, Npi, lam4pik)

compute $\int_0^{qsmax}dq_s qs^2\int_{-1}^{1}d\cos\theta \tilde{\partial_k}\mathrm{Im}V(q_0)$.

`dkVImintqs` only contains $V(q_0)$, for $-q_0$, we have $\int d\cos\theta V(q_0)=\int d\cos\theta V(-q_0)$,
so we need an extra $2$ at somewhere.

`dkVImintqs` doesn't have any delta function contribution, we include the type-1 delta function in `VImintqs_delta1` separately.

# Arguments
- `qsmax`: we integrate $q_s$ from $0$ to $k$, `qsmax` will set to `k` when we do the integration $dk'$, it should be distinguished from $k'$
- `m`: mass square, it will be $m(k')$ when we do the integration $dk'$.
- `lam4pik`: $\lambda_{4\pi}$, it will be $\lambda_{4\pi}(k')$ when we do the integration $dk'$ .
"""
function dkVImintqs(p0, ps, q0, qsmax, k, m, T, Npi, lam4pik)
    lam4pik^2 *
    (2 + Npi) *
    π *
    (
        3 * (
            dkF1Allintqs(p0 - q0, ps, qsmax, k, m, T)
            +
            dkF1Allintqs(p0 + q0, ps, qsmax, k, m, T)
            +
            dkF2Allintqs(p0 - q0, ps, qsmax, k, m, T)
            +
            dkF2Allintqs(p0 + q0, ps, qsmax, k, m, T)
        )
    )
end

function dkVImintqs_gpu(p0, ps, q0, qsmax, k, m, T, Npi, lam4pik)
    lam4pik^2 *
    (2 + Npi) *
    π *
    (
        3 * (
            dkF1Allintqs_gpu(p0 - q0, ps, qsmax, k, m, T)
            +
            dkF1Allintqs_gpu(p0 + q0, ps, qsmax, k, m, T)
            +
            dkF2Allintqs_gpu(p0 - q0, ps, qsmax, k, m, T)
            +
            dkF2Allintqs_gpu(p0 + q0, ps, qsmax, k, m, T)
        )
    )
end

function dkVImintqs_F4(p0, ps, q0, qsmax, k, m, T, Npi, lam4pik)
    lam4pik^2 *
    (2 + Npi) *
    π *
    (
        3 * (
             dkF4Allintqs(p0 - q0, ps, qsmax, k, m, T)
             +
             dkF4Allintqs(p0 + q0, ps, qsmax, k, m, T)
        )
    )
end

function dkVImintqs_zeroMomentum(p0, ps, q0, qsmax, k, m, T, Npi, lam4pik)
    lam4pik^2 * (2 + Npi) * π * (Npi + 2) * 2 / 3 *
    qsmax^3 *
    (dkF1All(k, m, T) + dkF2All(k, m, T))
end

@doc raw"""
    VImintqs(p0, ps, k, T, Npi,UVScale,mfun::Function,lamfun::Function)

compute $\int_0^{k}dq_s qs^2\int_{-1}^{1}d\cos\theta \mathrm{Im}V(q_0,k)$.
In our code, we perform integration over `kprim`, `q0` & `qs` does not involved,
so `qs=k`, `q0=Epi(k, mfun(k))`.

`VImintqs` don't have any delta function contribution, we include the type-1 delta function in `VImintqs_delta1` separately.

# Arguments
- `mfun::Function`: $m^2(k)$, input from zero momentum result
- `lampifun::Function`: $\lambda_{4\pi}(k)$, input from zero momentum result.
"""
function VImintqs(p0, ps, k, T, Npi,IRScale,UVScale, mfun, lamfun)
    -hquadrature(
        kprim -> dkVImintqs(
            p0,
            ps,
            Epi(k, mfun(k)),
            k,
            kprim,
            mfun(kprim),
            T,
            Npi,
            lamfun(kprim),
        ),
        k,
        UVScale,
        rtol = 1e-4,
        atol = 1e-4,
        initdiv=200,
    )[1]
end

function VImintqs(p0, ps, q0,k, T, Npi,IRScale,UVScale, mfun, lamfun)
    -hquadrature(
        kprim -> dkVImintqs(
            p0,
            ps,
            q0,
            k,
            kprim,
            mfun(kprim),
            T,
            Npi,
            lamfun(kprim),
        ),
        k,
        UVScale,
        rtol = 1e-4,
        atol = 1e-4,
        initdiv=200,
    )[1]
end

function VImintqs_zeroMomentum(p0, ps, k, T, Npi,IRScale,UVScale, mfun, lamfun)
    -hquadrature(
        kprim -> dkVImintqs_zeroMomentum(
            p0,
            ps,
            Epi(k, mfun(k)),
            k,
            kprim,
            mfun(kprim),
            T,
            Npi,
            lamfun(kprim),
        ),
        k,
        UVScale,
        rtol = 1e-2,
        atol = 1e-2,
    )[1]
end

@doc raw"""
    VImintqs_delta1(p0, ps, k, T, Npi,IRScale,UVScale, mfun::Spline1D, lamfun::Spline1D)

compute type-1 delta function contribution in $\mathrm{Im}V_k$, we do the triple integral analytically.
"""
function VImintqs_delta1(p0, ps, k, T, Npi,IRScale,UVScale, mfun, lamfun)
    -(2 + Npi) *
    π *
    3 *
    (
        deltasumkAll(p0 + Epi(k, mfun(k)), ps, k, T, Npi,IRScale, UVScale, mfun, lamfun)
        +
        deltasumkAll(p0 - Epi(k, mfun(k)), ps, k, T, Npi,IRScale, UVScale, mfun, lamfun)
    )
end

@doc raw"""
    Coeffgamm2(k, T, mfun::Spline1D)

The coefficient before the two point function flow.
"""
function Coeffgamm2(k, T, mfun)
    (
        k * (
            -(coth(Epi(k, mfun(k)) / (2 * T)) / Epi(k, mfun(k))^3) -
            csch(Epi(k, mfun(k)) / (2 * T))^2 / (2 * T * Epi(k, mfun(k))^2)
        )
    ) / (16 * pi^2)
end

function propImintqs(p0, ps, T, IRScale, UVScale, Npi, mfun, lamfun)
    return -hquadrature(
        k ->
            2 *
            VImintqs(p0, ps, k, T, Npi, IRScale, UVScale, mfun, lamfun) *
            Coeffgamm2(k, T, mfun),
        IRScale,
        UVScale,
        rtol = 1e-4,
        atol = 1e-4,
        maxevals = 2000
    )[1]
end

function propImintqs_zeroMomentum(p0, ps, T, IRScale, UVScale, Npi, mfun, lamfun)
    return -hquadrature(
        k ->
            2 *
            VImintqs_zeroMomentum(p0, ps, k, T, Npi, IRScale, UVScale, mfun, lamfun) *
            Coeffgamm2(k, T, mfun),
        IRScale,
        UVScale,
        rtol = 1e-4,
        atol = 1e-4,
        maxevals = 2000
    )[1]
end

function propImintqs_delta1(p0, ps, T, IRScale, UVScale, Npi, mfun, lamfun)
    return -hquadrature(
        k ->
            2 *
            VImintqs_delta1(p0, ps, k, T, Npi, IRScale, UVScale, mfun, lamfun) *
            Coeffgamm2(k, T, mfun),
        IRScale,
        UVScale,
        rtol = 1e-6,
        atol = 1e-6,
        maxevals = 2000,
    )[1]
end

function propImintqs_All(p0, ps, T, IRScale, UVScale, Npi, mfun, lamfun)
    return propImintqs(p0, ps, T, IRScale, UVScale, Npi, mfun, lamfun) +
           propImintqs_delta1(p0, ps, T, IRScale, UVScale, Npi, mfun, lamfun)+propImintqs_zeroMomentum(p0, ps, T, IRScale, UVScale, Npi, mfun, lamfun)
end

function fastpropImintqs(
    p0,
    ps,
    T,
    IRScale,
    UVScale,
    Npi,
    mfun,
    lamfun;
    kwarg...,
)
    return -hcubature(
        x ->
            (x[1] - UVScale) *
            2 *
            dkVImintqs(
                p0,
                ps,
                Epi(x[1], mfun(x[1])),
                x[1],
                x[2] * (x[1] - UVScale) + UVScale,
                mfun(x[2] * (x[1] - UVScale) + UVScale),
                T,
                Npi,
                lamfun(x[2] * (x[1] - UVScale) + UVScale),
            ) *
            Coeffgamm2(x[1], T, mfun),
        [IRScale, 0.0],
        [UVScale, 1.0];
        kwarg...,
    )[1]
end

function fastpropImintqs_F4(
    p0,
    ps,
    T,
    IRScale,
    UVScale,
    Npi,
    mfun,
    lamfun;
    kwarg...,
)
    return -hcubature(
        x ->
            (x[1] - UVScale) *
            2 *
            dkVImintqs_F4(
                p0,
                ps,
                Epi(x[1], mfun(x[1])),
                x[1],
                x[2] * (x[1] - UVScale) + UVScale,
                mfun(x[2] * (x[1] - UVScale) + UVScale),
                T,
                Npi,
                lamfun(x[2] * (x[1] - UVScale) + UVScale),
            ) *
            Coeffgamm2(x[1], T, mfun),
        [IRScale, 0.0],
        [UVScale, 1.0];
        kwarg...,
    )[1]
end

function fastpropImintqs_zeroMomentum(
    p0,
    ps,
    T,
    IRScale,
    UVScale,
    Npi,
    mfun,
    lamfun;
    kwarg...,
)
    return -hcubature(
        x ->
            (x[1] - UVScale) *
            2 *
            dkVImintqs_zeroMomentum(
                p0,
                ps,
                Epi(x[1], mfun(x[1])),
                x[1],
                x[2] * (x[1] - UVScale) + UVScale,
                mfun(x[2] * (x[1] - UVScale) + UVScale),
                T,
                Npi,
                lamfun(x[2] * (x[1] - UVScale) + UVScale),
            ) *
            Coeffgamm2(x[1], T, mfun),
        [IRScale, 0.0],
        [UVScale, 1.0];
        kwarg...,
    )[1]
end

function fastpropImintqs_All(p0, ps, T, IRScale, UVScale, Npi, mfun, lamfun;kwarg...)
    return fastpropImintqs(p0, ps, T, IRScale, UVScale, Npi, mfun, lamfun;kwarg...) +
           propImintqs_delta1(p0, ps, T, IRScale, UVScale, Npi, mfun, lamfun)+fastpropImintqs_F4(
               p0,
               ps,
               T,
               IRScale,
               UVScale,
               Npi,
               mfun,
               lamfun;
               kwarg...,
           )
end

function propImintqs_F4simple(p0, ps, T, IRScale, UVScale, Npi, mfun, lamfun)
    if p0^2<=mfun(IRScale)
        return 0.0
    elseif p0^2>mfun(IRScale)
        return (5*(2+Npi)*(-mfun(IRScale)+p0^2)^3*csch(p0/(2*T))^4*lamfun(sqrt(-mfun(IRScale)+p0^2))^2*sinh(p0/T)*(p0+T*sinh(p0/T)))/(4096*p0^4*pi^2*pi*T)
    end
end

#integrate delta function in F1 , we integrate out qs, cos(θ) and k' this will be used in Im part calculation
function deltasumImprop1(p0, ps, T,IRScale,UVScale, Npi, mfun, lamfun)
    #find the location of k0 where p0==2Epi(k,m)
    deltaf(x) = 3 * Epi(x, mfun(x)) - p0
    if deltaf(IRScale) * deltaf(UVScale) >= 0
        return 0.0
    else
        k0 = find_zero(deltaf, (IRScale, UVScale))
        # println(" k0=", k0)
        # k0 should lies between k~Λ, when k>k0 you will get 0
        if IRScale > k0
            return 0.0
            #δ function only appears in p<2k so we have the following division
        elseif IRScale <= k0
            if ps < k0
                return (
                    3 *
                    k0^2 *
                    (2 + Npi) *
                    coth(p0 / (6 * T)) *
                    csch(p0 / (6 * T))^2 *
                    lamfun(k0)^2 *
                    (
                        3075 * k0^4 - 1900 * k0^2 * ps^2 +
                        137 * ps^4 +
                        60 * ps^2 * (-20 * k0^2 + ps^2) * log(k0 / ps)
                    ) *
                    sign(2 * k0 + derivative(mfun, k0)) *
                    (p0 + 3 * T * sinh(p0 / (3 * T)))
                ) / (
                    25600 * p0^4 * pi^3 * T * (2 * k0 + derivative(mfun, k0))^2
                )
            elseif k0 < ps <= 2 * k0
                return (
                    3 *
                    k0^2 *
                    (2 + Npi) *
                    coth(p0 / (6 * T)) *
                    csch(p0 / (6 * T))^2 *
                    lamfun(k0)^2 *
                    (
                        (3 * k0 - ps) * (
                            132 * k0^4 + 3099 * k0^3 * ps - 2047 * k0^2 * ps^2 -
                            9 * k0 * ps^3 + 137 * ps^4
                        ) +
                        60 *
                        (2 * k0 - ps)^3 *
                        (4 * k0^2 + 6 * k0 * ps + ps^2) *
                        log(k0 / (2 * k0 - ps))
                    ) *
                    sign(2 * k0 + derivative(mfun, k0)) *
                    (p0 + 3 * T * sinh(p0 / (3 * T)))
                ) / (
                    51200 *
                    p0^4 *
                    pi^3 *
                    ps *
                    T *
                    (2 * k0 + derivative(mfun, k0))^2
                )
            elseif 2*k0 < ps <= 3 * k0
                return (
                    -3 *
                    k0^2 *
                    (2 + Npi) *
                    coth(p0 / (6 * T)) *
                    csch(p0 / (6 * T))^2 *
                    lamfun(k0)^2 *
                    (
                        -(
                            (3 * k0 - ps) * (
                                132 * k0^4 + 3099 * k0^3 * ps -
                                2047 * k0^2 * ps^2 - 9 * k0 * ps^3 + 137 * ps^4
                            )
                        ) +
                        60 *
                        (2 * k0 - ps)^3 *
                        (4 * k0^2 + 6 * k0 * ps + ps^2) *
                        log(-2 + ps / k0)
                    ) *
                    sign(2 * k0 + derivative(mfun, k0)) *
                    (p0 + 3 * T * sinh(p0 / (3 * T)))
                ) / (
                    51200 *
                    p0^4 *
                    pi^3 *
                    ps *
                    T *
                    (2 * k0 + derivative(mfun, k0))^2
                )
            elseif ps > 3 * k0
                return 0.0
            end
        end
    end
end

function deltasumImprop2(p0, ps, T,IRScale,UVScale, Npi, mfun, lamfun)
    #find the location of k0 where p0==2Epi(k,m)
    deltaf(x) = Epi(x, mfun(x)) - p0
    if deltaf(IRScale) * deltaf(UVScale) >= 0
        return 0.0
    else
        k0 = find_zero(deltaf, (IRScale, UVScale))
        # println(" k0=", k0)
        # k0 should lies between k~Λ, when k>k0 you will get 0
        if IRScale > k0
            return 0.0
            #δ function only appears in p<2k so we have the following division
        elseif IRScale <= k0
            if ps < k0
                return (k0^2*(2+Npi)*csch(p0/(2*T))^4*lamfun(k0)^2*(3075*k0^4-1900*k0^2*ps^2+137*ps^4+60*ps^2*(-20*k0^2+ps^2)*log(k0/ps))*sign(2*k0+derivative(mfun,k0))*sinh(p0/T)*(p0+T*sinh(p0/T)))/(153600*p0^4*pi^3*T*(2*k0+derivative(mfun,k0))^2)
            elseif k0 < ps <= 2 * k0
                return (k0^2*(2+Npi)*csch(p0/(2*T))^4*lamfun(k0)^2*((3*k0-ps)*(132*k0^4+3099*k0^3*ps-2047*k0^2*ps^2-9*k0*ps^3+137*ps^4)+60*(2*k0-ps)^3*(4*k0^2+6*k0*ps+ps^2)*log(k0/(2*k0-ps)))*sign(2*k0+derivative(mfun,k0))*sinh(p0/T)*(p0+T*sinh(p0/T)))/(307200*p0^4*pi^3*ps*T*(2*k0+derivative(mfun,k0))^2)
            elseif 2*k0 < ps <= 3 * k0
                return (k0^2*(2+Npi)*csch(p0/(2*T))^4*lamfun(k0)^2*((3*k0-ps)*(132*k0^4+3099*k0^3*ps-2047*k0^2*ps^2-9*k0*ps^3+137*ps^4)+60*(2*k0-ps)^3*(4*k0^2+6*k0*ps+ps^2)*log(k0/(2*k0-ps)))*sign(2*k0+derivative(mfun,k0))*sinh(p0/T)*(p0+T*sinh(p0/T)))/(307200*p0^4*pi^3*ps*T*(2*k0+derivative(mfun,k0))^2)
            elseif ps > 3 * k0
                return 0.0
            end
        end
    end
end


text:
# load dependencies
using Temporal, Indicators, Base.Dates
charts = false

# define backtest setup
from_date = today() - Year(10)
thru_date = today()
trade_qty = 10.0
initial_balance = 1e6
trade_field = :Open
close_field = :Settle

# define indicator parameters
ind_fun = Indicators.mama
ind_vals(x::TS) = hl2(x)
ind_args = (:fastlimit, 0.5,
            :slowlimit, 0.05)
ind_args_rng = (:fastlimit, 0.01:0.01:0.99,
                :slowlimit, 0.01:0.01:0.99)

# download historical data
X = quandl("CHRIS/CME_CL1", from=string(from_date), thru=string(thru_date))
ind_out = ind_fun(ind_vals(X); ind_args)

# initialize backtest variables
N = size(X,1)
go_long = crossover(ind_out[:MAMA].values, ind_out[:FAMA].values)
go_short = crossunder(ind_out[:MAMA].values, ind_out[:FAMA].values)
pnl = zeros(Float64, N)
pos = zeros(Float64, N)
in_price = NaN
trade_price = X[:,trade_field].values
close_price = X[:,close_field].values

# run trading logic for base case (default indicator args)
@inbounds for i in 2:N
    if go_long[i-1]
        pos[i] = trade_qty
        pnl[i] = pos[i] * (close_price[i] - trade_price[i])
    elseif go_short[i-1]
        pos[i] = -trade_qty
        pnl[i] = pos[i] * (close_price[i] - trade_price[i])
    else
        pos[i] = pos[i-1]
        pnl[i] = pos[i] * (close_price[i]-close_price[i-1])
    end
end
# summarize results from backtest
summary_ts = TS([Temporal.ohlc(X).values ind_out.values go_long go_short pos pnl cumsum(pnl)],
                X.index,
                vcat(Temporal.ohlc(X).fields, [:MAMA, :FAMA, :LongSignal, :ShortSignal, :Position, :PNL, :CumulativePNL]))

# visualize backtest results
using Plots
ℓ = @layout [ a; b{0.33h} ]
plotlyjs()
plot(plot(summary_ts[:,vcat(close_field,ind_out.fields)], color=[:black :magenta :green]),
     plot(summary_ts[:,end], color=:orange, fill=(0,:orange), fillalpha=0.5),
     layout = ℓ)

ps = ParameterSet([:fastlimit, :slowlimit], [0.5, 0.05])
ps.arg_ranges = [0.01:0.01:0.99, 0.01:0.01:0.99]
params = get_run_params(ps)

#TODO: make this more sophisticated
n_runs = get_n_runs(ps)
cum_pnl = zeros(Float64, n_runs)
@inbounds for j in 1:n_runs
    println("Iteration $j/$n_runs ($(j/n_runs*100)%)")
    ind_out = ind_fun(ind_vals(X); params[j]...)
    # initialize backtest variables
    N = size(X,1)
    go_long = crossover(ind_out[:MAMA].values, ind_out[:FAMA].values)
    go_short = crossunder(ind_out[:MAMA].values, ind_out[:FAMA].values)
    pnl = zeros(Float64, N)
    pos = zeros(Float64, N)
    in_price = NaN
    trade_price = X[:,trade_field].values
    close_price = X[:,close_field].values
    @inbounds for i in 2:N
        if go_long[i-1]
            pos[i] = trade_qty
            pnl[i] = pos[i] * (close_price[i] - trade_price[i])
        elseif go_short[i-1]
            pos[i] = -trade_qty
            pnl[i] = pos[i] * (close_price[i] - trade_price[i])
        else
            pos[i] = pos[i-1]
            pnl[i] = pos[i] * (close_price[i]-close_price[i-1])
        end
    end
    cum_pnl[j] = sum(pnl)
end

ext: jl


text:
using PackageCompiler
create_sysimage(:Plots, sysimage_path="sys_plots.so", precompile_execution_file="precompile_plots.jl")


text:
Base.:(==)(d1::BDays252, d2::BDays252) = d1.hc == d2.hc
Base.hash(d::BDays252) = 1 + hash(d.hc)

daycount(conv::BDays252, date_start::Date, date_end::Date) = bdayscount(conv.hc, date_start, date_end)
daycount(::Actual360, date_start::Date, date_end::Date) = Int(Dates.value(date_end - date_start))
daycount(::Actual365, date_start::Date, date_end::Date) = Int(Dates.value(date_end - date_start))

advancedays(conv::BDays252, date_start::Date, daycount::Int) = advancebdays(conv.hc, date_start, daycount)
advancedays(::Actual360, date_start::Date, daycount::Int) = date_start + Day(daycount)
advancedays(::Actual365, date_start::Date, daycount::Int) = date_start + Day(daycount)

function advancedays(conv::DayCountConvention, date_start::Date, daycount_vec::Vector{Int})
    l = length(daycount_vec)
    result = Vector{Date}(undef, l)
    for i in 1:l
        @inbounds result[i] = advancedays(conv, date_start, daycount_vec[i])
    end
    return result
end

advancedays(curve::AbstractIRCurve, daycount) = advancedays(curve_get_daycount(curve), curve_get_date(curve), daycount)

daysperyear(::BDays252) = 252
daysperyear(::Actual360) = 360
daysperyear(::Actual365) = 365

yearfraction(conv::DayCountConvention, date_start::Date, date_end::Date) = daycount(conv, date_start, date_end) / daysperyear(conv)
yearfraction(curve::AbstractIRCurve, maturity::Date) = yearfraction(curve_get_daycount(curve), curve_get_date(curve), maturity)

function days_to_maturity(curve::AbstractIRCurve, maturity::Date)
    d = daycount(curve_get_daycount(curve), curve_get_date(curve), maturity)
    @assert d >= 0 "Maturity date $(maturity) should be greater than curve observation date $(curve_get_date(curve))"
    return d
end


text:
function face_sampling(n::Int,dim::Int, tol::T) where T

	X1 = [collect(range(-tol,stop=tol, length=n)) for i in 1:dim]
	X = Iterators.product(X1 ...)
	sp = vec(collect.(X))
	sp = [sp[i] for i in 1:Int(length(sp)) if any(sp[i] .== tol) || any(sp[i] .== -tol) || sp[i] == zeros(dim)]
	return unique(sp)

end

function uniform_sampling(n::Int,dim::Int)
	X1 = [collect(range(-1,stop=1, length=n)) for i in 1:dim]
	X = Iterators.product(X1 ...)
	sp = vec(collect.(X))::Vector{Vector{Float64}}
	return unique(sp)
end

# Samples only the center point
function center_sampling(n::Int, dim::Int)
	return [zeros(dim)]
end

function sampleX(sampling_type::Symbol, ns::Int, dim::Int, T)

	if sampling_type == :uniform
		return uniform_sampling(ns, dim)
	end

	if sampling_type == :face
		return face_sampling(ns, dim, convert(T,0.98))
	end

	if sampling_type == :random
		return LHC_sampling(ns, dim)
	end

	if sampling_type == :center
		return center_sampling(ns, dim)
	end

	error("Unknown sampling method type. Please leave default for uniform sampling or select either uniform or face or random.")
end

# Auxiliary function for sampling inputs and disturbances
function sampleUW(n::Int, dim::Int, ub::Vector{T}, lb::Vector{T}) where T
	X1 = [collect(range(lb[i],stop=ub[i], length=n)) for i in 1:dim]
	X = Iterators.product(X1 ...)
	sp = vec(collect.(X))::Vector{Vector{T}}
	return unique(sp)
end

function sample_points!(y::Vector{Vector{T}}, bx::AbstractHyperrectangle, Xsp::Vector{Vector{T}}) where T
	@simd for i in 1:length(Xsp)
		@inbounds y[i] .= Xsp[i] .* radius_hyperrectangle(bx) .+ LS.center(bx) # [sp[i][1]*box.radius[1]+box.center[1], sp[i][2]*box.radius[2]+box.center[2]]
	end
end


text:
export correctLayerDims

function correctLayerDims(orderGraph::Array,layersDims)
  for j=1:size(orderGraph,1)
    if sum(orderGraph[:,j])!=0
      layersDims[j,1]=0
    end
    for i=1:size(orderGraph,1)
      layersDims[j,1]+=layersDims[i,2]*orderGraph[i,j]
    end
  end
  layersDims
end


text:
# Autogenerated wrapper script for Leptonica_jll for powerpc64le-linux-gnu
export convertfilestops, convertfilestopdf, convertsegfilestops, xtractprotos, fileinfo, converttops, converttopdf, convertformat, imagetops, liblept, convertsegfilestopdf

using Giflib_jll
using JpegTurbo_jll
using libpng_jll
using Libtiff_jll
using libwebp_jll
using Zlib_jll
using OpenJpeg_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"

# Relative path to `convertfilestops`
const convertfilestops_splitpath = ["bin", "convertfilestops"]

# This will be filled out by __init__() for all products, as it must be done at runtime
convertfilestops_path = ""

# convertfilestops-specific global declaration
function convertfilestops(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(convertfilestops_path)
    end
end

# Relative path to `convertfilestopdf`
const convertfilestopdf_splitpath = ["bin", "convertfilestopdf"]

# This will be filled out by __init__() for all products, as it must be done at runtime
convertfilestopdf_path = ""

# convertfilestopdf-specific global declaration
function convertfilestopdf(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(convertfilestopdf_path)
    end
end

# Relative path to `convertsegfilestops`
const convertsegfilestops_splitpath = ["bin", "convertsegfilestops"]

# This will be filled out by __init__() for all products, as it must be done at runtime
convertsegfilestops_path = ""

# convertsegfilestops-specific global declaration
function convertsegfilestops(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(convertsegfilestops_path)
    end
end

# Relative path to `xtractprotos`
const xtractprotos_splitpath = ["bin", "xtractprotos"]

# This will be filled out by __init__() for all products, as it must be done at runtime
xtractprotos_path = ""

# xtractprotos-specific global declaration
function xtractprotos(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(xtractprotos_path)
    end
end

# Relative path to `fileinfo`
const fileinfo_splitpath = ["bin", "fileinfo"]

# This will be filled out by __init__() for all products, as it must be done at runtime
fileinfo_path = ""

# fileinfo-specific global declaration
function fileinfo(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(fileinfo_path)
    end
end

# Relative path to `converttops`
const converttops_splitpath = ["bin", "converttops"]

# This will be filled out by __init__() for all products, as it must be done at runtime
converttops_path = ""

# converttops-specific global declaration
function converttops(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(converttops_path)
    end
end

# Relative path to `converttopdf`
const converttopdf_splitpath = ["bin", "converttopdf"]

# This will be filled out by __init__() for all products, as it must be done at runtime
converttopdf_path = ""

# converttopdf-specific global declaration
function converttopdf(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(converttopdf_path)
    end
end

# Relative path to `convertformat`
const convertformat_splitpath = ["bin", "convertformat"]

# This will be filled out by __init__() for all products, as it must be done at runtime
convertformat_path = ""

# convertformat-specific global declaration
function convertformat(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(convertformat_path)
    end
end

# Relative path to `imagetops`
const imagetops_splitpath = ["bin", "imagetops"]

# This will be filled out by __init__() for all products, as it must be done at runtime
imagetops_path = ""

# imagetops-specific global declaration
function imagetops(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(imagetops_path)
    end
end

# Relative path to `liblept`
const liblept_splitpath = ["lib", "liblept.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
liblept_path = ""

# liblept-specific global declaration
# This will be filled out by __init__()
liblept_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const liblept = "liblept.so.5"

# Relative path to `convertsegfilestopdf`
const convertsegfilestopdf_splitpath = ["bin", "convertsegfilestopdf"]

# This will be filled out by __init__() for all products, as it must be done at runtime
convertsegfilestopdf_path = ""

# convertsegfilestopdf-specific global declaration
function convertsegfilestopdf(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(convertsegfilestopdf_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"Leptonica")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (Giflib_jll.PATH_list, JpegTurbo_jll.PATH_list, libpng_jll.PATH_list, Libtiff_jll.PATH_list, libwebp_jll.PATH_list, Zlib_jll.PATH_list, OpenJpeg_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (Giflib_jll.LIBPATH_list, JpegTurbo_jll.LIBPATH_list, libpng_jll.LIBPATH_list, Libtiff_jll.LIBPATH_list, libwebp_jll.LIBPATH_list, Zlib_jll.LIBPATH_list, OpenJpeg_jll.LIBPATH_list,))

    global convertfilestops_path = normpath(joinpath(artifact_dir, convertfilestops_splitpath...))

    push!(PATH_list, dirname(convertfilestops_path))
    global convertfilestopdf_path = normpath(joinpath(artifact_dir, convertfilestopdf_splitpath...))

    push!(PATH_list, dirname(convertfilestopdf_path))
    global convertsegfilestops_path = normpath(joinpath(artifact_dir, convertsegfilestops_splitpath...))

    push!(PATH_list, dirname(convertsegfilestops_path))
    global xtractprotos_path = normpath(joinpath(artifact_dir, xtractprotos_splitpath...))

    push!(PATH_list, dirname(xtractprotos_path))
    global fileinfo_path = normpath(joinpath(artifact_dir, fileinfo_splitpath...))

    push!(PATH_list, dirname(fileinfo_path))
    global converttops_path = normpath(joinpath(artifact_dir, converttops_splitpath...))

    push!(PATH_list, dirname(converttops_path))
    global converttopdf_path = normpath(joinpath(artifact_dir, converttopdf_splitpath...))

    push!(PATH_list, dirname(converttopdf_path))
    global convertformat_path = normpath(joinpath(artifact_dir, convertformat_splitpath...))

    push!(PATH_list, dirname(convertformat_path))
    global imagetops_path = normpath(joinpath(artifact_dir, imagetops_splitpath...))

    push!(PATH_list, dirname(imagetops_path))
    global liblept_path = normpath(joinpath(artifact_dir, liblept_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global liblept_handle = dlopen(liblept_path)
    push!(LIBPATH_list, dirname(liblept_path))

    global convertsegfilestopdf_path = normpath(joinpath(artifact_dir, convertsegfilestopdf_splitpath...))

    push!(PATH_list, dirname(convertsegfilestopdf_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
using Nanomsg
using Base.Test

const nn = Nanomsg

@testset "PAIR" begin
    s1 = nn.socket(nn.AF_SP, nn.PAIR)
    s2 = nn.socket(nn.AF_SP, nn.PAIR)
    nn.bind(s1, "inproc://test")
    nn.connect(s2, "inproc://test")
    msg = "hello, nanomsg"
    nn.send(s1, msg)
    buf = Vector{UInt8}(1024)
    sz = nn.recv(s2, buf)
    nn.close(s1)
    nn.close(s2)
    @test sz == sizeof(msg)
    @test String(buf[1:sz]) == msg
end

@testset "PUBSUB" begin
    pub = nn.socket(nn.AF_SP, nn.PUB)
    sub1 = nn.socket(nn.AF_SP, nn.SUB)
    sub2 = nn.socket(nn.AF_SP, nn.SUB)
    nn.bind(pub, "inproc://test")
    nn.connect(sub1, "inproc://test")
    nn.connect(sub2, "inproc://test")
    nn.setsockopt(sub1, nn.SUB, nn.SUB_SUBSCRIBE, "")
    nn.setsockopt(sub2, nn.SUB, nn.SUB_SUBSCRIBE, "")
    msg = "hello, nanomsg"
    nn.send(pub, msg)
    buf1 = Vector{UInt8}(1024)
    buf2 = Vector{UInt8}(1024)
    sz1 = nn.recv(sub1, buf1)
    sz2 = nn.recv(sub2, buf2)
    nn.close(pub)
    nn.close(sub1)
    nn.close(sub2)
    @test sz1 == sz2 == sizeof(msg)
    @test String(buf1[1:sz1]) == String(buf2[1:sz2]) == msg
end


text:
# =============================================================================
# CALCULATRICE SCIENTIFIQUE
# =============================================================================

# opérateurs  classiques : +-*/\

1 + 3
3 / 7
3 \ 7 # ≈ 7 / 3
1e2
im # nombre complexe i
im^2
# Les fonctions usuelles : exp, log, log10, abs, sqrt, cbrt, sign

log10(100_000)
log(exp(2))
1 / (1 + exp(1))
# Les fonctions trigonométriques usuelles : sin, cos, tan, cot, sec, csc, sinh, cosh, tanh, coth, sech, csch, asin, acos, atan...

sin(π+√2) # Code Is Math !!  tapez "\pi<tab>" et "\sqrt<tab>2"

# =============================================================================
# VALEURS (tout ce qui est manipulable en Julia) & TYPES (Int32, Int64, Int128, Float32, Float64, Complex32, Complex64, Complex128...)
# =============================================================================

# types et collections

3 + 5 # nombre entier stockée en 32/64 bits en fonction de l'architecture machine
typeof(3 + 5)
3. * 4. # nombre en virgule flottante stockée en 32/64 bits en fonction de l'architecture machine
typeof(3. * 4.)
typeof(3. * 4)
typeof(Float64)
π # tapez "\pi<tab>"
typeof(π)
Float64(π)
Float32(π)
typeof(1 + 2im)
[1,2,3]
typeof([1,2,3])
1:5 # itérateur
typeof(1:5)
collect(1:5)
0:0.1:1
typeof(0:0.1:1)
collect(0:0.1:1)
[4. 5. 6.;7. 8. 9.]
(1, 2.0, "test") # tuple
typeof((1, 2.0, "test"))

# Caractères et chaînes de caractères

'A'
typeof('A')
"hello"

typeof("hello")

"racine carrée de 2 : $(sqrt(2))"

" √3 : $(√3)"

"Hello" * " John Doe"

# Boléens et opérateurs de comparaison | "==" | "!=" | > | ">=" | "<" | "<=" | et de conjonction | "&" | "|" | "!"|

1 > 2
typeof(1 > 2)

# supertype et sous-types

supertype(Int64)
subtypes(Integer)

# type utilisateur

struct MyComplex
    _real::Float64
    _imag::Float64
end
typeof(MyComplex)
MyComplex(1,2)
typeof(MyComplex(1, 2))

# =============================================================================
# REGLE DE NOMMAGE(variables, fonctions/méthodes, modules, types, macros)
# =============================================================================
#=
- noms des variables  en minuscule, la séparation des mots dans une variable se fait à l'aide d'un underscore ('_') ("snake case")
- noms de type et de modules commencent par une lettre majuscule, les majuscules séparant les différents mots du nom (exemple "MonModule") ("camel case")
- noms des fonctions et macros en minuscule sans underscores.
- une fonction qui modifie ses arguments s'écrit avec ! à la fin. =#

# =============================================================================
# VARIABLES (étiquettes associées aux valeurs)
# =============================================================================

# Les variables contenant des nombres scalaires

α = 6 # variable latex-unicode  names — tapez "\alpha<tab>"
typeof(α)
β = 2α + 1.5
typeof(β)
γ = α + β * im
typeof(γ)
δ = Float64(10)
δ = MyComplex(1, 2)
δ._real

# Les vecteurs (tableau de nombres scalaires ou non de dim=1)

v = [1, 17, 32, 15]
@show v
size(v)
n = length(v)
eltype(v)
v[2] # accès à l'élément d'index n°2 (indice de départ = 1)
for i in 1:n
	println("v[$i] = $(v[i])")
end
for elt in v
	println(elt)
end
v[2:3] # slicing
v' # transposée de x
v2 = Float64[1, 2, 3, 4]
v + v2
v' * v2
sin.(v) # opérateur . de broadcast ( sin de chaque terme)
v .* v2 # opérateur . de broadcast (ici multiplication terme à terme)
v3 = ones(Int64, 3)
v4 = []
typeof(v4)
push!(v4,1)     # => [1]
push!(v4,2.)     # => [1,2.]
push!(v4,"test")     # => [1,2, "test"]
append!(v4,v4)

# Les matrices (tableau de nombres scalaires ou non dim=2 )

A = [1. 2 3
     4 5 6
     7 8 10]
ni, nj = size(A)
length(A)
eltype(A)
A[2,3]
A[end,end]
A[2,:]
A[:,2]
A[2:3,2:3]
A[1:end]
for j in 1:nj
    for i in 1:ni
        println("A[$i,$j] = $(A[i,j])")
    end
end
for elt in A
	println(elt)
end
As = similar(A)
B = fill(2.0, (3, 3))
C = A + B
C = A * B
sin.(A)
AA = [A  A] # concaténation horizontale par bloc (hcat)
AAAA = [AA;AA] # concaténation verticale par bloc (vcat)
C = [ 1 2 3; 4 5 6; 7 8 9]
C = [ u * v for u = 0:0.1:1, v = 0:0.1:1]

# =============================================================================
# Programmation
# =============================================================================

for i = 1:5
    println("i = $i")
end
for i = 1:1000
    println(i)
    if i >= 5
       break
    end
end
i = 1
if i == 1
    println("i = 1")
end
i = 2
if i != 1
    println("i ≠ 1")
else
    println("i = 1")
end
i = 0
if i == 1
    println("i = 1")
elseif i > 1
    println("i > 1")
else
    println("i < 1")
end
i=0
while i<5
    global i+=1  # i =i+1
    println("i = $i")
end
i=0
while i<10
    global i+=1
    if i <= 5
    println("i = $i")
    else
        break
    end
end

# =============================================================================
# Algèbre linéaire
# Déterminant, valeurs et vecteurs propres...
# =============================================================================

using LinearAlgebra

A = [1. 2. 3.; 4. 5. 6.; 7. 8. 10.]
diag(A)
tr(A)
sum(A)
det(A)
F = eigen(A)
F.values
F.vectors
b = [3, 2, 1] # a vector
x = A \ b # résolution du système linéaire Ax=b
xx = inv(A) * b # résolution du système linéaire Ax=b
x - xx

# =============================================================================
# Fonctions, Méthodes, Compilation, Spécialisation, Multiple Dispatch
# =============================================================================

# Version longue

function fun_name(arg1, arg2)
    tmp = arg1 + arg2
    return tmp * arg1          # <- le mot-clé return est optionnel, la dernière valeur calculée est retournée
end

# Appel

fun_name(2, 3)

# Version courte

fun_name(arg1, arg2) = arg1 * (arg1 + arg2)

# Appel

fun_name(2, 3)

# Multiple Dispatch
# différentes implémentations (méthodes) de la fonction f

function f end
f(x)          = println("Generic method with 1 argument : $x")
f(x::Integer) = println("Specific method with 1 argument (Integer) : $x")
f(x, y)       = println("Generic method with 2 arguments: $x, $y")
f(x::Real, y) = println("Specific method with 2 arguments (Real, Any): $x, $y")

methods(f)

f(π)
f(1)
f(im, 2)
f(π, [1, 2, 3])

@which f(im, im)

using BenchmarkTools

function my_sum(X)
    acc = zero(eltype(X)) # zero du type de l'élément de l'argument
    for x in X
        acc += x
    end
    acc
end

V = rand(1_000_000)
@time my_sum(V) # Compilation + Exécution de la fonction  my_sum avec argument de type vecteur (spécialisation)
@time my_sum(V) # Exécution de la fonction  my_sum
@time my_sum(V)

@code_native debuginfo = :none my_sum(V)

UR = (1:1_000_000)
@time my_sum(UR) # # Compilation + Exécution de la fonction  my_sum avec argument de type UnitRange (spécialisation)
@time my_sum(UR)
@time my_sum(UR)

@code_native debuginfo = :none my_sum(UR) # code natif différent du précédent  avec même implémnetation => spécialisation

@btime my_sum(UR)
@btime Base.sum(UR) # pas de magie dans Julia !!

function my_sum(X::Array{Float64,2})
    acc = zero(eltype(X))
    ni, nj = size(X)
    for j in 1:nj
        for i in 1:ni
            acc += X[i,j]
        end
    end
    acc
end
A = rand(Float64, 1_000, 1_000)
methods(my_sum)
@time my_sum(A) # type est connu, donc pas de perte de temps pour le branchement
@time my_sum(A)
@time my_sum(A)
@btime my_sum(A)
@which my_sum(A)

@code_native debuginfo = :none my_sum(A)

# =============================================================================
# Graphique avec Plots
#  =============================================================================

using Plots
x=0:2π/100:2π
y=sin.(x)
theme(:dark)
plot(x,y,title="Mon premier graphique avec Plots")
y=cos.(x)
plot!(x,y)

# =============================================================================
# Sciences des données
# https://en.wikibooks.org/wiki/Introducing_Julia/DataFrames
# =============================================================================

# using Pkg; Pkg.add("CSV"); Pkg.add("HTTP");Pkg.add("DataFrames")
using CSV
using HTTP
using DataFrames

df = CSV.File(HTTP.get("https://raw.githubusercontent.com/nassarhuda/easy_data/master/programming_languages.csv").body) |> DataFrame
@show df

# =============================================================================
# Calcul parallèle sur GPU NVidia
# =============================================================================

using BenchmarkTools

mcpu = rand(2^12, 2^12)
@btime mcpu*mcpu

# using Pkg; Pkg.add("CUDA")
using CUDA

mgpu = cu(mcpu)
@btime  mgpu*mgpu

# =============================================================================
# Calcul symbolique
# =============================================================================

# using Pkg; Pkg.add("SymPy")
using SymPy

E, S, L, θ = symbols("E S L θ",real=true)
ke = E * S / L * Sym[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0 ]
Re=Sym[cos(θ) sin(θ) 0 0; -sin(θ) cos(θ) 0 0; 0 0 cos(θ) sin(θ); 0 0 -sin(θ) cos(θ)]
Ke=Re'*ke*Re

# using Pkg; Pkg.add("Symbolics");Pkg.add("SymbolicUtils")
using Symbolics

@variables Eₑ Sₑ Lₑ θₑ
kₑ = Eₑ * Sₑ / Lₑ * [1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0 ]
Rₑ = [cos(θₑ) sin(θₑ) 0 0; -sin(θₑ) cos(θₑ) 0 0; 0 0 cos(θₑ) sin(θₑ); 0 0 -sin(θₑ) cos(θₑ)]
Kₑ = Rₑ' * kₑ * Rₑ


text:
"""
# Hooke's atom

Find the total energy and density of the
[Hooke's atom](https://en.wikipedia.org/wiki/Hooke's_atom) using LDA DFT.
"""

using AtomEnergyLevels
using SpecialFunctions, PyPlot
pygui(true);

x = -30.0:0.1:5;
r  = exp.(x);

# exact density of the ground state, see
# S. Kais et al // Density functionals and dimensional
# renormalization for an exactly solvable model, JCP 99, 417 (1993);
# http://dx.doi.org/10.1063/1.465765

N² = 1 / (π^(3/2) * (8 + 5 * sqrt(π)));
ρₑ = @. 2N² * exp(-1/2 * r^2) * (sqrt(π/2) * (7/4 + 1/4 * r^2 + (r + 1/r) * erf(r/sqrt(2))) + exp(-1/2 * r^2));

xalpha = lda(2, x, conf = c"[He]", Vex = r -> 1/8 * r^2,
             xc = xc_func(xc_type_lda, false, ρ -> LDA_X(ρ, α = 0.798), ρ -> 0));
svwn5  = lda(2, x, conf = c"[He]", Vex = r -> 1/8 * r^2);

begin
    title("Hooke's atom density")
    xlabel("r, Bohr")
    ylabel("ρ(r)")

    ax1 = PyPlot.axes()
    ax1.set_xlim([0,5])

    plot(r, ρₑ, label = "Exact")
    plot(r, xalpha.density ./ r, label = "Xα, α=0.798")
    plot(r, svwn5.density ./ r, label = "LDA")

    legend(loc = "upper right", fancybox = "true")
    grid("on")
end


text:
using MinFEM
using Test

tests = ["poisson_test", "elasticity_test", "semilinear_test", "gradient_test"]
for t in tests
  include("$(t).jl")
end


text:
function pequal(p1::DynamicParticle{T}, p2::DynamicParticle{T}) where {T<:Any}
    return all([
        p1.x == p2.x,
        p1.ℓw == p2.ℓw,
        p1.ℓdens === p2.ℓdens,
        p1.cacheℓdens == p2.cacheℓdens,
        p1.τ == p2.τ])
end

function pequal(p1::StaticParticle{T}, p2::StaticParticle{T}) where {T<:Any}
    return all([
        p1.x == p2.x,
        p1.ℓw == p2.ℓw])
end


text:
@testset "runtestString" begin
    println("Testing : runtestString")

    include("../src/bfgaHelloWorld.jl")
    using .bfgaHelloWorld

    include("../src/Bf.jl")
    using .BfInterpreter

    goal = bfgaHelloWorld.getTarget()
    target_length = length(goal)

    @time model = test_serial()
    prems = model.refSet[length(model.refSet)]
    println("Generation : $(model.params.currentGeneration)")
    println("& fitness : $(model.specific_fitness.fitness(prems, model.instructionsSet))")
    #bfcode :: Array{Float64,1}  = prems.dna
    println("code : $(prems.program) ") #$(join(genesToBfInstructions(bfcode), "") ) ")

    #bft = bfType(bfcode)

    output, _ = execute(prems.program, model.instructionsSet)
    #output2 = execute(bft, input2)

    @test length(output) >= target_length
    @test collect(goal) == output[1:target_length]

    # do the test_parallel after ...
    #test_parallel()
end


text:
abstract type IRMethod <: MethodType end

####################################################################
# BM25
####################################################################

mutable struct BM25 <: IRMethod
    k1::Float64
    b::Float64
end

computeDocFreq(method::BM25, docFreq::DocFreq) = @. log(1.0 + (docFreq.nDocs - docFreq.doc_freq + 0.5) / (docFreq.doc_freq + 0.5))
compute_tf(method::BM25, tf, doc_len, doc_avg) = (tf * (k1 + 1.0)) / (tf + k1 * (1.0 - b + b * doc_len / doc_avg))
compute_term_sim(method::BM25, term_idf, q_weight, c_weight, doc_len, doc_avg) = compute_tf(method::BM25, c_weight, doc_len, doc_avg) * term_idf
generate_bow_vector(method::BM25, report::Report, df_vec::Vector{Float64}, vocab::Vocab, ngram::Int64) = transform_BOW(report, vocab::Vocab, ngram::Int64)

normalize_vec(method::BM25, bow::Vector{Pair{UInt32,Float64}}, df_vec::Vector{Float64})= bow

compute_final_similarity!(partial_similarities::Vector{Float64}, method::BM25, query::Report, candidates::Vector{Report}, df_vec::Vector{Float64}) = partial_similarities

####################################################################
# TF-IDF
####################################################################

@enum tf_scheme TF_RAW = 1  TF_SQRT = 2 TF_LOG = 3
@enum idf_scheme IDF_STD = 1 IDF_SMOOTH = 2 IDF_SMOOTH_PLUS = 3

function get_tf_scheme(tf_scheme_nm::String)
    if tf_scheme_nm == "raw"
        return TF_RAW
    elseif tf_scheme_nm == "log"
        return TF_LOG
    elseif tf_scheme_nm == "sqrt"
        return TF_SQRT
    end

    return TF_RAW
end

function get_idf_scheme(idf_scheme_nm::String)
    if idf_scheme_nm == "std"
        return IDF_STD
    elseif idf_scheme_nm == "smooth"
        return IDF_SMOOTH
    elseif idf_scheme_nm == "smooth_plus"
        return IDF_SMOOTH_PLUS
    end

    return IDF_SMOOTH
end

mutable struct TFIDF <: IRMethod
    tf_scheme::tf_scheme
    idf_scheme::idf_scheme
    len_norm::Vector{Float64}
    without_len_norm::Bool
    add_query_freq::Bool
end

function TFIDF(tf_scheme_nm::String, idf_scheme_nm::String, without_len_norm=false, add_query_freq=false)
    tf_scheme = get_tf_scheme(tf_scheme_nm)
    idf_scheme = get_idf_scheme(idf_scheme_nm)

    @info("TF: $(tf_scheme); IDF: $(idf_scheme)")

    if without_len_norm
        len_norm = [1.0]
    else
        len_norm = Vector{Float64}(undef, 12000)

        for doc_len in 1:12000
            len_norm[doc_len] = 1.0 / sqrt(doc_len)
        end

    end

    return TFIDF(tf_scheme, idf_scheme, len_norm, without_len_norm, add_query_freq)
end

tf_sqrt(tf) = sqrt(tf)
tf_log(tf) = log(1 + tf)

idf_std(docFreq::DocFreq) = @. log(docFreq.nDocs / docFreq.doc_freq)
idf_smooth(docFreq::DocFreq) = @. log((docFreq.nDocs + 1.0) / (docFreq.doc_freq + 1.0))
idf_smooth_plus(docFreq::DocFreq) = idf_smooth(docFreq) .+ 1.0

@inline function computeDocFreq(method::TFIDF, docFreq::DocFreq)
    if method.idf_scheme == IDF_SMOOTH
        return idf_smooth(docFreq)
    elseif method.idf_scheme == IDF_STD
        return idf_std(docFreq)
    elseif method.idf_scheme == IDF_SMOOTH_PLUS
        return idf_smooth_plus(docFreq)
    else
        return idf_smooth(docFreq)
    end
end

@inline function compute_tf(method::TFIDF, raw_tf, doc_len, doc_avg)
    len_norm = method.without_len_norm ? method.len_norm[1] : method.len_norm[doc_len]

    if method.tf_scheme == TF_SQRT
        return len_norm * tf_sqrt(raw_tf)
    elseif method.tf_scheme == TF_LOG
        return len_norm * tf_log(raw_tf)
    else
        return len_norm * raw_tf
    end
end

function similarity(method::TFIDF, q_bow, c_bow, idf_vec, cand_len)
    q_idx = 1
    c_idx = 1
    sim = 0.0

    q_len = length(q_bow) + 1
    c_len = length(c_bow) + 1

    @inbounds while q_idx != q_len && c_idx != c_len
        q_func_id, q_tf = q_bow[q_idx]
        c_func_id, c_tf = c_bow[c_idx]

        if q_func_id == c_func_id
            # Query and candidate contain the same function. Align positions starting from the left
            term_idf = idf_vec[q_func_id]
            sim += compute_term_sim(method, term_idf, q_tf, c_tf, cand_len, -0.0)

            q_idx +=1
            c_idx +=1
        elseif q_func_id > c_func_id
            # c_func_id does not exist in the query
            c_idx += 1
        else
            # q_func_id does not exist in the candidate
            q_idx += 1
        end
    end

    return sim
end

@inline function compute_term_sim(method::TFIDF, term_idf, q_weight, c_weight, doc_len, doc_avg)
    sim = compute_tf(method, c_weight, doc_len, doc_avg) * term_idf
    # println("==> sim=$(sim); tf=$(compute_tf(method, c_weight, doc_len, doc_avg)); idf=$(term_idf)")

    if method.add_query_freq
        sim *= compute_tf(method, q_weight, 1, doc_avg)
    end

    return sim
end

function generate_bow_vector(method::TFIDF, report::Report, df_vec::Vector{Float64}, vocab::Vocab, ngram::Int64)
    return transform_BOW(report, vocab::Vocab, ngram::Int64)
end
normalize_vec(method::TFIDF, bow::Vector{Pair{UInt32,Float64}}, df_vec::Vector{Float64})= bow
compute_final_similarity!(partial_similarities::Vector{Float64}, method::TFIDF, query::Report, candidates::Vector{Report}, df_vec::Vector{Float64}) = partial_similarities

####################################################################
# Cosine Simarity + TFIDF
####################################################################
using LinearAlgebra

mutable struct CosineTFIDF <: IRMethod
    tfidf::TFIDF
    precompute_unit_vec::Bool
end

function CosineTFIDF(tf_scheme_nm::String, idf_scheme_nm::String, precompute_unit_vec=false)
    tf_idf = TFIDF(tf_scheme_nm, idf_scheme_nm, true, false)

    return CosineTFIDF(tf_idf, precompute_unit_vec)
end

computeDocFreq(method::CosineTFIDF, docFreq::DocFreq) = computeDocFreq(method.tfidf, docFreq)

function similarity(method::CosineTFIDF, q_bow, c_bow, idf_vec, cand_len)
    q_idx = 1
    c_idx = 1
    sim = 0.0

    q_len = length(q_bow) + 1
    c_len = length(c_bow) + 1
    squared_sum_q_weight = 0.0
    squared_sum_c_weight = 0.0

    @inbounds while q_idx != q_len && c_idx != c_len
        q_func_id, q_tf = q_bow[q_idx]
        c_func_id, c_tf = c_bow[c_idx]

        if q_func_id == c_func_id
            # Query and candidate contain the same function. Align positions starting from the left
            q_weight = compute_tf(method.tfidf, q_tf, 0, 0.0) * idf_vec[q_func_id]
            c_weight = compute_tf(method.tfidf, c_tf, 0, 0.0) * idf_vec[c_func_id]

            sim+= q_weight * c_weight

            squared_sum_q_weight += q_weight^2
            squared_sum_c_weight += c_weight^2

            q_idx +=1
            c_idx +=1
        elseif q_func_id > c_func_id
            # c_func_id does not exist in the query
            c_weight = compute_tf(method.tfidf, c_tf, 0, 0.0) * idf_vec[c_func_id]
            squared_sum_c_weight += c_weight^2

            c_idx += 1
        else
            # q_func_id does not exist in the candidate
            q_weight = compute_tf(method.tfidf, q_tf, 0, 0.0) * idf_vec[q_func_id]
            squared_sum_q_weight += q_weight^2

            q_idx += 1
        end
    end

    @inbounds for j=q_idx:length(q_bow)
        func_id, tf = q_bow[j]
        weight = compute_tf(method.tfidf, tf, 0, 0.0) * idf_vec[func_id]
        squared_sum_q_weight += weight^2
    end

    @inbounds for j=c_idx:length(c_bow)
        func_id, tf = c_bow[j]
        weight = compute_tf(method.tfidf, tf, 0, 0.0) * idf_vec[func_id]
        squared_sum_c_weight += weight^2
    end

    return sim/(sqrt(squared_sum_q_weight) * sqrt(squared_sum_c_weight))
end

@inline function compute_term_sim(method::CosineTFIDF, term_idf, q_weight, c_weight, doc_len, doc_avg)
    if method.precompute_unit_vec
        return q_weight * c_weight
    else
        qtf = compute_tf(method.tfidf, q_weight, doc_len, doc_avg)
        ctf = compute_tf(method.tfidf, c_weight, doc_len, doc_avg)

        # println("===> $(qtf * term_idf^2 * ctf) $(qtf) $(term_idf) $(ctf)")
        return qtf * term_idf^2 * ctf
    end
end

function normalize_vec(method::CosineTFIDF, bow::Vector{Pair{UInt32,Float64}}, df_vec::Vector{Float64})
    if !method.precompute_unit_vec
        return bow
    end

    # Compute unit vector
    weights = Vector{Float64}(undef, length(bow))

    # println("Report: $(report.id)")
    @inbounds for idx in eachindex(weights)
        # println("idf= $(df_vec[bow[idx].first]); tf=$(compute_tf(method.tfidf, bow[idx].second, 0, 0.0))")
        df = df_vec[bow[idx].first] > 0.0 ? df_vec[bow[idx].first] : 1.0e-100
        weights[idx] = df * compute_tf(method.tfidf, bow[idx].second, 0, 0.0)
    end

    magnitude = norm(weights)
    # println("magnitude=$(magnitude)")
    bow_norm = Vector{Pair{UInt32,Float64}}(undef, length(bow))

    for idx in eachindex(weights)
        @inbounds bow_norm[idx] = Pair{UInt32,Float64}(bow[idx].first, weights[idx] / magnitude)
    end
    # println(bow_norm)

    return bow_norm
end

function magnitude(method::CosineTFIDF, report::Report, df_vec::Vector{Float64})
    magnitude = 0.0
    bow = report.bow

    # println("# $(report.id)")
    @inbounds for idx in eachindex(bow)
        df = df_vec[bow[idx].first] > 0.0 ? df_vec[bow[idx].first] : 1.0e-100
        magnitude+= (df * compute_tf(method.tfidf, bow[idx].second, 0, 0.0))^2
    end

    return sqrt(magnitude)
end

@inline function compute_final_similarity!(partial_similarities::Vector{Float64}, method::CosineTFIDF, query::Report, candidates::Vector{Report}, df_vec::Vector{Float64})
    # Compute magnitude of query
    if method.precompute_unit_vec
        return partial_similarities
    end

    q_magnitude = magnitude(method, query, df_vec)

    Threads.@threads for idx in eachindex(partial_similarities)
        @inbounds if partial_similarities[idx] != 0.0
            c_magnitude = magnitude(method, candidates[idx], df_vec)
            partial_similarities[idx] /=  c_magnitude * q_magnitude
        end
    end

    return partial_similarities
end

max_stars_repo_name: irving-muller/FaST

max_issues_repo_name: irving-muller/FaST

max_forks_repo_name: irving-muller/FaST


text:
# multiple dispatch

# ----------------------------------------------------------------------
# setup

mutable struct Position
    x::Int
    y::Int
end

struct Size
    width::Int
    height::Int
end

# A thing should have a position and size
abstract type Thing end
position(t::Thing) = t.position
size(t::Thing) = t.size
shape(t::Thing) = :unknown

# Type of weapons
@enum Weapon Laser Missile

# Spaceship
struct Spaceship <: Thing
    position::Position
    size::Size
    weapon::Weapon
end
shape(s::Spaceship) = :saucer

# Asteroid
struct Asteroid <: Thing
    position::Position
    size::Size
end

# ----------------------------------------------------------------------
# dispatch to the narrowest type

s1 = Spaceship(Position(0,0), Size(30,5), Missile);
s2 = Spaceship(Position(10,0), Size(30,5), Laser);
a1 = Asteroid(Position(20,0), Size(20,20));
a2 = Asteroid(Position(0,20), Size(20,20));

#=
julia> s1 = Spaceship(Position(0,0), Size(30,5), Missile);

julia> s2 = Spaceship(Position(10,0), Size(30,5), Laser);

julia> a1 = Asteroid(Position(20,0), Size(20,20));

julia> a2 = Asteroid(Position(0,20), Size(20,20));

julia> position(s1), size(s1), shape(s1)
(Position(0, 0), Size(30, 5), :saucer)

julia> position(a1), size(a1), shape(a1)
(Position(20, 0), Size(20, 20), :unknown)
=#

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

struct Rectangle
    top::Int
    left::Int
    bottom::Int
    right::Int
    # return two upper-left and lower-right points of the rectangle
    Rectangle(p::Position, s::Size) =
        new(p.y+s.height, p.x, p.y, p.x+s.width)
end

# check if the two rectangles (A & B) overlap
function overlap(A::Rectangle, B::Rectangle)
    return A.left < B.right && A.right > B.left &&
        A.top > B.bottom && A.bottom < B.top
end

# using abstract type Thing
function collide(A::Thing, B::Thing)
    println("Checking collision of thing vs. thing")
    rectA = Rectangle(position(A), size(A))
    rectB = Rectangle(position(B), size(B))
    return overlap(rectA, rectB)
end

#=
julia> collide(s1, s2);
Checking collision of thing vs. thing

julia> collide(a1, a2);
Checking collision of thing vs. thing

julia> collide(s1, a1);
Checking collision of thing vs. thing

julia> collide(a1, s1);
Checking collision of thing vs. thing
=#

# more specific
function collide(A::Spaceship, B::Spaceship)
    println("Checking collision of spaceship vs. spaceship")
    return true   # just a test
end

#=
julia> collide(s1, s2);
Checking collision of spaceship vs. spaceship

julia> collide(a1, a2);
Checking collision of thing vs. thing

julia> collide(s1, a1);
Checking collision of thing vs. thing

julia> collide(a1, s1);
Checking collision of thing vs. thing
=#

# ----------------------------------------------------------------------
# Ambiguity

function collide(A::Asteroid, B::Thing)
    println("Checking collision of asteroid vs. thing")
    return true
end

function collide(A::Thing, B::Asteroid)
    println("Checking collision of thing vs. asteroid")
    return false
end

collide(a1, s1);
collide(s1, a1);
collide(a1, a2);

#= REPL
julia> collide(a1, s1);
Checking collision of asteroid vs. thing

julia> collide(s1, a1);
Checking collision of thing vs. asteroid

julia> collide(a1, a2);
ERROR: MethodError: collide(::Asteroid, ::Asteroid) is ambiguous. Candidates:
  collide(A::Thing, B::Asteroid) in Main at REPL[23]:2
  collide(A::Asteroid, B::Thing) in Main at REPL[22]:5
Possible fix, define
  collide(::Asteroid, ::Asteroid)
=#

function collide(A::Asteroid, B::Asteroid)
    println("Checking collision of asteroid vs. asteroid")
    return true  # just a test
end
collide(a1, a2);

#= REPL
julia> collide(a1, a2);
Checking collision of asteroid vs. asteroid
=#

# ----------------------------------------------------------------------
# detecting ambiguity

using Test

module Foo
    foo(x, y) = 1
    foo(x::Integer, y) = 2
    foo(x, y::Integer) = 3
end

detect_ambiguities(Main.Foo)

#= REPL
julia> using Test

julia> module Foo
            foo(x, y) = 1
            foo(x::Integer, y) = 2
            foo(x, y::Integer) = 3
       end
Main.Foo

julia> detect_ambiguities(Main.Foo)
1-element Array{Tuple{Method,Method},1}:
 (foo(x::Integer, y) in Main.Foo at REPL[25]:3, foo(x, y::Integer) in Main.Foo at REPL[25]:4)
=#

# -- fixed
module Foo
    foo(x, y) = 1
    foo(x::Integer, y) = 2
    foo(x, y::Integer) = 3
    foo(x::Integer, y::Integer) = 4
end

detect_ambiguities(Main.Foo)

#=
julia> module Foo
            foo(x, y) = 1
            foo(x::Integer, y) = 2
            foo(x, y::Integer) = 3
            foo(x::Integer, y::Integer) = 4
        end
Main.Foo

julia> detect_ambiguities(Main.Foo)
0-element Array{Tuple{Method,Method},1}
=#

# -- multiple modules
module Foo2
    foo(x, y) = 1
    foo(x::Integer, y) = 2
    foo(x, y::Integer) = 3
end

module Foo4
    import Main.Foo2
    Foo2.foo(x::Integer, y::Integer) = 4
end

detect_ambiguities(Main.Foo2, Main.Foo4)

#=
julia> module Foo2
       foo(x, y) = 1
       foo(x::Integer, y) = 2
       foo(x, y::Integer) = 3
       end
Main.Foo2

julia> module Foo4
       import Main.Foo2
       Foo2.foo(x::Integer, y::Integer) = 4
       end
Main.Foo4

julia> detect_ambiguities(Main.Foo2, Main.Foo4)
0-element Array{Tuple{Method,Method},1}
=#

# ----------------------------------------------------------------------
# dynamic dispatch

# randomly pick two things and check
function check_randomly(things)
    for i in 1:5
        two = rand(things, 2)
        collide(two...)
    end
end

check_randomly([s1, s2, a1, a2])

#= REPL
julia> check_randomly([s1, s2, a1, a2])
Checking collision of spaceship vs. spaceship
Checking collision of thing vs. asteroid
Checking collision of asteroid vs. thing
Checking collision of asteroid vs. thing
Checking collision of asteroid vs. asteroid
=#


text:
using Gadfly, Distributions

plot(x=rand(Poisson(20), 1000), Scale.x_discrete, Geom.histogram)


text:
# --------------------------------------------------------------------------- #
# Circumvent current limitations of Knet for cat(3, ...) with 4-dim arrays

# Copied from issue #198 of Knet.jl
function setidx!(a::KnetArray, v, I...)
    crange = CartesianRange(to_indices(a, I))
    linind = [sub2ind(size(a), t.I...) for t in crange]
    setindex!(a, v, vec(linind))
end

# This can only be done efficiently with knet-Arrays in case of
# batchsize = 1, when permutedims is a no-op
function cat3(a, b)

  if size(a, 4) == 1
    k = size(a)[1:2]
    sa, sb = size(a, 3), size(b, 3)

    ar, br = reshape(a, prod(k), sa), reshape(b, prod(k), sb)
    cr = reshape(hcat(ar, br), k..., sa + sb, 1)

    return cr

  else
    ap, bp = permutedims(a, (1,2,4,3)), permutedims(b, (1,2,4,3))
    sa, sb = size(ap, 4), size(bp, 4)

    k = size(ap)[1:3]
    ap, bp = reshape(ap, prod(k), sa), reshape(bp, prod(k), sb)

    cp = hcat(ap, bp)
    cp = reshape(cp, k..., sa + sb)

    return permutedims(cp, (1, 2, 4, 3))
  end

end

cat3(a :: Array{Float32, 4}, b :: Array{Float32, 4}) = cat(3, a, b)

# --------------------------------------------------------------------------- #
# Primitives for multires networks

# Padding for convolution such that dimensions are preserved
pad(w) = floor(Int, size(w, 1) / 2)

# Convolution and relu activation
rconv(w, x) = relu.(conv4(w, x, padding=pad(w)))

# Convolution, batch normalization, and relu activation
rbconv(wb, m, wc, x) = relu.(batchnorm(conv4(wc, x, padding=pad(wc)), m, wb))

# Up-convolution
uconv(w, x) = deconv4(w, x, stride=2, padding=1)

# Up-convolution and batch normalization
ubconv(wb, m, wc, x) = batchnorm(deconv4(wc, x, stride=2, padding=1), m, wb)

# Up-convolution and stacking
uconv(w, x, y) = cat3(x, deconv4(w, y, stride=2, padding=1))

# Up-convolution, batch normalization, and stacking
ubconv(wb, m, wc, x, y) =
    cat3(x, batchnorm(deconv4(wc, y, stride=2, padding=1), m, wb))

# Default initialization for convolution kernel
wr(i, o, size = (3, 3)) = gaussian(Float32, size..., i, o, std=0.060)

# Default initialization for up-convolution kernel
wu(i, o) = bilinear(Float32, 2, 2, o, i)

# Knet parameters for batch normalization
bn(c) = bnparams(Float32, c)

# Check if a number is odd
odd(x::Integer) = (x % 2 != 0)

max_stars_repo_name: tscode/DCellC

max_issues_repo_name: tscode/DCellC