text:
"""
    waveguide_debug

A module with an implementation of the waveguide functionality
using MATLAB callbacks.
Used as rerefence implementation to check that Julia implementation works.
"""
module waveguide_debug
push!(LOAD_PATH, string(@__DIR__,"/../"))	# looks for modules in the current directory
push!(LOAD_PATH, string(@__DIR__,"/../../"))	# looks for modules in the current directory

using MATLAB

using Statistics
using Random
using LinearAlgebra
using NonlinearEigenproblems
using NonlinearEigenproblems.NEPCore
using NonlinearEigenproblems.NEPTypes
using NonlinearEigenproblems.LinSolvers
using NonlinearEigenproblems.NEPSolver
using NonlinearEigenproblems.Gallery
using GalleryWaveguide

export matlab_debug_WEP_FD
export matlab_debug_full_matrix_WEP_FD_SPMF
export debug_sqrt_schur
export fft_debug_mateq
export debug_sqrt_derivative
export matlab_debug_Schur_WEP_FD
export debug_Mlincomb_FD_WEP
export debug_Sylvester_SMW_WEP
export matlab_debug_eigval_comp_WEP_FD_and_SPMF
export debug_WEP_FD_preconditioner
export debug_eigval_comp_WEP_FD

#OBS: Explicit imports of non-exported functions
import GalleryWaveguide.generate_wavenumber_fd
import GalleryWaveguide.generate_fd_interior_mat
import GalleryWaveguide.generate_fd_boundary_mat
import GalleryWaveguide.generate_R_matvecs
import GalleryWaveguide.generate_S_function
import GalleryWaveguide.sqrt_schur_pos_imag
import GalleryWaveguide.sqrt_derivative
import GalleryWaveguide.SchurMatVec
import GalleryWaveguide.solve_wg_sylvester_fft!
import GalleryWaveguide.generate_smw_matrix
import GalleryWaveguide.solve_smw
import GalleryWaveguide.construct_WEP_schur_complement

########### SOME REFERENCE IMPLEMENTATIONS ################
###########################################################

# Generate P-matrix
# Is the lower right part of the system matrix, from the DtN maps Jarlebring-(1.5)(1.6) and Ringh-(2.4)(2.8)
function generate_P_matrix(nz::Integer, hx, Km, Kp)

    R, Rinv = generate_R_matvecs(nz::Integer)
    p = (nz-1)/2;

    # Constants from the problem
    d0 = -3/(2*hx);
    a = ones(ComplexF64,nz);
    b = 4*pi*1im * (-p:p);
    cM = Km^2 .- 4*pi^2 * ((-p:p).^2);
    cP = Kp^2 .- 4*pi^2 * ((-p:p).^2);

    function betaM(γ)
        return a*γ^2 + b*γ + cM
    end
    function betaP(γ)
        return a*γ^2 + b*γ + cP
    end

    signM = 1im*sign.(imag(betaM(-1-1im))); # OBS! LEFT HALF-PLANE!
    signP = 1im*sign.(imag(betaP(-1-1im))); # OBS! LEFT HALF-PLANE!

    function sM(γ::Number)
        return signM.*sqrt.(betaM(γ)) .+d0;
    end
    function sP(γ::Number)
        return signP.*sqrt.(betaP(γ)) .+d0;
    end

    function p_sM(γ)
        return signM.*(2*a*γ+b)./(2*sqrt(a*γ^2+b*γ+cM));
    end
    function p_sP(γ)
        return signP.*(2*a*γ+b)./(2*sqrt(a*γ^2+b*γ+cP));
    end

    # BUILD THE FOURTH BLOCK P
    function P(γ,x::Union{Vector{ComplexF64}, Vector{Float64}})
        return vec( [R(Rinv(x[1:Int(end/2)]) .* sM(γ));
                     R(Rinv(x[Int(end/2)+1:end]) .* sP(γ))  ])
    end

    # BUILD THE DERIVATIVE OF P
    function p_P(γ,x::Union{Vector{ComplexF64}, Vector{Float64}})
        return vec( [R(Rinv(x[1:Int(end/2)]) .* p_sM(γ));
                     R(Rinv(x[Int(end/2)+1:end]) .* p_sP(γ))  ])
    end

    return P, p_P
end

# Generate a function for mat-vecs with P^{-1}-matrix
# P is the lower right part of the system matrix, from the DtN maps Jarlebring-(1.5)(1.6) and Ringh-(2.4)(2.8)
function generate_Pinv_matrix(nz::Integer, hx, Km, Kp)

   R, Rinv = generate_R_matvecs(nz::Integer)
   p = (nz-1)/2;

   # Constants from the problem
   d0 = -3/(2*hx);
   a = ones(ComplexF64,nz);
   b = 4*pi*1im * (-p:p);
   cM = Km^2 .- 4*pi^2 * ((-p:p).^2);
   cP = Kp^2 .- 4*pi^2 * ((-p:p).^2);

   function betaM(γ)
       return a*γ^2 + b*γ + cM
   end
   function betaP(γ)
       return a*γ^2 + b*γ + cP
   end

   function sM(γ::Number)
       bbeta = betaM(γ)
       return 1im*sign.(imag(bbeta)).*sqrt.(bbeta) .+ d0
   end
   function sP(γ::Number)
       bbeta = betaP(γ)
       return 1im*sign.(imag(bbeta)).*sqrt.(bbeta) .+ d0
   end

   # BUILD THE INVERSE OF THE FOURTH BLOCK P
   function Pinv(γ,x::Union{Vector{ComplexF64}, Vector{Float64}})
       return vec(  [R(Rinv(x[1:Int(end/2)]) ./ sM(γ));
                     R(Rinv(x[Int(end/2)+1:end]) ./ sP(γ))  ]  )
   end

   return Pinv
end

###########################################################
# Compute the matrix square root
# (only reference implementation, see sqrt_schur_pos_imag)
function sqrt_schur(A::AbstractMatrix)
    n = size(A,1);
    (T, Q) = schur(complex(A))
    U = zeros(ComplexF64,n,n);
    for i = 1:n
        U[i,i] = sqrt(T[i,i])
    end
    for j = 2:n
        for i = (j-1):-1:1
            temp = zero(ComplexF64)
            for k = (i+1):(j-1)
                temp += U[i,k]*U[k,j]
            end
            U[i,j] = (T[i,j] - temp)/(U[i,i]+U[j,j])
        end
    end
    return Q*U*Q'
end

######################## DEBUG ############################
###########################################################

# Test the basic generated FD matrices against MATLAB code
function matlab_debug_WEP_FD(nx::Integer, nz::Integer, delta::Number)
    println("\n\n--- Debugging Matrices FD against MATLAB ---\n")
    if(nx > 200 || nz > 200)
        @warn "This debug is 'naive' and might be slow for the discretization used."
    end

    #The error observed likely comes from difference in linspace-implementation.
    γ = -rand() - 1im*rand()
    gamma = γ

    for waveguide = ["TAUSCH", "JARLEBRING"]
        println("\n")
        println("Testing waveguide: ", waveguide)

        K, hx, hz, Km, Kp = generate_wavenumber_fd( nx, nz, waveguide, delta)
        Dxx, Dzz, Dz = generate_fd_interior_mat( nx, nz, hx, hz)
        C1, C2T = generate_fd_boundary_mat( nx, nz, hx, hz)
        P, p_P = generate_P_matrix(nz, hx, Km, Kp)
        Pinv = generate_Pinv_matrix(nz, hx, Km, Kp)
        P_j = zeros(ComplexF64, 2*nz,2*nz)
        Iz = Matrix(1.0*I, 2*nz, 2*nz)
        for i = 1:2*nz
            P_j[:,i] = P(γ, Iz[:,i])
        end

        R, Rinv = generate_R_matvecs(nz)
        S = generate_S_function(nz, hx, Km, Kp)
        P_j2 = zeros(ComplexF64, 2*nz,2*nz)
        D1 = zeros(ComplexF64, nz,nz)
        D2 = zeros(ComplexF64, nz,nz)
        for j = 1:nz
            D1[j,j] = S(reshape([γ],1,1),j)[1]
        end
        for j = 1:nz
            D2[j,j] = S(reshape([γ],1,1),j+nz)[1]
        end
        Iz = Matrix(1.0*I, nz, nz);
        for j = 1:nz
            P_j2[1:nz,j] = R(D1*Rinv(Iz[:,j]))
            P_j2[(nz+1):(2*nz),j+nz] = R(D2*Rinv(Iz[j,:]))
        end

        if waveguide == "JARLEBRING"
            waveguide_str = "CHALLENGE"
        else
            waveguide_str = waveguide
        end

        println("  -- Matlab printouts start --")
        WEP_path = string(@__DIR__, "/WEP-ref")
        mat"""
            addpath($WEP_path)
            nxx = double($nx);
            nzz = double($nz);
            gamma_m = double($gamma);
            options = [];
            options.delta = $delta;
            options.wg = $waveguide_str;
            matlab_nep = nep_wg_generator(nxx, nzz, options);

            $P_m = NaN(2*nzz, 2*nzz);
            Iz = eye(2*nzz);
            for i = 1:2*nzz
              $P_m(:,i) = matlab_nep.P(gamma_m, Iz(:,i));
            end
            $C1_m = matlab_nep.C1;
            $C2T_m = matlab_nep.C2T;
            $K_m = matlab_nep.K;
            $hx_m = matlab_nep.hx;
            $hz_m = matlab_nep.hz;
        """
        println("  -- Matlab printouts end --")

        println("Difference hx_m - hx = ", abs(hx_m-hx))
        println("Relative difference (hx_m - hx)/hx = ", abs(hx_m-hx)/abs(hx))
        println("Difference hz_m - hz = ", abs(hz_m-hz))
        println("Difference K_m  -K = ", opnorm(K_m-K))
        println("Difference C1_m - C1 = ", opnorm(Matrix(C1_m-C1)))
        println("Relative difference opnorm(C1_m - C1)/opnorm(C1) = ", opnorm(Matrix(C1_m-C1))/opnorm(Matrix(C1)))
        println("Difference C2T_m - C2T = ", opnorm(Matrix(C2T_m-C2T)))
        println("Relative difference opnorm(C2T-m - C2T)/opnorm(C2T) = ", opnorm(Matrix(C2T_m-C2T))/opnorm(Matrix(C2T)))
        println("Difference P_m(γ) - P(γ) = ", opnorm(P_m-P_j))
        println("Relative difference opnorm(P_m(γ) - P(γ))/opnorm(P(γ)) = ", opnorm(P_m-P_j)/opnorm(P_j))
        println("Difference P_m(γ) - P_2(γ) = ", opnorm(P_m-P_j2))
        println("Relative difference opnorm(P_m(γ) - P_2(γ))/opnorm(P_2(γ)) = ", opnorm(P_m-P_j2)/opnorm(Matrix(P_j2)))
        v = rand(ComplexF64,2*nz)
        println("Relative difference norm(v - Pinv_j((γ))*P_j(γ)*v)/norm(v) = ", norm(v - Pinv(γ,P(γ,v)) )/norm(v))
    end
    println("\n--- End Matrices FD against MATLAB ---\n")
end

###########################################################
# Test the full generated system-matrix against MATLAB code
function matlab_debug_full_matrix_WEP_FD_SPMF(nx::Integer, nz::Integer, delta::Number)
    println("\n\n--- Debugging Full Matrix FD SPMF against MATLAB ---\n")
    if(nx > 40 || nz > 40)
        @warn "This debug is 'naive' and might be slow for the discretization used."
    end

    γ = -rand() - 1im*rand()
    gamma = γ

    for waveguide = ["TAUSCH", "JARLEBRING"]
        println("\n")
        println("Testing full matrix M for waveguide: ", waveguide)

        nep_j = nep_gallery(WEP, nx = nx, nz = nz, benchmark_problem = waveguide, neptype = "SpmF", delta = delta)
        M_j = compute_Mder(nep_j,γ)

        if waveguide == "JARLEBRING"
            waveguide_str = "CHALLENGE"
        else
            waveguide_str = waveguide
        end

        println("  -- Matlab printouts start --")
        WEP_path = string(@__DIR__, "/WEP-ref")
        mat"""
            addpath($WEP_path)
            nxx = double($nx);
            nzz = double($nz);
            options = [];
            options.delta = double($delta);
            options.wg = $waveguide_str;
            matlab_nep = nep_wg_generator(nxx, nzz, options);

            Ixz = eye(nxx*nzz+2*nzz);
            $M_m = NaN(nxx*nzz+2*nzz, nxx*nzz+2*nzz);
            for i = 1:(nxx*nzz+2*nzz)
              $M_m(:,i) = matlab_nep.M($gamma, Ixz(:,i));
            end
        """
        println("  -- Matlab printouts end --")

        println("Difference M_m(γ) - M(γ) = ", opnorm(Matrix(M_m-M_j)))
        println("Relative difference opnorm(M_m(γ) - M(γ))/opnorm(M(γ)) = ", opnorm(Matrix(M_m-M_j))/opnorm(Matrix(M_j)))
    end
    println("\n--- End Full Matrix FD SPMF against MATLAB ---\n")
end

###########################################################
# Test the square-root implementations. Both reference and negative-imaginary-part.
function debug_sqrt_schur(n::Integer)
    println("\n\n--- Debugging square root implementations ---\n")
    A = rand(n,n) + 0.1im*rand(n,n);
    sqrtA = sqrt(A);
    sqrtA2 = sqrt_schur(A);

    println("Relative error between sqrtm and Schur-fact: ", opnorm(sqrtA-sqrtA2)/opnorm(sqrtA))
    println("Relative error between Schur-fact² and A: ", opnorm(A - sqrtA2^2)/opnorm(A))

    sqrtA3 = sqrt_schur_pos_imag(A);
    println("Relative error between Schur-fact-pos-imag² and A: ", opnorm(A - sqrtA3^2)/opnorm(A))
    v,_ = eigen(sqrtA3)
    test_var = zeros(n)
    TOL = 1e-15;
    for i = 1:n
        test_var[i] = (sign(imag(v[i])) > 0) || abs(imag(v[i])) < TOL
    end
    if sum(test_var) == n
        println("All eigenvalues have negative real part or absolut value smaller than ", TOL)
    else
        println("Eigenvalues with negative real part and absolut value larger than ", TOL, " :")
        for i = 1:n
            if (test_var[i] == 0)
                println("   ", v[i])
            end
        end
    end
    println("\n--- End square root implementations ---\n")
end

###########################################################
# Test computation of derivative <d> of sqrt(ax^2 + bx + c) in <x>
function debug_sqrt_derivative()
    println("\n\n--- Debugging derivatives of square root of polynomials against MATLAB symbolic ---\n")
    a = 2*rand()
    b = 2*pi*rand(ComplexF64)
    c = 1.67*rand()
    d_vec = [0 1 2 3 4 11 19 20 21 22 30 35 45 60] #Factorial for Int64 overflows at 21!
    x = 25*rand()

        WEP_path = string(@__DIR__, "/WEP-ref")
        println("  -- Matlab printouts start --")
        mat"""
            addpath($WEP_path)
            $der_val = sqrt_derivative_test($a,$b,$c, $d_vec, $x);
        """
        println("  -- Matlab printouts end --")

    julia_der_vec = sqrt_derivative(a,b,c, maximum(d_vec), x)
    for i = 1:size(d_vec,2)
        d = d_vec[i]
        println("Derivative number d = ", d)
        println("    MATLAB symbolic = ", der_val[i])
        julia_der = julia_der_vec[d+1]
        println("    Implemented recursion = ", julia_der)
        println("    Relative error = ", abs(der_val[i]-julia_der)/abs(julia_der))
    end

    println("\n--- End derivatives of square root of polynomials against MATLAB symbolic ---\n")
end

###########################################################
# Test the full generated native-WEP system-matrix against the SPMF-matrix (Mlincomb)
function debug_Mlincomb_FD_WEP(nx::Integer, nz::Integer, delta::Number)
    println("\n\n--- Debugging Mlincomb (Full Matrix FD WEP-native-format against SPMF) ---\n")
    if(nx > 40 || nz > 40)
        @warn "This debug is 'naive' and might be slow for the discretization used."
    end

    γ = -rand() - 1im*rand()
    gamma = γ
    n = nx*nz+2*nz
    V = rand(ComplexF64, n, 4)

    for waveguide = ["TAUSCH", "JARLEBRING"]
        println("\n")
        println("Testing full matrix M for waveguide: ", waveguide)

        nep_j = nep_gallery(WEP, nx = nx, nz = nz, benchmark_problem = waveguide, neptype = "WeP", delta = delta)
        nep_SPMF = nep_gallery(WEP, nx = nx, nz = nz, benchmark_problem = waveguide, neptype = "SpmF", delta = delta)

        for d = [0 1 2 5]
            M_j = zeros(ComplexF64, n, n)
            for i = 1:n
                V = zeros(ComplexF64, n)
                V[i] = 1
                M_j[:,i] += compute_Mlincomb(nep_j, γ, V, [1], d)
            end

            M_SPMF = zeros(ComplexF64, n, n)
            for i = 1:n
                V = zeros(ComplexF64, n)
                V[i] = 1
                M_SPMF[:,i] += compute_Mlincomb(nep_SPMF, γ, V, [1], d)
            end
            println("Derivative d = ", d)
            println("    Difference M_SPMF(γ) - M(γ) = ", opnorm(Matrix(M_SPMF-M_j)))
            println("    Relative difference opnorm(M_m(γ) - M(γ))/opnorm(M(γ)) = ", opnorm(Matrix(M_SPMF-M_j))/opnorm(Matrix(M_j)))

            v_j = compute_Mlincomb(nep_j, γ, V)
            v_SPMF = compute_Mlincomb(nep_SPMF, γ, V)
            println("    Relative difference on random set of 4 vectors = ", norm(v_j - v_SPMF)/norm(v_j))
        end

    end
    println("\n--- End Mlincomb (Full Matrix FD WEP-native-format against SPMF) ---\n")
end

###########################################################
# Test the Schur-complement of native-WEP against MATLAB code
function matlab_debug_Schur_WEP_FD(nx::Integer, nz::Integer, delta::Number)
    println("\n\n--- Debugging Schur-complement of native-WEP against MATLAB ---\n")
    if(nx > 50 || nz > 50)
        @warn "This debug is 'naive' and might be slow for the discretization used."
    end

    γ = -rand() - 1im*rand()
    gamma = γ
    n = nx*nz;

    for waveguide = ["TAUSCH", "JARLEBRING"]
        println("\n")
        println("Testing Schur-complement for waveguide: ", waveguide)

        nep_j = nep_gallery(WEP, nx = nx, nz = nz, benchmark_problem = waveguide, neptype = "WEp", delta = delta)
        Schur_fun = SchurMatVec(nep_j, γ)
        Schur_j = zeros(ComplexF64, n, n)
        for i = 1:n
            V = zeros(ComplexF64, n)
            V[i] = 1
            Schur_j[:,i] += Schur_fun* V
        end

        Schur_jj = Matrix(construct_WEP_schur_complement(nep_j, γ))

        if waveguide == "JARLEBRING"
            waveguide_str = "CHALLENGE"
        else
            waveguide_str = waveguide
        end

        println("  -- Matlab printouts start --")
        WEP_path = string(@__DIR__, "/WEP-ref")
        mat"""
            addpath($WEP_path)
            nxx = double($nx);
            nzz = double($nz);
            options = [];
            options.delta = double($delta);
            options.wg = $waveguide_str;
            matlab_nep = nep_wg_generator(nxx, nzz, options);

            Ixz = eye(nxx*nzz);
            $Schur_m = NaN(nxx*nzz, nxx*nzz);
            for i = 1:(nxx*nzz)
              $Schur_m(:,i) = matlab_nep.S($gamma, Ixz(:,i));
            end
         """
        println("  -- Matlab printouts end --")

        norm_diff = opnorm(Schur_m-Schur_j)
        println("Difference Schur_m(γ) - Schur_j(γ) = ", norm_diff)
        println("Relative difference opnorm(Schur_m(γ) - Schur_j(γ))/opnorm(Schur_j(γ)) = ", norm_diff/opnorm(Schur_j))
        norm_diff = opnorm(Schur_m-Schur_jj)
        println("Difference Schur_m(γ) - Schur_jj(γ) = ", norm_diff)
        println("Relative difference opnorm(Schur_m(γ) - Schur_jj(γ))/opnorm(Schur_jj(γ)) = ", norm_diff/opnorm(Schur_jj))
    end
    println("\n--- End Schur-complement of native-WEP against MATLAB ---\n")
end

###########################################################
# Test the FFT-based Sylvester solver against naive solver
function fft_debug_mateq(nx::Integer, nz::Integer, delta::Number)
    println("\n\n--- Debugging FFT-Sylvester ---\n")
    γ = -rand(ComplexF64)
    gamma = γ
    C = rand(ComplexF64, nz, nx);
    waveguide = "JARLEBRING"

        if waveguide == "JARLEBRING"
            waveguide_str = "CHALLENGE"
        else
            waveguide_str = waveguide
        end
        println("  -- Matlab printouts start --")
        WEP_path = string(@__DIR__, "/WEP-ref")
         mat"""
            addpath($WEP_path)
            nxx = double($nx);
            nzz = double($nz);
            options = [];
            options.delta = double($delta);
            options.wg = $waveguide_str;
            nep = nep_wg_generator(nxx, nzz, options);

            CC = double($C);
            sigma = double($gamma);

            kk = mean(nep.K(:));
            K = nep.K - kk;

            $X_m = fft_wg( CC, sigma, kk, nep.hx, nep.hz );
         """
        println("  -- Matlab printouts end --")

    K, hx, hz, Km, Kp = generate_wavenumber_fd( nx, nz, waveguide, delta)
    Dxx, Dzz, Dz = generate_fd_interior_mat( nx, nz, hx, hz)

    k_bar = mean(K)

    A = Matrix(Dzz + 2*γ*Dz + (γ^2+k_bar)*I)
    B = complex(Matrix(Dxx))

    println("\nBuilt-in Sylvester solver (X_jj)")
    X_jj = @time sylvester(A,B,-C)
    println("Relative residual norm = ", opnorm(A*X_jj+X_jj*B-C)/opnorm(C))

    println("FFT-based Sylvester solver for WG (X_j)")
    X_j::Matrix{ComplexF64} = copy(C)
    @time solve_wg_sylvester_fft!( X_j, γ, k_bar, hx, hz )
    println("Relative residual norm = ", opnorm(A*X_j+X_j*B-C)/opnorm(C))

    println("MATLAB implemented FFT-based Sylvester solver for WG")
    println("Relative residual norm = ", opnorm(A*X_m+X_m*B-C)/opnorm(C))

    println("\nRelative difference opnorm(X_m - X_j)/opnorm(X_j) = ", opnorm(X_m - X_j)/opnorm(X_j))
    println("Relative difference opnorm(X_m - X_jj)/opnorm(X_jj) = ", opnorm(X_m - X_jj)/opnorm(X_jj))
    println("Relative difference opnorm(X_j - X_jj)/opnorm(X_j) = ", opnorm(X_j - X_jj)/opnorm(X_j))

    println("\n--- End FFT-Sylvester ---\n")
end

###########################################################
# Test the Sylvester SMW
function debug_Sylvester_SMW_WEP(nx::Integer, nz::Integer, delta::Number, N::Integer)
    println("\n\n--- Debugging Sylvester SMW for WEP ---\n")

    γ = -rand() - 1im*rand()
    gamma = γ
    σ = γ

    C = rand(Float64, nz, nx)

    for waveguide = ["TAUSCH", "JARLEBRING"]
        println("\n")
        println("Testing Sylvester SMW for waveguide: ", waveguide)

        nep = nep_gallery(WEP, nx = nx, nz = nz, benchmark_problem = waveguide, neptype = "wEp", delta = delta)

        println("  Generate SMW-matrix")
        M_j = @time generate_smw_matrix(nep, N, σ)
        M_jj = Matrix(M_j)

        println("  Solve SMW-problem")
        X_j = @time solve_smw(nep, M_j, C, σ)

        if waveguide == "JARLEBRING"
            waveguide_str = "CHALLENGE"
        else
            waveguide_str = waveguide
        end

        println("  -- Matlab printouts start --")
        WEP_path = string(@__DIR__, "/WEP-ref")
        mat"""
            addpath($WEP_path)
            nxx = double($nx);
            nzz = double($nz);
            nn = nzz;
            NN = double($N);
            CC = double($C);
            option = [];
            options.delta = double($delta);
            options.wg = $waveguide_str;
            nep = nep_wg_generator(nxx, nzz, options);

            sigma = double($gamma);

            %OBS: The below (MATLAB) code is copied from 'main.m' in the public code
            % CONSTRUCT THE PRECONDITIONER
            kk = mean(nep.K(:));
            K = nep.K - kk;

            Linv=@(X) fft_wg( X, sigma, kk, nep.hx, nep.hz );
            Pm_inv=@(x) -nep.Pm_inv(sigma, x);              %OBS! The minus sign!
            Pp_inv=@(x) -nep.Pp_inv(sigma, x);              %OBS! The minus sign!

            dd1 = nep.d1/nep.hx^2;
            dd2 = nep.d2/nep.hx^2;
            $M_m = generate_smw_matrix( nn, NN, Linv, dd1, dd2, Pm_inv, Pp_inv, K, false );

            $X_m = solve_smw( $M_m, CC, Linv, dd1, dd2, Pm_inv, Pp_inv, K );
        """
        println("  -- Matlab printouts end --")

        println("    Difference SMW_matrix_m - SMW_matrix_j = ", opnorm(M_m - M_jj))
        println("    Relative difference opnorm(SMW_matrix_m - SMW_matrix_j)/opnorm(SMW_matrix_j) = ", opnorm(M_m - M_jj)/opnorm(M_jj))
        println("\n    X is the solution to a Sylvester-SMW system")
        println("    Difference X_m - X_j = ", opnorm(X_m - X_j))
        println("    Relative difference opnorm(X_m - X_j)/opnorm(X_j) = ", opnorm(X_m - X_j)/opnorm(X_j))

    end
    println("\n--- End Debugging Sylvester SMW for WEP ---\n")
end

###########################################################
# Test convergence for the preconditioner
function debug_WEP_FD_preconditioner(delta::Number)
    println("\n\n--- Debugging WEP_FD preconditioner ---\n")

    nz = 3*3*5
    nx = nz + 4

    γ = -rand(ComplexF64)
    gamma = γ

    for waveguide = ["TAUSCH", "JARLEBRING"]
        println("\n")
        println("Preconditioner in GMRES for waveguide: ", waveguide)

        nep = nep_gallery(WEP, nx = nx, nz = nz, benchmark_problem = waveguide, neptype = "WeP", delta = delta)

        bb = rand(ComplexF64, nx*nz)
        precond = wep_generate_preconditioner(nep, 45, γ)
        Schur_fun = SchurMatVec(nep, γ)
        bbb = ldiv!(precond, copy(Schur_fun*bb))
        println("    Preconditioner * (Schur complement * b): Relative residual norm = ", norm(bbb - bb)/norm(bb))
        bb = rand(ComplexF64, nx*nz)
        bbb = ldiv!(precond, copy(bb))
        bbb = Schur_fun * bbb
        println("    Schur complement * (Preconditioner * b): Relative residual norm = ", norm(bbb - bb)/norm(bb))

        b = rand(ComplexF64, nx*nz+2*nz)
        for N = [1, 3, 9, 15, 45]
            println("    Testing for n = ", nz, " and N = ", N)
            precond = @time wep_generate_preconditioner(nep, N, γ)

            gmres_kwargs = ((:maxiter,100), (:restart,100), (:log,true), (:Pl,precond), (:tol, 1e-13), (:verbose,true))
            linsolvercreator=GalleryWaveguide.WEPLinSolverCreator(solver_type=:gmres,kwargs=gmres_kwargs)
            wep_solver = create_linsolver(linsolvercreator, nep, γ)
            x = lin_solve(wep_solver, b)

            println("    Relative residual norm = ", norm(compute_Mlincomb(nep, gamma, x) - b)/norm(b))

        end

    end
    println("\n--- End WEP_FD preconditioner  ---\n")
end

###########################################################
# Test to compute eigenvalues and eigenvectors. Test preconditioned RESINV as in Ringh for native-WEP and MATLAB code, and pre-factorized-matrix RESINV for SPMF.
function matlab_debug_eigval_comp_WEP_FD_and_SPMF(nz::Integer, N::Integer, delta::Number)
    println("\n\n--- Debugging eigenvalue computations for WEP_FD and SPMF against MATLAB ---\n")

    nx = nz + 4;

    for waveguide = ["TAUSCH", "JARLEBRING"]
        println("\n")
        println("Testing eigenvalue computations for waveguide: ", waveguide)

        nep_j_WEPFD =    nep_gallery(WEP, nx = nx, nz = nz, benchmark_problem = waveguide, neptype = "WeP", delta = delta)
        nep_j_SPMF  =    nep_gallery(WEP, nx = nx, nz = nz, benchmark_problem = waveguide, neptype = "SpmF", delta = delta)
        nep_j_SPMF_pre = nep_gallery(WEP, nx = nx, nz = nz, benchmark_problem = waveguide, neptype = "SpmF_prE", delta = delta)

        if waveguide == "JARLEBRING"
            waveguide_str = "CHALLENGE"
            γ = -3-3.5im
        else
            waveguide_str = waveguide
            γ = -0.015-5.1im
        end
        gamma = γ

        println("    Generating preconditioner")
        precond = @time wep_generate_preconditioner(nep_j_WEPFD, N, γ)
        gmres_kwargs = ((:maxiter,100), (:restart,100), (:log,false), (:Pl,precond), (:tol, 1e-13))

        println("    Compute for WEP_FD")
        eigval_j_WEPFD = NaN
        eigvec_j_WEPFD = NaN
        try
            eigval_j_WEPFD, eigvec_j_WEPFD = @time resinv(nep_j_WEPFD, logger=1, λ=γ, maxit = 30,
                                                          tol = 1e-10, v=ones(ComplexF64,nx*nz+2*nz), c=zeros(ComplexF64,nx*nz+2*nz),
                                                          linsolvercreator=GalleryWaveguide.WEPLinSolverCreator(solver_type=:gmres,kwargs=gmres_kwargs)
                                                          );
        catch err
            # Only catch NoConvergence
            isa(err, NoConvergenceException) || rethrow(err)
            println("No convergence because: ", err.msg)
            # still access the approximations
            eigval_j_WEPFD = err.λ
            eigvec_j_WEPFD = err.v
        end

        println("    Compute for SPMF")
        eigval_j_SPMF = NaN
        eigvec_j_SPMF = NaN
        try
            eigval_j_SPMF, eigvec_j_SPMF = @time resinv(nep_j_SPMF, logger=1, λ=γ, maxit = 30, tol = 1e-10, v=ones(ComplexF64,nx*nz+2*nz), c=zeros(ComplexF64,nx*nz+2*nz))
        catch err
            # Only catch NoConvergence
            isa(err, NoConvergenceException) || rethrow(err)
            println("No convergence because: ", err.msg)
            # still access the approximations
            eigval_j_SPMF = err.λ
            eigvec_j_SPMF = err.v
        end

        println("    Compute for SPMF (pre-Schur-fact)")
        eigval_j_SPMF_pre = NaN
        eigvec_j_SPMF_pre = NaN
        try
            eigval_j_SPMF_pre, eigvec_j_SPMF_pre = @time resinv(nep_j_SPMF_pre, logger=1, λ=γ, maxit = 30, tol = 1e-10, v=ones(ComplexF64,nx*nz+2*nz), c=zeros(ComplexF64,nx*nz+2*nz))
        catch err
            # Only catch NoConvergence
            isa(err, NoConvergenceException) || rethrow(err)
            println("No convergence because: ", err.msg)
            # still access the approximations
            eigval_j_SPMF_pre = err.λ
            eigvec_j_SPMF_pre = err.v
        end

        println("  -- Matlab printouts start --")
        WEP_path = string(@__DIR__, "/WEP-ref")
        mat"""
            addpath($WEP_path)

            nzz = double($nz);
            lambda = double($gamma);
            NN = double($N);

            [$eigval_m, $eigvec_m] = main_func_WEP(nzz, NN, lambda, $delta, $waveguide_str)

        """
        println("  -- Matlab printouts end --")

        println("    Difference between WEP_FD and SPMF computed eigenvalue = ", abs(eigval_j_WEPFD - eigval_j_SPMF))
        println("    Relative difference between WEP_FD and SPMF computed eigenvalue = ", abs(eigval_j_WEPFD - eigval_j_SPMF)/abs(eigval_j_WEPFD))
        println("    Difference between WEP_FD and SPMF computed eigenvectors = ", norm(eigvec_j_WEPFD/eigvec_j_WEPFD[1] - eigvec_j_SPMF/eigvec_j_SPMF[1]))
        println("    Relative difference between WEP_FD and SPMF computed eigenvalue = ", norm(eigvec_j_WEPFD/eigvec_j_WEPFD[1] - eigvec_j_SPMF/eigvec_j_SPMF[1])/norm(eigvec_j_WEPFD/eigvec_j_WEPFD[1]))
        println("")
        println("    Difference between WEP_FD and MATLAB computed eigenvalue = ", abs(eigval_j_WEPFD - eigval_m))
        println("    Relative difference between WEP_FD and MATLAB computed eigenvalue = ", abs(eigval_j_WEPFD - eigval_m)/abs(eigval_m))
        println("    Difference between WEP_FD and MATLAB computed eigenvectors = ", norm(eigvec_j_WEPFD/eigvec_j_WEPFD[1] - eigvec_m/eigvec_m[1]))
        println("    Relative difference between WEP_FD and MATLAB computed eigenvalue = ", norm(eigvec_j_WEPFD/eigvec_j_WEPFD[1] - eigvec_m/eigvec_m[1])/norm(eigvec_m/eigvec_m[1]))
        println("")
        println("    Difference between MATLAB and SPMF computed eigenvalue = ", abs(eigval_m - eigval_j_SPMF))
        println("    Relative difference between MATLAB and SPMF computed eigenvalue = ", abs(eigval_m - eigval_j_SPMF)/abs(eigval_m))
        println("    Difference between MATLAB and SPMF computed eigenvectors = ", norm(eigvec_m/eigvec_m[1] - eigvec_j_SPMF/eigvec_j_SPMF[1]))
        println("    Relative difference between MATLAB and SPMF computed eigenvalue = ", norm(eigvec_m/eigvec_m[1] - eigvec_j_SPMF/eigvec_j_SPMF[1])/norm(eigvec_m/eigvec_m[1]))
        println("")
        println("    Difference between SPMF and SPMF-pre computed eigenvalue = ", abs(eigval_j_SPMF_pre - eigval_j_SPMF))
        println("    Relative difference between SPMF and SPMF-pre computed eigenvalue = ", abs(eigval_j_SPMF_pre - eigval_j_SPMF)/abs(eigval_j_SPMF))
        println("    Difference between SPMF and SPMF-pre computed eigenvectors = ", norm(eigvec_j_SPMF_pre/eigvec_j_SPMF_pre[1] - eigvec_j_SPMF/eigvec_j_SPMF[1]))
        println("    Relative difference between SPMF and SPMF-pre computed eigenvalue = ", norm(eigvec_j_SPMF_pre/eigvec_j_SPMF_pre[1] - eigvec_j_SPMF/eigvec_j_SPMF[1])/norm(eigvec_j_SPMF/eigvec_j_SPMF[1]))
        println("")
        println("    WEP_FD converged to the eigenvalue = ", eigval_j_WEPFD)

    end
    println("\n--- End eigenvalue computations of WEP_FD and SPMF against MATLAB ---\n")
end

end


text:
# Implements a pre-computed table of recursion coefficients

"""
    struct LegendreNormCoeff{N<:AbstractLegendreNorm,T<:Real} <: AbstractLegendreNorm

Precomputed recursion relation coefficients for the normalization `N` and value type
`T`.

# Example
```jldoctest
julia> LegendreNormCoeff{LegendreSphereNorm,Float64}(1)
LegendreNormCoeff{LegendreSphereNorm,Float64} for lmax = 1, mmax = 1 with coefficients:
    μ: [0.0, 1.22474]
    ν: [1.73205, 2.23607]
    α: [0.0 0.0; 1.73205 0.0]
    β: [0.0 0.0; -0.0 0.0]
```
"""
struct LegendreNormCoeff{N<:AbstractLegendreNorm,T<:Real} <: AbstractLegendreNorm
    μ::Vector{T}
    ν::Vector{T}
    α::Matrix{T}
    β::Matrix{T}

    function LegendreNormCoeff{N,T}(lmax::Integer, mmax::Integer) where {N,T}
        (lmax ≥ 0) || throw(DomainError(lmax, "lmax must be positive"))
        (0 ≤ mmax ≤ lmax) || throw(DomainError(mmax, "mmax must be bounded in 0 to $lmax"))

        μ = zeros(T, mmax+1)
        ν = zeros(T, mmax+1)
        α = zeros(T, lmax+1, mmax+1)
        β = zeros(T, lmax+1, mmax+1)

        @inbounds for m in 0:mmax
            μ[m+1] = m == 0 ? zero(T) : coeff_μ(N(), T, m)
            # N.B.: Need access to mmax+1 but will never use m==0 term, so offset storage
            ν[m+1] = coeff_ν(N(), T, m+1)

            for l in (m+1):lmax
                α[l+1,m+1] = coeff_α(N(), T, l, m)
                β[l+1,m+1] = coeff_β(N(), T, l, m)
            end
        end

        return new(μ, ν, α, β)
    end

    function LegendreNormCoeff{N,T1}(norm::LegendreNormCoeff{N,T2}) where {N,T1,T2}
        return new(convert(Vector{T1}, norm.μ),
                   convert(Vector{T1}, norm.ν),
                   convert(Matrix{T1}, norm.α),
                   convert(Matrix{T1}, norm.β))
    end
end

LegendreNormCoeff{N,T}(lmax::Integer) where {N,T} = LegendreNormCoeff{N,T}(lmax, lmax)

Base.convert(::Type{LegendreNormCoeff{N,T}}, norm::LegendreNormCoeff{N,T}) where {N,T} =
    norm
Base.convert(::Type{LegendreNormCoeff{N,T}}, norm::LegendreNormCoeff{N,S}) where {N,T,S} =
    LegendreNormCoeff{N,T}(norm)

Base.eltype(::Type{LegendreNormCoeff{N,T}}) where {N,T} = T

# Improve printing somewhat
Base.show(io::IO, norm::LegendreNormCoeff{N,T}) where {N,T} =
    print(io, LegendreNormCoeff, "{$N,$T}")
function Base.show(io::IO, ::MIME"text/plain", N::LegendreNormCoeff)
    lmax,mmax = size(N.α) .- 1
    println(io, N, " for lmax = $lmax, mmax = $mmax with coefficients:")
    io′ = IOContext(io, :compact => true)
    println(io′, "    μ: ", N.μ)
    println(io′, "    ν: ", N.ν)
    println(io′, "    α: ", N.α)
    println(io′, "    β: ", N.β)
end

# Implements the legendre interface for the normalization

@inline function
initcond(::LegendreNormCoeff{N}, ::Type{T}) where {N<:AbstractLegendreNorm, T}
    return initcond(N(), T)
end

@propagate_inbounds function
coeff_μ(norm::LegendreNormCoeff, ::Type{T}, l::Integer) where T
    return norm.μ[l+1]
end

@propagate_inbounds function
coeff_ν(norm::LegendreNormCoeff, ::Type{T}, l::Integer) where T
    # N.B.: Storage is offset by 1 compared to other arrays. See constructor.
    return norm.ν[l]
end

@propagate_inbounds function
coeff_α(norm::LegendreNormCoeff, ::Type{T}, l::Integer, m::Integer) where T
    return norm.α[l+1,m+1]
end

@propagate_inbounds function
coeff_β(norm::LegendreNormCoeff, ::Type{T}, l::Integer, m::Integer) where T
    return norm.β[l+1,m+1]
end

function boundscheck_hook(norm::LegendreNormCoeff, lmax, mmax)
    @noinline _throw(α, lmax, mmax) = throw(BoundsError(α, (lmax+1, mmax+1)))
    lmax′,mmax′ = size(norm.α)
    (lmax < lmax′ && mmax < mmax′) || _throw(norm.α, lmax, mmax)
    nothing
end

ext: jl


text:
export bpfcall

@generated function _bpfcall(::Val{cmd}, ::Type{rettype}, ::Type{argtypes}, args::Vararg{Any}) where {cmd,rettype,argtypes}
    Context() do ctx
        T_ret = convert(LLVMType, rettype; ctx)
        T_jlptr = convert(LLVMType, Ptr{Cvoid}; ctx)
        T_ptr_i8 = LLVM.PointerType(LLVM.Int8Type(ctx))

        outer_args = filter(arg->!(arg <: RTMap), args)
        T_outer_args = LLVMType[convert(LLVMType, arg; ctx) for arg in outer_args]

        llvm_f, _ = create_function(T_ret, T_outer_args)
        mod = LLVM.parent(llvm_f)

        outer_args_ex = Expr[]
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            inner_args = LLVM.Value[]
            T_inner_args = LLVMType[]
            jlidx = 1
            pidx = 1
            for (idx, arg) in enumerate(argtypes.parameters)
                if arg <: RTMap
                    map_gv = _genmap!(mod, arg, ctx)
                    push!(inner_args, map_gv)
                    push!(T_inner_args, llvmtype(map_gv))
                else
                    llvm_arg = parameters(llvm_f)[pidx]
                    T_arg = T_outer_args[pidx]
                    if arg <: Ptr
                        llvm_arg = inttoptr!(builder, llvm_arg, T_ptr_i8)
                        T_arg = llvmtype(llvm_arg)
                    end
                    push!(inner_args, llvm_arg)
                    push!(T_inner_args, T_arg)
                    push!(outer_args_ex, Expr(:ref, :args, jlidx))
                    pidx += 1
                end
                jlidx += 1
            end

            ft = LLVM.FunctionType(T_ret, T_inner_args)
            ftp = LLVM.PointerType(ft)
            f = inttoptr!(builder, ConstantInt(cmd; ctx), ftp)
            value = call!(builder, f, inner_args)
            ret!(builder, value)
        end
        #outer_args_ex = Expr(:tuple, outer_args_ex...)
        call_function(llvm_f, rettype, Base.to_tuple_type(outer_args), outer_args_ex...)
    end
end
@inline bpfcall(cmd::API.BPFHelper, RT, AT, args...) =
    _bpfcall(Val(Int(cmd)), RT, AT, args...)
@inline bpfcall(cmd::API.BPFHelper, RT) = _bpfcall(Val(Int(cmd)), RT, Tuple{})
@inline bpfcall(cmd::API.BPFHelper) = _bpfcall(Val(Int(cmd)), Cvoid, Tuple{})


text:
    using DifferentialEquations
    using OrdinaryDiffEq
    using DiffEqBase
    using Sundials
    using ODEInterfaceDiffEq
    using Turing
    using MCMCChains

        # du -> Derivative
        # u -> State at time t
        # p-> paramter vector
        # t-> time as tuple (init, end)

        function Hes1ExampleODE!(du,u,p,t)
            m, p1, p2 = u;

            k1, P0, v, h, inps = p;

        du[1] = dm = -0.03*m + (1/(1+(p2/P0)^h))*inps;
        du[2] = dp1 = -0.03*p1 + v*m - k1*p1;
        du[3] = dp2 = -0.03*p2 + k1*p1;

        end

        # p -> parameter vector (plus inducer vector at the end)
        # I -> initial Y0 vector

        function Hes1ExampleSteadyState(p,I)

            m, p1, p2 = I;
               return(I)

        end

        # ts -> time vector going from t=0 to t=end every 1
        # p -> parameter vector
        # sp -> vector with switching times (t=0 and t=end needs to be included). Length of vector will be length(inducer steps)+1
        # inputs -> Vecor with the inducers (if more than 1 inducer in the system, the shape will be a1,b1,c1,..., a2,b2,c2,... aend,bend,cend,... where each leter represents a different inducer)
        # ivss -> Initial Y0 vector
        # pre -> Vector with the inducers in the ON. It can be empty

        @model function Hes1Example_TuringModel(data,prob1)

            p1 ~ Truncated(Distributions.Normal(0.05, 0.025), 0.0, 0.1)
            p2 ~ Truncated(Distributions.Normal(1.0, 0.5), 0.0, 2.0)
            p3 ~ Truncated(Distributions.Normal(0.05, 0.025), 0.0, 0.1)
            p4 ~ Truncated(Distributions.Normal(5.5, 2.25), 1.0, 10.0)

            p = [p1,p2,p3,p4,];

            Neq = 3;

            for exp in 1:1

                ts=collect(0.:1:data["finalTime"][exp])
                sp = data["switchT"][exp];
                inputs = data["uInd"][exp];
                ivss = data["y0"][exp];
                pre=data["preInd"][exp];

                maxtime = length(ts);
                Nsp = length(sp);
                Nevents = length(sp)-1;

                if typeof(p) == Array{Float64, 1}
                    if pre != []
                        pSS = vcat(p, pre);
                    else
                        pSS = p;
                    end
                else
                    if pre != []
                        pSS = vcat([p[i].value for i in 1:4], pre);
                    else
                        pSS = [p[i].value for i in 1:4];
                    end
                end

                final = zeros(maxtime, Neq);
                finalBay2 = Array{Any,2}(undef, 1, maxtime)

                y_al = Hes1ExampleSteadyState(pSS,ivss) # Calculation of initial guesses for steady state
                y0 = y_al;
                initialV = y0;
                i = 1;

                for q in collect(1:Nevents)

                    lts = length(ts[(sp[q]+1):sp[q+1]+1]);  # General way to define the number of elements in each event series

                    Tevent = ts[(sp[q]+1):sp[q+1]+1];  # General way to extract the times of each event
                    I = inputs[i:(i+(0))];
                    pSte = vcat(p, I);

                    if q == 1
                        if typeof(p) == Array{Float64, 1}
                            prob = ODEProblem(Hes1ExampleODE!,initialV,(ts[(sp[q]+1)],ts[sp[q+1]+1]),pSte);
                            part1 = DifferentialEquations.solve(prob, DifferentialEquations.Tsit5(),reltol=1.0e-9,abstol=1.0e-9,saveat=1);
                        else
                            prob = remake(prob1, p=pSte)
                            prob = remake(prob, u0=initialV)
                            prob = remake(prob, tspan = convert.(Int, (ts[(sp[q]+1)],ts[sp[q+1]+1])))
                            part1 = DifferentialEquations.solve(prob,DifferentialEquations.Tsit5(),reltol=1.0e-9,abstol=1.0e-9,saveat=1)
                        end
                    else
                        if typeof(p) == Array{Float64, 1}
                            prob = ODEProblem(Hes1ExampleODE!,initialV,(ts[(sp[q]+1)],ts[sp[q+1]+1]),pSte);
                            part1 = DifferentialEquations.solve(prob, DifferentialEquations.Tsit5(),reltol=1.0e-9,abstol=1.0e-9,saveat=1);
                        else
                            prob = remake(prob1, p=pSte)
                            prob = remake(prob, u0=initialV)
                            prob = remake(prob, tspan = (ts[(sp[q]+1)],ts[sp[q+1]+1]))
                            part1 = DifferentialEquations.solve(prob,DifferentialEquations.Tsit5(),reltol=1.0e-9,abstol=1.0e-9,saveat=1)
                        end
                    end

                    if typeof(p) == Array{Float64, 1}
                        initialV = part1[lts];
                    else
                        initialV = [part1[:,lts][h].value for h in 1:3];
                    end

                    i+=1;

                    if typeof(p) == Array{Float64, 1}
                        for d in collect((sp[q]+1):(sp[q]+lts))
                            final[d,:] = part1[d-sp[q]];
                        end
                    else
                        for d in collect((sp[q]+1):(sp[q]+lts))
                            finalBay2[1,d] = part1[1, d-sp[q]];

                        end
                    end

                end

                for ob in 1:1
                    if typeof(p) == Array{Float64, 1}
                        obs = final[convert.(Int,round.(data["tsamps"][exp])).+1, [1]];
                    else
                        obs = finalBay2[ob, convert.(Int,round.(data["tsamps"][exp])).+1];
                    end

                    for dat in 1:length(obs)
                        data["DataMean"][ob][dat] ~ Normal(obs[dat], data["DataError"][ob][1][dat])
                    end
                end

            end

        end


text:
module Dashlorikeet
using Dash

const resources_path = realpath(joinpath( @__DIR__, "..", "deps"))
const version = "0.0.1"

include("dashlorikeet.jl")

function __init__()
    DashBase.register_package(
        DashBase.ResourcePkg(
            "dashlorikeet",
            resources_path,
            version = version,
            [
                DashBase.Resource(
    relative_package_path = "dashlorikeet.min.js",
    external_url = "https://unpkg.com/dashlorikeet@0.0.1/dashlorikeet/dashlorikeet.min.js",
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dashlorikeet.min.js.map",
    external_url = "https://unpkg.com/dashlorikeet@0.0.1/dashlorikeet/dashlorikeet.min.js.map",
    dynamic = true,
    async = nothing,
    type = :js
)
            ]
        )

    )
end
end


text:
include("visualize.jl")

q = [[0.0, 0.5 + 0.5 * sin(1.0 * π * (t - 1) / 10.0), 0.0, 0.5] for t = 1:11]

vis = Visualizer()
render(vis)
visualize!(vis, q, 0.1)


text:
"""
Composition-Rejection Direct Method (DirectCR), implementation combining
features from the original article:
*A constant-time kinetic Monte Carlo algorithm for simulation of large biochemical reaction networks*,
by A. Slepoy, A.P. Thompson and S.J. Plimpton, J. Chem. Phys, 128, 205101 (2008).
and
*Efficient Formulations for Exact Stochastic Simulation of Chemical Systems*,
by S. Mauch and M. Stalzer, ACM Trans. Comp. Biol. and Bioinf., 8, No. 1, 27-35 (2010).
"""

const MINJUMPRATE = 2.0^exponent(1e-12)

mutable struct DirectCRJumpAggregation{T,S,F1,F2,RNG,DEPGR,U<:PriorityTable,W<:Function} <: AbstractSSAJumpAggregator
    next_jump::Int
    next_jump_time::T
    end_time::T
    cur_rates::Vector{T}
    sum_rate::T
    ma_jumps::S
    rates::F1
    affects!::F2
    save_positions::Tuple{Bool,Bool}
    rng::RNG
    dep_gr::DEPGR
    minrate::T
    maxrate::T   # initial maxrate only, table can increase beyond it!
    rt::U
    ratetogroup::W
  end

function DirectCRJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T,
                                      maj::S, rs::F1, affs!::F2, sps::Tuple{Bool,Bool},
                                      rng::RNG; num_specs, dep_graph=nothing,
                                      minrate=convert(T,MINJUMPRATE), maxrate=convert(T,Inf),
                                      kwargs...) where {T,S,F1,F2,RNG}

    # a dependency graph is needed and must be provided if there are constant rate jumps
    if dep_graph === nothing
        if (get_num_majumps(maj) == 0) || !isempty(rs)
            error("To use ConstantRateJumps with the DirectCR algorithm a dependency graph must be supplied.")
        else
            dg = make_dependency_graph(num_specs, maj)
        end
    else
        dg = dep_graph
    end

    # make sure each jump depends on itself
    add_self_dependencies!(dg)

    # mapping from jump rate to group id
    minexponent = exponent(minrate)

    # use the largest power of two that is <= the passed in minrate
    minrate = 2.0^minexponent
    ratetogroup = rate -> priortogid(rate, minexponent)

    # construct an empty initial priority table -- we'll overwrite this in init anyways...
    rt = PriorityTable{T,Int,Int,typeof(ratetogroup)}(minrate, 2*minrate,
                                                        Vector{PriorityGroup{T,Vector{Int}}}(),
                                                        Vector{T}(), zero(T),
                                                        Vector{Tuple{Int,Int}}(), ratetogroup)

    DirectCRJumpAggregation{T,S,F1,F2,RNG,typeof(dg),typeof(rt),typeof(ratetogroup)}(
                                            nj, njt, et, crs, sr, maj, rs, affs!, sps, rng,
                                            dg, minrate, maxrate, rt, ratetogroup)
end

########### The following routines should be templates for all SSAs ###########

# condition for jump to occur
@inline function (p::DirectCRJumpAggregation)(u, t, integrator)
    p.next_jump_time == t
end

# executing jump at the next jump time
function (p::DirectCRJumpAggregation)(integrator)
    execute_jumps!(p, integrator, integrator.u, integrator.p, integrator.t)
    generate_jumps!(p, integrator, integrator.u, integrator.p, integrator.t)
    register_next_jump_time!(integrator, p, integrator.t)
    nothing
end

# setting up a new simulation
function (p::DirectCRJumpAggregation)(dj, u, t, integrator) # initialize
    initialize!(p, integrator, u, integrator.p, t)
    register_next_jump_time!(integrator, p, t)
    nothing
end

############################# Required Functions ##############################

# creating the JumpAggregation structure (function wrapper-based constant jumps)
function aggregate(aggregator::DirectCR, u, p, t, end_time, constant_jumps,
                   ma_jumps, save_positions, rng; kwargs...)

    # handle constant jumps using function wrappers
    rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps)

    build_jump_aggregation(DirectCRJumpAggregation, u, p, t, end_time, ma_jumps,
                           rates, affects!, save_positions, rng; num_specs=length(u), kwargs...)
end

# set up a new simulation and calculate the first jump / jump time
function initialize!(p::DirectCRJumpAggregation, integrator, u, params, t)

    # initialize rates
    fill_rates_and_sum!(p, u, params, t)

    # if no maxrate was set, use largest initial rate (pad by 2 for an extra group)
    isinf(p.maxrate) && (p.maxrate = 2*maximum(p.cur_rates))

    # setup PriorityTable
    p.rt   = PriorityTable(p.ratetogroup, p.cur_rates, p.minrate, p.maxrate)

    generate_jumps!(p, integrator, u, params, t)
    nothing
end

# execute one jump, changing the system state
function execute_jumps!(p::DirectCRJumpAggregation, integrator, u, params, t)
    # execute jump
    u = update_state!(p, integrator, u)

    # update current jump rates
    update_dependent_rates!(p, u, params, t)
    nothing
end

# calculate the next jump / jump time
function generate_jumps!(p::DirectCRJumpAggregation, integrator, u, params, t)
    @fastmath p.next_jump_time = t + calc_next_jump!(p, u, params, t)
    nothing
end

######################## SSA specific helper routines #########################

# searches down the rate list for the next reaction
@fastmath function calc_next_jump!(p::DirectCRJumpAggregation, u, params, t)

    # next jump type
    p.next_jump = sample(p.rt, p.cur_rates, p.rng)

    # return time to next jump
    randexp(p.rng) / p.sum_rate
end

# recalculate jump rates for jumps that depend on the just executed jump
# requires dependency graph
function update_dependent_rates!(p::DirectCRJumpAggregation, u, params, t)
    @unpack cur_rates, rates, ma_jumps, rt = p
    @inbounds dep_rxs = p.dep_gr[p.next_jump]
    num_majumps = get_num_majumps(ma_jumps)

    @inbounds for rx in dep_rxs
        oldrate = cur_rates[rx]

        # update rate
        if rx <= num_majumps
            newrate = evalrxrate(u, rx, ma_jumps)
        else
            newrate = rates[rx-num_majumps](u, params, t)
        end
        cur_rates[rx] = newrate

        # update table
        update!(rt, rx, oldrate, newrate)
    end

    p.sum_rate = groupsum(rt)
    nothing
  end


text:
using RDatasets, DataArrays, Gadfly

plot(dataset("car", "SLID"), y=:Wages, color=:Language,
     Geom.histogram(position=:dodge, orientation=:horizontal))


text:
function ERI_2e2c(BS::BasisSet, T::DataType = Float64)

    # Pre allocate output
    out = zeros(T, BS.nbas, BS.nbas)

    # Pre compute a list of angular momentum numbers (l) for each shell
    lvals = [Libcint.CINTcgtos_spheric(i-1, BS.lc_bas) for i = 1:BS.nshells]
    Lmax = maximum(lvals)

    # Offset list for each shell, used to map shell index to AO index
    ao_offset = [sum(lvals[1:(i-1)]) for i = 1:BS.nshells]

    buf_arrays = [zeros(Cdouble, Lmax^4) for _ = 1:Threads.nthreads()]

    @sync for i in 1:BS.nshells
        Threads.@spawn begin
            @inbounds begin
                Li = lvals[i]
                buf = buf_arrays[Threads.threadid()]
                ioff = ao_offset[i]
                for j in i:BS.nshells
                    Lj = lvals[j]
                    joff = ao_offset[j]

                    # Call libcint
                    cint2c2e_sph!(buf, Cint.([i-1,j-1]), BS.lc_atoms, BS.natoms, BS.lc_bas, BS.nbas, BS.lc_env)

                    # Loop through shell block and save unique elements
                    for js = 1:Lj
                        J = joff + js
                        for is = 1:Li
                            I = ioff + is
                            J < I ? break : nothing
                            out[I,J] = buf[is + Li*(js-1)]
                            out[J,I] = out[I,J]
                        end
                    end
                end
            end #inbounds
        end #spawn
    end #sync
    return out
end


text:
using Test
using PolynomialBasis
# using Revise
using CutCellDG
include("../useful_routines.jl")

L, W = 1.0, 1.0
levelsetbasis = LagrangeTensorProductBasis(2, 1)
dim, nf = size(interpolation_points(levelsetbasis))

mesh = CutCellDG.CGMesh([0.0, 0.0], [L, W], [1, 1], nf)
levelset = CutCellDG.LevelSet(
    x -> plane_distance_function(x, [1.0, 0.0], [0.75, 0.0]),
    mesh,
    levelsetbasis,
)
CutCellDG.load_coefficients!(levelset, 1)
cellmap = CutCellDG.cell_map(mesh, 1)

xq = [0.0, 0.0]
x0 = [0.75, 0.5]
X0 = CutCellDG.inverse(cellmap, x0)
invjac = CutCellDG.inverse_jacobian(mesh)
func = CutCellDG.interpolater(levelset)

grad(X) = CutCellDG.spatial_gradient(func, X, invjac)
hess(X) = CutCellDG.spatial_hessian(func, X, invjac)
invfunc(x) = func(CutCellDG.inverse(cellmap, x))

g0 = grad(X0)
l0 = (xq - x0)' * g0 / (g0' * g0)

p0 = func(X0)
∇p0 = grad(X0)
∇2p0 = hess(X0)
gf = CutCellDG.saye_newton_gradient(x0, l0, xq, p0, ∇p0)
hf = CutCellDG.saye_newton_hessian(l0, ∇p0, ∇2p0)

x1, l1 = CutCellDG.step_saye_newton_iterate(x0, l0, gf, hf, 2.0)
testcp = [0.75, 0.0]
@test allapprox(x1, testcp)

X1 = CutCellDG.inverse(cellmap, x1)
p1 = func(X1)
∇p1 = grad(X1)
∇2p1 = hess(X1)
gf = CutCellDG.saye_newton_gradient(x1, l1, xq, p1, ∇p1)
hf = CutCellDG.saye_newton_hessian(l1, ∇p1, ∇2p1)

x2, l2 = CutCellDG.step_saye_newton_iterate(x1, l1, gf, hf, 2.0)
@test allapprox(x2, testcp)
################################################################################

x0 = [0.75, 0.5]
X0 = CutCellDG.inverse(cellmap, x0)
p0 = func(X0)
∇p0 = grad(X0)
∇2p0 = hess(X0)

x1, l1 = CutCellDG.step_chopp_iterate(x0, xq, p0, ∇p0, 0.5)
@test allapprox(x1, testcp)

X1 = CutCellDG.inverse(cellmap, x1)
p1 = func(X1)
∇p1 = grad(X1)
∇2p0 = hess(X1)

x2, l2 = CutCellDG.step_chopp_iterate(x1, xq, p1, ∇p1, 0.5)
@test allapprox(x2, testcp)

################################################################################
tol = 1e3eps()
x1, flag = CutCellDG.spatial_closest_point(
    xq,
    x0,
    func,
    cellmap,
    tol,
    5.0,
    20,
    1e4eps(),
)
@test allapprox(x1, testcp)
@test flag

################################################################################
mesh = CutCellDG.CGMesh([0.0, 0.0], [L, W], [2, 3], nf)
normal = [1.0, 0.0]
xI = [0.75, 0.0]
levelset = CutCellDG.LevelSet(
    x -> plane_distance_function(x, normal, xI),
    mesh,
    levelsetbasis,
)
cutmesh = CutCellDG.CutMesh(mesh, levelset)

refseedpoints, seedcellids = CutCellDG.seed_zero_levelset(2, levelset, cutmesh)
seedpoints = CutCellDG.map_to_spatial(refseedpoints, seedcellids, mesh)

querypoints = [
    0.0 0.5 0.8 0.9 1.2 -0.8
    0.1 0.8 0.3 0.5 0.4 0.2
]

dx = CutCellDG.element_size(mesh)
tol = 1e3eps()
boundingradius = 1.5maximum(dx)
closestpoints, closestcellids, flags =
    CutCellDG.closest_points_on_zero_levelset(
        querypoints,
        seedpoints,
        seedcellids,
        levelset,
        tol,
        boundingradius,
    )

testcp = closest_point_on_plane(querypoints, normal, xI)
testcellids = [4, 6, 4, 5, 5, 4]
@test all(flags)
@test allapprox(testcp, closestpoints)
@test allequal(closestcellids, testcellids)
################################################################################

################################################################################
L, W = 1.0, 1.0
levelsetbasis = HermiteTensorProductBasis(2)

dim, nf = size(interpolation_points(levelsetbasis))

mesh = CutCellDG.CGMesh([0.0, 0.0], [L, W], [1, 1], nf)
levelset = CutCellDG.LevelSet(
    x -> plane_distance_function(x, [1.0, 0.0], [0.75, 0.0]),
    mesh,
    levelsetbasis,
)
CutCellDG.load_coefficients!(levelset, 1)
cellmap = CutCellDG.cell_map(mesh, 1)

xq = [0.0, 0.0]
x0 = [0.75, 0.5]
X0 = CutCellDG.inverse(cellmap, x0)
invjac = CutCellDG.inverse_jacobian(mesh)
func = CutCellDG.interpolater(levelset)

grad(X) = CutCellDG.spatial_gradient(func, X, invjac)
hess(X) = CutCellDG.spatial_hessian(func, X, invjac)
invfunc(x) = func(CutCellDG.inverse(cellmap, x))

g0 = grad(X0)
l0 = (xq - x0)' * g0 / (g0' * g0)

p0 = func(X0)
∇p0 = grad(X0)
∇2p0 = hess(X0)
gf = CutCellDG.saye_newton_gradient(x0, l0, xq, p0, ∇p0)
hf = CutCellDG.saye_newton_hessian(l0, ∇p0, ∇2p0)

x1, l1 = CutCellDG.step_saye_newton_iterate(x0, l0, gf, hf, 2.0)
testcp = [0.75, 0.0]
@test allapprox(x1, testcp, 1e3eps())

X1 = CutCellDG.inverse(cellmap, x1)
p1 = func(X1)
∇p1 = grad(X1)
∇2p1 = hess(X1)
gf = CutCellDG.saye_newton_gradient(x1, l1, xq, p1, ∇p1)
hf = CutCellDG.saye_newton_hessian(l1, ∇p1, ∇2p1)

x2, l2 = CutCellDG.step_saye_newton_iterate(x1, l1, gf, hf, 2.0)
@test allapprox(x2, testcp, 1e4eps())
################################################################################

################################################################################
x0 = [0.75, 0.5]
X0 = CutCellDG.inverse(cellmap, x0)
p0 = func(X0)
∇p0 = grad(X0)
∇2p0 = hess(X0)

x1, l1 = CutCellDG.step_chopp_iterate(x0, xq, p0, ∇p0, 0.5)
@test allapprox(x1, testcp, 1e4eps())

X1 = CutCellDG.inverse(cellmap, x1)
p1 = func(X1)
∇p1 = grad(X1)
∇2p0 = hess(X1)

x2, l2 = CutCellDG.step_chopp_iterate(x1, xq, p1, ∇p1, 0.5)
@test allapprox(x2, testcp, 1e4eps())
################################################################################

################################################################################
tol = 1e3eps()
x1, flag = CutCellDG.spatial_closest_point(
    xq,
    x0,
    func,
    cellmap,
    tol,
    5.0,
    20,
    1e4eps(),
)
@test allapprox(x1, testcp, 1e4eps())
@test flag
################################################################################

################################################################################
mesh = CutCellDG.CGMesh([0.0, 0.0], [L, W], [2, 3], nf)
normal = [1.0, 0.0]
xI = [0.75, 0.0]
levelset = CutCellDG.LevelSet(
    x -> plane_distance_function(x, normal, xI),
    mesh,
    levelsetbasis,
)
cutmesh = CutCellDG.CutMesh(mesh, levelset)

refseedpoints, seedcellids = CutCellDG.seed_zero_levelset(2, levelset, cutmesh)
seedpoints = CutCellDG.map_to_spatial(refseedpoints, seedcellids, mesh)

querypoints = [
    0.0 0.5 0.8 0.9 1.2 -0.8
    0.1 0.8 0.3 0.5 0.4 0.2
]

dx = CutCellDG.element_size(mesh)
tol = 1e3eps()
boundingradius = 1.5maximum(dx)
closestpoints, closestcellids, flags =
    CutCellDG.closest_points_on_zero_levelset(
        querypoints,
        seedpoints,
        seedcellids,
        levelset,
        tol,
        boundingradius,
    )

testcp = closest_point_on_plane(querypoints, normal, xI)
testcellids = [4, 6, 4, 5, 5, 4]
@test all(flags)
@test allapprox(testcp, closestpoints)
@test allequal(closestcellids, testcellids)


text:
@testset "Base models: known capacities" begin
    opt = ECOS.Optimizer

    g = MetaDiGraph()
    add_vertex!(g, :name, "A")
    add_vertex!(g, :name, "B")
    add_vertex!(g, :name, "C")
    add_vertex!(g, :name, "D")
    add_edge!(g, 1, 2, :capacity, 2.0)
    add_edge!(g, 1, 3, :capacity, 2.0)
    add_edge!(g, 2, 4, :capacity, 2.0)
    add_edge!(g, 3, 4, :capacity, 2.0)

    k = SimpleDiGraph(4)
    add_edge!(k, 1, 4)
    d = Edge(1, 4)
    dm = Dict{Edge{Int}, Float64}(d => 4.0)

    @testset "Base case" begin
        m = Model(opt)
        mt = ModelType(Load(), MinimumTotal(), FlowFormulation(), false, CuttingPlane(), ObliviousUncertainty(), UncertainDemand())
        rd = RoutingData(g, k, opt, mt)
        rm = basic_routing_model_unitary(m, rd, FlowFormulation())

        @test_throws ErrorException basic_routing_model_unitary(m, rd, 42)
        @test_throws ErrorException basic_routing_model_unscaled(m, rd, 42)

        @test_throws ErrorException total_flow_in_edge(rm, Edge(1, 2), 42)
        @test_throws ErrorException total_flow_in_edge(rm, Edge(1, 2), dm, 42)
    end

    @testset "Flow" begin
        # Only useful parameter: FlowFormulation.
        mt = ModelType(Load(), MinimumTotal(), FlowFormulation(), false, CuttingPlane(), ObliviousUncertainty(), UncertainDemand())
        rd = RoutingData(g, k, opt, mt, traffic_matrix=dm)

        function __check_feasible_knowncapa_flow(rm::RoutingModel)
            @test rm.routing !== nothing
            @test rm.constraints_source_in !== nothing && length(rm.constraints_source_in) >= 0
            @test rm.constraints_source_out !== nothing && length(rm.constraints_source_out) > 0
            @test rm.constraints_target_in !== nothing && length(rm.constraints_target_in) > 0
            @test rm.constraints_target_out !== nothing && length(rm.constraints_target_out) >= 0
            @test rm.constraints_balance !== nothing && length(rm.constraints_balance) > 0
            @test rm.constraints_convexity === nothing || length(rm.constraints_convexity) == 0

            # Check that the model has a feasible solution.
            set_silent(rm.model)
            optimize!(rm.model)
            @test termination_status(rm.model) == MOI.OPTIMAL

            @test value(rm.routing[d, Edge(1, 2)]) == value(rm.routing[d, Edge(2, 4)])
            @test value(rm.routing[d, Edge(1, 3)]) == value(rm.routing[d, Edge(3, 4)])

            # Check that the model correctly enforces the connectivity.
            for e in [Edge(1, 2), Edge(2, 4), Edge(1, 3), Edge(3, 4)]
                @objective(rm.model, Max, rm.routing[d, e])
                optimize!(rm.model)
                @test termination_status(rm.model) == MOI.OPTIMAL

                @test value(rm.routing[d, Edge(1, 2)]) ≈ value(rm.routing[d, Edge(2, 4)]) atol=1.0e-5
                @test value(rm.routing[d, Edge(1, 3)]) ≈ value(rm.routing[d, Edge(3, 4)]) atol=1.0e-5
            end
        end

        function __check_feasible_knowncapa_withcapa_flow(rm::RoutingModel)
            # Capacity helpers.
            @test rm.constraints_capacity === nothing || length(rm.constraints_capacity) == 0
            capacity_constraints(rm, dm)
            @test rm.constraints_capacity !== nothing && length(rm.constraints_capacity) > 0

            @objective(rm.model, Max, rm.routing[d, Edge(1, 2)])
            optimize!(rm.model)
            @test termination_status(rm.model) == MOI.OPTIMAL

            @test value(rm.routing[d, Edge(1, 2)]) ≈ value(rm.routing[d, Edge(2, 4)]) atol=1.0e-5
            @test value(rm.routing[d, Edge(1, 3)]) ≈ value(rm.routing[d, Edge(3, 4)]) atol=1.0e-5

            @test value(total_flow_in_edge(rm, Edge(1, 2), dm)) ≈ 2
            @test value(total_flow_in_edge(rm, Edge(1, 3), dm)) ≈ 2

            if rm.scaled_flows == UnitaryFlows
                @test_throws AssertionError value(total_flow_in_edge(rm, Edge(1, 2))) ≈ 2
            else
                @test value(total_flow_in_edge(rm, Edge(1, 2))) ≈ 2
            end

            # Incompatible with other capacity constraints.
            @test_throws AssertionError mu_capacity_constraints(rm, dm)
        end

        @testset "Unitary flows" begin
            m = Model(opt)
            rm = basic_routing_model_unitary(m, rd, FlowFormulation())

            @test rm.data == rd
            @test rm.model == m
            @test rm.scaled_flows == UnitaryFlows

            __check_feasible_knowncapa_flow(rm)
            __check_feasible_knowncapa_withcapa_flow(rm)

            # Helpers.
            @test total_flow_in_edge(rm, Edge(1, 2), dm) == 4 * rm.routing[d, Edge(1, 2)]
        end

        @testset "Unscaled flows" begin
            m = Model(opt)
            rm = basic_routing_model_unscaled(m, rd, dm, FlowFormulation())

            @test rm.data == rd
            @test rm.model == m
            @test rm.scaled_flows == UnscaledFlows

            __check_feasible_knowncapa_flow(rm)
            __check_feasible_knowncapa_withcapa_flow(rm)

            # Helpers.
            @test total_flow_in_edge(rm, Edge(1, 2), dm) == rm.routing[d, Edge(1, 2)]
        end

        @testset "µ-capacited unitary flows" begin
            m = Model(opt)
            rm = basic_routing_model_unitary(m, rd, FlowFormulation())
            mu_capacity_constraints(rm, dm)
            @objective(m, Min, rm.mu)

            __check_feasible_knowncapa_flow(rm)

            # Helpers.
            @test total_flow_in_edge(rm, Edge(1, 2), dm) == 4 * rm.routing[d, Edge(1, 2)]

            # Incompatible with other capacity constraints.
            @test_throws AssertionError capacity_constraints(rm, dm)
            @test_throws AssertionError mu_capacity_constraints(rm)
        end
    end

    @testset "Path" begin
        # Only useful parameter: PathFormulation.
        mt = ModelType(Load(), MinimumTotal(), PathFormulation(), false, CuttingPlane(), ObliviousUncertainty(), UncertainDemand())
        rd = RoutingData(g, k, opt, mt, traffic_matrix=dm)

        function __check_feasible_knowncapa_path(rm::RoutingModel)
            @test rm.routing !== nothing
            @test rm.constraints_source_in === nothing || length(rm.constraints_source_in) == 0
            @test rm.constraints_source_out === nothing || length(rm.constraints_source_out) == 0
            @test rm.constraints_target_in === nothing || length(rm.constraints_target_in) == 0
            @test rm.constraints_target_out === nothing || length(rm.constraints_target_out) == 0
            @test rm.constraints_balance === nothing || length(rm.constraints_balance) == 0
            @test rm.constraints_convexity !== nothing && length(rm.constraints_convexity) > 0

            # Check that the model has a feasible solution.
            set_silent(rm.model)
            optimize!(rm.model)
            @test termination_status(rm.model) == MOI.OPTIMAL

            if rm.scaled_flows == UnitaryFlows
                @test 0 <= value(rm.routing[d, 1]) <= 1
                @test 0 <= value(rm.routing[d, 2]) <= 1
            else
                @test rm.scaled_flows == UnscaledFlows
                @test value(rm.routing[d, 1]) >= 0
                @test value(rm.routing[d, 2]) >= 0
            end

            # Check that the model correctly enforces the connectivity.
            for p in [1, 2]
                @objective(rm.model, Max, rm.routing[d, p])
                optimize!(rm.model)
                @test termination_status(rm.model) == MOI.OPTIMAL

                @test value(rm.routing[d, 1]) + value(rm.routing[d, 2]) ≈ ((rm.scaled_flows == UnitaryFlows) ? 1 : 4)
            end

            # Capacity helpers.
            @test rm.constraints_capacity === nothing || length(rm.constraints_capacity) == 0
            capacity_constraints(rm, dm)
            @test rm.constraints_capacity !== nothing && length(rm.constraints_capacity) > 0

            @objective(rm.model, Max, rm.routing[d, 1])
            optimize!(rm.model)
            @test termination_status(rm.model) == MOI.OPTIMAL

            @test value(total_flow_in_edge(rm, Edge(1, 2), dm)) ≈ 2
            @test value(total_flow_in_edge(rm, Edge(1, 3), dm)) ≈ 2

            if rm.scaled_flows == UnitaryFlows
                @test_throws AssertionError value(total_flow_in_edge(rm, Edge(1, 2))) ≈ 2
            else
                @test value(total_flow_in_edge(rm, Edge(1, 2))) ≈ 2
            end

            # Helpers.
            tf12 = total_flow_in_edge(rm, Edge(1, 2), dm)
            if rm.scaled_flows == UnitaryFlows
                @test tf12 == 4 * rm.routing[d, 1]
            else
                @test rm.scaled_flows == UnscaledFlows
                @test tf12 == rm.routing[d, 1]
            end
        end

        @testset "Unitary flows" begin
            m = Model(opt)
            rm = basic_routing_model_unitary(m, rd, PathFormulation())

            @test rm.data == rd
            @test rm.model == m
            @test rm.scaled_flows == UnitaryFlows

            __check_feasible_knowncapa_path(rm)
        end

        @testset "Unscaled flows" begin
            m = Model(opt)
            rm = basic_routing_model_unscaled(m, rd, dm, PathFormulation())

            @test rm.data == rd
            @test rm.model == m
            @test rm.scaled_flows == UnscaledFlows

            __check_feasible_knowncapa_path(rm)
        end
    end
end


text:
next_prime = # https://projecteuler.net/problem=3
# The prime factors of 13195 are 5, 7, 13 and 29.
# What is the largest prime factor of the number 600851475143 ?

function find_next_prime(start, stop, primes)
    candidates = []
    for p in primes
        candidates = ones(Int, stop - start)
        println(candidates)
        i = (start + p) % p
        while i < size(candidates, 1)
            println(i)
            candidates[i+1] = 0
            i += p
        end
    end
    println("-------")
    println(candidates)

    for i in 1:(stop - start + 1)
        if candidates[i] == 1
            return start + i - 1
        end
    end
end

primes = [2, 3, 5, 7]
next_prime = find_next_prime(10, 15, primes)
println(next_prime)

max_stars_repo_name: the-gigi/project-euler

max_issues_repo_name: the-gigi/project-euler

max_forks_repo_name: the-gigi/project-euler


text:
using CompScienceMeshes
using BEAST

using Test

o, x, y, z = euclidianbasis(3)
tet = simplex(x,y,z,o)
ctr = center(tet)

T = Float64
iref = BEAST.LagrangeRefSpace{T,1,4,4}()
oref = BEAST.NDLCCRefSpace{T}()
ishp = BEAST.Shape{T}(1, 3, 1.0)

output = BEAST.gradient(iref, ishp, tet)

oval = point(0,0,0)
for oshp in output
    global oval += oref(ctr)[oshp.refid].value * oshp.coeff
end

@test oval ≈ point(0,0,1)


text:
#=
testODE:
- Julia version:
- Author: ymocquar
- Date: 2019-11-13
=#
include("../src/twoscales_pure_ab.jl")
using DifferentialEquations
using LinearAlgebra
using Plots
using Random

function twoscales_solve(par_u0::PrepareU0, order, t, nb)

    parGen = PrepareTwoScalePureAB(nb, t, order, par_u0)

    return twoscales_pure_ab(parGen, only_end = true)

end

function getmindif(tab::Vector{Vector{BigFloat}})
    nmmin = Inf
    ret = (1, 2)
    for i = 1:size(tab, 1)
        for j = (i+1):size(tab, 1)
            nm = norm(tab[i] - tab[j], Inf)
            if nm < nmmin
                nmmin = nm
                ret = i, j
            end
        end
    end
    return ret, min(nmmin, 1.1)
end

function fctMain(n_tau)

    #   u0 =[big"0.12345678", big"0.1209182736", big"0.1290582671", big"0.1239681094" ]
    seed = 12996
    Random.seed!(seed)
    u0 = rand(BigFloat, 4)
    println("seed = $seed")
    #   tab_eps = zeros(BigFloat,7)
    #   tab_eps= [big"0.5", big"0.2", big"0.1",big"0.05", big"0.02", big"0.01",big"0.005", big"0.002", big"0.001",big"0.0005", big"0.0002", big"0.0001"]
    tab_eps = [
        big"0.1",
        big"0.01",
        big"0.001",
        big"0.0001",
        big"0.00001",
        big"0.000001",
        big"0.0000001",
        big"0.00000001",
        big"0.000000001",
    ]
    # epsilon=big"0.1"
    # for i=1:7
    #     tab_eps[i] = epsilon
    #     epsilon /= 100
    # end
    nbmaxtest = 11
    t_max = big"1.0"
    A = [0 0 1 0; 0 0 0 0; -1 0 0 0; 0 0 0 0]
    tabtabsol = Vector{Vector{Vector{BigFloat}}}()
    for order = 2:14
        ordprep = order + 2
        y = ones(Float64, nbmaxtest, size(tab_eps, 1))
        y_big = ones(BigFloat, nbmaxtest, size(tab_eps, 1))
        x = zeros(Float64, nbmaxtest)
        ind = 1
        for epsilon in tab_eps
            #        fct = u -> [ 1.11u[2]^2+0.8247u[4]^3+0.12647u[3]^4,
            #       0.4356789141u[1]-0.87u[3]*u[4],
            #       1.9898985/(1+1.1237u[4]^3), 0.97u[1]*u[3]+0.8111u[2]^2 ]
            fct =
                (u, p, t) -> [
                    u[3]^3 + sin(t^2),
                    u[4] + u[1]^2 + 1 / (1 + exp(t)) - 0.35,
                    u[1] * u[2] * u[4] * (1 + t),
                    u[2] - u[2]^2 + 1 / (1 + u[4]^2) + u[1]^2 + cos(exp(t)),
                ]
            parphi = PreparePhi(epsilon, n_tau, A, fct)
            println("prepareU0 eps=$epsilon n_tau=$n_tau")
            nb = 10
            @time par_u0 = PrepareU0(parphi, ordprep, u0)
            @time pargen = PrepareTwoScalesPureAB(nb * 2^nbmaxtest, t_max, order, par_u0)
            @time solref = twoscales_pure_ab(pargen, only_end = true)
            if size(tabtabsol, 1) < ind
                push!(tabtabsol, Vector{Vector{BigFloat}}())
            end
            tabsol = tabtabsol[ind]
            push!(tabsol, solref)
            res_gen = Array{Array{BigFloat,1},1}(undef, nbmaxtest)
            indref = 1
            eps_v = convert(Float32, epsilon)
            println("epsilon = $eps_v solref=$solref")
            indc = 1
            labels = Array{String,2}(undef, 1, ind)
            while indc <= nbmaxtest
                @time pargen = PrepareTwoScalesPureAB(nb, t_max, order, par_u0)
                @time sol = twoscales_pure_ab(pargen, only_end = true)
                push!(tabsol, sol)
                res_gen[indc] = sol
                diff = solref - sol
                println("tabsol=$tabsol")
                (a, b), nm = getmindif(tabsol)
                if a != indref
                    println("New solref !!!! a=$a, b=$b nm=$nm")
                    indref = a
                    solref = tabsol[a]
                    for i = 1:indc
                        nm2 = min(norm(res_gen[i] - solref, Inf), 1.1)
                        y[i, ind] = nm2 == 0 ? nm : nm2
                        y_big[i, ind] = nm2 == 0 ? nm : nm2
                    end
                else
                    diff = solref - sol
                    y[indc, ind] = min(norm(diff, Inf), 1.1)
                    y_big[indc, ind] = min(norm(diff, Inf), 1.1)
                end
                println("solref=$solref")
                println("nb=$nb sol=$sol")
                x[indc] = 1.0 / nb
                println("epsilon=$epsilon result=$y")
                println("epsilon=$epsilon reslog2=$(log2.(y))")
                println("epsilon=$epsilon result=$y_big")
                println("epsilon=$epsilon reslog2=$(log2.(y_big))")
                nb *= 2
                indc += 1
            end
            if ind == size(tab_eps, 1)
                for i = 1:ind
                    labels[1, i] = " epsilon,order=$(convert(Float32,tab_eps[i])),$order "
                end
                gr()
                p = Plots.plot(
                    x,
                    view(y, :, 1:ind),
                    xlabel = "delta t",
                    xaxis = :log,
                    ylabel = "error",
                    yaxis = :log,
                    legend = :bottomright,
                    label = labels,
                    marker = 2,
                )

                prec_v = precision(BigFloat)
                Plots.savefig(
                    p,
                    "out/r20_$(prec_v)_$(eps_v)_$(order)_$(ordprep)_$(n_tau)_epsilon_fct.pdf",
                )
            end
            ind += 1
        end
    end
end

# testODESolverEps()

# for i=3:9
#     fctMain(2^i)
# end
# setprecision(512)
fctMain(32)