text:
import PackageCompiler, TOML

if length(ARGS) < 1 || length(ARGS) > 2
    println("Usage: julia $PROGRAM_FILE target_dir [major|minor]")
    println()
    println("where:")
    println("    target_dir is the directory to use to create the library bundle")
    println("    [major|minor] is the (optional) compatibility version (default: major).")
    println("                  Use 'minor' if you use new/non-backwards-compatible functionality.")
    println()
    println("[major|minor] is only useful on OSX.")
    exit(1)
end

const build_dir = @__DIR__
const target_dir = ARGS[1]
const project_toml = realpath(joinpath(build_dir, "..", "Project.toml"))
const version = VersionNumber(TOML.parsefile(project_toml)["version"])

const compatibility = length(ARGS) == 2 ? ARGS[2] : "major"

PackageCompiler.create_library(".", target_dir;
                            lib_name="cg",
                            precompile_execution_file=[joinpath(build_dir, "generate_precompile.jl")],
                            precompile_statements_file=[joinpath(build_dir, "additional_precompile.jl")],
                            incremental=false,
                            filter_stdlibs=true,
                            header_files = [joinpath(build_dir, "cg.h")],
                            force=true,
                            version=version,
                            compat_level=compatibility,
                        )

max_stars_repo_name: simonbyrne/libcg

max_issues_repo_name: simonbyrne/libcg

max_forks_repo_name: simonbyrne/libcg


text:
# This file is part of the TaylorSeries.jl Julia package, MIT license
#
# Luis Benet & David P. Sanders
# UNAM
#
# MIT Expat license
#

# Arithmetic operations: +, -, *, /

## Equality ##
for T in (:Taylor1, :TaylorN)

    @eval begin
        ==(a::$T{T}, b::$T{S}) where {T<:Number, S<:Number} = ==(promote(a,b)...)

        function ==(a::$T{T}, b::$T{T}) where {T<:Number}
            if a.order != b.order
                a, b = fixorder(a, b)
            end
            return a.coeffs == b.coeffs
        end
    end
end

function ==(a::HomogeneousPolynomial, b::HomogeneousPolynomial)
    a.order == b.order && return a.coeffs == b.coeffs
    return iszero(a.coeffs) && iszero(b.coeffs)
end

for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
    @eval iszero(a::$T) = iszero(a.coeffs)
end

## zero and one ##
for T in (:Taylor1, :TaylorN), f in (:zero, :one)
    @eval ($f)(a::$T) = $T(($f)(a[0]), a.order)
end

function zero(a::HomogeneousPolynomial{T}) where {T<:Number}
    v = zero.(a.coeffs)
    return HomogeneousPolynomial(v, a.order)
end

function zeros(::HomogeneousPolynomial{T}, order::Int) where {T<:Number}
    order == 0 && return [HomogeneousPolynomial([zero(T)], 0)]
    v = Array{HomogeneousPolynomial{T}}(undef, order+1)
    @simd for ord in eachindex(v)
        @inbounds v[ord] = HomogeneousPolynomial([zero(T)], ord-1)
    end
    return v
end

zeros(::Type{HomogeneousPolynomial{T}}, order::Int) where {T<:Number} =
    zeros( HomogeneousPolynomial([zero(T)], 0), order)

function one(a::HomogeneousPolynomial{T}) where {T<:Number}
    v = one.(a.coeffs)
    return HomogeneousPolynomial(v, a.order)
end

function ones(a::HomogeneousPolynomial{T}, order::Int) where {T<:Number}
    order == 0 && return [HomogeneousPolynomial([one(a[1])], 0)]
    v = Array{HomogeneousPolynomial{T}}(undef, order+1)
    @simd for ord in eachindex(v)
        @inbounds num_coeffs = size_table[ord]
        @inbounds v[ord] = HomogeneousPolynomial(ones(T, num_coeffs), ord-1)
    end
    return v
end

ones(::Type{HomogeneousPolynomial{T}}, order::Int) where {T<:Number} =
    ones( HomogeneousPolynomial([one(T)], 0), order)

## Addition and substraction ##
for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))

    for T in (:Taylor1, :TaylorN)
        @eval begin
            ($f)(a::$T{T}, b::$T{S}) where {T<:Number, S<:Number} =
                $f(promote(a,b)...)

            function $f(a::$T{T}, b::$T{T}) where {T<:Number}
                if a.order != b.order
                    a, b = fixorder(a, b)
                end
                v = similar(a.coeffs)
                @__dot__ v = $f(a.coeffs, b.coeffs)
                return $T(v, a.order)
            end

            function $f(a::$T)
                v = similar(a.coeffs)
                @__dot__ v = $f(a.coeffs)
                return $T(v, a.order)
            end

            ($f)(a::$T{T}, b::S) where {T<:Number, S<:Number} =
                $f(promote(a,b)...)

            function $f(a::$T{T}, b::T) where {T<:Number}
                coeffs = copy(a.coeffs)
                @inbounds coeffs[1] = $f(a[0], b)
                return $T(coeffs, a.order)
            end

            ($f)(b::S, a::$T{T}) where {T<:Number, S<:Number} =
                $f(promote(b,a)...)

            function $f(b::T, a::$T{T}) where {T<:Number}
                coeffs = similar(a.coeffs)
                @__dot__ coeffs = ($f)(a.coeffs)
                @inbounds coeffs[1] = $f(b, a[0])
                return $T(coeffs, a.order)
            end

            ## add! and subst! ##
            function ($fc)(v::$T{T}, a::$T{T}, k::Int) where {T}
                @inbounds v[k] = ($f)(a[k])
                return nothing
            end
            function ($fc)(v::$T, a::NumberNotSeries, k::Int)
                @inbounds v[k] = k==0 ? ($f)(zero(v[0]),a) : zero(v[k])
                return nothing
            end
            function ($fc)(v::$T, a::$T, b::$T, k::Int)
                @inbounds v[k] = ($f)(a[k], b[k])
                return nothing
            end
            function ($fc)(v::$T, a::$T, b::NumberNotSeries, k::Int)
                @inbounds v[k] = k==0 ? ($f)(a[0], b) : ($f)(a[k], zero(b))
                return nothing
            end
            function ($fc)(v::$T, a::NumberNotSeries, b::$T, k::Int)
                @inbounds v[k] = k==0 ? ($f)(a, b[0]) : ($f)(zero(a), b[k])
                return nothing
            end
        end
    end

    @eval begin
        ($f)(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{S}) where
            {T<:NumberNotSeriesN,S<:NumberNotSeriesN} = $f(promote(a,b)...)

        function $f(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{T}) where
                T<:NumberNotSeriesN
            @assert a.order == b.order
            v = similar(a.coeffs)
            @__dot__ v = $f(a.coeffs, b.coeffs)
            return HomogeneousPolynomial(v, a.order)
        end

        function $f(a::HomogeneousPolynomial)
            v = similar(a.coeffs)
            @__dot__ v = $f(a.coeffs)
            return HomogeneousPolynomial(v, a.order)
        end

        function ($f)(a::TaylorN{Taylor1{T}}, b::Taylor1{S}) where
                {T<:NumberNotSeries, S<:NumberNotSeries}
            @inbounds aux = $f(a[0][1], b)
            R = eltype(aux)
            coeffs = Array{HomogeneousPolynomial{Taylor1{R}}}(undef, a.order+1)
            coeffs .= a.coeffs
            @inbounds coeffs[1] = aux
            return TaylorN(coeffs, a.order)
        end

        function ($f)(b::Taylor1{S}, a::TaylorN{Taylor1{T}}) where
                {T<:NumberNotSeries,S<:NumberNotSeries}
            @inbounds aux = $f(b, a[0][1])
            R = eltype(aux)
            coeffs = Array{HomogeneousPolynomial{Taylor1{R}}}(undef, a.order+1)
            @__dot__ coeffs = $f(a.coeffs)
            @inbounds coeffs[1] = aux
            return TaylorN(coeffs, a.order)
        end

        function ($f)(a::Taylor1{TaylorN{T}}, b::TaylorN{S}) where
                {T<:NumberNotSeries,S<:NumberNotSeries}
            @inbounds aux = $f(a[0], b)
            R = eltype(aux)
            coeffs = Array{TaylorN{R}}(undef, a.order+1)
            coeffs .= a.coeffs
            @inbounds coeffs[1] = aux
            return Taylor1(coeffs, a.order)
        end

        function ($f)(b::TaylorN{S}, a::Taylor1{TaylorN{T}}) where
                {T<:NumberNotSeries,S<:NumberNotSeries}
            @inbounds aux = $f(b, a[0])
            R = eltype(aux)
            coeffs = Array{TaylorN{R}}(undef, a.order+1)
            @__dot__ coeffs = $f(a.coeffs)
            @inbounds coeffs[1] = aux
            return Taylor1(coeffs, a.order)
        end
    end
end

+(a::Taylor1{T}, b::TaylorN{S}) where {T<:NumberNotSeries,S<:NumberNotSeries} =
    +(promote(a,b)...)
-(a::Taylor1{T}, b::TaylorN{S}) where {T<:NumberNotSeries,S<:NumberNotSeries} =
    -(promote(a,b)...)

## Multiplication ##
for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)

    @eval begin
        *(a::Bool, b::$T) = *(convert(Int, a), b)

        *(a::$T, b::Bool) = b * a

        function *(a::T, b::$T) where {T<:NumberNotSeries}
            @inbounds aux = a * b.coeffs[1]
            v = Array{typeof(aux)}(undef, length(b.coeffs))
            @__dot__ v = a * b.coeffs
            return $T(v, b.order)
        end

        *(b::$T, a::T) where {T<:NumberNotSeries} = a * b
    end
end

for T in (:HomogeneousPolynomial, :TaylorN)

    @eval begin
        function *(a::Taylor1{T}, b::$T{Taylor1{S}}) where
                {T<:NumberNotSeries, S<:NumberNotSeries}

            @inbounds aux = a * b.coeffs[1]
            R = typeof(aux)
            coeffs = Array{R}(undef, length(b.coeffs))
            @__dot__ coeffs = a * b.coeffs
            return $T(coeffs, b.order)
        end

        *(b::$T{Taylor1{R}}, a::Taylor1{T}) where
            {T<:NumberNotSeries,R<:NumberNotSeries} = a * b

        function *(a::$T{T}, b::Taylor1{$T{S}}) where {T<:NumberNotSeries, S<:NumberNotSeries}
            @inbounds aux = a * b[0]
            R = typeof(aux)
            coeffs = Array{R}(undef, length(b.coeffs))
            @__dot__ coeffs = a * b.coeffs
            return Taylor1(coeffs, b.order)
        end

        *(b::Taylor1{$T{S}}, a::$T{T}) where
            {T<:NumberNotSeries,S<:NumberNotSeries} = a * b
    end
end

for (T, W) in ((:Taylor1, :Number), (:TaylorN, :NumberNotSeriesN))
    @eval *(a::$T{T}, b::$T{S}) where {T<:$W, S<:$W} = *(promote(a,b)...)

    @eval function *(a::$T{T}, b::$T{T}) where {T<:$W}
        if a.order != b.order
            a, b = fixorder(a, b)
        end
        c = $T(zero(a[0]), a.order)
        for ord in eachindex(c)
            mul!(c, a, b, ord) # updates c[ord]
        end
        return c
    end
end

*(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{S}) where
    {T<:NumberNotSeriesN,S<:NumberNotSeriesN} = *(promote(a,b)...)

function *(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{T}) where
        {T<:NumberNotSeriesN}

    order = a.order + b.order

    order > get_order() && return HomogeneousPolynomial(zero(T), get_order())

    res = HomogeneousPolynomial(zero(T), order)
    mul!(res, a, b)
    return res
end

# Internal multiplication functions
for T in (:Taylor1, :TaylorN)
    @eval @inline function mul!(c::$T{T}, a::$T{T}, b::$T{T}, k::Int) where {T}

        if $T == Taylor1
            @inbounds c[k] = a[0] * b[k]
        else
            @inbounds mul!(c[k], a[0], b[k])
        end
        @inbounds for i = 1:k
            if $T == Taylor1
                c[k] += a[i] * b[k-i]
            else
                mul!(c[k], a[i], b[k-i])
            end
        end

        return nothing
    end

    @eval @inline function mul!(v::$T, a::$T, b::NumberNotSeries, k::Int)
        @inbounds v[k] = a[k] * b
        return nothing
    end
    @eval @inline function mul!(v::$T, a::NumberNotSeries, b::$T, k::Int)
        @inbounds v[k] = a * b[k]
        return nothing
    end
end

@doc doc"""
    mul!(c, a, b, k::Int) --> nothing

Update the `k`-th expansion coefficient `c[k]` of `c = a * b`,
where all `c`, `a`, and `b` are either `Taylor1` or `TaylorN`.

The coefficients are given by

```math
c_k = \sum_{j=0}^k a_j b_{k-j}.
```

""" mul!

"""
    mul!(c, a, b) --> nothing

Return `c = a*b` with no allocation; all arguments are `HomogeneousPolynomial`.

"""
@inline function mul!(c::HomogeneousPolynomial, a::HomogeneousPolynomial,
        b::HomogeneousPolynomial)

    (iszero(b) || iszero(a)) && return nothing

    @inbounds num_coeffs_a = size_table[a.order+1]
    @inbounds num_coeffs_b = size_table[b.order+1]

    @inbounds posTb = pos_table[c.order+1]

    @inbounds indTa = index_table[a.order+1]
    @inbounds indTb = index_table[b.order+1]

    @inbounds for na in 1:num_coeffs_a
        ca = a[na]
        # iszero(ca) && continue
        inda = indTa[na]

        @inbounds for nb in 1:num_coeffs_b
            cb = b[nb]
            # iszero(cb) && continue
            indb = indTb[nb]

            pos = posTb[inda + indb]
            c[pos] += ca * cb
        end
    end

    return nothing
end

## Division ##
function /(a::Taylor1{Rational{T}}, b::S) where {T<:Integer, S<:NumberNotSeries}
    R = typeof( a[0] // b)
    v = Array{R}(undef, a.order+1)
    @__dot__ v = a.coeffs // b
    return Taylor1(v, a.order)
end

for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)

    @eval function /(a::$T{T}, b::S) where {T<:NumberNotSeries, S<:NumberNotSeries}
        R = promote_type(T,S)
        return convert($T{R}, a) * inv(convert(R, b))
    end

    @eval /(a::$T, b::T) where {T<:NumberNotSeries} = a * inv(b)
end

for T in (:HomogeneousPolynomial, :TaylorN)
    @eval function /(b::$T{Taylor1{S}}, a::Taylor1{T}) where
            {T<:NumberNotSeries,S<:NumberNotSeries}
        @inbounds aux = b.coeffs[1] / a
        R = typeof(aux)
        coeffs = Array{R}(undef, length(b.coeffs))
        @__dot__ coeffs = b.coeffs / a
        return $T(coeffs, b.order)
    end

    @eval function /(b::Taylor1{$T{S}}, a::$T{T}) where
            {T<:NumberNotSeries,S<:NumberNotSeries}
        @inbounds aux = b[0] / a
        R = typeof(aux)
        coeffs = Array{R}(undef, length(b.coeffs))
        @__dot__ coeffs = b.coeffs / a
        return Taylor1(coeffs, b.order)
    end
end

/(a::Taylor1{T}, b::Taylor1{S}) where {T<:Number, S<:Number} = /(promote(a,b)...)

function /(a::Taylor1{T}, b::Taylor1{T}) where {T<:Number}
    if a.order != b.order
        a, b = fixorder(a, b)
    end

    # order and coefficient of first factorized term
    ordfact, cdivfact = divfactorization(a, b)

    c = Taylor1(cdivfact, a.order)
    for ord in eachindex(c)
        div!(c, a, b, ord) # updates c[ord]
    end

    return c
end

/(a::TaylorN{T}, b::TaylorN{S}) where
    {T<:NumberNotSeriesN, S<:NumberNotSeriesN} = /(promote(a,b)...)

function /(a::TaylorN{T}, b::TaylorN{T}) where {T<:NumberNotSeriesN}
    @assert !iszero(constant_term(b))

    if a.order != b.order
        a, b = fixorder(a, b)
    end

    # first coefficient
    @inbounds cdivfact = a[0] / constant_term(b)
    c = TaylorN(cdivfact, a.order)
    for ord in eachindex(c)
        div!(c, a, b, ord) # updates c[ord]
    end

    return c
end

@inline function divfactorization(a1::Taylor1, b1::Taylor1)
    # order of first factorized term; a1 and b1 assumed to be of the same order
    a1nz = findfirst(a1)
    b1nz = findfirst(b1)
    a1nz = a1nz ≥ 0 ? a1nz : a1.order
    b1nz = b1nz ≥ 0 ? b1nz : a1.order
    ordfact = min(a1nz, b1nz)
    cdivfact = a1[ordfact] / b1[ordfact]

    # Is the polynomial factorizable?
    iszero(b1[ordfact]) && throw( ArgumentError(
        """Division does not define a Taylor1 polynomial;
        order k=$(ordfact) => coeff[$(ordfact)]=$(cdivfact).""") )

    return ordfact, cdivfact
end

## TODO: Implement factorization (divfactorization) for TaylorN polynomials

# Homogeneous coefficient for the division
@doc doc"""
    div!(c, a, b, k::Int)

Compute the `k-th` expansion coefficient `c[k]` of `c = a / b`,
where all `c`, `a` and `b` are either `Taylor1` or `TaylorN`.

The coefficients are given by

```math
c_k =  \frac{1}{b_0} \big(a_k - \sum_{j=0}^{k-1} c_j b_{k-j}\big).
```

For `Taylor1` polynomials, a similar formula is implemented which
exploits `k_0`, the order of the first non-zero coefficient of `a`.
""" div!

@inline function div!(c::Taylor1, a::Taylor1, b::Taylor1, k::Int)

    # order and coefficient of first factorized term
    ordfact, cdivfact = divfactorization(a, b)
    if k == 0
        @inbounds c[0] = cdivfact
        return nothing
    end

    imin = max(0, k+ordfact-b.order)
    @inbounds c[k] = c[imin] * b[k+ordfact-imin]
    @inbounds for i = imin+1:k-1
        c[k] += c[i] * b[k+ordfact-i]
    end
    if k+ordfact ≤ b.order
        @inbounds c[k] = (a[k+ordfact]-c[k]) / b[ordfact]
    else
        @inbounds c[k] = - c[k] / b[ordfact]
    end
    return nothing
end

@inline function div!(v::Taylor1, a::Taylor1, b::NumberNotSeries, k::Int)
    @inbounds v[k] = a[k] / b
    return nothing
end

div!(v::Taylor1, b::NumberNotSeries, a::Taylor1, k::Int) =
    div!(v::Taylor1, Taylor1(b, a.order), a, k)

@inline function div!(c::TaylorN, a::TaylorN, b::TaylorN, k::Int)
    if k==0
        @inbounds c[0] = a[0] / constant_term(b)
        return nothing
    end

    @inbounds for i = 0:k-1
        mul!(c[k], c[i], b[k-i])
    end
    @inbounds c[k] = (a[k] - c[k]) / constant_term(b)
    return nothing
end

"""
    mul!(Y, A, B)

Multiply A*B and save the result in Y.
"""
function mul!(y::Vector{Taylor1{T}},
        a::Union{Matrix{T},SparseMatrixCSC{T}},
        b::Vector{Taylor1{T}}) where {T<:Number}

    n, k = size(a)
    @assert (length(y)== n && length(b)== k)

    # determine the maximal order of b
    # order = maximum([b1.order for b1 in b])
    order = maximum(get_order.(b))

    # Use matrices of coefficients (of proper size) and mul!
    # B = zeros(T, k, order+1)
    B = Array{T}(undef, k, order+1)
    B = zero.(B)
    for i = 1:k
        @inbounds ord = b[i].order
        @inbounds for j = 1:ord+1
            B[i,j] = b[i][j-1]
        end
    end
    Y = Array{T}(undef, n, order+1)
    mul!(Y, a, B)
    @inbounds for i = 1:n
        y[i] = Taylor1( collect(Y[i,:]), order)
    end

    return y
end

# Adapted from (Julia v1.2) stdlib/v1.2/LinearAlgebra/src/dense.jl#721-734,
# licensed under MIT "Expat".
# Specialize a method of `inv` for Matrix{Taylor1{T}}. Simply, avoid pivoting,
# since the polynomial field is not an ordered one.
# function Base.inv(A::StridedMatrix{Taylor1{T}}) where T
#     checksquare(A)
#     S = Taylor1{typeof((one(T)*zero(T) + one(T)*zero(T))/one(T))}
#     AA = convert(AbstractArray{S}, A)
#     if istriu(AA)
#         Ai = triu!(parent(inv(UpperTriangular(AA))))
#     elseif istril(AA)
#         Ai = tril!(parent(inv(LowerTriangular(AA))))
#     else
#         # Do not use pivoting !!
#         Ai = inv!(lu(AA, Val(false)))
#         Ai = convert(typeof(parent(Ai)), Ai)
#     end
#     return Ai
# end

# Adapted from (Julia v1.2) stdlib/v1.2/LinearAlgebra/src/lu.jl#240-253
# and (Julia v1.4.0-dev) stdlib/LinearAlgebra/v1.4/src/lu.jl#270-274,
# licensed under MIT "Expat".
# Specialize a method of `lu` for Matrix{Taylor1{T}}, which avoids pivoting,
# since the polynomial field is not an ordered one.
# We can't assume an ordered field so we first try without pivoting
function lu(A::AbstractMatrix{Taylor1{T}}; check::Bool = true) where T
    S = Taylor1{lutype(T)}
    F = lu!(copy_oftype(A, S), Val(false); check = false)
    if issuccess(F)
        return F
    else
        return lu!(copy_oftype(A, S), Val(true); check = check)
    end
end

ext: jl


text:
@testset "interpolate" begin
    @testset "basic" begin
        @system SInterpolate(Controller) begin
            m => ([1 => 10, 2 => 20, 3 => 30]) ~ interpolate
            n(m) ~ interpolate(reverse)
            a(m) => m(2.5) ~ track
            b(n) => n(25) ~ track
       end
       s = instance(SInterpolate)
       @test s.a' == 25
       @test s.b' == 2.5
    end

    @testset "matrix" begin
        @system SInterpolateMatrix(Controller) begin
            m => ([1 10; 2 20; 3 30]) ~ interpolate
            n(m) ~ interpolate(reverse)
            a(m) => m(2.5) ~ track
            b(n) => n(25) ~ track
       end
       s = instance(SInterpolateMatrix)
       @test s.a' == 25
       @test s.b' == 2.5
    end

    @testset "config" begin
        @system SInterpolateConfig(Controller) begin
            m => ([1 => 10]) ~ interpolate(parameter)
            n(m) ~ interpolate(reverse)
            a(m) => m(2.5) ~ track
            b(n) => n(25) ~ track
       end
       o = SInterpolateConfig => :m => [1 => 10, 2 => 20, 3 => 30]
       s = instance(SInterpolateConfig; config=o)
       @test s.a' == 25
       @test s.b' == 2.5
    end

    @testset "unit" begin
        @system SInterpolateUnit(Controller) begin
            m => ([1 => 10, 2 => 20, 3 => 30]) ~ interpolate(u"s", knotunit=u"m")
            n(m) ~ interpolate(u"m", reverse)
            a(m) => m(2.5u"m") ~ track(u"s")
            b(n) => n(25u"s") ~ track(u"m")
       end
       s = instance(SInterpolateUnit)
       @test s.a' == 25u"s"
       @test s.b' == 2.5u"m"
    end
end


text:
@doc raw"
    einsum(::EinCode{ixs, iy}, xs, size_dict) where {ixs, iy}

return the tensor that results from contracting the tensors `xs` according
to their indices `ixs`, where all indices that do not appear in the output `iy` are
summed over.
The result is permuted according to `out`.

- `ixs` - tuple of tuples of index-labels of the input-tensors `xs`

- `iy` - tuple of index-labels of the output-tensor

- `xs` - tuple of tensors

- `size_dict` - `IndexSize`-object that maps index-labels to their sizes

# example

```jldoctest; setup = :(using OMEinsum)
julia> a, b = rand(2,2), rand(2,2);

julia> einsum(EinCode((('i','j'),('j','k')),('i','k')), (a, b)) ≈ a * b
true

julia> einsum(EinCode((('i','j'),('j','k')),('k','i')), (a, b)) ≈ permutedims(a * b, (2,1))
true
```
"
@generated function einsum(code::EinCode{ixs, iy}, xs, size_dict) where {ixs, iy}
    rule = match_rule(ixs, iy)
    :(einsum($rule, code, xs, size_dict))
end

function einsum(::Tr, ::EinCode, xs, size_dict)
    @debug "Tr" size.(xs)
    asarray(tr(xs[1]), xs[1])
end

using TensorOperations

function einsum(::PTrace, ::EinCode{ixs,iy}, xs::NTuple{<:Any, AbstractArray{<:BlasFloat}}, size_dict) where {ixs, iy}
    @debug "PTrace tensortrace" ixs => iy size.(xs)
    asarray(tensortrace(xs[1], ixs[1], iy), xs[1])
end

# note that dispatching to Hadamard if some `ixs` are permuted has inferior
# performance compared to the fallback
function einsum(::Hadamard, ::EinCode{ixs, iy}, xs, size_dict) where {ixs, iy}
    @debug "Hadamard" ixs => iy size.(xs)
    asarray(broadcast(*, xs...), xs[1])
end

function einsum(::PairWise, ::EinCode{ixs, iy}, xs::NTuple{NT,AbstractArray}, size_dict) where {ixs, iy, NT}
    @debug "PairWise optcontract" ixs => iy size.(xs)
    optcontract(ixs, xs, iy)
end

function einsum(sm::Sum, code::EinCode{ixs, iy}, xs, size_dict) where {ixs, iy}
    dims = (findall(i -> i ∉ iy, ixs[1])...,)
    (ix1,) = ixs
    ix1f = filter!(i -> i in iy, collect(ix1))
    perm = map(i -> findfirst(==(i), ix1f), iy)
    res = dropdims(sum(xs[1], dims=dims), dims=dims)
    if perm == iy
        @debug "Sum" ixs => iy size.(xs)
        res
    else
        @debug "Sum permutedims" ixs => iy size.(xs) perm
        tensorpermute(res, perm)
    end
end

@generated function einsum(sm::MatMul, code::EinCode{ixs, iy}, xs, size_dict) where {ixs, iy}
    ix1, ix2 = ixs
    l = ifelse(ix1[1] in ix2, ix1[1], ix1[2])
    if ix1[2] == l && ix2[1] == l
        if iy == (ix1[1], ix2[2])
            return :(@debug "MatMul1 (ij,jk -> ik)" ixs => iy size.(xs);
                xs[1] * xs[2]
            )
        else
            return :(@debug "MatMul2 (ij,jk -> ki) permutedims(ik)" ixs => iy size.(xs);
                permutedims(xs[1] * xs[2])
            )
        end
    elseif ix1[1] == l && ix2[1] == l
        if iy == (ix1[2], ix2[2])
            return :(@debug "MatMul3 (ji,jk -> ik) transpose(ji)" ixs => iy size.(xs);
                transpose(xs[1]) * xs[2]
            )
        else
            return :(@debug "MatMul4 (ji,jk -> ki) transpose(jk)" ixs => iy size.(xs);
                transpose(xs[2]) * xs[1]
            )
        end
    elseif ix1[2] == l && ix2[2] == l
        if iy == (ix1[1], ix2[1])
            return :(@debug "MatMul5 (ij,kj -> ik) transpose(kj)" ixs => iy size.(xs);
                xs[1] * transpose(xs[2])
            )
        else
            return :(@debug "MatMul6 (ij,kj -> ki) transpose(ij)" ixs => iy size.(xs);
                xs[2] * transpose(xs[1])
            )
        end
    else #ix1[1] == l && ix2[2] == l
        if iy == (ix1[2], ix2[1])
            return :(@debug "MatMul7 (ji,kj -> ik) permutedims(ki)" ixs => iy size.(xs);
                permutedims(xs[2] * xs[1])
            )
        else
            return :(@debug "MatMul8 (ji,kj -> ki)" ixs => iy size.(xs);
                xs[2] * xs[1]
            )
        end
    end
end

function einsum(::Permutedims, code::EinCode{ixs, iy}, xs, size_dict) where {ixs, iy}
    (ix,) = ixs
    (x,) = xs
    perm = map(i -> findfirst(==(i), ix), iy)
    @debug "Permutedims" ix => iy size(xs[1]) perm
    return tensorpermute(x, perm)
end

function einsum(::Identity, ::EinCode{ixs, iy}, xs, size_dict) where {ixs, iy}
    xs[1]
end

# the fallback
function einsum(::DefaultRule, code::EinCode{ixs, iy}, xs, size_dict) where {ixs, iy}
    @debug "DefaultRule loop_einsum" ixs => iy size.(xs)
    loop_einsum(code, xs, size_dict)
end

function einsum(::PTrace, code::EinCode{ixs, iy}, xs::NTuple{NT, Any}, size_dict) where {ixs, iy, NT}
    @debug "PTrace loop_einsum" ixs => iy size.(xs)
    loop_einsum(code, xs, size_dict)
end

function einsum(::BatchedContract, code::EinCode{ixs, iy}, xs::NTuple{NT, Any}, size_dict) where {ixs, iy, NT}
    @debug "BatchedContract loop_einsum" ixs => iy size.(xs)
    loop_einsum(code, xs, size_dict)
end

@generated function _preprocess_dupindices(::Val{ix}, x) where {ix}
    if length(tunique(ix)) != length(ix)
        iy = [l for l in ix if count(==(l), ix) == 1]
        :(($(Val((iy...,))), einsum($(EinCode((ix,), (iy...,))), (x,), get_size_dict(($ix,), (x,)))))
    else
        :(($(Val(ix)), x))
    end
end

@generated function einsum(::BatchedContract, ::EinCode{ixs,iy}, xs::NTuple{<:Any, AbstractArray{<:BlasFloat}}, size_dict) where {ixs, iy}
    quote
        ixs1, xs1 = _preprocess_dupindices($(Val(ixs[1])), xs[1])
        ixs2, xs2 = _preprocess_dupindices($(Val(ixs[2])), xs[2])
        @debug "BatchedContract" ixs => iy Tuple(ixs1) Tuple(ixs2) size.((xs1, xs2))
        batched_contract(ixs1, xs1, ixs2, xs2, $(Val(iy)))
    end
end


text:
# 1D diffusion problem

# Packages and inclusions
using ModelingToolkit,DiffEqOperators,LinearAlgebra,Test,OrdinaryDiffEq, DomainSets
using ModelingToolkit: Differential

# Tests
@testset "Dt(u(t,x)) ~ Dxx(u(t,x)), 0 ~ Dxx(v(t,x)) + sin(x), Dirichlet BCs" begin
    # Method of Manufactured Solutions
    u_exact = (x,t) -> exp.(-t) * cos.(x)
    v_exact = (x,t) -> exp.(-t) * sin.(x)

    # Parameters, variables, and derivatives
    @parameters t x
    @variables u(..) v(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dxx = Dx^2

    # 1D PDE and boundary conditions
    eqs = [Dt(u(t,x)) ~ Dxx(u(t,x)),
           0 ~ Dxx(v(t,x)) + exp(-t)*sin(x)]
    bcs = [u(0,x) ~ cos(x),
           v(0,x) ~ sin(x),
           u(t,0) ~ exp(-t),
           Dx(u(t,1)) ~ -exp(-t) * sin(1),
           Dx(v(t,0)) ~ exp(-t),
           v(t,1) ~ exp(-t) * sin(1)]

    # Space and time domains
    domains = [t ∈ Interval(0.0,1.0),
               x ∈ Interval(0.0,1.0)]

    # PDE system
    @named pdesys = PDESystem(eqs,bcs,domains,[t,x],[u(t,x),v(t,x)])

    # Method of lines discretization
    l = 20
    dx = range(0.0,1.0,length=l)
    order = 2
    discretization = MOLFiniteDifference([x=>dx],t)

    # Convert the PDE problem into an ODE problem
    prob = discretize(pdesys,discretization)

    # Solve ODE problem
    sol = solve(prob,Rodas4(),saveat=0.1)

    x_sol = dx
    t_sol = sol.t

    @variables u[1:l](..) v[1:l](..)

    # Test against exact solution
    for i in 1:l
        @test all(isapprox.(u_exact(x_sol[i], t_sol), sol[u[i](t)], atol=0.01))
        @test all(isapprox.(v_exact(x_sol[i], t_sol), sol[v[i](t)], atol=0.01))
    end
end


text:
using DiffPrecTest
using Statistics, StatsBase, LinearAlgebra
using SparseArrays
using ProximalBase, CoordinateDescent, CovSel
using Random, Distributions
using JLD

@show gethostname()

pArr = [25, 100, 200]
elemArr = [(5,5), (8, 7), (50, 25)]
n = 300
est      = Array{Any}(undef, 9)   # number of methods

# rep   = parse(Int,ARGS[1])
# ip    = parse(Int,ARGS[2])
# iElem = parse(Int,ARGS[3])
# dir = ARGS[4]

rep   = 1
ip    = 1
iElem = 1
dir = "."

Random.seed!(134)

p = pArr[ip]
# generate model
Ω = Matrix{Float64}(I, p, p)
ρ = 0.9
for k=1:p-1
    for l=1:p-k
        Ω[l  , l+k] = ρ^k
        Ω[l+k, l  ] = ρ^k
    end
end
d = Vector{Float64}(undef, p)
rand!(Uniform(0.5, 2.5), d)
d .= sqrt.(d)
D = Diagonal(d)
Σ = inv(Symmetric(D * Ω * D))

dist_X = MvNormal(convert(Matrix, Σ))
dist_Y = MvNormal(convert(Matrix, Σ))

# generate data
Random.seed!(134 + rep)
X = rand(dist_X, n)'
Y = rand(dist_Y, n)'

Sx = Symmetric( X'X / n )
Sy = Symmetric( Y'Y / n )

@show (ip, iElem)
supportOracle =  Vector{Int64}()
ri, ci = elemArr[iElem]

@time est[1], eΔ, esuppΔ, eω, esuppω, eSupport = DiffPrecTest.estimate(SymmetricNormal(), X, Y, ri, ci; Sx=Sx, Sy=Sy)
@time est[2], _, _, _, _, _ = DiffPrecTest.estimate(AsymmetricNormal(), X, Y, ri, ci; Sx=Sx, Sy=Sy, Δ=eΔ, suppΔ=esuppΔ)
@time est[3] = DiffPrecTest.estimate(SeparateNormal(), X, Y, ri, ci; Sx=Sx, Sy=Sy)
@time est[4] = DiffPrecTest.estimate(SymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, eSupport)
@time est[5] = DiffPrecTest.estimate(AsymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, eSupport)

@time est[6] = DiffPrecTest.estimate(SymmetricOracleNormal(), Sx, Sy, X, Y, ri, ci, supportOracle)
@time est[7] = DiffPrecTest.estimate(AsymmetricOracleNormal(), Sx, Sy, X, Y, ri, ci, supportOracle)
@time est[8] = DiffPrecTest.estimate(SymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, supportOracle)
@time est[9] = DiffPrecTest.estimate(AsymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, supportOracle)

@save "$(dir)/res_$(ip)_$(iElem)_$(rep).jld" est


text:
"""
# module LUP

# Examples

```jldoctest
julia>
```
"""
module LUP

end


text:
@doc doc"""
    default_estimator(model, updatemethod!)

Determines estimator to be used when `param["Estimator"]` is not set.
"""
default_estimator(model::Model, update) = simple_estimator
default_estimator(model::QuantumXXZ, update) = improved_estimator
default_estimator(model::Union{Ising, Potts}, update) = ifelse(update==SW_update!, improved_estimator, simple_estimator)

function simple_estimator(model, param::Parameter, _=nothing)
    p = convert_parameter(model,param)
    return simple_estimator(model, p..., nothing)
end

@doc """
    simple_estimator(model::Ising, T::Real, Js::AbstractArray)
    simple_estimator(model::Potts, T::Real, Js::AbstractArray)

Returns the following observables as `Dict{String, Any}`

# Observables
- `"Energy"`
    - Energy per spin (site)
- `"Energy^2"`
- `"Magnetization"`
    - Total magnetization per spin (order paremeter)
- `"|Magnetization|"`
- `"Magnetization^2"`
- `"Magnetization^4"`
"""
function simple_estimator(model::Ising, T::Real, Js::AbstractArray, _=nothing)
    nsites = numsites(model)
    nbonds = numbonds(model)

    M = mean(model.spins)
    E = 0.0
    @inbounds for b in bonds(model)
        s1, s2 = source(b), target(b)
        E += ifelse(model.spins[s1] == model.spins[s2], -1.0, 1.0) * Js[bondtype(b)]
    end
    E /= nsites

    res = Measurement()
    res["Magnetization"] = M
    res["|Magnetization|"] = abs(M)
    res["Magnetization^2"] = M^2
    res["Magnetization^4"] = M^4
    res["Energy"] = E
    res["Energy^2"] = E^2
    return res
end

function simple_estimator(model::Potts, T::Real, Js::AbstractArray, _=nothing)
    nsites = numsites(model)
    nbonds = numbonds(model)
    invQ = 1.0/model.Q

    M = 0.0
    @inbounds for s in 1:nsites
        M += ifelse(model.spins[s]==1, 1.0-invQ, -invQ)
    end
    M /= nsites
    E = 0.0
    @inbounds for b in bonds(model)
        s1, s2 = source(b), target(b)
        E -= ifelse(model.spins[s1] == model.spins[s2], 1.0, 0.0) * Js[bondtype(b)]
    end
    E /= nsites

    res = Measurement()
    res["Magnetization"] = M
    res["|Magnetization|"] = abs(M)
    res["Magnetization^2"] = M^2
    res["Magnetization^4"] = M^4
    res["Energy"] = E
    res["Energy^2"] = E^2
    return res
end

@doc """
    simple_estimator(model::Clock, T::Real, Js::AbstractArray)
    simple_estimator(model::XY, T::Real, Js::AbstractArray)

Returns the following observables as `Dict{String, Any}`

# Observables
- `"Energy"`
    - Energy per spin (site)
- `"Energy^2"`
- `"|Magnetization|"`
    - Absolute value of total magnetization per spin (order paremeter)
- `"|Magnetization|^2"`
- `"|Magnetization|^4"`
- `"Magnetization x"`
    - x component of total magnetization per spin (order paremeter)
- `"|Magnetization x|"`
- `"Magnetization x^2"`
- `"Magnetization x^4"`
- `"Magnetization y"`
    - y component of total magnetization per spin (order paremeter)
- `"|Magnetization y|"`
- `"Magnetization y^2"`
- `"Magnetization y^4"`
- `"Helicity Modulus x"`
- `"Helicity Modulus y"`
"""
function simple_estimator(model::Clock, T::Real, Js::AbstractArray, _=nothing)
    nsites = numsites(model)
    nbonds = numbonds(model)
    D = dim(model)
    invN = 1.0/nsites
    beta = 1.0/T

    M = zeros(2)
    E = 0.0
    U1 = zeros(2)
    U2 = zeros(2)

    @inbounds for s in 1:nsites
        M[1] += model.cosines[model.spins[s]]
        M[2] += model.sines[model.spins[s]]
    end

    @inbounds for b in bonds(model)
        i, j = source(b), target(b)
        dir = bonddirection(b)
        dt = mod1(model.spins[j] - model.spins[i], model.Q)
        E -= model.cosines[dt] * Js[bondtype(b)]

        for d in 1:D
            U1[d] += model.cosines[dt] * dir[d]^2
            U2[d] += model.sines[dt] * dir[d]
        end
    end
    for d in 1:2
        M[d] *= invN
        U1[d] -= beta * U2[d]^2
        U1[d] *= invN
    end
    E /= nsites

    res = Measurement()
    for (i,c) in enumerate(["x","y"])
        res["Magnetization $c"] = M[i]
        res["|Magnetization $c|"] = abs(M[i])
        res["Magnetization $c^2"] = M[i]^2
        res["Magnetization $c^4"] = M[i]^4
        res["Helicity Modulus $c"] = U1[i]
    end
    M2 = sum(abs2, M)
    res["|Magnetization|"] = sqrt(M2)
    res["|Magnetization|^2"] = M2
    res["|Magnetization|^4"] = M2^2
    res["Energy"] = E
    res["Energy^2"] = E^2

    return res
end

function simple_estimator(model::XY, T::Real, Js::AbstractArray, _=nothing)
    nsites = numsites(model)
    nbonds = numbonds(model)
    D = dim(model)
    invN = 1.0/nsites
    beta = 1.0/T

    M = zeros(2)
    E = 0.0
    U1 = zeros(2)
    U2 = zeros(2)

    @inbounds for s in 1:nsites
        M[1] += cospi(2model.spins[s])
        M[2] += sinpi(2model.spins[s])
    end

    @inbounds for b in bonds(model)
        i, j = source(b), target(b)
        dir = bonddirection(b)
        dt = mod(model.spins[j] - model.spins[i] + 2.0, 1.0)
        dt = ifelse(0.5 <= dt, dt - 1.0, dt)
        E -= cospi(2dt) * Js[bondtype(b)]

        for d in 1:D
            U1[d] += cospi(2dt) * dir[d]^2
            U2[d] += sinpi(2dt) * dir[d]
        end
    end
    for d in 1:2
        M[d] *= invN
        U1[d] -= beta * U2[d]^2
        U1[d] *= invN
    end
    E /= nsites

    res = Measurement()
    for (i,c) in enumerate(["x","y"])
        res["Magnetization $c"] = M[i]
        res["|Magnetization $c|"] = abs(M[i])
        res["Magnetization $c^2"] = M[i]^2
        res["Magnetization $c^4"] = M[i]^4
        res["Helicity Modulus $c"] = U1[i]
    end
    M2 = sum(abs2, M)
    res["|Magnetization|"] = sqrt(M2)
    res["|Magnetization|^2"] = M2
    res["|Magnetization|^4"] = M2^2
    res["Energy"] = E
    res["Energy^2"] = E^2

    return res
end


text:
module SyntaxHighlighter

using Tokenize
using Tokenize.Tokens
import Tokenize.Tokens: Token, kind, exactkind, iskeyword, untokenize

using Crayons

import OhMyREPL: add_pass!, PASS_HANDLER

mutable struct ColorScheme
    symbol::Crayon
    comment::Crayon
    string::Crayon
    call::Crayon
    op::Crayon
    keyword::Crayon
    text::Crayon
    function_def::Crayon
    error::Crayon
    argdef::Crayon
    _macro::Crayon
    number::Crayon
end

symbol!(cs, crayon::Crayon) = cs.symbol = crayon
comment!(cs, crayon::Crayon) = cs.comment = crayon
string!(cs, crayon::Crayon) = cs.string = crayon
call!(cs, crayon::Crayon) = cs.call = crayon
op!(cs, crayon::Crayon) = cs.op = crayon
keyword!(cs, crayon::Crayon) = cs.keyword = crayon
text!(cs, crayon::Crayon) = cs.text = crayon
function_def!(cs, crayon::Crayon) = cs.function_def = crayon
error!(cs, crayon::Crayon) = cs.error = crayon
argdef!(cs, crayon::Crayon) = cs.argdef = crayon
macro!(cs, crayon::Crayon) = cs._macro = crayon
number!(cs, crayon::Crayon) = cs.number = crayon

function Base.show(io::IO, cs::ColorScheme)
    for n in fieldnames(ColorScheme)
        crayon = getfield(cs, n)
        print(io, crayon, "█ ", inv(crayon))
    end
    print(io, Crayon(foreground = :default))
end

ColorScheme() = ColorScheme([Crayon() for _ in 1:length(fieldnames(ColorScheme))]...)

include("colorschemes.jl")

mutable struct SyntaxHighlighterSettings
    active::ColorScheme
    schemes::Dict{String, ColorScheme}
end

function Base.show(io::IO, sh::SyntaxHighlighterSettings)
    first = true
    l = maximum(x->length(x), keys(sh.schemes))
    for (k, v) in sort(collect(sh.schemes), by = x -> x[1])
        if !first
            print(io, "\n\n")
        end
        first = false
        print(io, rpad(k, l+1, " "))
        print(io, v)
    end
    println(io)
end

function SyntaxHighlighterSettings()
    def = _create_juliadefault()
    d = Dict("JuliaDefault" => def)
    SyntaxHighlighterSettings(def, d)
end

SYNTAX_HIGHLIGHTER_SETTINGS = SyntaxHighlighterSettings()

add!(sh::SyntaxHighlighterSettings, name::String, scheme::ColorScheme) = sh.schemes[name] = scheme
add!(name::String, scheme::ColorScheme) = add!(SYNTAX_HIGHLIGHTER_SETTINGS, name, scheme)
activate!(sh::SyntaxHighlighterSettings, name::String) = sh.active = sh.schemes[name]

add!(SYNTAX_HIGHLIGHTER_SETTINGS, "Monokai24bit", _create_monokai_24())
add!(SYNTAX_HIGHLIGHTER_SETTINGS, "Monokai256", _create_monokai_256())
add!(SYNTAX_HIGHLIGHTER_SETTINGS, "Monokai16", _create_monokai())
add!(SYNTAX_HIGHLIGHTER_SETTINGS, "BoxyMonokai256", _create_boxymonokai_256())
add!(SYNTAX_HIGHLIGHTER_SETTINGS, "TomorrowNightBright", _create_tomorrow_night_bright())
add!(SYNTAX_HIGHLIGHTER_SETTINGS, "TomorrowNightBright24bit", _create_tomorrow_night_bright_24())
add!(SYNTAX_HIGHLIGHTER_SETTINGS, "Tomorrow24bit", _create_tomorrow_24())
add!(SYNTAX_HIGHLIGHTER_SETTINGS, "Tomorrow", _create_tomorrow_256())
add!(SYNTAX_HIGHLIGHTER_SETTINGS, "Distinguished", _create_distinguished())
add!(SYNTAX_HIGHLIGHTER_SETTINGS, "OneDark", _create_onedark())
add!(SYNTAX_HIGHLIGHTER_SETTINGS, "Base16MaterialDarker", _create_base16_material_darker())
# Added by default
# add!(SYNTAX_HIGHLIGHTER_SETTINGS, "JuliaDefault", _create_juliadefault())

if !Sys.iswindows()
    activate!(SYNTAX_HIGHLIGHTER_SETTINGS, "Monokai256")
else
    activate!(SYNTAX_HIGHLIGHTER_SETTINGS, "Monokai16")
end
add_pass!(PASS_HANDLER, "SyntaxHighlighter", SYNTAX_HIGHLIGHTER_SETTINGS, false)

function (highlighter::SyntaxHighlighterSettings)(crayons::Vector{Crayon}, tokens::Vector{Token}, ::Int)
    cscheme = highlighter.active
    prev_t = Tokens.Token()
    pprev_t = Tokens.Token()
    for (i, t) in enumerate(tokens)
        # a::x
        if exactkind(prev_t) == Tokens.DECLARATION
            crayons[i-1] = cscheme.argdef
            crayons[i] = cscheme.argdef
        # :foo
        elseif kind(t) == Tokens.IDENTIFIER && exactkind(prev_t) == Tokens.COLON &&
               kind(pprev_t) ∉ (Tokens.INTEGER, Tokens.FLOAT, Tokens.IDENTIFIER, Tokens.RPAREN)
            crayons[i-1] = cscheme.symbol
            crayons[i] = cscheme.symbol
        # function
        elseif iskeyword(kind(t))
            if kind(t) == Tokens.TRUE || kind(t) == Tokens.FALSE
                crayons[i] = cscheme.symbol
            else
                crayons[i] = cscheme.keyword
            end
        # "foo"
        elseif kind(t) == Tokens.STRING || kind(t) == Tokens.TRIPLE_STRING || kind(t) == Tokens.CHAR || kind(t) == Tokens.CMD || kind(t) == Tokens.TRIPLE_CMD
            crayons[i] = cscheme.string
        # * -
        elseif Tokens.isoperator(kind(t)) || exactkind(t) == Tokens.TRUE || exactkind(t) == Tokens.FALSE
            crayons[i] = cscheme.op
        # #= foo =#
        elseif kind(t) == Tokens.COMMENT
            crayons[i] = cscheme.comment
        # f(...)
        elseif kind(t) == Tokens.LPAREN
            if kind(prev_t) == Tokens.IDENTIFIER && !(i > 2 && exactkind(tokens[i-2]) == Tokens.AT_SIGN)
                crayons[i-1] = cscheme.call
            elseif exactkind(prev_t) == Tokens.DOT && kind(pprev_t) == Tokens.IDENTIFIER
                crayons[i-1] = cscheme.call
                crayons[i-2] = cscheme.call
            end
             # function f(...)
            if i > 3 && kind(tokens[i-2]) == Tokens.WHITESPACE && exactkind(tokens[i-3]) == Tokens.FUNCTION
                crayons[i-1] = cscheme.function_def
            end
        # @fdsafds
        elseif kind(t) == Tokens.IDENTIFIER && exactkind(prev_t) == Tokens.AT_SIGN
            crayons[i-1] = cscheme._macro
            crayons[i] = cscheme._macro
        # 2] = 32.32
        elseif kind(t) ∈ (Tokens.INTEGER, Tokens.BIN_INT, Tokens.OCT_INT, Tokens.HEX_INT, Tokens.FLOAT) || (kind(t) == Tokens.IDENTIFIER && untokenize(t) == "NaN")
            crayons[i] = cscheme.number
        elseif kind(t) == Tokens.WHITESPACE
            crayons[i] = Crayon()
        else
            crayons[i] = cscheme.text
        end
        pprev_t = prev_t
        prev_t = t
    end
    return
end

end


text:
using Plots
include("SplineModel.jl")

function nrbplot(nurbs::NURBS, subd::Array)#; light="off", cmap="summer")
    # light = lowercase(light)
    # if !isa(light, String)
    #     error("Parameter for light must be a string (on or off)")
    # end
    #
    # if !isa(cmap, String)
    #     error("Parameter for cmap must be a string")
    # end
    # if (light!="off") || (light=="on")
    #     shade = true
    # else
    #     shade = false
    # end
    subd .+=  1
    #plot the curve or surface
    if length(nurbs.knots)==2
        knt = nurbs.knots
        order = nurbs.order
        evalPtsU = LinRange(knt[1][order[1]],knt[1][end-order[1]+1], subd[1])
        evalPtsV = LinRange(knt[2][order[2]],knt[2][end-order[2]+1], subd[2])
        p = nrbeval(nurbs, [collect(evalPtsU), collect(evalPtsV)])
        plt = Plots.surface!(p[1,:,:], p[2,:,:], p[3,:,:], fc = :viridis)
    elseif length(nurbs.knots)==3
        bnd = nrbextract(nurbs)
        nrbplot(bnd[1], subd[2:3])
        nrbplot(bnd[2], subd[2:3])
        nrbplot(bnd[3], subd[[1,3]])
        nrbplot(bnd[4], subd[[1,3]])
        nrbplot(bnd[5], subd[1:2])
        plt = nrbplot(bnd[6], subd[1:2])
    elseif length(nurbs.knots)==1
        order = nurbs.order
        knt = nurbs.knots
        evalPts = LinRange(knt[1][order[1]],knt[1][end-order[1]+1], subd[1])
        p = nrbeval(nurbs, [collect(evalPts)])
        if any(nurbs.coefs[3,:].!=0)
            plt = Plots.plot!(p[1,:],p[2,:],p[3,:],lw=2, aspect_ratio=:equal)
        else
            plt = Plots.plot!(p[1,:],p[2,:],lw=2, aspect_ratio=:equal)
        end
        return plt
    end
end

function nrbkntplot(nurbs::NURBS, numDim::String="auto")
    nsub = 10

    if length(nurbs.knots)==2 #plot a NURBS surface
        plt = nrbplot(nurbs, [nsub, nsub])

        #plot the knots
        knt1 = unique(nurbs.knots[1][nurbs.order[1]:end-nurbs.order[1]+1])
        knt2 = unique(nurbs.knots[2][nurbs.order[2]:end-nurbs.order[2]+1])
        p1 = nrbeval(nurbs, [knt1, collect(LinRange(knt2[1], knt2[end], nsub))])
        p2 = nrbeval(nurbs, [collect(LinRange(knt1[1], knt1[end], nsub)), knt2])
        if numDim=="auto" && any(nurbs.coefs[3,:,:].!=0)
            numDim="3"
        elseif numDim=="auto" && !any(nurbs.coefs[3,:,:].!=0)
            numDim="2"
        end
        if numDim=="3"
            #surface in a 3D space
            for ii = 1:length(knt1)
                plt = Plots.plot!(p1[1,ii,:], p1[2,ii,:], p1[3,ii,:], linewidth = 2,
                         color=:black, legend=false)
            end
            for ii = 1:length(knt2)
                plt = Plots.plot!(p2[1,:,ii], p2[2,:,ii], p2[3,:,ii], linewidth = 2,
                        color=:black, legend=false)
            end
        else
            #plane surface
            for ii = 1:length(knt1)
                plt = Plots.plot!(p1[1,ii,:], p1[2,ii,:], color=:black, legend=false)
            end
            for ii = 1:length(knt2)
                plt = Plots.plot!(p2[1,ii,:], p2[2,ii,:], color=:black, legend=false)
            end
        end
    elseif length(nurbs.knots)==3 #plot a NURBS volume
        if numDim=="auto"
            numDim="3"
        end
        # Plot the boundaries
        bnd = nrbextract(nurbs)
        for iface = 1:6
            plt = nrbkntplot(bnd[iface], numDim)
        end
    else #plot a NURBS curve
        nrbplot(nurbs, [nsub])
        if numDim=="auto" && any(nurbs.coefs[3,:].!=0)
            numDim="3"
        elseif numDim=="auto" && !any(nurbs.coefs[3,:].!=0)
            numDim="2"
        end
        # And plot the knots
        order = nurbs.order[1]
        p = nrbeval(nurbs, unique(nurbs.knots[1][order:end-order+1]))
        if numDim=="3"
            plt = Plots.scatter!(p[1,:], p[2,:], p[3,:], color=:red, markersize=1,
             shape=:xcross, legend=false)
        else
            plt = Plots.scatter!(p[1,:], p[2,:], color=:red, shape=:xcross,
                    legend=false)
        end
    end
    return plt
end

function nrbctrlplot(nurbs::NURBS, numDim::String="auto")
    if length(nurbs.knots)==3
        nsub = 100
        plt = nrbplot(nurbs, [nsub, nsub, nsub])
        # Plot the control points
        coefs = broadcast(/, nurbs.coefs[1:3,:,:,:], nurbs.coefs[[4],:,:,:])
        coefs = reshape(coefs, 3, :)
        @show coefs
        readline(stdin)
        plt = Plots.scatter!(coefs[1,:], coefs[2,:], coefs[3,:], color=:red,
                markersize=2, legend = false)
        # Plot the control net
        for ii = 1:size(nurbs.coefs, 2)
            for jj = 1:size(nurbs.coefs, 3)
                coefs = reshape(nurbs.coefs[1:3,ii,jj,:], 3, :)
                weights = reshape(nurbs.coefs[[4],ii,jj,:], :)
                xval = coefs[1,:]./weights
                yval = coefs[2,:]./weights
                zval = coefs[3,:]./weights
                plt = Plots.plot!(xval, yval, zval, color = :black,
                    linewidth = 2, linestyle = :dash)
            end
            for kk = 1:size(nurbs.coefs, 4)
                coefs = reshape(nurbs.coefs[1:3,ii,:,kk], 3, :)
                weights = reshape(nurbs.coefs[4,ii,:,kk], :)
                xval = coefs[1,:]./weights
                yval = coefs[2,:]./weights
                zval = coefs[3,:]./weights
                plt = Plots.plot!(xval, yval, zval, color = :black,
                    linewidth = 2, linestyle = :dash)
            end
        end
        for jj = 1:size(nurbs.coefs, 3)
            for kk = 1:size(nurbs.coefs, 4)
                coefs = reshape(nurbs.coefs[1:3, :, jj, kk], 3, :)
                weights = reshape(nurbs.coefs[4, :, jj, kk], :)
                xval = coefs[1,:]./weights
                yval = coefs[2,:]./weights
                zval = coefs[3,:]./weights
                plt = Plots.plot!(xval, yval, zval, color = :black,
                    linewidth = 2, linestyle = :dash)
            end
        end
    elseif length(nurbs.knots)==2
        # Plot a NURBS surface
        nsub = 100
        nrbplot(nurbs, [nsub, nsub])

        # And plot the control net
        for ii = 1:size(nurbs.coefs, 2)
            coefs = reshape(nurbs.coefs[1:3, ii, :], 3, :)
            weights = reshape(nurbs.coefs[4, ii, :], :)
            xval = coefs[1,:]./weights
            yval = coefs[2,:]./weights
            zval = coefs[3,:]./weights
            plt = Plots.plot!(xval, yval, zval, color = :black,
                linewidth = 2, linestyle = :dash)
            plt = Plots.scatter!(xval, yval, zval, color=:red,
                    markersize=2, legend = false)
        end
        for jj = 1:size(nurbs.coefs, 3)
            coefs = reshape(nurbs.coefs[1:3, :, jj], 3, :)
            weights = reshape(nurbs.coefs[4, :, jj], :)
            xval = coefs[1,:]./weights
            yval = coefs[2,:]./weights
            zval = coefs[3,:]./weights
            plt = Plots.plot!(xval, yval, zval, color = :black,
                linewidth = 2, linestyle = :dash)
            plt = Plots.scatter!(xval, yval, zval, color=:red,
                    markersize=2, legend = false)
        end
    else
        # Plot a NURBS curve
        nsub = 1000
        nrbplot(nurbs, [nsub])
        if numDim=="auto" && any(nurbs.coefs[3,:].!=0)
            numDim="3"
        elseif numDim=="auto" && !any(nurbs.coefs[3,:].!=0)
            numDim="2"
        end
        # And plot the control polygon
        coefs = nurbs.coefs[1:3, :]
        weights = nurbs.coefs[4,:]
        xval = coefs[1,:]./weights
        yval = coefs[2,:]./weights
        zval = coefs[3,:]./weights
        if numDim=="3"
            plt = Plots.plot!(xval, yval, zval, color = :black,
                linewidth = 2, linestyle = :dash)
            plt = Plots.scatter!(xval, yval, zval, color=:red, markersize=2,
                legend=false)
        else
            plt = Plots.plot!(xval, yval, color = :black,
                    linewidth = 2, linestyle = :dash)
            plt = Plots.scatter!(xval, yval, color=:red, markersize=5,
                    legend=false)
        end
    end
    return plt
end


text:
using NewtonMethod2
using Test

@testset "NewtonMethod2.jl" begin

### 2 & 3 & 4
f(x)=x^2+2*x-3
f_prime(x)=2*x+2
x0=2.0
t2=newtonroot(f, f_prime; x0, tol=1E-7, maxiter=1000)[1]
t3=newtonroot2(f; x0, tol=1E-7, maxiter=1000)[1]
t4=newtonroot(f, f_prime; x0=BigFloat(x0), tol=1E-7, maxiter=1000)[1]
t5=newtonroot2(f; x0=BigFloat(x0), tol=1E-7, maxiter=1000)[1]
@test t2 ≈ 1.0
@test t3 ≈ 1.0
@test t4 ≈ 1.0
@test t5 ≈ 1.0

f(x)=x^2+2*x-8
f_prime(x)=2*x+2
x0=1.0
t2=newtonroot(f, f_prime; x0, tol=1E-7, maxiter=1000)[1]
t3=newtonroot2(f; x0, tol=1E-7, maxiter=1000)[1]
t4=newtonroot(f, f_prime; x0=BigFloat(x0), tol=1E-7, maxiter=1000)[1]
t5=newtonroot2(f; x0=BigFloat(x0), tol=1E-7, maxiter=1000)[1]
@test t2 ≈ 2.0
@test t3 ≈ 2.0
@test t4 ≈ 2.0
@test t5 ≈ 2.0

f(x)=x-2
f_prime(x)=1
x0=0.0
t2=newtonroot(f, f_prime; x0, tol=1E-7, maxiter=1000)[1]
t3=newtonroot2(f; x0, tol=1E-7, maxiter=1000)[1]
t4=newtonroot(f, f_prime; x0=BigFloat(x0), tol=1E-7, maxiter=1000)[1]
t5=newtonroot2(f; x0=BigFloat(x0), tol=1E-7, maxiter=1000)[1]
@test t2 ≈ 2.0
@test t3 ≈ 2.0
@test t4 ≈ 2.0
@test t5 ≈ 2.0

### 5
f(x)=2+x^2
f_prime(x)=2*x
x0=1.0
t2=newtonroot(f, f_prime; x0, tol=1E-7, maxiter=1000)[1]
t3=newtonroot2(f; x0, tol=1E-7, maxiter=1000)[1]
@test t2 == nothing
@test t3 == nothing

### 6
f(x)=x^2+2*x-3
f_prime(x)=2*x+2
x0=2.0
t2=newtonroot(f, f_prime; x0, tol=1E-7, maxiter=3)[1]
t3=newtonroot2(f; x0, tol=1E-7, maxiter=3)[1]
@test t2 == nothing
@test t3 == nothing

### 7
f(x)=x^2+2*x-3
f_prime(x)=2*x+2
x0=2.0
t2=newtonroot(f, f_prime; x0, tol=1E-7, maxiter=1000)[2]
t3=newtonroot2(f; x0, tol=1E-7, maxiter=1000)[2]
@test t2 <= 1E-7
@test t3 <= 1E-7 ###

end


text:
export ReservoirTrajectory

using MacroTools: @forward
using Random

mutable struct ReservoirTrajectory{B,R<:AbstractRNG} <: AbstractTrajectory
    buffer::B
    n::Int
    capacity::Int
    rng::R
end

@forward ReservoirTrajectory.buffer Base.keys, Base.haskey, Base.getindex

Base.length(x::ReservoirTrajectory) = length(x.buffer[1])

function ReservoirTrajectory(
    capacity,
    kw::Pair{Symbol,DataType}...;
    n = 0,
    rng = Random.GLOBAL_RNG,
)
    buffer = Trajectory(; (s => Vector{t}() for (s, t) in kw)...)
    ReservoirTrajectory(buffer, n, capacity, rng)
end

function Base.push!(b::ReservoirTrajectory; kw...)
    b.n += 1
    if b.n <= b.capacity
        push!(b.buffer; kw...)
    else
        i = rand(b.rng, 1:b.n)
        if i <= b.capacity
            for (k, v) in kw
                b.buffer[k][i] = v
            end
        end
    end
end


text:
# Newton.

function optimize(f::Function, g::Function, H::Function,
                  x0::Vector{T}; kmax::Int64 = 5000, δ::Float64 = 1e-6) where {T<:Real}
    k::Int64 = 0
    δ *= δ
    x = x0
    n::Int64 = length(x)
    dfx = ones(n)
    d2fx = Matrix{Float64}(I, n, n)
    while norm(dfx) > δ && k < kmax
        dfx = g(x)
        d2fx = H(x)
        x -= d2fx \ dfx
        k += 1
    end
    return x, k
end


text:
# A rudimentary XML parser

export parse_elems

namespace{ns}(el::Elem{ns}) = ns
tag{ns, tag}(el::Elem{ns, tag}) = tag

function make_tag(str)
    # This function will get strings such as /tag>
    # or tag attr="val"...>

    if first(str) == '/'
        # this is a closing tag.
        ms = match(r"^/([a-zA-Z\-_0-9]+)>", str)
        if ms === nothing
            error("Invalid closing tag, $str")
        end
        tag = ms.captures[1]
        return Elem(:closing_tag, Symbol(tag)) # Hack alert!
    else
        # this is an opening tag
        ms = match(r"^([a-zA-Z\-_0-9]+)\s*(.*)>", str)
        if ms === nothing
            error("Invalid opening tag, $str")
        end
        tag = ms.captures[1]
        attrs = ms.captures[2]
        if attrs != ""
            attr_ms = matchall(r"([^\s=]*\s*=\s*[^\s=]*)", attrs)
            attr = Any[split(s, "=") for s in attr_ms]
            kvs = Any[(k, strip(v, Set("\"' "))) for (k, v) in attr]
            return Elem(Symbol(tag)) & kvs
        else
            return Elem(Symbol(tag))
        end
    end
end

# Given a string, returns a vector of elements
function parse_elems(str::AbstractString, ns=:xhtml)

    node_stack = Node[Elem(:wrapper, :dummy)] # hack again

    i = start(str)

    while !done(str, i)
        c, i = next(str, i)
        if c == '<'
            j = search(str, '>', i)
            if j===0
                error("Tag at $i never closes?")
            end
            el = make_tag(str[i:j])
            if namespace(el) !== :closing_tag
                # A tag just opened, push it to the stack
                push!(node_stack, el)
            else
                # A tag just closed, check to see if it was on top
                top_el = pop!(node_stack)
                if tag(top_el) === tag(el)
                    # yes, you just closed the tag that you were supposed to
                    # add it to the previous elem in the stack
                    parent = pop!(node_stack)
                    push!(node_stack, parent << top_el)
                else
                    error("Wrong closing tag at $i")
                end
            end
        else
            i = prevind(str, i)
            j = search(str, '<', i)
            if j==0
                substr = str[i:end]
                el = pop!(node_stack)
                push!(node_stack, el << TextNode(str[i:end]))
                break
            else
                el = pop!(node_stack)
                j = prevind(str, j)
                push!(node_stack, el << TextNode(str[i:j]))
            end
        end
        i = j
        c, i = next(str, i)
    end

    if length(node_stack) !== 1
        error("Looks like the HTML is malformed. Is there an un closed tag?")
    else
        el = pop!(node_stack)
        return children(el)
    end
end


text:
using Base.Test
using D3GraphPlot

# write your own tests here
@test 1 == 2


text:
export Add, Activation, Concatenate, Permute
export CausalMeanPool, CausalMaxPool, CausalMinPool

mutable struct Add{F}
    fs::F
end

Add(fs...) = Add(fs)

Flux.@functor Add

(m::Add)(x) = .+(m.fs.(x)...)

mutable struct Activation{F}
    f::F
end

(m::Activation)(x) = m.f.(x)

mutable struct Concatenate{F}
    fs::F
    dims::Int
end

Concatenate(fs...; dims = 1) = Concatenate(fs, dims)

Flux.@functor Concatenate

call(f, x) = f(x)

function (m::Concatenate)(xs)
    if xs isa Tuple
        ys = [f(x) for (f, x) in zip(m.fs, xs)]
    else
        ys = [f(xs) for f in m.fs]
    end
    if m.dims == 1
        reduce(vcat, ys)
    elseif m.dims == 2
        reduce(hcat, ys)
    else
        cat(ys...; dims = m.dims)
    end
end

mutable struct LeftPad{N}
    pad::NTuple{N, Int}
end

(m::LeftPad)(x) = add_padding(x, Tuple(reduce(vcat, collect.(zip(m.pad, zero.(m.pad))))))

mutable struct RightPad{N}
    pad::NTuple{N, Int}
end

(m::RightPad)(x) = add_padding(x, Tuple(reduce(vcat, collect.(zip(zero.(m.pad), m.pad)))))

mutable struct Permute{N}
    dims::NTuple{N, Int}
end

(m::Permute)(x) = permutedims(x, m.dims)

mutable struct FixPoolEdge{N}
    k::NTuple{N, Int}
end

function (m::FixPoolEdge)(x)
    k = min.(m.k, size(x)[1:length(m.k)])
    λ = prod(m.k) ./ Float32[prod(I.I) for I in CartesianIndices(k)]
    λ = Flux.CUDA.Adapt.adapt(typeof(unwrap(x)), λ)
    [λ .* x[UnitRange.(1, k)..., :, :]; x[UnitRange.(k .+ 1, size(x, 1))..., :, :]]
end

CausalMeanPool(k) = Chain(LeftPad(k .- 1), MeanPool(k, stride = one.(k)), FixPoolEdge(k))
CausalMaxPool(k) = Chain(LeftPad(k .- 1), MaxPool(k, stride = one.(k)))
CausalMinPool(k) = Chain(LeftPad(k .- 1), x -> -x, MaxPool(k, stride = one.(k)), x -> -x)


text:
# This file was generated, do not modify it. # hide
v = categorical(["AA", "BB", missing, "AA", "BB", "CC"]);


text:
#
#--------------------------------------#
export linearFEM
#--------------------------------------#
"""
 1D linear FEM operators in (0,1)
"""
function linearFEM(N)

	n  = N + 1
	dz = 1/N
    z  = linspace(0,1,n)

    dn = ones(n)
    dN = ones(N)

    A = Tridiagonal(-1*dN,2*dn ,-1*dN) / dz # == D'*B*D
    B = Tridiagonal(1dN/6,2dn/3,1dN/6) * dz
    C = Tridiagonal(-dN/2,0dn  ,dN/2)        # == B * D
    D = Tridiagonal( 0dN ,-1dn ,1dN)   / dz  # attempt second order?

    A[1,1]     = 0.5 * A[1,1]
    A[end,end] = 0.5 * A[end,end]

    B[1,1]     = 0.5 * B[1,1]
    B[end,end] = 0.5 * B[end,end]

    w = B * ones(n)

return z,w,A,B,C,D
end


text:
using Colab
using Test
using HTTP

@testset "Colab.jl" begin
    @test HTTP.get("https://raw.githubusercontent.com/Chandu-4444/Colab.jl/main/src/assets/Config_Notebook.json").status == 200
end


text:
using Primes

function factors(n)
    f = [one(n)]
    for (p,e) in factor(n)
        f = reduce(vcat, [f*p^j for j in 1:e], init=f)
    end
    return length(f) == 1 ? [one(n), n] : sort!(f)
end

function multorder(a, m)
    gcd(a,m) == 1 || error("$a and $m are not coprime")
    res = one(m)
    for (p,e) in factor(m)
        m = p^e
        t = div(m, p) * (p-1)
        for f in factors(t)
            if powermod(a, f, m) == 1
                res = lcm(res, f)
                break
            end
        end
    end
    res
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
include("InstanceLoader.jl")
include("VRP_PLNE.jl")
include("LSP_PLNE.jl")
include("ResolvePlne.jl")
include("PDI_exact_resolution.jl")
include("VRP_Heuristic.jl")
include("Helper.jl")
include("PDI_heuristic_resolution.jl")

include("GraphHelper.jl")

include("BranchAndCut.jl")
include("allFiles.jl")

#Instances A

INSTANCE_PATH = "../PRP_instances/A_014_#ABS1_15_1.prp"
#INSTANCE_PATH = "../PRP_instances/A_050_ABS12_50_3.prp"
#INSTANCE_PATH = "../PRP_instances/A_050_ABS14_50_1.prp"
#INSTANCE_PATH = "../PRP_instances/A_100_ABS5_100_4.prp"

#INSTANCE_PATH = "/Users/davidpinaud/Desktop/Projet_Probleme-de-Production-et-Distribution-Integre/PRP_instances/A_014_ABS1_15_1.prp"
#INSTANCE_PATH="/Users/davidpinaud/GitHub/Projet_Probleme-de-Production-et-Distribution-Integre/PRP_instances/A_050_ABS14_50_1.prp"

#Instances B
#INSTANCE_PATH = "../PRP_instances/B_200_instance20.prp"

function testGenerateGraph()

	params, nodes, demands, costs = readPRP(INSTANCE_PATH)

	g = generateGraphComplet(nodes)

	generateGraphPDF(".", "graphe.pdf", g)

end

function testLSP(solve=false, verbose=3)

	params, nodes, demands, costs = readPRP(INSTANCE_PATH)

	model = createLSP(params, nodes, demands, costs)

	if solve
		resolvePlne(model, verbose)
	end

end

function testVRP_MTZ(solve=false, t=1, verbose=1)

	params, nodes, demands, costs = readPRP(INSTANCE_PATH)

	model = createVRP_MTZ(params, nodes, demands, costs, t)

	if solve
		resolvePlne(model, verbose)
	end

end

function testPDI_Bard_Nananukul(solve=false, verbose=1)

	params, nodes, demands, costs = readPRP(INSTANCE_PATH)

	model = createPDI_Bard_Nananukul_compacte(params, nodes, demands, costs)
	if solve
		resolvePlne(model, verbose)
	end
	return model

end

function testPDI_heuristique(resoudreVRPwithHeuristic=true, nbMaxIte=10)

	lsp_model, params, nodes, demands, costs, SC, fonctionObjInitial=initialisation_PDI_heuristique(INSTANCE_PATH)
	lsp_model, circuits = PDI_heuristique(lsp_model, params, nodes, demands, costs, SC, fonctionObjInitial, nbMaxIte, resoudreVRPwithHeuristic)

end

#=
t: choisir le temps sur lequel tester
choice: Choisir l'heuristique à utiliser
metaChoice: Choisir la métaheuristique à utiliser (0 si aucune)

showCircuits: Montrer ou non les circuits des heuristiques/métaheuristiques

useLSP: Si vrai, appliquer d'abord le LSP pour obtenir les vrai valeur de VRP

showMTZ:
	0: Ne pas montrer le résultat de MTZ
	1: Montrer le résultat de MTZ
	2: Montrer le circuit de MTZ
	3: Ne montrer que MTZ et pas l'heuristique

	ATTENTION: MTZ sur des instances de plus de 14 risque de prendre très longtemps

heuristicExtraParam:
	Paramètre en plus de ceux de base pour l'heuristique (Exemple l'angle pour sectorielle)

=#
function testHeuristicVRP(;t=1, choice=1, metaChoice=0, showCircuits=false, useLSP=false, showMTZ=0, heuristicExtraParam=[], savePath="")
	allHeuristic = [[binPacking, "Bin packing"], [clark_wright, "Clark-Wright"], [sectorielle, "Sectorielle"]]
	allMetaheuristic = [[TSPBoost, "TSP Local search"], [mixMetaheuristic, "Mix Metaheuristic"]]

	params, nodes, demands, costs = readPRP(INSTANCE_PATH)

	if useLSP

		tToVRP = getTrueVRP(params, nodes, demands, costs)
		params, nodes, demands, costs = tToVRP[t]

	end

	if showMTZ > 0

		modelVRP = createVRP_MTZ(params, nodes, demands, costs, t)
		resolvePlne(modelVRP, 0)

		println("\nMTZ objective value: ", objective_value(modelVRP))

		if showMTZ > 1
			circuits = vrpToCircuit(modelVRP, params)
			totalCost = getCircuitsCost(circuits, costs)

			println("MTZ circuits: ", circuits)
			println("MTZ circuits cost: ", totalCost)

			if showMTZ > 2
				return
			end

		end

		println()

	end

	heuristic = allHeuristic[choice]

	if choice < 3
		heuristicExtraParam = []
	end
	circuits = heuristic[1](params, nodes, demands, costs, t, heuristicExtraParam...)
	totalCost = getCircuitsCost(circuits, costs)

	if showCircuits
		println("$(heuristic[2]) circuits: ", circuits)
	end

	println("$(heuristic[2]) circuits cost: ", totalCost)

	if metaChoice > 0
		println()

		metaheuristic = allMetaheuristic[metaChoice]

		circuits = metaheuristic[1](circuits, params, costs)
		totalCost = getCircuitsCost(circuits, costs)

		if showCircuits
			println("$(metaheuristic[2]) circuits: ", circuits)
		end

		println("$(metaheuristic[2]) circuits cost: ", totalCost)

	end

	if length(savePath) > 0
		saveCircuits(params, circuits, savePath)
	end

end

function testPDI_Boudia(solve=false, verbose=1)

	params, nodes, demands, costs = readPRP(INSTANCE_PATH)

	model = createPDI_Boudia(params, nodes, demands, costs)
	if solve
		resolvePlne(model, verbose)
	end

end

#testGenerateGraph()
#testLSP(true)
#testVRP_MTZ(true)

function testBranchAndCutPDI(filePath)

	params, nodes, demands, costs = readPRP(INSTANCE_PATH)

	model = BranchAndCutPDI(params, nodes, demands, costs)

	println("optimum = ", objective_value(model))

	allCircuits = PDItoCircuits(model, params, nodes, demands, costs)

	saveMultiCircuits(params, allCircuits, filePath)

end

function testCompareExact(logPath)

	open(logPath, "a") do io

		t = 1

		path_to_file = "../PRP_instances/"

		for file in allFiles

			params, nodes, demands, costs = readPRP(INSTANCE_PATH)

			model, totalExactTime = @timed BranchAndCutPDI(params, nodes, demands, costs)

			exactValue = objective_value(model)

			resoudreVRPwithHeuristic=true
			nbMaxIte=10

			lsp_model, params, nodes, demands, costs, SC, fonctionObjInitial=initialisation_PDI_heuristique(INSTANCE_PATH)
			(lsp_model, circuits), totalTimeHeuristic = @timed PDI_heuristique(lsp_model, params, nodes, demands, costs, SC, fonctionObjInitial, nbMaxIte, resoudreVRPwithHeuristic)

			heuristicValue = objective_value(lsp_model)

			println("Exacte: ", exactValue)
			println("Heuristique: ", heuristicValue)

			println(io,"$file,$(exactValue),$(totalExactTime),$(heuristicValue),$(totalTimeHeuristic)")

		end #Loop file

	end #Open

end

function testLogCompareHeuristic(logPath; heuristicExtraParam=[30])
	allHeuristic = [[binPacking, "Bin packing"], [clark_wright, "Clark-Wright"], [sectorielle, "Sectorielle"]]
	allMetaheuristic = [[TSPBoost, "TSP Local search"], [mixMetaheuristic, "Mix Metaheuristic"]]

	# modelVRP = createVRP_MTZ(params, nodes, demands, costs, t)
	# resolvePlne(modelVRP, 0)

	# objective_value(modelVRP)

	open(logPath, "a") do io

		t = 1

		path_to_file = "../PRP_instances/"

		for file in allFiles

			params, nodes, demands, costs = readPRP(path_to_file*file)

			rez = Vector{Float64}(undef, length(allHeuristic) * length(allMetaheuristic))

			for (i, heuristic) in enumerate(allHeuristic)

				heuristicExtraParamCurrent = heuristicExtraParam

				if i < 3
					heuristicExtraParamCurrent = []
				end
				circuits = heuristic[1](params, nodes, demands, costs, t, heuristicExtraParamCurrent...)

				for (j, metaheuristic) in enumerate(allMetaheuristic)

					circuits = metaheuristic[1](circuits, params, costs)
					totalCost = getCircuitsCost(circuits, costs)

					#println(heuristic[2], " ", metaheuristic[2], " = ", totalCost)
					rez[i + length(allHeuristic) * (j - 1)] = totalCost

				end

			end

			println(io,"$file,$(rez[1]),$(rez[2]),$(rez[3]),$(rez[4]),$(rez[5]),$(rez[6])")

		end #Loop file

	end #Open
end

function testLogCompareMTZ_Heuristic(logPath)

	heuristic = clark_wright
	metaheuristic = mixMetaheuristic

	# modelVRP = createVRP_MTZ(params, nodes, demands, costs, t)
	# resolvePlne(modelVRP, 0)

	# objective_value(modelVRP)

	open(logPath, "a") do io

		t = 1

		path_to_file = "../PRP_instances/"

		for file in allFiles

			params, nodes, demands, costs = readPRP(path_to_file*file)

			modelVRP = createVRP_MTZ(params, nodes, demands, costs, t)
			modelVRP, totalMTZtime = @timed resolvePlne(modelVRP, 0)

			mtzVal = objective_value(modelVRP)

			totalTimeHeuristic = 0

			circuits, elapsedTime = @timed heuristic(params, nodes, demands, costs, t)
			totalTimeHeuristic += elapsedTime

			circuits, elapsedTime = @timed metaheuristic(circuits, params, costs)
			totalTimeHeuristic += elapsedTime

			totalCost = getCircuitsCost(circuits, costs)

			println(io,"$file,$(mtzVal),$(totalMTZtime),$(totalCost),$(totalTimeHeuristic)")

		end #Loop file

	end #Open
end

testCompareExact("pdiCompare.csv")

#testLogCompareHeuristic("metaheuristic.csv", heuristicExtraParam=[10])

#testLogCompareMTZ_Heuristic("mtz_meta.csv")

#Nouvelle fonction qui réunis toutes les heuristique, le MTZ et le LSP

#testHeuristicVRP(t=3, choice=2, metaChoice=1, showCircuits=false, useLSP=false, heuristicExtraParam=[10], showMTZ=0, savePath="../Save/test.png")

#testHeuristicVRP(t=4, choice=2, metaChoice=1, showCircuits=false, useLSP=true, heuristicExtraParam=[10], showMTZ=0, savePath="../Save/test.png")

#testPDI_heuristique()
#testPDI_Boudia(true)
#testPDI_Bard_Nananukul(true,2)

#testBranchAndCutPDI("./Save/BAC_David.png")


text:
"""
Data pre- / post-processing for Slack problems. Takes a ProfileData object and prepares it for use in GP.
"""

"""
get_predictors_targets(vavg::Array, problem::Slack_KPP)
# Description
    Returns x and y, the scaled predictors and target predictions from which to extract the training and verification data for temperature profiles.
    Scales the predictors and targets using min-max scaling based on the initial temperature profile from the les simulation.

    model( predictors ) = targets
    model( KPP(T[i])  ) = T[i] - KPP(T[i])

# Arguments
- `vavg`: (Array)               Nt-length array of D-length vectors. Data from which to extract x and y, the predictors and corresponding predictions.
-  `problem`: (Slack_KPP)    Slack_KPP object associated with the data (output of get_problem)

# *--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*
# | Slack_KPP                                     |
# |                                                  |
# |   predictor            target                    |
# |   KPP(T[i]) --model--> T[i] - KPP(T[i])          |
# |                                                  |
# |   model output                                   |
# |   model(KPP(T[i]))                               |
# |                                                  |
# |   predicted T[i] = model(KPP(T[i])) + KPP(T[i])  |
# *--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*

# *--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*
# | Slack_TKE                                     |
# |                                                  |
# |   predictor            target                    |
# |   TKE(T[i]) --model--> T[i] - TKE(T[i])          |
# |                                                  |
# |   model output                                   |
# |   model(TKE(T[i]))                               |
# |                                                  |
# |   predicted T[i] = model(TKE(T[i])) + TKE(T[i])  |
# *--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*
"""
function get_predictors_targets(vavg::Array, problem::Union{Slack_KPP, Slack_TKE})

    # scale according to initial temperature profile vavg[1]
    scaling = min_max_scaling(vavg[1][end]-vavg[1][1], minimum(vavg[1]))
    vavg = [scale(vec, scaling) for vec in vavg]

    # store the scaling as an attribute in the problem for later use
    problem.scaling = scaling

    # scale kpp data with same scaling
    predictors = [scale(vec, scaling) for vec in problem.physics_data] # kpp(T[i])

    targets = vavg .- predictors # residual
    return (predictors, targets)
end

"""
# Description
Takes in a predictor, T[i], the GP prediction on the predictor, G(T[i]), and a Slack_KPP object.
Returns the predicted temperature profile, T[i+1], computed from T[i] and G(T[i]) by

           prediction = model(predictor) # residual
     predicted T[i+1] = model(KPP(T[i])) + KPP(T[i])

# Arguments
- `predictor`: (Array)            T[i], the predictor for a temperature profile
- `prediction`: (Array)           model(T[i), the prediction for a temperature profile
- `problem`: (Slack_KPP)
"""
function postprocess_prediction(predictor, prediction, problem::Union{Slack_KPP, Slack_TKE})

    #unscale predictor
    predictor = [unscale(vec, problem.scaling) for vec in predictor] # kpp data

    # predictor = problem.kpp_data

    #unscale prediction
    prediction = [(vec * problem.scaling.Δv) for vec in prediction] # residual

    return predictor .+ prediction

    # return [unscale(a .+ b, problem.scaling) for a in predictor, b in prediction]
end

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
for N = 1:48
    run(`julia --threads $N run_X_threads.jl`)
    run(`mkdir N$N`)

    for f in readdir(@__DIR__)
        if occursin(".out", f)
            mv(f, "N$N/$f")
        end
    end
end


text:
mach_test = BaseMachine(   0.0, #R
                       0.2995, #Xd_p
                       0.588, #eq_p
                       100.0)  #MVABase

mach2_test = OneDOneQMachine(0.0, #R
                      0.8, #Xd
                       0.2, #Xq
                       0.1813, #Xd_p
                       0.05, #Xq_p
                       5.89, #Td0_p
                       0.6, #Tq0_p
                       100.0)   #MVABase

shaft_test = SingleMass(3.148, #H
                       2.0) #D

shaft_test2 = SingleMass(3.01, #H
                      2.0) #D

pss_test = PSSFixed(0.0)
tg_test = TGFixed(1.0)
avr_test = AVRFixed(1.0)
avr_test2 = AVRFixed(0.26854)
#avr_test3 = AVRFixed(0.3337)
avr_test3 = AVRFixed(0.34992)

load_test = PowerLoad("Bus3", true, nodes_test[1], PowerSystems.ConstantPower, 0.0, 0.0, 0.0, 0.0)
load_test2 = PowerLoad("Bus3", true, nodes_test2[3], PowerSystems.ConstantPower, 0.6, 0.05, 0.6, 0.05)

load_test3a = PowerLoad("Bus3", true, nodes_test2[2], PowerSystems.ConstantPower, 0.3, 0.05, 0.6, 0.05)
load_test3b = PowerLoad("Bus3", true, nodes_test2[3], PowerSystems.ConstantPower, 0.3, 0.05, 0.6, 0.05)

Gen_test = DynGenerator(1, #Number
                 :TestGen,
                 nodes_test[2], #bus
                 1.0, # ω_ref,
                 1.0, #V_ref
                 1.0, #P_ref
                 mach_test,
                 shaft_test,
                 avr_test, #avr
                 tg_test, #tg
                 pss_test)

Gen_test2 = DynGenerator(1, #Number
                :TestGen,
                nodes_test[2], #bus
                1.0, # ω_ref,
                1.0, #V_ref
                1.0, #P_ref
                mach2_test,
                shaft_test,
                avr_test2, #avr
                tg_test, #tg
                pss_test)

Gen_test3 = DynGenerator(1, #Number
              :TestGen,
              nodes_test[2], #bus
              1.0, # ω_ref,
              1.0, #V_ref
              0.5, #P_ref
              mach2_test,
              shaft_test,
              avr_test3, #avr
              tg_test, #tg
              pss_test)

OMIB_test =  DynamicSystem(nodes_test, branch_test, [Gen_test], [load_test], 100.0, 60.0)
OMIB_test2 =  DynamicSystem(nodes_test, branch_test, [Gen_test2], [load_test], 100.0, 60.0)
OMIB_test_load = DynamicSystem(nodes_test2, branch_test2, [Gen_test3], [load_test2], 100.0, 60.0)
OMIB_test_load2 = DynamicSystem(nodes_test2, branch_test2, [Gen_test3], [load_test3a, load_test3b], 100.0, 60.0)


text:
abstract type AbstractOptimizer end
abstract type AbstractConstrainedOptimizer <: AbstractOptimizer end
abstract type ZerothOrderOptimizer <: AbstractOptimizer end
abstract type FirstOrderOptimizer  <: AbstractOptimizer end
abstract type SecondOrderOptimizer <: AbstractOptimizer end
abstract type UnivariateOptimizer  <: AbstractOptimizer end

abstract type AbstractOptimizerState end
abstract type ZerothOrderState <: AbstractOptimizerState end

"""
Configurable options with defaults (values 0 and NaN indicate unlimited):
```
x_abstol::Real = 0.0,
x_reltol::Real = 0.0,
f_abstol::Real = 0.0,
f_reltol::Real = 0.0,
g_abstol::Real = 1e-8,
g_reltol::Real = 1e-8,
outer_x_abstol::Real = 0.0,
outer_x_reltol::Real = 0.0,
outer_f_abstol::Real = 0.0,
outer_f_reltol::Real = 0.0,
outer_g_abstol::Real = 1e-8,
outer_g_reltol::Real = 1e-8,
f_calls_limit::Int = 0,
g_calls_limit::Int = 0,
h_calls_limit::Int = 0,
allow_f_increases::Bool = true,
allow_outer_f_increases::Bool = true,
successive_f_tol::Int = 1,
iterations::Int = 1_000,
outer_iterations::Int = 1000,
store_trace::Bool = false,
show_trace::Bool = false,
extended_trace::Bool = false,
show_every::Int = 1,
callback = nothing,
time_limit = NaN
```
See http://julianlsolvers.github.io/Optim.jl/stable/#user/config/
"""
struct Options{T, TCallback}
    x_abstol::T
    x_reltol::T
    f_abstol::T
    f_reltol::T
    g_abstol::T
    g_reltol::T
    outer_x_abstol::T
    outer_x_reltol::T
    outer_f_abstol::T
    outer_f_reltol::T
    outer_g_abstol::T
    outer_g_reltol::T
    f_calls_limit::Int
    g_calls_limit::Int
    h_calls_limit::Int
    allow_f_increases::Bool
    allow_outer_f_increases::Bool
    successive_f_tol::Int
    iterations::Int
    outer_iterations::Int
    store_trace::Bool
    trace_simplex::Bool
    show_trace::Bool
    extended_trace::Bool
    show_every::Int
    callback::TCallback
    time_limit::Float64
end

function Options(;
        x_tol = nothing,
        f_tol = nothing,
        g_tol = nothing,
        x_abstol::Real = 0.0,
        x_reltol::Real = 0.0,
        f_abstol::Real = 0.0,
        f_reltol::Real = 0.0,
        g_abstol::Real = 1e-8,
        g_reltol::Real = 1e-8,
        outer_x_tol = 0.0,
        outer_f_tol = 0.0,
        outer_g_tol = nothing,
        outer_x_abstol::Real = 0.0,
        outer_x_reltol::Real = 0.0,
        outer_f_abstol::Real = 0.0,
        outer_f_reltol::Real = 0.0,
        outer_g_abstol::Real = 1e-8,
        outer_g_reltol::Real = 1e-8,
        f_calls_limit::Int = 0,
        g_calls_limit::Int = 0,
        h_calls_limit::Int = 0,
        allow_f_increases::Bool = true,
        allow_outer_f_increases::Bool = true,
        successive_f_tol::Int = 1,
        iterations::Int = 1_000,
        outer_iterations::Int = 1000,
        store_trace::Bool = false,
        trace_simplex::Bool = false,
        show_trace::Bool = false,
        extended_trace::Bool = false,
        show_every::Int = 1,
        callback = nothing,
        time_limit = NaN)
    show_every = show_every > 0 ? show_every : 1
    #if extended_trace && callback == nothing
    #    show_trace = true
    #end
    if !(x_tol == nothing)
        x_abstol = x_tol
    end
    if !(g_tol == nothing)
        g_abstol = g_tol
    end
    if !(f_tol == nothing)
        f_reltol = f_tol
    end
    if !(outer_x_tol == nothing)
        outer_x_abstol = outer_x_tol
    end
    if !(outer_g_tol == nothing)
        outer_g_abstol = outer_g_tol
    end
    if !(outer_f_tol == nothing)
        outer_f_reltol = outer_f_tol
    end
    Options(promote(x_abstol, x_reltol, f_abstol, f_reltol, g_abstol, g_reltol, outer_x_abstol, outer_x_reltol, outer_f_abstol, outer_f_reltol, outer_g_abstol, outer_g_reltol)..., f_calls_limit, g_calls_limit, h_calls_limit,
        allow_f_increases, allow_outer_f_increases, successive_f_tol, Int(iterations), Int(outer_iterations), store_trace, trace_simplex, show_trace, extended_trace,
        Int(show_every), callback, Float64(time_limit))
end

function Base.show(io::IO, o::Optim.Options)
    for k in fieldnames(typeof(o))
        v = getfield(o, k)
        if v isa Nothing
            @printf io "%24s = %s\n" k "nothing"
        else
            @printf io "%24s = %s\n" k v
        end
    end
end

function print_header(options::Options)
    if options.show_trace
        @printf "Iter     Function value   Gradient norm \n"
    end
end

function print_header(method::AbstractOptimizer)
        @printf "Iter     Function value   Gradient norm \n"
end

struct OptimizationState{Tf<:Real, T <: AbstractOptimizer}
    iteration::Int
    value::Tf
    g_norm::Tf
    metadata::Dict
end

const OptimizationTrace{Tf, T} = Vector{OptimizationState{Tf, T}}

abstract type OptimizationResults end

mutable struct MultivariateOptimizationResults{O, T, Tx, Tc, Tf, M, Tls, Tsb} <: OptimizationResults
    method::O
    initial_x::Tx
    minimizer::Tx
    minimum::Tf
    iterations::Int
    iteration_converged::Bool
    x_converged::Bool
    x_abstol::T
    x_reltol::T
    x_abschange::Tc
    x_relchange::Tc
    f_converged::Bool
    f_abstol::T
    f_reltol::T
    f_abschange::Tc
    f_relchange::Tc
    g_converged::Bool
    g_abstol::T
    g_residual::Tc
    f_increased::Bool
    trace::M
    f_calls::Int
    g_calls::Int
    h_calls::Int
    ls_success::Tls
    time_limit::Float64
    time_run::Float64
    stopped_by::Tsb
end

# pick_best_x and pick_best_f are used to pick the minimizer if we stopped because
# f increased and we didn't allow it
pick_best_x(f_increased, state) = f_increased ? state.x_previous : state.x
pick_best_f(f_increased, state, d) = f_increased ? state.f_x_previous : value(d)

function Base.show(io::IO, t::OptimizationState)
    @printf io "%6d   %14e   %14e\n" t.iteration t.value t.g_norm
    if !isempty(t.metadata)
        for (key, value) in t.metadata
            @printf io " * %s: %s\n" key value
        end
    end
    return
end

function Base.show(io::IO, tr::OptimizationTrace)
    @printf io "Iter     Function value   Gradient norm \n"
    @printf io "------   --------------   --------------\n"
    for state in tr
        show(io, state)
    end
    return
end

function Base.show(io::IO, r::MultivariateOptimizationResults)
    take = Iterators.take
    failure_string = "failure"
    if iteration_limit_reached(r)
        failure_string *= " (reached maximum number of iterations)"
    end
    if f_increased(r) && !iteration_limit_reached(r)
        failure_string *= " (objective increased between iterations)"
    end
    if isa(r.ls_success, Bool) && !r.ls_success
        failure_string *= " (line search failed)"
    end
    if time_run(r) > time_limit(r)
        failure_string *= " (exceeded time limit of $(time_limit(r)))"
    end

    @printf io " * Status: %s\n\n" converged(r) ? "success" : failure_string

    @printf io " * Candidate solution\n"
    @printf io "    Final objective value:     %e\n" minimum(r)
    @printf io "\n"

    @printf io " * Found with\n"
    @printf io "    Algorithm:     %s\n" summary(r)

    @printf io "\n"
    @printf io " * Convergence measures\n"
    if isa(r.method, NelderMead)
        @printf io "    √(Σ(yᵢ-ȳ)²)/n %s %.1e\n" g_converged(r) ? "≤" : "≰" g_tol(r)
    else
        @printf io "    |x - x'|               = %.2e %s %.1e\n"  x_abschange(r) x_abschange(r)<=x_abstol(r) ? "≤" : "≰" x_abstol(r)
        @printf io "    |x - x'|/|x'|          = %.2e %s %.1e\n"  x_relchange(r) x_relchange(r)<=x_reltol(r) ? "≤" : "≰" x_reltol(r)
        @printf io "    |f(x) - f(x')|         = %.2e %s %.1e\n"  f_abschange(r) f_abschange(r)<=f_abstol(r) ? "≤" : "≰" f_abstol(r)
        @printf io "    |f(x) - f(x')|/|f(x')| = %.2e %s %.1e\n"  f_relchange(r) f_relchange(r)<=f_reltol(r) ? "≤" : "≰" f_reltol(r)
        @printf io "    |g(x)|                 = %.2e %s %.1e\n" g_residual(r) g_residual(r)<=g_tol(r) ?  "≤" : "≰" g_tol(r)
    end

    @printf io "\n"

    @printf io " * Work counters\n"
    @printf io "    Seconds run:   %d  (vs limit %d)\n" time_run(r) isnan(time_limit(r)) ? Inf : time_limit(r)
    @printf io "    Iterations:    %d\n" iterations(r)
    @printf io "    f(x) calls:    %d\n" f_calls(r)
    if !(isa(r.method, NelderMead) || isa(r.method, SimulatedAnnealing))
        @printf io "    ∇f(x) calls:   %d\n" g_calls(r)
    end
    if isa(r.method, Newton) || isa(r.method, NewtonTrustRegion)
        @printf io "    ∇²f(x) calls:  %d\n" h_calls(r)
    end
    return
end

function Base.append!(a::MultivariateOptimizationResults, b::MultivariateOptimizationResults)
    a.iterations += iterations(b)
    a.minimizer = minimizer(b)
    a.minimum = minimum(b)
    a.iteration_converged = iteration_limit_reached(b)
    a.x_converged = x_converged(b)
    a.f_converged = f_converged(b)
    a.g_converged = g_converged(b)
    append!(a.trace, b.trace)
    a.f_calls += f_calls(b)
    a.g_calls += g_calls(b)
end


text:
using LinearAlgebra
using StaticArrays
using OffsetArrays
using Plots
using StiffKinetic
using Printf

Nx = 2^8 + 1
x  = range(0.0, 4.0, length=Nx)
dx = step(x)

uL, uR = 4.0, 2.0
U0 = uL*Float64.(0.0 .<= x .<= 2.0) + uR*Float64.(2.0 .< x .<= 4.0)
V0 = zero(U0)

p = u -> u
f = u -> u # /!\ we may use |f'| = 1 in the schemes

pb = RelaxPb(U0,V0,f,p)
pbs = [pb]

ε = 2^-1

T = 0.25
dt = (2^-4) * dx^2

dt_anim = 2^-6 * T
nt_anim = Int(T/dt_anim)

nt_pf = Int(dt_anim/dt)

anim = @animate for t in 0 : dt_anim : T
    str_t = @sprintf "t = %1.3f" t
    plot(x, pbs[1].U.data, title = str_t, ylims=(1.8, 4.2))
    # plot!(x, pbs[1].V.data, title = "t = $t", ylims=(-5, 4.2))
    for _ in 1 : nt_pf
        imex_bdf4!(pbs, dt, dx, ε)
    end
end

gif(anim, "riemann.gif", fps=24)

max_stars_repo_name: tremelow/StiffKinetic

max_issues_repo_name: tremelow/StiffKinetic

max_forks_repo_name: tremelow/StiffKinetic