<jupyter_start><jupyter_text>Spectrum Slicing Test Suite DemoWe will demo the spectrum slicing test suite by generating test problems from a given matrix and visualizing the expected performance of two different rational filters in FEAST.First we load the `SpectrumSlicingTestSuite` and other required modules:<jupyter_code>import SpectrumSlicingTestSuite
import MatrixDepot
import BenchmarkProfiles
using Measures
using LinearAlgebra
using Plots<jupyter_output>┌ Info: Precompiling SpectrumSlicingTestSuite [9d626e9e-037c-4cf1-92ae-c6e6dcb749aa]
└ @ Base loading.jl:1342<jupyter_text>The rational filtersWe use rational filter functions taken from [here](https://arxiv.org/abs/1704.03255)<jupyter_code>zGauss = [
 -0.998055+0.0623361im
 -0.949425+0.313993im 
 -0.73489+0.678186im 
 -0.284168+0.958775im 
];
aGauss = [
 0.0252579-0.00157755im
 0.0527835-0.0174565im 
 0.057635-0.0531879im 
 0.0257658-0.086933im  
];

zSLiSe = [
 -0.999716+0.0100714im
 -0.985299+0.0834839im
 -0.890636+0.30265im  
 -0.435696+0.697878im 
];
aSLiSe = [
  0.00522242-0.000328736im
  0.0197875-0.00532048im 
  0.0532178-0.0321933im  
  0.0536992-0.121132im   
];<jupyter_output><empty_output><jupyter_text>Now, we plot the two functions on the interval `[-1,1]`<jupyter_code>x = range(-1.5, 1.5, length=1001)
f = function(z,a) SpectrumSlicingTestSuite.phi(z,a,[-1,1],x) |> real end
plot(x, f(zGauss,aGauss), label="Gauss")
plot!(x,f(zSLiSe,aSLiSe), label="SLiSe")<jupyter_output><empty_output><jupyter_text>Obtaining the Test IntervalsFirst, we get an example matrix from the UoF collection:<jupyter_code>try
# MatrixDepot.update()
MatrixDepot.matrixdepot("PARSEC/Si2", :get)
catch
end

const B = MatrixDepot.matrixdepot("PARSEC/Si2")<jupyter_output><empty_output><jupyter_text>Additionally, we require some parameters:- `n` the size of the matrix- `nexPercentage` the search size increment for the subspace iteration. In FEAST parlance the size of the subspace, `M0`, is `nexPercentage` times the number of eigenvalues inside each interval - `nvec` the number of Lanczos procedures for the Density of States. - `V`, the starting vectors for the Lanczos procedures.<jupyter_code>const n = size(B)[1];
const nexPercentage = 0.5;
const nvec = 10;
V = (randn(n,nvec) + randn(n,nvec) * 1im) / sqrt(2);<jupyter_output><empty_output><jupyter_text>Now to the generation of the intervals. `W` is the spectrum of the matrix `B`, mapped to `[-1,1]`.<jupyter_code>intvls,W = SpectrumSlicingTestSuite.testIntervals(B,V; l=0.02, u=0.04);<jupyter_output><empty_output><jupyter_text>We get some basic information of generated intervals<jupyter_code>@info "number of generated intervals: " size(intvls)[1]<jupyter_output>┌ Info: number of generated intervals: 
│   (size(intvls))[1] = 398
└ @ Main In[7]:1<jupyter_text>We plot the endoints of all the intervals, the eigenvalues `W`, and one of the Density of States that are used to generate the intervals<jupyter_code>p = SpectrumSlicingTestSuite.kpm(Diagonal(W),V,25);

x = range(-1, 1, length=1001);
plot( x, p.(x) , label="Density of States", title="Si2")
scatter!(W, 0.0 .*ones(size(W)...), label="Eigenvalues")
intervalEdge = union( map( x-> x[1], intvls), map( x-> x[2], intvls))
scatter!(intervalEdge, 0.5.*ones(size(intervalEdge)...), label="Interval Edges")<jupyter_output><empty_output><jupyter_text>Performance of Rational Filters in FEASTGet the FEAST convergence ratios for these intervals<jupyter_code>convergenceRatioGauss = SpectrumSlicingTestSuite.convergenceratio(zGauss,aGauss,W,intvls,nexPercentage);
convergenceRatioSLiSe = SpectrumSlicingTestSuite.convergenceratio(zSLiSe,aSLiSe,W,intvls,nexPercentage);<jupyter_output><empty_output><jupyter_text>Plot the convergence rates as a performance profile, see [here](https://link.springer.com/article/10.1007%2Fs101070100263) and [here](https://dl.acm.org/citation.cfm?id=2950048)<jupyter_code>BenchmarkProfiles.performance_profile(
    BenchmarkProfiles.PlotsBackend(),
    [convergenceRatioGauss convergenceRatioSLiSe],
    ["Gauss","SLiSe"],
    logscale=false, title="Si2")<jupyter_output><empty_output><jupyter_text>Filter performance for additional matricesIn order to obtain a more well-rounded comparison of the two filters, we apply the above process to multiple matrices from the [SuiteSparse Matrix Collection](https://sparse.tamu.edu/)<jupyter_code>function getRatios( s, nexPercentage = 0.5, nvec = 10 )
    try
        # MatrixDepot.update()
        MatrixDepot.matrixdepot(s)
    catch
    end

    B = MatrixDepot.matrixdepot(s);
    n = size(B)[1];
    V = (randn(n,nvec) + randn(n,nvec) * 1im) / sqrt(2);
    intvls,W = SpectrumSlicingTestSuite.testIntervals(B,V;l=0.02, u=0.04);
    convergenceRatioGauss = SpectrumSlicingTestSuite.convergenceratio(zGauss,aGauss,W,intvls,nexPercentage);
    convergenceRatioSLiSe = SpectrumSlicingTestSuite.convergenceratio(zSLiSe,aSLiSe,W,intvls,nexPercentage);
    (convergenceRatioGauss,convergenceRatioSLiSe)
end

problems1 = [
    "PARSEC/Na5"
    "PARSEC/SiH4"
    "PARSEC/SiNa"
    "PARSEC/benzene"
]

ratios1 = map(getRatios,problems1);<jupyter_output><empty_output><jupyter_text>Shown below is a plot generated from four matrices from the PARSEC group. These matrices are obtained from electronic structure calculations. Note that the plots looks very similar to the performance profile above.<jupyter_code>plots = map(ratios1) do x
    BenchmarkProfiles.performance_profile(
        BenchmarkProfiles.PlotsBackend(),
        [x[1] x[2]],
        ["",""],
        logscale=false)
end

plot(plots..., 
    legend=false, 
    xlabel="", 
    ylabel="", 
    margin=5mm,
    xlims = (1,6),
    title=reshape(problems1,1,length(problems1))
)<jupyter_output><empty_output><jupyter_text>Now we present performance profiles for four matrices from the VDOL group, which are optimal control problems. Here the `SLiSe` filter performs significantly better than the Gauss filter.<jupyter_code>problems2 = [
    "VDOL/spaceStation_9"
    "VDOL/freeFlyingRobot_10"
    "VDOL/lowThrust_3"
    "VDOL/lowThrust_2"
]

ratios2 = map(getRatios,problems2)
plots = map(ratios2) do x
    BenchmarkProfiles.performance_profile(
        BenchmarkProfiles.PlotsBackend(),
        [x[1] x[2]],
        ["",""],
        logscale=false)
end

Plots.plot(plots..., 
    legend=false,
    xlims = (1,50),
    xlabel="",
    ylabel="",
    margin=5mm,
    title=reshape(problems2,1,length(problems2))
)<jupyter_output><empty_output><jupyter_text>We can combine the data from all nine matrices into a single performance profile:<jupyter_code>ratioGauss = convergenceRatioGauss;
ratioSLiSe = convergenceRatioSLiSe;

for x in ratios1
   append!(ratioGauss,x[1]) 
   append!(ratioSLiSe,x[2]) 
end
for x in ratios2
   append!(ratioGauss,x[1]) 
   append!(ratioSLiSe,x[2]) 
end

p = BenchmarkProfiles.performance_profile(
        BenchmarkProfiles.PlotsBackend(),
        [ratioGauss ratioSLiSe],
        ["",""],
        logscale=false);

Plots.plot(p, xlims = (1,50), title="Combined Matrices")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1206087327549424, 'fraction_numerical': 0.057744275352012514, 'mean_word_length': 5.43956043956044, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 58, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '517345', 'license': 'MIT', 'n_tokens_mistral': 2808, 'n_tokens_neox': 2550, 'n_words': 535, 'path': 'doc/notebooks/demo.ipynb', 'repo_name': 'brunowu/SpectrumSlicingTestSuite.jl', 'size': '7031'}