Laplacian spectral analysis¶

from __future__ import annotations

import numpy as np
from scipy.sparse import linalg as splng
import pandas as pd
import plotly.express as px

from dxtr import Cochain
from dxtr.complexes import sphere
from dxtr.operators import laplacian

sizes = ['small', 'large', 'massive'] 
spheres = {size: sphere(size=size, manifold=True) for size in sizes}

N = 15

eigvals_continuous = [-n*(n+1) for n in np.arange(N)]
multiplicity_continuous = [2*n+1 for n in np.arange(N)] 

eigenvalues, multiplicity= [], []

# Expected eigenvalues and multiplicity for the Laplacian on 0-simplices
eigenvalues.append(eigvals_continuous)
multiplicity.append(multiplicity_continuous)

# Expected eigenvalues and multiplicity for the Laplacian on 1-simplices
eigenvalues.append(eigvals_continuous[1:])
multiplicity.append([2*m for m in multiplicity_continuous[1:]])

# Expected eigenvalues and multiplicity for the Laplacian on 2-simplices
eigenvalues.append(eigvals_continuous)
multiplicity.append(multiplicity_continuous)

def error(values:np.ndarray[float], 
          expectations:np.ndarray[float], average:bool=True) -> float:
    """Computes the relative RMS error.
    
    Parameters
    ----------
    values
        The list of computed values to assess.
    expectations
        The list of expected values.
    average
        If True a single, averaged value is returned,
        else, an array is returned.
    
    Returns
    -------
        The seeked error value.

    Notes
    -----
    * The use of this function to estimate the precision of our framework is
      inspired by the Ben-Chen et al (2010) paper.
    """
    
    try:
        values = np.asarray(values)
        expectations = np.asarray(expectations)

        assert values.shape == expectations.shape, ('WARNING:'
            + f'values shape {values.shape} !='
            + f'expectations shape {expectations.shape}')
        
        res = (values - expectations) ** 2
        nrm = (expectations ** 2).sum()
        if nrm < 1e-10: nrm = 1
        
        if average:
            return np.sqrt(res.sum() / nrm)
        else:
            return np.sqrt(res) / np.sqrt(nrm)

    except AssertionError as msg:
        print(msg)
        return None

results = pd.DataFrame()
for size, mfld in spheres.items():
    for k in range(3):
        
        lpc = laplacian(mfld, k).values

        j=0
        for i, (xpct_ev, multi) in enumerate(zip(eigenvalues[k],
                                                multiplicity[k])):
            rnk = multi*[i]
            eig = multi*[xpct_ev]
            mlt = multi*[multi]
            odr = np.arange(j, j+multi)

            evals, evects = splng.eigs(lpc, k=multi, sigma=xpct_ev, which='LM')
            evals = evals.real
            evects = evects.real

            err = multi*[error(evals[:multi], eig)]
            j += multi

            # Wrapping results into a DataFrame
            res = pd.DataFrame({'Operator': 'Laplacian',
                                'Manifold': 'Spheres',
                                'Complex size': size,
                                'Degree': k,
                                'Value': evals,
                                'Order':odr,
                                'Error': err,
                                'Rank': np.array(rnk)+k,
                                'Multiplicity': mlt,
                                'Expectation': eig})

            results = pd.concat((results, res))

def eigenvalue_plot(data:pd.DataFrame, show_expectation:bool=True,
                    hover_data:bool=True,
                    **kwargs) -> plotly.graph_objs._figure.Figure:
    """Plots the eigenvalues of a differential operator.

    Parameters
    ----------
    data :
        The pandas.DataFrame containing all the relevant information to plot.
    show_expectation :
        If True  shows the expected values and
        a confidence interval around them.
    hover_data:
        Display some info when hovering points if True.

    Other parameters
    ----------------
    confidence_interval : float
        The width of error band to plot around the expected values.
        Should be given as `.xx`, e.g. for 5% write `0.05`.
    operator_name : str
        The name of the operator to display and use for recording.
    manifold_name : str
        The name of the manifold to display and use for recording.
    title : str
        The name of the graph.

    Returns
    -------
    The Figure to display.

    Notes
    -----
    The DataFrame to pass as argument should contain the following columns:
    * 'Rank': int, the eigenvalue rank.
    * 'Value': float, the computed eigenvalues.
    * 'Multiplicity': int, the eigenvalue multiplicity.
    * 'Error: float, the relative RMS error between
      computed and expected eigenvalues.
    * 'Complex size': str, the size/name of the SimplicialManifold
      the eigenvalues have been computed on.
    """
    if hover_data:
        hd = {'Complex size': False,
              'Multiplicity': True,
              'Value':':.2f',
              'Rank': True,
              'Expectation': True,
              'Error': ':.2e'}
    else:
        hd = None
    
    fig = px.scatter(data, x='Rank', y='Value', 
                     color='Complex size',
                     color_discrete_sequence=px.colors.qualitative.G10,
                     labels={'Value': 'Eigenvalues',
                             'Complex size': 'Sphere sizes'},
                     hover_name='Complex size',
                     hover_data=hd)

    # Expectation
    if show_expectation:
        err = kwargs.get('Confidence interval', .01)

        rank = data['Rank'].unique()
        expectation = data['Expectation'].unique()
        theo = px.line(x=rank, y=expectation)
        uplim = px.line(x=rank, y=(1+err)*expectation)
        lolim = px.line(x=rank, y=(1-err)*expectation)

        fig.add_trace(theo.data[0])
        fig.add_trace(uplim.data[0])
        fig.add_trace(lolim.data[0])

    # Cosmetics
    fig['data'][-3]['line']['color']='rgba(0,0,0,.25)'
    fig['data'][-3]['showlegend']= True
    fig['data'][-3]['name'] = 'Theoretical expectation'
    fig['data'][-2]['line']['color'] ='rgba(0,0,0,.25)'
    fig['data'][-2]['line']['dash'] ='dot'
    fig['data'][-1]['line']['color'] ='rgba(0,0,0,.25)'
    fig['data'][-1]['line']['dash'] ='dot'
    fig['data'][-1]['showlegend'] = True
    fig['data'][-1]['name'] = f'{err:.0%} Confidence interval'

    if 'Operator' in data.columns:
        default_operator_name = data['Operator'].unique()[0]
    else:
        default_operator_name = 'Unknown operator'
    
    if 'Manifold' in data.columns:
        default_manifold_name = data['Manifold'].unique()[0]
    else:
        default_manifold_name = 'Unknown operator'

    operator = kwargs.get('operator_name', default_operator_name)
    manifold = kwargs.get('manifold_name', default_manifold_name)

    fig.update_layout(template='plotly_white',
                      title=kwargs.get('title', 
                                       f'<b>{operator}</b> | {manifold}'),
                      showlegend=kwargs.get('showlegend', True))


    return fig