text:
using FillArrays, LinearAlgebra, SparseArrays, Random, Test

import FillArrays: AbstractFill

@testset "fill array constructors and convert" begin
    for (Typ, funcs) in ((:Zeros, :zeros), (:Ones, :ones))
        @eval begin
            @test_throws BoundsError $Typ((-1,5))
            @test $Typ(5) isa AbstractVector{Float64}
            @test $Typ(5,5) isa AbstractMatrix{Float64}
            @test $Typ(5) == $Typ((5,))
            @test $Typ(5,5) == $Typ((5,5))
            @test eltype($Typ(5,5)) == Float64

            for T in (Int, Float64)
                Z = $Typ{T}(5)
                @test eltype(Z) == T
                @test Array(Z) == $funcs(T,5)
                @test Array{T}(Z) == $funcs(T,5)
                @test Array{T,1}(Z) == $funcs(T,5)

                @test convert(AbstractArray,Z) ≡ Z
                @test convert(AbstractArray{T},Z) ≡ AbstractArray{T}(Z) ≡ Z
                @test convert(AbstractVector{T},Z) ≡ AbstractVector{T}(Z) ≡ Z

                @test $Typ{T,1}(2ones(T,5)) == Z
                @test $Typ{T}(2ones(T,5)) == Z
                @test $Typ(2ones(T,5)) == Z

                Z = $Typ{T}(5, 5)
                @test eltype(Z) == T
                @test Array(Z) == $funcs(T,5,5)
                @test Array{T}(Z) == $funcs(T,5,5)
                @test Array{T,2}(Z) == $funcs(T,5,5)

                @test convert(AbstractArray,Z) ≡ convert(AbstractFill,Z) ≡ Z
                @test convert(AbstractArray{T},Z) ≡ convert(AbstractFill{T},Z) ≡ AbstractArray{T}(Z) ≡ Z
                @test convert(AbstractMatrix{T},Z) ≡ convert(AbstractFill{T,2},Z) ≡ AbstractMatrix{T}(Z) ≡ Z

                @test $Typ{T,2}(2ones(T,5,5)) == Z
                @test $Typ{T}(2ones(T,5,5)) == Z
                @test $Typ(2ones(T,5,5)) == Z

                @test AbstractArray{Float32}(Z) == $Typ{Float32}(5,5)
                @test AbstractArray{Float32,2}(Z) == $Typ{Float32}(5,5)
            end
        end
    end

    @test_throws BoundsError Fill(1,(-1,5))
    @test Fill(1.0,5) isa AbstractVector{Float64}
    @test Fill(1.0,5,5) isa AbstractMatrix{Float64}
    @test Fill(1,5) == Fill(1,(5,))
    @test Fill(1,5,5) == Fill(1,(5,5))
    @test eltype(Fill(1.0,5,5)) == Float64

    @test Matrix{Float64}(Zeros{ComplexF64}(10,10)) == zeros(10,10)
    @test_throws InexactError Matrix{Float64}(Fill(1.0+1.0im,10,10))

    for T in (Int, Float64)
        F = Fill{T}(one(T), 5)

        @test eltype(F) == T
        @test Array(F) == fill(one(T),5)
        @test Array{T}(F) == fill(one(T),5)
        @test Array{T,1}(F) == fill(one(T),5)

        F = Fill{T}(one(T), 5, 5)
        @test eltype(F) == T
        @test Array(F) == fill(one(T),5,5)
        @test Array{T}(F) == fill(one(T),5,5)
        @test Array{T,2}(F) == fill(one(T),5,5)

        @test convert(AbstractArray,F) ≡ F
        @test convert(AbstractArray{T},F) ≡ AbstractArray{T}(F) ≡ F
        @test convert(AbstractMatrix{T},F) ≡ AbstractMatrix{T}(F) ≡ F

        @test convert(AbstractArray{Float32},F) == AbstractArray{Float32}(F) ==
                Fill{Float32}(one(Float32),5,5)
        @test convert(AbstractMatrix{Float32},F) == AbstractMatrix{Float32}(F) ==
                Fill{Float32}(one(Float32),5,5)
    end

    @test_throws BoundsError Eye((-1,5))
    @test Eye(5,5) isa AbstractMatrix{Float64}
    @test Eye(5,5) == Eye(5) == Eye((5,5)) == Eye{Float64}(5) == Eye{Float64}(5, 5)
    @test eltype(Eye(5,5)) == Float64

    for T in (Int, Float64)
        E = Eye{T}(5, 5)
        M = Matrix{T}(I, 5, 5)

        @test eltype(E) == T
        @test Array(E) == M
        @test Array{T}(E) == M
        @test Array{T,2}(E) == M

        @test convert(AbstractArray,E) === E
        @test convert(AbstractArray{T},E) === E
        @test convert(AbstractMatrix{T},E) === E

        @test AbstractArray{Float32}(E) == Eye{Float32}(5,5)
        @test AbstractArray{Float32,2}(E) == Eye{Float32}(5,5)

        @test Eye{T}(ones(T, 5, 5)) == E
        @test Eye(ones(T, 5, 5)) == E
    end

    @testset "Bool should change type" begin
        x = Fill(true,5)
        y = x + x
        @test y isa Fill{Int,1}
        @test y[1] == 2

        x = Ones{Bool}(5)
        y = x + x
        @test y isa Fill{Int,1}
        @test y[1] == 2
        @test x + Zeros{Bool}(5) ≡ x
        @test x - Zeros{Bool}(5) ≡ x
        @test Zeros{Bool}(5) + x ≡ x
        @test -x ≡ Fill(-1,5)
    end
end

# Check that all pair-wise combinations of + / - elements of As and Bs yield the correct
# type, and produce numerically correct results.
function test_addition_and_subtraction(As, Bs, Tout::Type)
    for A in As, B in Bs
        @test A + B isa Tout{promote_type(eltype(A), eltype(B))}
        @test Array(A + B) == Array(A) + Array(B)

        @test A - B isa Tout{promote_type(eltype(A), eltype(B))}
        @test Array(A - B) == Array(A) - Array(B)

        @test B + A isa Tout{promote_type(eltype(B), eltype(A))}
        @test Array(B + A) == Array(B) + Array(A)

        @test B - A isa Tout{promote_type(eltype(B), eltype(A))}
        @test Array(B - A) == Array(B) - Array(A)
    end
end

# Check that all permutations of + / - throw a `DimensionMismatch` exception.
function test_addition_and_subtraction_dim_mismatch(a, b)
    @test_throws DimensionMismatch a + b
    @test_throws DimensionMismatch a - b
    @test_throws DimensionMismatch b + a
    @test_throws DimensionMismatch b - a
end

@testset "FillArray addition and subtraction" begin
    test_addition_and_subtraction_dim_mismatch(Zeros(5), Zeros(6))
    test_addition_and_subtraction_dim_mismatch(Zeros(5), Zeros{Int}(6))
    test_addition_and_subtraction_dim_mismatch(Zeros(5), Zeros(6,6))
    test_addition_and_subtraction_dim_mismatch(Zeros(5), Zeros{Int}(6,5))

    # Construct FillArray for repeated use.
    rng = MersenneTwister(123456)
    A_fill, B_fill = Fill(randn(rng, Float64), 5), Fill(4, 5)

    # Unary +/- constructs a new FillArray.
    @test +A_fill === A_fill
    @test -A_fill === Fill(-A_fill.value, 5)

    # FillArray +/- FillArray should construct a new FillArray.
    test_addition_and_subtraction([A_fill, B_fill], [A_fill, B_fill], Fill)
    test_addition_and_subtraction_dim_mismatch(A_fill, Fill(randn(rng), 5, 2))

    # FillArray + Array (etc) should construct a new Array using `getindex`.
    A_dense, B_dense = randn(rng, 5), [5, 4, 3, 2, 1]
    test_addition_and_subtraction([A_fill, B_fill], [A_dense, B_dense], Array)
    test_addition_and_subtraction_dim_mismatch(A_fill, randn(rng, 5, 2))

    # FillArray + StepLenRange / UnitRange (etc) should yield an AbstractRange.
    A_ur, B_ur = 1.0:5.0, 6:10
    test_addition_and_subtraction([A_fill, B_fill], (A_ur, B_ur), AbstractRange)
    test_addition_and_subtraction_dim_mismatch(A_fill, 1.0:6.0)
    test_addition_and_subtraction_dim_mismatch(A_fill, 5:10)
end

@testset "Other matrix types" begin
    @test Diagonal(Zeros(5)) == Diagonal(zeros(5))

    @test Diagonal(Zeros(8,5)) == Diagonal(zeros(5))
    @test convert(Diagonal, Zeros(5,5)) == Diagonal(zeros(5))
    @test_throws BoundsError convert(Diagonal, Zeros(8,5))

    @test convert(Diagonal{Int}, Zeros(5,5)) == Diagonal(zeros(Int,5))
    @test_throws BoundsError convert(Diagonal{Int}, Zeros(8,5))

    @test Diagonal(Eye(8,5)) == Diagonal(ones(5))
    @test convert(Diagonal, Eye(5)) == Diagonal(ones(5))
    @test_throws BoundsError convert(Diagonal, Eye(8,5))

    @test convert(Diagonal{Int}, Eye(5)) == Diagonal(ones(Int,5))
    @test_throws BoundsError convert(Diagonal{Int}, Eye(8,5))
end

@testset "Sparse vectors and matrices" begin
    @test SparseVector(Zeros(5)) ==
            SparseVector{Float64}(Zeros(5)) ==
            SparseVector{Float64,Int}(Zeros(5)) ==
            convert(AbstractSparseArray,Zeros(5)) ==
            convert(AbstractSparseVector,Zeros(5)) ==
            convert(AbstractSparseArray{Float64},Zeros(5)) ==
            convert(AbstractSparseVector{Float64},Zeros(5)) ==
            convert(AbstractSparseVector{Float64,Int},Zeros(5)) ==
            spzeros(5)

    for (Mat, SMat) in ((Zeros(5,5), spzeros(5,5)), (Zeros(6,5), spzeros(6,5)),
                        (Eye(5), sparse(I,5,5)), (Eye(6,5), sparse(I,6,5)))
        @test SparseMatrixCSC(Mat) ==
                SparseMatrixCSC{Float64}(Mat) ==
                SparseMatrixCSC{Float64,Int}(Mat) ==
                convert(AbstractSparseArray,Mat) ==
                convert(AbstractSparseMatrix,Mat) ==
                convert(AbstractSparseArray{Float64},Mat) ==
                convert(AbstractSparseArray{Float64,Int},Mat) ==
                convert(AbstractSparseMatrix{Float64},Mat) ==
                convert(AbstractSparseMatrix{Float64,Int},Mat) ==
                SMat
    end
end

@testset "Rank" begin
    @test rank(Zeros(5,4)) == 0
    @test rank(Ones(5,4)) == 1
    @test rank(Fill(2,5,4)) == 1
    @test rank(Fill(0,5,4)) == 0
    @test rank(Eye(2,5)) == rank(Eye(5,2)) == rank(Eye(2)) == 2
end

@testset "BigInt indices" begin
    for A in (Zeros(BigInt(100)), Ones(BigInt(100)), Fill(2, BigInt(100)))
        @test length(A) isa BigInt
        @test axes(A) == tuple(Base.OneTo{BigInt}(BigInt(100)))
        @test size(A) isa Tuple{BigInt}
    end
    let A = Eye(BigInt(100))
        @test length(A) isa BigInt
        @test axes(A) == tuple(Base.OneTo{BigInt}(BigInt(100)),Base.OneTo{BigInt}(BigInt(100)))
        @test size(A) isa Tuple{BigInt,BigInt}
    end
    for A in (Zeros(BigInt(10), 10), Ones(BigInt(10), 10), Fill(2.0, (BigInt(10), 10)), Eye(BigInt(10), 10))
        @test size(A) isa Tuple{BigInt,Int}
    end

end

@testset "IndexStyle" begin
    @test IndexStyle(Zeros(5,5)) == IndexStyle(typeof(Zeros(5,5))) == IndexLinear()
end

@testset "Identities" begin
    @test Zeros(3,4) * randn(4,5) === randn(3,4) * Zeros(4,5) === Zeros(3, 5)
    @test_throws DimensionMismatch randn(3,4) * Zeros(3, 3)
    @test eltype(Zeros{Int}(3,4) * fill(1, 4, 5)) == Int
    @test eltype(Zeros{Int}(3,4) * fill(3.4, 4, 5)) == Float64
    @test Zeros(3, 4) * randn(4) == Zeros(3, 4) * Zeros(4) == Zeros(3)
    @test Zeros(3, 4) * Zeros(4, 5) === Zeros(3, 5)

    @test [1,2,3]*Zeros(1) ≡ Zeros(3)
    @test [1,2,3]*Zeros(1,3) ≡ Zeros(3,3)
    @test_throws DimensionMismatch [1,2,3]*Zeros(3)

    # Check multiplication by Adjoint vectors works as expected.
    @test randn(4, 3)' * Zeros(4) === Zeros(3)
    @test randn(4)' * Zeros(4) === zero(Float64)
    @test [1, 2, 3]' * Zeros{Int}(3) === zero(Int)
    @test_throws DimensionMismatch randn(4)' * Zeros(3)

    # Check multiplication by Transpose-d vectors works as expected.
    @test transpose(randn(4, 3)) * Zeros(4) === Zeros(3)
    @test transpose(randn(4)) * Zeros(4) === zero(Float64)
    @test transpose([1, 2, 3]) * Zeros{Int}(3) === zero(Int)
    @test_throws DimensionMismatch transpose(randn(4)) * Zeros(3)

    @test +(Zeros{Float64}(3, 5)) === Zeros{Float64}(3, 5)
    @test -(Zeros{Float32}(5, 2)) === Zeros{Float32}(5, 2)

    # `Zeros` are closed under addition and subtraction (both unary and binary).
    z1, z2 = Zeros{Float64}(4), Zeros{Int}(4)
    @test +(z1) === z1
    @test -(z1) === z1

    test_addition_and_subtraction([z1, z2], [z1, z2], Zeros)
    test_addition_and_subtraction_dim_mismatch(z1, Zeros{Float64}(4, 2))

    # `Zeros` +/- `Fill`s should yield `Fills`.
    fill1, fill2 = Fill(5.0, 4), Fill(5, 4)
    test_addition_and_subtraction([z1, z2], [fill1, fill2], Fill)
    test_addition_and_subtraction_dim_mismatch(z1, Fill(5, 5))

    X = randn(3, 5)
    for op in [+, -]

        # Addition / subtraction with same eltypes.
        @test op(Zeros(6, 4), Zeros(6, 4)) === Zeros(6, 4)
        @test_throws DimensionMismatch op(X, Zeros(4, 6))
        @test eltype(op(Zeros(3, 5), X)) == Float64

        # Different eltypes, the other way around.
        @test op(X, Zeros{Float32}(3, 5)) isa Matrix{Float64}
        @test !(op(X, Zeros{Float32}(3, 5)) === X)
        @test op(X, Zeros{Float32}(3, 5)) == X
        @test !(op(X, Zeros{ComplexF64}(3, 5)) === X)
        @test op(X, Zeros{ComplexF64}(3, 5)) == X

        # Addition / subtraction of Zeros.
        @test eltype(op(Zeros{Float64}(4, 5), Zeros{Int}(4, 5))) == Float64
        @test eltype(op(Zeros{Int}(5, 4), Zeros{Float32}(5, 4))) == Float32
        @test op(Zeros{Float64}(4, 5), Zeros{Int}(4, 5)) isa Zeros{Float64}
        @test op(Zeros{Float64}(4, 5), Zeros{Int}(4, 5)) === Zeros{Float64}(4, 5)
    end

    # Zeros +/- dense where + / - have different results.
    @test +(Zeros(3, 5), X) == X && +(X, Zeros(3, 5)) == X
    @test !(Zeros(3, 5) + X === X) && !(X + Zeros(3, 5) === X)
    @test -(Zeros(3, 5), X) == -X

    # Addition with different eltypes.
    @test +(Zeros{Float32}(3, 5), X) isa Matrix{Float64}
    @test !(+(Zeros{Float32}(3, 5), X) === X)
    @test +(Zeros{Float32}(3, 5), X) == X
    @test !(+(Zeros{ComplexF64}(3, 5), X) === X)
    @test +(Zeros{ComplexF64}(3, 5), X) == X

    # Subtraction with different eltypes.
    @test -(Zeros{Float32}(3, 5), X) isa Matrix{Float64}
    @test -(Zeros{Float32}(3, 5), X) == -X
    @test -(Zeros{ComplexF64}(3, 5), X) == -X

    # Tests for ranges.
    X = randn(5)
    @test !(Zeros(5) + X === X)
    @test Zeros{Int}(5) + (1:5) === (1:5) && (1:5) + Zeros{Int}(5) === (1:5)
    @test Zeros(5) + (1:5) === (1.0:1.0:5.0) && (1:5) + Zeros(5) === (1.0:1.0:5.0)
    @test (1:5) - Zeros{Int}(5) === (1:5)
    @test Zeros{Int}(5) - (1:5) === -1:-1:-5
    @test Zeros(5) - (1:5) === -1.0:-1.0:-5.0
end

@testset "maximum/minimum" begin
    @test maximum(Fill(1, 1_000_000_000)) == minimum(Fill(1, 1_000_000_000)) == 1
end

@testset "Cumsum" begin
    @test sum(Fill(3,10)) ≡ 30
    @test cumsum(Fill(3,10)) ≡ 3:3:30

    @test sum(Ones(10)) ≡ 10.0
    @test cumsum(Ones(10)) ≡ 1.0:10.0

    @test sum(Ones{Int}(10)) ≡ 10
    @test cumsum(Ones{Int}(10)) ≡ Base.OneTo(10)

    @test sum(Zeros(10)) ≡ 0.0
    @test cumsum(Zeros(10)) ≡ Zeros(10)

    @test sum(Zeros{Int}(10)) ≡ 0
    @test cumsum(Zeros{Int}(10)) ≡ Zeros{Int}(10)

    @test cumsum(Zeros{Bool}(10)) ≡ Zeros{Bool}(10)
    @test cumsum(Ones{Bool}(10)) ≡ Base.OneTo{Int}(10)
    @test cumsum(Fill(true,10)) ≡ 1:1:10
end

@testset "Broadcast" begin
    x = Fill(5,5)
    @test (.+)(x) ≡ x
    @test (.-)(x) ≡ -x
    @test exp.(x) ≡ Fill(exp(5),5)
    @test x .+ 1 ≡ Fill(6,5)
    @test x .+ x ≡ Fill(10,5)
    @test x .+ Ones(5) ≡ Fill(6.0,5)
    f = (x,y) -> cos(x*y)
    @test f.(x, Ones(5)) ≡ Fill(f(5,1.0),5)

    y = Ones(5,5)
    @test (.+)(y) ≡ Fill(1.0,5,5)
    @test (.-)(y) ≡ Fill(-1.0,5,5)
    @test exp.(y) ≡ Fill(exp(1),5,5)
    @test y .+ 1 ≡ Fill(2.0,5,5)
    @test y .+ y ≡ Fill(2.0,5,5)
    @test y .* y ≡ y ./ y ≡ y .\ y ≡ y

    @test Zeros{Int}(5) .+ Zeros(5) isa Zeros{Float64}
end

@testset "Sub-arrays" begin
    A = Fill(3.0,5)
    @test A[1:3] ≡ Fill(3.0,3)
    @test A[1:3,1:1] ≡ Fill(3.0,3,1)
    @test_broken A[1:3,2] ≡ Zeros{Int}(3)
    @test_throws BoundsError A[1:26]
    @test A[[true, false, true, false, false]] ≡ Fill(3.0, 2)
    A = Fill(3.0, 2, 2)
    @test A[[true true; true false]] ≡ Fill(3.0, 3)
    @test_throws DimensionMismatch A[[true, false]]

    A = Ones{Int}(5,5)
    @test A[1:3] ≡ Ones{Int}(3)
    @test A[1:3,1:2] ≡ Ones{Int}(3,2)
    @test_broken A[1:3,2] ≡ Ones{Int}(3)
    @test_throws BoundsError A[1:26]
    A = Ones{Int}(2,2)
    @test A[[true false; true false]] ≡ Ones{Int}(2)
    @test A[[true, false, true, false]] ≡ Ones{Int}(2)
    @test_throws DimensionMismatch A[[true false false; true false false]]

    A = Zeros{Int}(5,5)
    @test A[1:3] ≡ Zeros{Int}(3)
    @test A[1:3,1:2] ≡ Zeros{Int}(3,2)
    @test_broken A[1:3,2] ≡ Zeros{Int}(3)
    @test_throws BoundsError A[1:26]
    A = Zeros{Int}(2,2)
    @test A[[true false; true false]] ≡ Zeros{Int}(2)
    @test A[[true, false, true, false]] ≡ Zeros{Int}(2)
    @test_throws DimensionMismatch A[[true false false; true false false]]
end


text:
# ##############################################################################
#
# State space representation of Max-Plus linear systems.
# Creation of Max-Plus dynamical linear systems in implicit state form:
#    X(n) = D.X(n) + A.X(n-1) + B.U(n),
#    Y(n) = C.X(n)
#
# ##############################################################################

# Utility function: return an unified size for a matrix or a vector
size2(A::SparseMatrixCSC) = (size(A, 1), size(A, 2))
size2(A::Array) = (size(A, 1), size(A, 2))
size2(A::Vector) = (size(A, 1), size(A, 2))

# ==============================================================================
# Implicit dynamic linear Max-Plus system.

"""
    MPSysLin(A::AbsMatMP, B::AbsMatMP, C::AbsMatMP [, D::AbsMatMP, [x0::AbsMatMP]])

Structure for state space representation of Max-Plus linear systems.
Creation of max-plus dynamical linear systems in implicit state form:
```
    X(n) = D ⨂ X(n) ⨁ A ⨂ X(n-1) ⨁ B ⨂ U(n),
    Y(n) = C ⨂ X(n)
```

`MPSysLin` is the equivalent to ScicosLab function `mpsyslin`. It accepts dense
or sparse Max-Plus arrays. `D` and `x0` can be omited: they will be filled
with Max-Plus zeros `ε`.

# Arguments
- `A` shall be squared.
- `B` shall has its row dimension in accordance with dimensions of `A`.
- `C` shall has its column dimension in accordance with dimensions of `A`.
- `D` shall has its dimension in accordance with dimensions of `A`.
- `x0` shall has its dimensions in accordance with dimensions of `A`.

Sizes are checked and an error is thrown during the compilation.

Elements of the system can access directly ie. `S.A`, `S.B`, `S.C`, `S.D`,
`S.x0`.

# Examples
```julia-repl
julia> S = MPSysLin(MP([1.0 2 3;  4 5 6;  7 8 9]), # Max-Plus matrix A
                    MP([    0.0;      0;      0]), # Max-Plus column vector B
                    MP([    0.0       0       0]), # Max-Plus row vector C
                    mpeye(3, 3),                   # Max-Plus matrix D
                    full(zeros(MP, 3, 1)))         # Max-Plus column vector x0

Implicit dynamic linear Max-Plus system:
  x(n) = D*x(n) + A*x(n-1) + B*u(n)
  y(n) = C*x(n)
  x(0) = x0

with:
D = 3×3 Max-Plus dense matrix:
  0   .   .
  .   0   .
  .   .   0

A = 3×3 Max-Plus dense matrix:
  1   2   3
  4   5   6
  7   8   9

B = 3-element Max-Plus vector:
  0
  0
  0

C = 1×3 Max-Plus dense matrix:
  0   0   0

x0 = 3-element Max-Plus vector:
  .
  .
  .

julia> typeof(S)
MPSysLin

julia> S.A
3×3 Max-Plus dense matrix:
  1   2   3
  4   5   6
  7   8   9
```
"""
struct MPSysLin
    A::AbsMatMP
    B::AbsMatMP
    C::AbsMatMP
    D::AbsMatMP
    x0::AbsMatMP

    # Constructor with Max-Plus matrices: A, B, C, D and x0
    function MPSysLin(A::AbsMatMP, B::AbsMatMP, C::AbsMatMP, D::AbsMatMP, x0::AbsMatMP)
        # show(stdout, (size2(A),size2(B),size2(C),size2(D),size2(x0)))

        (ma,na) = size2(A)
        (ma != na) &&
            error("Matrix A shall be squared")

        (mb,nb) = size2(B)
        (mb != na) &&
            error("The row dimension of B: $(mb) is not in accordance with dimensions of A: $(na)")

        (mc,nc) = size2(C)
        (nc != na) &&
            error("The column dimension of C: $(nc) is not in accordance with dimensions of A: $(na)")

        (mx0,nx0) = size2(x0)
        (mx0 != na) &&
            error("Dimensions of x0: $(mx0) do not in accordance with dimensions of A: $(na)")

        (md,nd) = size2(D)
        ((md != na) || (nd != na)) &&
            error("The column dimension of D: $(md)×$(nd) is not in accordance with dimensions of A: $(na)")

        new(A, B, C, D, x0)
    end

    # Constructor with Max-Plus matrices: A, B, C, D (and implicit x0 set as zeros)
    function MPSysLin(A::AbsMatMP, B::AbsMatMP, C::AbsMatMP, D::AbsMatMP)
        na = size(A,2)
        MPSysLin(A, B, C, D, spzeros(MP, na, 1))
    end

    # Constructor with Max-Plus matrices: A, B, C (and implicit D and x0 set as zeros)
    function MPSysLin(A::AbsMatMP, B::AbsMatMP, C::AbsMatMP)
        na = size(A,2)
        MPSysLin(A, B, C, spzeros(MP, na, na), spzeros(MP, na, 1))
    end
end # MPSysLin

# ==============================================================================

"""
    Base.:(==)(x::MPSysLin, y::MPSysLin)

Test whether two Max-Plus linear system are equal.
```
"""
function Base.:(==)(x::MPSysLin, y::MPSysLin)
    (x.A == y.A) && (x.B == y.B) && (x.C == y.C) && (x.D == y.D) && (x.x0 == y.x0)
end

# ==============================================================================

"""
    mpshow(io::IO, S::MPSysLin)

Base function for displaying a Max-Plus linear system.

# Examples
```julia-repl
julia>
```
"""
function mpshow(io::IO, S::MPSysLin)
    (@printf io "Implicit dynamic linear Max-Plus system:\n")
    (@printf io "  x(n) = D*x(n) + A*x(n-1) + B*u(n)\n  y(n) = C*x(n)\n  x(0) = x0\n\nwith:")
    (@printf io "\nD = "); mpshow(io, S.D)
    (@printf io "\nA = "); mpshow(io, S.A)
    (@printf io "\nB = "); mpshow(io, S.B)
    (@printf io "\nC = "); mpshow(io, S.C)
    (@printf io "\nx0 = "); mpshow(io, S.x0)
end

"""
    show(io::IO, S::MPSysLin)

Display a Max-Plus linear system.

# Examples
```julia-repl
julia>
```
"""
Base.show(io::IO, S::MPSysLin) = mpshow(io, S)

"""
    show(io::IO, ::MIME"text/plain", S::MPSysLin)

Display a Max-Plus linear system.

# Examples
```julia-repl
julia>
```
"""
Base.show(io::IO, ::MIME"text/plain", S::MPSysLin) = mpshow(io, S)

# ==============================================================================

"""
    LaTeX(S::MPSysLin)

Return the LaTeX code as string from the given Max-Plus linear system.

# Examples
```julia-repl
julia>
```
"""
function LaTeX(S::MPSysLin)
    "\\left\\{\\begin{array}{lcl}\nx_n & = & " *
    LaTeX(S.D) * " x_n \\oplus " *
    LaTeX(S.A) * " x_{n-1} \\oplus " *
    LaTeX(S.B) * " u_n\\\\ y_n & = & " *
    LaTeX(S.C) * " x_n\\\\ x_0 & = & " *
    LaTeX(S.x0) * "\\end{array}\\right."
end

"""
    LaTeX(io::IO, S::MPSysLin)

Display the LaTeX code from the given Max-Plus linear system.

# Examples
```julia-repl
julia>
```
"""
function LaTeX(io::IO, S::MPSysLin)
    (@printf io "\\left\\{\\begin{array}{lcl}\n")
    (@printf io "x_n & = & ")
    LaTeX(io, S.D)
    (@printf io " x_n \\oplus ")
    LaTeX(io, S.A)
    (@printf io " x_{n-1} \\oplus ")
    LaTeX(io, S.B)
    (@printf io " u_n\\\\")
    (@printf io "y_n & = & ")
    LaTeX(io, S.C)
    (@printf io " x_n\\\\")
    (@printf io "x_0 & = & ")
    LaTeX(io, S.x0)
    (@printf io "\\end{array}\\right.")
end

Base.show(io::IO, ::MIME"text/latex", S::MPSysLin) = LaTeX(io, S)

# ==============================================================================
# Parallel composition.
# From ScicosLab file: %mpls_a_mpls.sci

"""
    Base.:(+)(y::MPSysLin, x::MPSysLin)

Parallel composition of two Max-Plus linear systems.

# Examples
```julia-repl
julia> S1 = MPSysLin(MP([1.0 2 3;  4 5 6;  7 8 9]),
                     MP([    1.0;      2;      3]),
                     MP([    4.0       5       6]),
                     mpeye(3, 3),
                     full(zeros(MP, 3, 1)));

julia> S2 = MPSysLin(MP([10.0 11 12;  13 14 15;  16 17 18]),
                     MP([    0.0;      0;      0]),
                     MP([    0.0       0       0]),
                     mpeye(3, 3),
                     full(zeros(MP, 3, 1)));

julia> S1 + S2
Implicit dynamic linear Max-Plus system:
  x(n) = D*x(n) + A*x(n-1) + B*u(n)
  y(n) = C*x(n)
  x(0) = x0

with:
D = 6×6 Max-Plus dense matrix:
  0   .   .   .   .   .
  .   0   .   .   .   .
  .   .   0   .   .   .
  .   .   .   0   .   .
  .   .   .   .   0   .
  .   .   .   .   .   0

A = 6×6 Max-Plus dense matrix:
  1   2   3    .    .    .
  4   5   6    .    .    .
  7   8   9    .    .    .
  .   .   .   10   11   12
  .   .   .   13   14   15
  .   .   .   16   17   18

B = 6-element Max-Plus vector:
  1
  2
  3
  0
  0
  0

C = 1×6 Max-Plus dense matrix:
  4   5   6   0   0   0

x0 = 6-element Max-Plus vector:
  .
  .
  .
  .
  .
  .
```
"""
function Base.:(+)(x::MPSysLin, y::MPSysLin)
    n1 = size(x.A, 1)
    n2 = size(y.A, 1)
    MPSysLin([x.A spzeros(MP, n1, n2); spzeros(MP, n2, n1) y.A],
             [x.B; y.B],
             [x.C y.C],
             [x.D spzeros(MP, n1, n2); spzeros(MP, n2, n1) y.D],
             [x.x0; y.x0])
end

# ==============================================================================
# Series composition.
# From ScicosLab file: %mpls_m_mpls.sci

"""
    Base.:(*)(y::MPSysLin, x::MPSysLin)

Series composition of two Max-Plus linear systems.

# Examples
```julia-repl
julia> S1 = MPSysLin(MP([1.0 2 3;  4 5 6;  7 8 9]),
                     MP([    1.0;      2;      3]),
                     MP([    4.0       5       6]),
                     mpeye(3, 3),
                     full(zeros(MP, 3, 1)));

julia> S2 = MPSysLin(MP([10.0 11 12;  13 14 15;  16 17 18]),
                     MP([    0.0;      0;      0]),
                     MP([    0.0       0       0]),
                     mpeye(3, 3),
                     full(zeros(MP, 3, 1)));

julia> S1 * S2
TODO
```
"""
function Base.:(*)(y::MPSysLin, x::MPSysLin)
    n1 = size(x.A, 1)
    n2 = size(y.A, 1)
    nb1 = size(x.B, 2)
    nc2 = size(y.C, 1)
    MPSysLin([x.A spzeros(MP, n1, n2); spzeros(MP, n2, n1) y.A],
             [x.B; spzeros(MP, n2, nb1)],
             [spzeros(MP, nc2, n1) y.C],
             [x.D spzeros(MP, n1, n2); y.B * x.C y.D],
             [x.x0; y.x0])
end

# ==============================================================================
# Diagonal composition.
# From ScicosLab file: %mpls_g_mpls.sci

"""
    Base.:(|)(y::MPSysLin, x::MPSysLin)

Diagonal composition of two Max-Plus linear systems.

# Examples
```julia-repl
julia> S1 = MPSysLin(MP([1.0 2 3;  4 5 6;  7 8 9]),
                     MP([    1.0;      2;      3]),
                     MP([    4.0       5       6]),
                     mpeye(3, 3),
                     full(zeros(MP, 3, 1)));

julia> S2 = MPSysLin(MP([10.0 11 12;  13 14 15;  16 17 18]),
                     MP([    0.0;      0;      0]),
                     MP([    0.0       0       0]),
                     mpeye(3, 3),
                     full(zeros(MP, 3, 1)));

julia> S1 | S2

Implicit dynamic linear Max-Plus system:
  x(n) = D*x(n) + A*x(n-1) + B*u(n)
  y(n) = C*x(n)
  x(0) = x0

with:
D = 6×6 Matrix{MP{Float64}}:
  0   .   .   .   .   .
  .   0   .   .   .   .
  .   .   0   .   .   .
  .   .   .   0   .   .
  .   .   .   .   0   .
  .   .   .   .   .   0

A = 6×6 Matrix{MP{Float64}}:
  1   2   3    .    .    .
  4   5   6    .    .    .
  7   8   9    .    .    .
  .   .   .   10   11   12
  .   .   .   13   14   15
  .   .   .   16   17   18

B = 6×2 Matrix{MP{Float64}}:
  1   .
  2   .
  3   .
  .   0
  .   0
  .   0

C = 2×6 Matrix{MP{Float64}}:
  4   5   6   .   .   .
  .   .   .   0   0   0

x0 = 6-element Vector{MP{Float64}}:
  .
  .
  .
  .
  .
  .
```
"""
function Base.:(|)(x::MPSysLin, y::MPSysLin)
    n1 = size(x.A, 1)
    n2 = size(y.A, 1)
    nb1 = size(x.B, 2)
    nb2 = size(y.B, 2)
    nc1 = size(x.C, 1)
    nc2 = size(y.C, 1)
    MPSysLin([x.A spzeros(MP, n1, n2); spzeros(MP, n2, n1) y.A],
             [x.B spzeros(MP, n1, nb2); spzeros(MP, n2, nb1) y.B],
             [x.C spzeros(MP, nc1, n2); spzeros(MP, nc2, n1) y.C],
             [x.D spzeros(MP, n1, n2); spzeros(MP, n2, n1) y.D],
             [x.x0; y.x0])
end

# ==============================================================================
# computes [S1;S2] that is : inputs in common, concatenation of outputs.
# From ScicosLab file: %mpls_f_mpls.sci

"""
    Base.:vcat(x::MPSysLin, y::MPSysLin)

Composition of two Max-Plus linear systems: inputs in common, concatenation of
outputs.

# Examples
```julia-repl
julia> S1 = MPSysLin(MP([1.0 2 3;  4 5 6;  7 8 9]),
                     MP([    1.0;      2;      3]),
                     MP([    4.0       5       6]),
                     mpeye(3, 3),
                     full(zeros(MP, 3, 1)));

julia> S2 = MPSysLin(MP([10.0 11 12;  13 14 15;  16 17 18]),
                     MP([    0.0;      0;      0]),
                     MP([    0.0       0       0]),
                     mpeye(3, 3),
                     full(zeros(MP, 3, 1)));

julia> [S1 ; S2]
Implicit dynamic linear Max-Plus system:
  x(n) = D*x(n) + A*x(n-1) + B*u(n)
  y(n) = C*x(n)
  x(0) = x0

with:
D = 6×6 Matrix{MP{Float64}}:
  0   .   .   .   .   .
  .   0   .   .   .   .
  .   .   0   .   .   .
  .   .   .   0   .   .
  .   .   .   .   0   .
  .   .   .   .   .   0

A = 6×6 Matrix{MP{Float64}}:
  1   2   3    .    .    .
  4   5   6    .    .    .
  7   8   9    .    .    .
  .   .   .   10   11   12
  .   .   .   13   14   15
  .   .   .   16   17   18

B = 6-element Vector{MP{Float64}}:
  1
  2
  3
  0
  0
  0

C = 2×6 Matrix{MP{Float64}}:
  4   5   6   .   .   .
  .   .   .   0   0   0

x0 = 6-element Vector{MP{Float64}}:
  .
  .
  .
  .
  .
  .
```
"""
function Base.:vcat(x::MPSysLin, y::MPSysLin)
    n1 = size(x.A, 1)
    n2 = size(y.A, 1)
    nc1 = size(x.C, 1)
    nc2 = size(y.C, 1)
    MPSysLin([x.A spzeros(MP, n1, n2); spzeros(MP, n2, n1) y.A],
             [x.B; y.B],
             [x.C spzeros(MP, nc1, n2); spzeros(MP, nc2, n1) y.C],
             [x.D spzeros(MP, n1, n2); spzeros(MP, n2, n1) y.D],
             [x.x0; y.x0])
end

# ==============================================================================
# computes [S1,S2] that is concatenation of inputs, addition of outputs
# From ScicosLab file: %mpls_c_mpls.sci

"""
    Base.:hcat(x::MPSysLin, y::MPSysLin)

Composition of two Max-Plus linear systems: concatenation of inputs, addition of
outputs.

# Examples
```julia-repl
julia> S1 = MPSysLin(MP([1.0 2 3;  4 5 6;  7 8 9]),
                     MP([    1.0;      2;      3]),
                     MP([    4.0       5       6]),
                     mpeye(3, 3),
                     full(zeros(MP, 3, 1)));

julia> S2 = MPSysLin(MP([10.0 11 12;  13 14 15;  16 17 18]),
                     MP([    0.0;      0;      0]),
                     MP([    0.0       0       0]),
                     mpeye(3, 3),
                     full(zeros(MP, 3, 1)));

julia> [S1 S2]

Implicit dynamic linear Max-Plus system:
  x(n) = D*x(n) + A*x(n-1) + B*u(n)
  y(n) = C*x(n)
  x(0) = x0

with:
D = 6×6 Matrix{MP{Float64}}:
  0   .   .   .   .   .
  .   0   .   .   .   .
  .   .   0   .   .   .
  .   .   .   0   .   .
  .   .   .   .   0   .
  .   .   .   .   .   0

A = 6×6 Matrix{MP{Float64}}:
  1   2   3    .    .    .
  4   5   6    .    .    .
  7   8   9    .    .    .
  .   .   .   10   11   12
  .   .   .   13   14   15
  .   .   .   16   17   18

B = 6×2 Matrix{MP{Float64}}:
  1   .
  2   .
  3   .
  .   0
  .   0
  .   0

C = 1×6 Matrix{MP{Float64}}:
  4   5   6   0   0   0

x0 = 6-element Vector{MP{Float64}}:
  .
  .
  .
  .
  .
  .
```
"""
function Base.:hcat(x::MPSysLin, y::MPSysLin)
    n1 = size(x.A, 1)
    n2 = size(y.A, 1)
    nb1 = size(x.B, 2)
    nb2 = size(y.B, 2)
    MPSysLin([x.A spzeros(MP, n1, n2); spzeros(MP, n2, n1) y.A],
             [x.B spzeros(MP, n1, nb2); spzeros(MP, n2, nb1) y.B],
             [x.C y.C],
             [x.D spzeros(MP, n1, n2); spzeros(MP, n2, n1) y.D],
             [x.x0; y.x0])
end

# ==============================================================================
# Feedback composition: computes star(S1*S2)*S1 in state-space form.
# From ScicosLab file: %mpls_v_mpls.sci

"""
    Base.:(/)(y::MPSysLin, x::MPSysLin)

Feedback composition: computes `mpstar(S1*S2)*S1` in state-space form.

# Examples
```julia-repl
julia> S1 = MPSysLin(MP([1.0 2 3;  4 5 6;  7 8 9]),
                     MP([    1.0;      2;      3]),
                     MP([    4.0       5       6]),
                     mpeye(3, 3),
                     full(zeros(MP, 3, 1)));

julia> S2 = MPSysLin(MP([10.0 11 12;  13 14 15;  16 17 18]),
                     MP([    0.0;      0;      0]),
                     MP([    0.0       0       0]),
                     mpeye(3, 3),
                     full(zeros(MP, 3, 1)));

julia> S1 / S2

Implicit dynamic linear Max-Plus system:
Implicit dynamic linear Max-Plus system:
  x(n) = D*x(n) + A*x(n-1) + B*u(n)
  y(n) = C*x(n)
  x(0) = x0

with:
D = 6×6 Matrix{MP{Float64}}:
  0   .   .   1   1   1
  .   0   .   2   2   2
  .   .   0   3   3   3
  4   5   6   0   .   .
  4   5   6   .   0   .
  4   5   6   .   .   0

A = 6×6 Matrix{MP{Float64}}:
  1   2   3    .    .    .
  4   5   6    .    .    .
  7   8   9    .    .    .
  .   .   .   10   11   12
  .   .   .   13   14   15
  .   .   .   16   17   18

B = 6-element Vector{MP{Float64}}:
  1
  2
  3
  .
  .
  .

C = 1×6 Matrix{MP{Float64}}:
  4   5   6   .   .   .

x0 = 6-element Vector{MP{Float64}}:
  .
  .
  .
  .
  .
  .
```
"""
function Base.:(/)(x::MPSysLin, y::MPSysLin)
    n1 = size(x.A, 1)
    n2 = size(y.A, 1)
    nb1 = size(x.B, 2)
    nc1 = size(y.C, 1)
    MPSysLin([x.A spzeros(MP, n1, n2); spzeros(MP, n2, n1) y.A],
             [x.B; spzeros(MP, n2, nb1)],
             [x.C spzeros(MP, nc1, n2)],
             [x.D x.B * y.C; y.B * x.C y.D],
             [x.x0; y.x0])
end

# ==============================================================================
# From ScicosLab file: mplssize.sci
function Base.:size(S::MPSysLin)
    [size(S.B, 1), size(S.B, 2), size(S.C, 1)]
end

# ==============================================================================
# From ScicosLab file: explicit.sci

"""
    mpexplicit(S::MPSysLin)

Convert to an explicit linear system.

# Examples
```julia-repl
julia> S = MPSysLin(MP([1.0 2; 3 4]),
                    MP([0.0; 0]),
                    MP([0.0 0]),
                    mpeye(Float64, 2,2))

Implicit dynamic linear Max-Plus system:
  x(n) = D*x(n) + A*x(n-1) + B*u(n)
  y(n) = C*x(n)
  x(0) = x0

with:
D = 2×2 Matrix{MP{Float64}}:
  0   .
  .   0

A = 2×2 Matrix{MP{Float64}}:
  1   2
  3   4

B = 2-element Vector{MP{Float64}}:
  0
  0

C = 1×2 Matrix{MP{Float64}}:
  0   0

x0 = 2-element Vector{MP{Float64}}:
  .
  .

julia> mpexplicit(S)

Implicit dynamic linear Max-Plus system:
  x(n) = D*x(n) + A*x(n-1) + B*u(n)
  y(n) = C*x(n)
  x(0) = x0

with:
D = 2×2 Matrix{MP{Float64}}:
  0   .
  .   0

A = 2×2 Matrix{MP{Float64}}:
  1   2
  3   4

B = 2-element Vector{MP{Float64}}:
  0
  0

C = 1×2 Matrix{MP{Float64}}:
  0   0

x0 = 2-element Vector{MP{Float64}}:
  .
  .
```
"""
function mpexplicit(S::MPSysLin)
    ds = mpstar(S.D)
    A = ds * S.A
    B = ds * S.B
    AC = [A; S.C]
    MPSysLin(AC, B, S.C,
             spzeros(MP, size(AC,1), size(AC,2)), S.x0)
end

# ==============================================================================

"""
    mpsimul(S::MPSysLin, u::AbsMatMP, history::Bool)

Compute states X of an autonomous linear max-plus system:
```
    x(n+1) = A ⨂ x(n)`  for n = 0 .. k
```

Where:
- `A` is a system matrix
- `x0` is a initial state vector,
- `k` is the number of iterations
- when history is set to true save all computed states, else return the last one.

# Arguments
- u: inputs to inject in the system
- history: if true then return a vector holding all results, else return the
  last result.

# Examples
```julia-repl
julia> S = MPSysLin(MP([1.0 2; 3 4]),
                    MP([0.0; 0]),
                    MP([0.0 0]),
                    mpeye(Float64, 2,2))

Implicit dynamic linear Max-Plus system:
  x(n) = D*x(n) + A*x(n-1) + B*u(n)
  y(n) = C*x(n)
  x(0) = x0

with:
D = 2×2 Matrix{MP{Float64}}:
  0   .
  .   0

A = 2×2 Matrix{MP{Float64}}:
  1   2
  3   4

B = 2-element Vector{MP{Float64}}:
  0
  0

C = 1×2 Matrix{MP{Float64}}:
  0   0

x0 = 2-element Vector{MP{Float64}}:
  .
  .

julia> mpsimul(S, MP(1:10))
1×10 Array{MP{Float64},2}:
 1  5  9  13  17  21  25  29  33  37

julia> mpsimul(S1, MP(1:10), history=false)
1×1 Array{MP{Float64},2}:
 37
```
"""
function mpsimul(S::MPSysLin, u::AbsMatMP, history::Bool)
    x = S.x0
    k = size(u, 1)
    if history
        Y = mpones(k, size(S.C, 1))
        for i = 1:k
            x = S.A * x + S.B * u[i,:]
            Y[i,:] = S.C * x
        end
        Y
    else
        for i = 1:k
            x = S.A * x + S.B * u[i,:]
        end
        Y = S.C * x
        Y[1,:]
    end
end

function mpsimul(S::MPSysLin, u::VecMP, history::Bool)
    mpsimul(S, map(x -> MP(x.λ), reshape(u, length(u), 1)), history)
end

ext: jl


text:
module ExtendableSparseTest
using ExtendableSparse
using SparseArrays
using Printf

function randtest!(A::AbstractSparseMatrix{Tv,Ti},xnnz::Int,nsplice::Int) where {Tv,Ti}
    m,n=size(A)
    S=spzeros(m,n)
    for isplice=1:nsplice
        for inz=1:xnnz
            i=rand((1:m))
            j=rand((1:n))
            a=1.0+rand(Float64)
            S[i,j]+=a
            A[i,j]+=a
            @assert(nnz(S)==nnz(A))
        end
        flush!(A)
        for j=1:n
            @assert(issorted(A.cscmatrix.rowval[A.cscmatrix.colptr[j]:A.cscmatrix.colptr[j+1]-1]))
        end
        @assert(nnz(S)==nnz(A))

        (I,J,V)=findnz(S)
        for inz=1:nnz(S)
            @assert(A[I[inz],J[inz]]==V[inz])
        end

        (I,J,V)=findnz(A)
        for inz=1:nnz(A)
            @assert(S[I[inz],J[inz]]==V[inz])
        end

    end
    return true
end

function check(;m=1000,n=1000,nnz=5000,nsplice=1)
    mat=ExtendableSparseMatrix{Float64,Int64}(m,n)
    return randtest!(mat,nnz,nsplice)
end

function benchmark(;n=10000,m=10000,nnz=50000)

    println("SparseMatrixCSC:")
    mat=spzeros(Float64,Int64,m,n)
    @time sprand!(mat,nnz)

    println("SparseMatrixExtension:")
    extmat=SparseMatrixExtension{Float64,Int64}(m,n)
    @time sprand!(extmat,nnz)

    println("ExtendableSparseMatrix:")
    xextmat=ExtendableSparseMatrix{Float64,Int64}(m,n)
    @time begin
        sprand!(xextmat,nnz)
        @inbounds flush!(xextmat)
    end
    return
    b
end

end


text:
module mod

const constant = nothing

function func()
end

myfunc() = nothing

abstract type MyType end

struct MyStruct <: MyType
    MyStruct() = new()
    method() = nothing
end

end


text:
"""
    CGDescent

Conjugate gradient method (Hager-Zhang version [W.Hager, H.Zhang // SIAM J. Optim (2006) Vol. 16, pp. 170-192])
"""
mutable struct CGDescent{T<:AbstractFloat,V<:AbstractVector{T}} <: CoreMethod
    x::V
    xpre::V
    g::V
    gpre::V
    gdiff::V
    dir::V
    y::T
    α::T
    α0::T
end

@inline fnval(M::CGDescent) = M.y
@inline gradientvec(M::CGDescent) = M.g
@inline argumentvec(M::CGDescent) = M.x
@inline step_origin(M::CGDescent) = M.xpre

function CGDescent(x::AbstractVector)
    F = float(eltype(x))
    CGDescent(similar(x, F),
              similar(x, F),
              similar(x, F),
              similar(x, F),
              similar(x, F),
              similar(x, F),
              zero(F),
              convert(F, 0.01),
              convert(F, 0.01)
             )
end

function __descent_dir!(M::CGDescent)
    d = M.dir
    y = M.gdiff
    dty = dot(d, y)
    if iszero(dty)
        # d' * y == 0 means we are at the first iteration
        β = dty
    else
        β = (dot(y, M.g) - 2 * dot(d, M.g) * dot(y, y) / dty) / dty
    end
    η = one(β) / 100
    β = max(β, -1 / (norm(d) * min(η, norm(M.g))))
    map!(d, d, M.g) do a, b
        muladd(a, β, -b)
    end
    return d
end

function init!(M::CGDescent{T}, optfn!, x0; reset, constrain_step = infstep) where {T}
    y, g = optfn!(x0, zero(T), x0)
    __update_gpre!(M, M.g)
    map!(-, M.dir, M.g)
    M.α = M.α0
    return
end

@inline reset!(::CGDescent) = nothing

function reset!(M::CGDescent, x0)
    if length(M.x) != length(x0)
        foreach((M.x, M.xpre, M.g, M.gpre, M.dir)) do v
            resize!(v, length(x0))
        end
    end
    return
end

function callfn!(M::CGDescent, fdf, x, α, d)
    __update_arg!(M, x, α, d)
    let x = argumentvec(M)
        (y, g) = fdf(x, gradientvec(M))
        __update_grad!(M, g)
        M.y = y
        return (y, g)
    end
end

function step!(M::CGDescent, optfn!; constrain_step = infstep)
    M.x, M.xpre = M.xpre, M.x
    map!(-, M.gdiff, M.g, M.gpre)

    d = __descent_dir!(M)
    xpre = M.xpre
    M.g, M.gpre = M.gpre, M.g
    ypre = M.y
    maxstep = constrain_step(xpre, d)
    α = strong_backtracking!(optfn!, xpre, d, ypre, M.gpre, α = M.α, αmax = maxstep, β = 0.01, σ = 0.1)
    fdiff = M.y - ypre
    if fdiff < 0
        M.α = 2 * fdiff / dot(d, M.gpre)
    end
    return α
end

function __update_arg!(M::CGDescent, x, α, d)
    map!(M.x, d, x) do a, b
        muladd(α, a, b)
    end
    return
end

@inline function __update_arg!(M::CGDescent, x)
    if x !== M.x
        copy!(M.x, x)
    end
    return
end

@inline function __update_xpre!(M::CGDescent, x)
    if x !== M.xpre
        copy!(M.xpre, x)
    end
    return
end

@inline function __update_grad!(M::CGDescent, g)
    if M.g !== g
        copy!(M.g, g)
    end
    return
end

@inline function __update_gpre!(M::CGDescent, g)
    if M.gpre !== g
        copy!(M.gpre, g)
    end
    return
end

@inline function __zero_dir!(M::CGDescent)
    fill!(M.dir, 0)
    return
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %%
using Distributions
using HypothesisTests
using StatsPlots
default(size=(500, 300))

function pval_t(null, X)
    μ = mean(null)
    n = length(X)
    T = (mean(X) - μ)/√(var(X)/n)
    2ccdf(TDist(n-1), abs(T))
end

function pval_wilcoxon(null, X)
    m = quantile(null, 0.5) # median of null
    W = HypothesisTests.SignedRankTest(X .- m)
    pvalue(W)
end

function plot_pvals(;
        pval = pval_wilcoxon,
        null = Gamma(4, 1),
        n = 100,
        L = 10^6
    )
    PVal = similar(zeros(), L)
    for i in 1:L
        X = rand(null, n)
        PVal[i] = pval(null, X)
    end
    realalpha = mean(PVal .< 0.05)
    histogram(PVal; norm=true, alpha=0.3, bin=0.001:0.025:1.001, label="")
    plot!(; xtick=0:0.05:1, xrotation=90, bottom_margin=3Plots.mm)
    title!("""
        pval = $pval,  n = $n
        null = $null
        P(p-value < 0.05) = $realalpha,  niters = $L""",
        titlefontsize=8)
end

function plot_both(;
        null = Gamma(4, 1),
        n = 100,
        L = 10^6
    )
    P = plot_pvals(; pval = pval_t, null, n, L)
    Q = plot_pvals(; pval = pval_wilcoxon, null, n, L)
    plot(P, Q; size=(500, 600), layout=(2, 1), tickfontsize=8)
end

# %%
plot(Gamma(4, 1); label="Gamma(4, 1)")

# %%
plot_both(; null = Gamma(4, 1), n = 30)

# %%
plot_both(; null = Gamma(4, 1), n = 100)

# %%
plot_both(null = Gamma(4, 1), n = 300, L=10^5)

# %%
plot_both(; null = Gamma(4, 1), n = 1000, L=10^5)

# %%
plot(Gamma(10, 1); label="Gamma(10, 1)")

# %%
plot_both(; null = Gamma(10, 1))

# %%
plot(Gamma(20, 1); label="Gamma(20, 1)")

# %%
plot_both(; null = Gamma(20, 1))

# %%
plot(Gamma(50, 1), 20, 80; label="Gamma(50, 1)")

# %%
plot_both(; null = Gamma(50, 1), n = 30)

# %%
plot_both(; null = Gamma(50, 1), n = 100)

# %%
plot_both(; null=Gamma(50, 1), n = 300, L=10^5)

# %%
plot_both(; null=Gamma(50, 1), n = 1000, L=10^5)

# %%
plot(Uniform(-1, 1), -3, 3; label="Uniform(-1, 1)")

# %%
plot_both(; null = Uniform(-1, 1))

# %%
null = TDist(3)
plot(null, -8, 8; label="TDist(3)")
plot!(Normal(mean(null), std(null)), -8, 8; label="normal approx.", ls=:dash)

# %%
plot_both(; null)

# %%
null = MixtureModel([Normal(), Normal(10, 1)], [0.95, 0.05])
plot(null, -5, 15; label="", title="$null", titlefontsize=8) |> display
plot_both(; null, n = 30)

# %%
plot_both(; null, n = 100)

# %%
plot_both(; null, n = 300, L = 10^5)

# %%
plot_both(; null, n = 1000, L = 10^5)

# %%
null = Exponential()
plot(null; label="", title="$null", titlefontsize=8) |> display
plot_both(; null, n = 30)

# %%
plot_both(; null, n = 100)

# %%

max_stars_repo_path: 0022/Wilcoxon test vs. t-test for asymmetric null.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
module FileTransportTests

using Thrift
using Test

function testfiletransport()
    println("\nTesting file transport...")

    fname = tempname()
    s1 = "Hello World"
    open(fname, "w") do f
        t = TFileTransport(f)
        p = TCompactProtocol(t)

        @test open(t) === nothing
        @test close(t) === nothing
        @test flush(t) === nothing
        @test isa(Thrift.rawio(t), IO)
        @test isopen(t)

        for i in 1:5
            write(p, s1)
        end
    end

    open(fname, "r") do f
        t = TFileTransport(f)
        p = TCompactProtocol(t)
        for i in 1:5
            s2 = read(p, typeof(s1))
            @test s2 == s1
        end
        println("passed.")
    end
    rm(fname)
end

testfiletransport()

end


text:
function CatKeyMul(A::Assoc,B::Assoc)
    if isa(Col(A)[1],String) && isa(Row(B)[1],String)
        AB = sortedintersect(A.col,B.row)
        A = A[:,AB]
        B = B[AB,:]
        rrr,ccc,vvv = findnz(Adj(A*B))
        ABVal = Array(Union{AbstractString,Number},length(rrr))
        for i in 1:length(rrr)
            r = rrr[i]
            c = ccc[i]
            ABvalList = sortedintersect(Col(A[r,:]),Row(B[:,c]))
            if length(ABvalList) > 0
                val = join(ABvalList,";")*";"
                ABVal[i] = val
            end
        end
        return Assoc(Row(A)[rrr],Col(B)[ccc],ABVal)
    else
        return A*B
    end

end


text:
#
# IRKCoefficients.jl file:
#   PolInterp
#   MuCoefficients!
#   HCoefficients!
#   EstimateCoeffs!

function PolInterp(X::AbstractVector{ctype}, Y::AbstractMatrix{ctype},
                   Z::AbstractVector{ctype}) where {ctype}
    N = length(X)
    M = length(Z)
    K = size(Y,1)
    if size(Y,2)!=N
        error("columns(Y) != length(X)")
    end
    pz = zeros(ctype,K,M)

    @inbounds begin
    for i=1:N
        lag = 1.
        for j=1:N
            if (j!=i)
                lag *= X[i]-X[j]
            end
        end
        lag = 1/lag
        for m=1:M
            liz = lag
            for j=1:N
                if (j!=i)
                    liz *= Z[m]-X[j]
                end
            end
            for k=1:K
                pz[k,m] += Y[k,i]*liz
            end
        end
    end
    end
    return pz
end

function MuCoefficients!(mu,T::Type{<:CompiledFloats})

    mu[1,1]=convert(T,0.5)

    mu[2,1]=convert(T,+1.0818949631055814971365081647359309e+00)
    mu[2,2]=convert(T,0.5)

    mu[3,1]=convert(T,+9.5995729622205494766003095439844678e-01)
    mu[3,2]=convert(T,+1.0869589243008327233290709646162480e+00)
    mu[3,3]=convert(T,0.5)

    mu[4,1]=convert(T,+1.0247213458032003748680445816450829e+00)
    mu[4,2]=convert(T,+9.5505887369737431186016905653386876e-01)
    mu[4,3]=convert(T,+1.0880938387323083134422138713913203e+00)
    mu[4,4]=convert(T,0.5)

    mu[5,1]=convert(T,+9.8302382676362890697311829123888390e-01)
    mu[5,2]=convert(T,+1.0287597754747493109782305570410685e+00)
    mu[5,3]=convert(T,+9.5383453518519996588326911440754302e-01)
    mu[5,4]=convert(T,+1.0883471611098277842507073806008045e+00)
    mu[5,5]=convert(T,0.5)

    mu[6,1]=convert(T,+1.0122259141132982060539425317219435e+00)
    mu[6,2]=convert(T,+9.7998287236359129082628958290257329e-01)
    mu[6,3]=convert(T,+1.0296038730649779374630125982121223e+00)
    mu[6,4]=convert(T,+9.5383453518519996588326911440754302e-01)
    mu[6,5]=convert(T,+1.0880938387323083134422138713913203e+00)
    mu[6,6]=convert(T,0.5)

    mu[7,1]=convert(T,+9.9125143323080263118822334698608777e-01)
    mu[7,2]=convert(T,+1.0140743558891669291459735166525994e+00)
    mu[7,3]=convert(T,+9.7998287236359129082628958290257329e-01)
    mu[7,4]=convert(T,+1.0287597754747493109782305570410685e+00)
    mu[7,5]=convert(T,+9.5505887369737431186016905653386876e-01)
    mu[7,6]=convert(T,+1.0869589243008327233290709646162480e+00)
    mu[7,7]=convert(T,0.5)

    mu[8,1]=convert(T,+1.0054828082532158826793409353214951e+00)
    mu[8,2]=convert(T,+9.9125143323080263118822334698608777e-01)
    mu[8,3]=convert(T,+1.0122259141132982060539425317219435e+00)
    mu[8,4]=convert(T,+9.8302382676362890697311829123888390e-01)
    mu[8,5]=convert(T,+1.0247213458032003748680445816450829e+00)
    mu[8,6]=convert(T,+9.5995729622205494766003095439844678e-01)
    mu[8,7]=convert(T,+1.0818949631055814971365081647359309e+00)
    mu[8,8]=convert(T,0.5)

    s=8
    @inbounds begin
    for i in 1:s
        for j in i+1:s
            mu[i,j] = 1 - mu[j,i]
        end
    end
    end

end

function MuCoefficients!(mu,T)

#          println("Mu: NO CompiledFloats !!!")

          mu[1,1]=parse(T,"0.5")

          mu[2,1]=parse(T,"+1.0818949631055814971365081647359309e+00")
          mu[2,2]=parse(T,"0.5")

          mu[3,1]=parse(T,"+9.5995729622205494766003095439844678e-01")
          mu[3,2]=parse(T,"+1.0869589243008327233290709646162480e+00")
          mu[3,3]=parse(T,"0.5")

          mu[4,1]=parse(T,"+1.0247213458032003748680445816450829e+00")
          mu[4,2]=parse(T,"+9.5505887369737431186016905653386876e-01")
          mu[4,3]=parse(T,"+1.0880938387323083134422138713913203e+00")
          mu[4,4]=parse(T,"0.5")

          mu[5,1]=parse(T,"+9.8302382676362890697311829123888390e-01")
          mu[5,2]=parse(T,"+1.0287597754747493109782305570410685e+00")
          mu[5,3]=parse(T,"+9.5383453518519996588326911440754302e-01")
          mu[5,4]=parse(T,"+1.0883471611098277842507073806008045e+00")
          mu[5,5]=parse(T,"0.5")

          mu[6,1]=parse(T,"+1.0122259141132982060539425317219435e+00")
          mu[6,2]=parse(T,"+9.7998287236359129082628958290257329e-01")
          mu[6,3]=parse(T,"+1.0296038730649779374630125982121223e+00")
          mu[6,4]=parse(T,"+9.5383453518519996588326911440754302e-01")
          mu[6,5]=parse(T,"+1.0880938387323083134422138713913203e+00")
          mu[6,6]=parse(T,"0.5")

          mu[7,1]=parse(T,"+9.9125143323080263118822334698608777e-01")
          mu[7,2]=parse(T,"+1.0140743558891669291459735166525994e+00")
          mu[7,3]=parse(T,"+9.7998287236359129082628958290257329e-01")
          mu[7,4]=parse(T,"+1.0287597754747493109782305570410685e+00")
          mu[7,5]=parse(T,"+9.5505887369737431186016905653386876e-01")
          mu[7,6]=parse(T,"+1.0869589243008327233290709646162480e+00")
          mu[7,7]=parse(T,"0.5")

          mu[8,1]=parse(T,"+1.0054828082532158826793409353214951e+00")
          mu[8,2]=parse(T,"+9.9125143323080263118822334698608777e-01")
          mu[8,3]=parse(T,"+1.0122259141132982060539425317219435e+00")
          mu[8,4]=parse(T,"+9.8302382676362890697311829123888390e-01")
          mu[8,5]=parse(T,"+1.0247213458032003748680445816450829e+00")
          mu[8,6]=parse(T,"+9.5995729622205494766003095439844678e-01")
          mu[8,7]=parse(T,"+1.0818949631055814971365081647359309e+00")
          mu[8,8]=parse(T,"0.5")

          s=8
          @inbounds begin
          for i in 1:s
              for j in i+1:s
                  mu[i,j] = 1 - mu[j,i]
              end
          end
          end

end

function HCoefficients!(mu,hc,hb,nu,h,hprev,T::Type{<:CompiledFloats})

          hb .= convert.(T,
                [ +5.0614268145188129576265677154981094e-02,
                  +1.1119051722668723527217799721312045e-01,
                  +1.5685332293894364366898110099330067e-01,
                  +1.8134189168918099148257522463859781e-01,
                  +1.8134189168918099148257522463859781e-01,
                  +1.5685332293894364366898110099330067e-01,
                  +1.1119051722668723527217799721312045e-01,
                  +5.0614268145188129576265677154981094e-02
               ])

          hc .= convert.(T,
                [ +1.9855071751231884158219565715263505e-02,
                  +1.0166676129318663020422303176208480e-01,
                  +2.3723379504183550709113047540537686e-01,
                  +4.0828267875217509753026192881990801e-01,
                  +5.9171732124782490246973807118009203e-01,
                  +7.6276620495816449290886952459462321e-01,
                  +8.9833323870681336979577696823791522e-01,
                  +9.8014492824876811584178043428473653e-01
          ])

  s = length(hb)

""" Interpolate coefficients """

   if (hprev==0.)
       nu.=zeros(T,s,s)
   else
       lambda=h/hprev
       X=vcat(-hc[end:-1:1],[0])
       Y=hcat(mu,zeros(T,s))
       Z=lambda*hc
       nu.=-(PolInterp(X,Y,Z))'
   end

""" hb """

  hb.=hb*h
  hb1 = (h-sum(hb[2:end-1]))/2
  hb[1] = hb1
  hb[end] = hb1

  hc.=hc*h

return

end

function HCoefficients!(mu,hc,hb,nu,h,hprev,T)

#          println("HCoefficients: NO Compiled")

          hb .= [parse(T,"+5.0614268145188129576265677154981094e-02"),
                 parse(T,"+1.1119051722668723527217799721312045e-01"),
                 parse(T,"+1.5685332293894364366898110099330067e-01"),
                 parse(T,"+1.8134189168918099148257522463859781e-01"),
                 parse(T,"+1.8134189168918099148257522463859781e-01"),
                 parse(T,"+1.5685332293894364366898110099330067e-01"),
                 parse(T,"+1.1119051722668723527217799721312045e-01"),
                 parse(T,"+5.0614268145188129576265677154981094e-02")
               ]

          hc .= [parse(T,"+1.9855071751231884158219565715263505e-02"),
                 parse(T,"+1.0166676129318663020422303176208480e-01"),
                 parse(T,"+2.3723379504183550709113047540537686e-01"),
                 parse(T,"+4.0828267875217509753026192881990801e-01"),
                 parse(T,"+5.9171732124782490246973807118009203e-01"),
                 parse(T,"+7.6276620495816449290886952459462321e-01"),
                 parse(T,"+8.9833323870681336979577696823791522e-01"),
                 parse(T,"+9.8014492824876811584178043428473653e-01")
          ]

  s = length(hb)

""" Interpolate coefficients """

   if (hprev==0.)
       nu.=zeros(T,s,s)
   else
       lambda=h/hprev
       X=vcat(-hc[end:-1:1],[0])
       Y=hcat(mu,zeros(T,s))
       Z=lambda*hc
       nu.=-(PolInterp(X,Y,Z))'
   end

""" hb """

  hb.=hb*h
  hb1 = (h-sum(hb[2:end-1]))/2
  hb[1] = hb1
  hb[end] = hb1

  hc.=hc*h

return

end

function EstimateCoeffs!(alpha, T::Type{<:CompiledFloats})

    c= convert.(T,
          [ +1.9855071751231884158219565715263505e-02,
            +1.0166676129318663020422303176208480e-01,
            +2.3723379504183550709113047540537686e-01,
            +4.0828267875217509753026192881990801e-01,
            +5.9171732124782490246973807118009203e-01,
            +7.6276620495816449290886952459462321e-01,
            +8.9833323870681336979577696823791522e-01,
            +9.8014492824876811584178043428473653e-01
    ])

    s = length(c)

    B=vcat(zeros(T,s-1),1)
    M=[(c[i]-1/2)^(k-1) for i in 1:s, k in 1:s]'
    alpha.=1000*M\B

return

end

function EstimateCoeffs!(alpha,T)

#    println("Estimate2: NO Compiled")

    c = [parse(T,"+1.9855071751231884158219565715263505e-02"),
          parse(T,"+1.0166676129318663020422303176208480e-01"),
          parse(T,"+2.3723379504183550709113047540537686e-01"),
          parse(T,"+4.0828267875217509753026192881990801e-01"),
          parse(T,"+5.9171732124782490246973807118009203e-01"),
          parse(T,"+7.6276620495816449290886952459462321e-01"),
          parse(T,"+8.9833323870681336979577696823791522e-01"),
          parse(T,"+9.8014492824876811584178043428473653e-01")
        ]

    s = length(c)

    B=vcat(zeros(T,s-1),1)
    M=[(c[i]-1/2)^(k-1) for i in 1:s, k in 1:s]'
    alpha.=1000*M\B

return

end


text:
using OrthoPoly
using Polynomials
using Test

function orthogonal_poly_test(F,Ps,num=[])
    if num == []
        num = ones(eltype(Ps[1]), length(Ps))
    end
    for n=0:length(Ps)-1
        f = F(n)
        p = Poly(Ps[n+1])/num[n+1]
        @test f == p
    end
end

println("Testing the first 11 Legendre polynomials")
orthogonal_poly_test(legendre, [[1],[0,1],[-1,0,3],[0,-3,0,5],[3,0,-30,0,35],[0,15,0,-70,0,63],[-5,0,105,0,-315,0,231],[0,-35,0,315,0,-693,0,429],[35,0,-1260,0,6930,0,-12012,0,6435],[0,315,0,-4620,0,18018,0,-25740,0,12155],[-63,0,3465,0,-30030,0,90090,0,-109395,0,46189]], [1 1 2 2 8 8 16 16 128 128 256])

println("Testing the first 16 Chebyshev polynomials of the first kind")
orthogonal_poly_test(chebyshevT, [[1],[0,1],[-1,0,2],[0,-3,0,4],[1,0,-8,0,8],[0,5,0,-20,0,16],[-1,0,18,0,-48,0,32],[0,-7,0,56,0,-112,0,64],[1,0,-32,0,160,0,-256,0,128],[0,9,0,-120,0,432,0,-576,0,256],[-1,0,50,0,-400,0,1120,0,-1280,0,512],[0,-11,0,220,0,-1232,0,2816,0,-2816,0,1024],[1,0,-72,0,840,0,-3584,0,6912,0,-6144,0,2048],[0,13,0,-364,0,2912,0,-9984,0,16640,0,-13312,0,4096],[-1,0,98,0,-1568,0,9408,0,-26880,0,39424,0,-28672,0,8192],[0,-15,0,560,0,-6048,0,28800,0,-70400,0,92160,0,-61440,0,16384]])

println("Testing the first 16 Chebyshev polynomials of the second kind")
orthogonal_poly_test(chebyshevU, [[1],[0,2],[-1,0,4],[0,-4,0,8],[1,0,-12,0,16],[0,6,0,-32,0,32],[-1,0,24,0,-80,0,64],[0,-8,0,80,0,-192,0,128],[1,0,-40,0,240,0,-448,0,256],[0,10,0,-160,0,672,0,-1024,0,512],[-1,0,60,0,-560,0,1792,0,-2304,0,1024],[0,-12,0,280,0,-1792,0,4608,0,-5120,0,2048],[1,0,-84,0,1120,0,-5376,0,11520,0,-11264,0,4096]])

println("Testing the first 11 Hermite polynomials")
orthogonal_poly_test(hermite_prob, [[1],[0,1],[-1,0,1],[0,-3,0,1],[3,0,-6,0,1],[0,15,0,-10,0,1],[-15,0,45,0,-15,0,1],[0,-105,0,105,0,-21,0,1],[105,0,-420,0,210,0,-28,0,1],[0,945,0,-1260,0,378,0,-36,0,1],[-945,0,4725,0,-3150,0,630,0,-45,0,1]])
orthogonal_poly_test(hermite, [[1],[0,2],[-2,0,4],[0,-12,0,8],[12,0,-48,0,16],[0,120,0,-160,0,32],[-120,0,720,0,-480,0,64],[0,-1680,0,3360,0,-1344,0,128],[1680,0,-13440,0,13440,0,-3584,0,256],[0,30240,0,-80640,0,48384,0,-9216,0,512],[-30240,0,302400,0,-403200,0,161280,0,-23040,0,1024]])

println("Testing the first 9 Laguerre polynomials")
orthogonal_poly_test(n -> laguerre(n, Poly([0,1//1])), [[1],[1,-1],[2,-4,1],[6,-18,9,-1],[24,-96,72,-16,1],[120,-600,600,-200,25,-1],[720,-4320,5400,-2400,450,-36,1],[5040,-35280,52920,-29400,7350,-882,49,-1],[40320,-322560,564480,-376320,117600,-18816,1568,-64,1]], [1 1 2 6 24 120 720 5040 40320]//1)

println("Testing the first 4 associated Laguerre polynomials")
for α=[1,1.0,2,2.0,2.1]
    orthogonal_poly_test(n -> laguerre_assoc(n,α,Poly([0,one(α)])),
                         [[1],[1+α,-1],[(α+2)*(α+1),-2(α+2),1],[(α+1)*(α+2)*(α+3),-3(α+2)*(α+3),3(α+3),-1]], [1.0 1 2 6])
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -12.331138821460632 <= q <= 11.287379993959185)

                     add_NL_constraint(m, :(log(1 + exp(0.16893740747119512 + 0.361142465903193*log(1 + exp(0.3759145594542668 + 0.23452173305395663*log(1 + exp(-0.8645963938709578 + -0.1467408759040536*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.7957413591360925*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + 0.6670298352855966*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + 0.192500882174786*log(1 + exp(-0.9609056399980664 + 0.06060702618594416*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + -0.49474851412125664*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.5046120388520428*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + -0.6554806276835787*log(1 + exp(-0.9973666264873402 + 0.1126936245130743*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.2338493567409592*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.8532702233572196*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))))) + 0.21914525995597334*log(1 + exp(-0.017859443401365915 + -0.7776638835685663*log(1 + exp(-0.8645963938709578 + -0.1467408759040536*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.7957413591360925*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + 0.6670298352855966*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + -0.3390165730930925*log(1 + exp(-0.9609056399980664 + 0.06060702618594416*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + -0.49474851412125664*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.5046120388520428*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + -0.16349638009140932*log(1 + exp(-0.9973666264873402 + 0.1126936245130743*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.2338493567409592*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.8532702233572196*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))))) + 0.7287835676612944*log(1 + exp(0.6615910141517634 + 0.18731902906517428*log(1 + exp(-0.8645963938709578 + -0.1467408759040536*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.7957413591360925*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + 0.6670298352855966*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + -0.853933007586007*log(1 + exp(-0.9609056399980664 + 0.06060702618594416*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + -0.49474851412125664*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.5046120388520428*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + -0.9912964912419922*log(1 + exp(-0.9973666264873402 + 0.1126936245130743*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.2338493567409592*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.8532702233572196*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))))))) + log(1 + exp(-0.9123133183860701 + 0.273513749523993*log(1 + exp(0.3759145594542668 + 0.23452173305395663*log(1 + exp(-0.8645963938709578 + -0.1467408759040536*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.7957413591360925*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + 0.6670298352855966*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + 0.192500882174786*log(1 + exp(-0.9609056399980664 + 0.06060702618594416*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + -0.49474851412125664*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.5046120388520428*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + -0.6554806276835787*log(1 + exp(-0.9973666264873402 + 0.1126936245130743*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.2338493567409592*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.8532702233572196*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))))) + 0.052599029592083735*log(1 + exp(-0.017859443401365915 + -0.7776638835685663*log(1 + exp(-0.8645963938709578 + -0.1467408759040536*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.7957413591360925*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + 0.6670298352855966*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + -0.3390165730930925*log(1 + exp(-0.9609056399980664 + 0.06060702618594416*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + -0.49474851412125664*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.5046120388520428*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + -0.16349638009140932*log(1 + exp(-0.9973666264873402 + 0.1126936245130743*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.2338493567409592*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.8532702233572196*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))))) + 0.2642412570560455*log(1 + exp(0.6615910141517634 + 0.18731902906517428*log(1 + exp(-0.8645963938709578 + -0.1467408759040536*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.7957413591360925*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + 0.6670298352855966*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + -0.853933007586007*log(1 + exp(-0.9609056399980664 + 0.06060702618594416*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + -0.49474851412125664*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.5046120388520428*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + -0.9912964912419922*log(1 + exp(-0.9973666264873402 + 0.1126936245130743*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.2338493567409592*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.8532702233572196*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))))))) + log(1 + exp(0.5514504616591211 + 0.27732859795053555*log(1 + exp(0.3759145594542668 + 0.23452173305395663*log(1 + exp(-0.8645963938709578 + -0.1467408759040536*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.7957413591360925*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + 0.6670298352855966*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + 0.192500882174786*log(1 + exp(-0.9609056399980664 + 0.06060702618594416*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + -0.49474851412125664*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.5046120388520428*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + -0.6554806276835787*log(1 + exp(-0.9973666264873402 + 0.1126936245130743*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.2338493567409592*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.8532702233572196*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))))) + -0.08702067567597149*log(1 + exp(-0.017859443401365915 + -0.7776638835685663*log(1 + exp(-0.8645963938709578 + -0.1467408759040536*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.7957413591360925*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + 0.6670298352855966*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + -0.3390165730930925*log(1 + exp(-0.9609056399980664 + 0.06060702618594416*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + -0.49474851412125664*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.5046120388520428*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + -0.16349638009140932*log(1 + exp(-0.9973666264873402 + 0.1126936245130743*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.2338493567409592*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.8532702233572196*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))))) + 0.40713382739201665*log(1 + exp(0.6615910141517634 + 0.18731902906517428*log(1 + exp(-0.8645963938709578 + -0.1467408759040536*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.7957413591360925*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + 0.6670298352855966*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + -0.853933007586007*log(1 + exp(-0.9609056399980664 + 0.06060702618594416*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + -0.49474851412125664*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.5046120388520428*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + -0.9912964912419922*log(1 + exp(-0.9973666264873402 + 0.1126936245130743*log(1 + exp(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3]))) + 0.2338493567409592*log(1 + exp(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3]))) + -0.8532702233572196*log(1 + exp(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# using Julia 0.6.2
iter = 100
println("$(iter) iterations for each n")
for n in (10,50,100,150,200,300,400,500,10000,20000,50000,100000)
	# construct the polynomial
	A = randn!(Array{Float64}(n+1)) # coefficients of the polynomial
	T = Array{Float64}(0)
	# run several times and find the mean time
	for i in 1:iter
		X = fill!(Array{Float64}(n+1), 1)
		x = Float64(rand()) # parameter of the polynomial
		tic()
		for j in 2:n+1
			X[j] = x * X[j-1]
		end
		P_x = dot(A, X)
		push!(T, toq())
	end
	Tmean = mean(T)
	println("n = $(n) \t\tTmean = $(Tmean)")
end

#using Julia 0.6.2
iter = 100
println("$(iter) iterations for each n")
for n in (10,50,100,150,200,300,400,500,10000,20000,50000,100000)
	# construct the polynomial
	A = randn!(Array{Float64}(n+1))
	T = Array{Float64}(0)
	# run several times and find the mean time
	for i in 1:iter
		x = Float64(rand())
		P_x = 0.
		tic()
		for j in 1:n+1
			P_x += A[j] * x^(j-1)
		end
		push!(T, toq())
	end
	Tmean = mean(T)
	println("n = $(n) \t\tTmean = $(Tmean)")
end

#using Julia 0.6.2
iter = 100
println("$(iter) iterations for each n")
for n in (10,50,100,150,200,300,400,500,10000,20000,50000,100000)
	# construct the polynomial
	A = randn!(Array{Float64}(n+1))
	T = Array{Float64}(0)
	# run several times and find the mean time
	for i in 1:iter
		x = Float64(rand())
		P_x = 0.
		xj = 1.
		tic()
		for j in 1:n+1
			P_x += xj * A[j]
			xj *= x
		end
		push!(T, toq())
	end
	Tmean = mean(T)
	println("n = $(n) \t\tTmean = $(Tmean)")
end

#using Julia 0.6.2
iter = 100
println("$(iter) iterations for each n")
for n in (10,50,100,150,200,300,400,500,10000,20000,50000,100000)
	# construct the polynomial
	A = randn!(Array{Float64}(n+1))
	T = Array{Float64}(0)
	# run several times and find the mean time
	for i in 1:iter
		x = Float64(rand())
		P_x = A[end]
		tic()
		for j in n:-1:1
			P_x = x * P_x + A[j]
		end
		push!(T, toq())
	end
	Tmean = mean(T)
	println("n = $(n) \t\tTmean = $(Tmean)")
end

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes