Operators¶

from dxtr import Cochain
from dxtr.cochains import random_cochain
from dxtr.operators import exterior_derivative, d

c = random_cochain('icosa') # Let's consider a random cochain defined on a icosahedral simplicial complex as example.

ext_deri_c = exterior_derivative(c)

print('Congratulations !')
print(f'  * You have just computed the {ext_deri_c}')
if isinstance(ext_deri_c, Cochain) & (ext_deri_c.complex == c.complex):
    print(f'  * The exterior derivative of a {c.dim}D cochain is a {ext_deri_c.dim}D one, defined on the same complex.')

dc = d(c)
if (dc.values == ext_deri_c.values).all():
    print('  * The alias `d()` yields the same values as the function `exterior_derivative()`.')

c  = random_cochain('icosa')

incicence_matrix = c.complex[c.dim].coboundary

if (d(c).values == incicence_matrix @ c.values).all():
    print(f'The exterior derivative of our {c.name} is properly computed.')

c = random_cochain('icosa', dim=2)
dc = d(c)

if c.dim == c.complex.dim:
    assert (dc.values == 0).all()
    print(dc.dim)

c = random_cochain('icosa', dual=True)
dc = d(c)
if dc.isdual:
    print(f'The exterior derivative of a dual {c.dim}-cochain is a dual {dc.dim}-cochain.')

from dxtr import SimplicialComplex

indices = [[0,1,2,3,4], [1,2,3,4,5], [4,5,6], [6,7]]
asc = SimplicialComplex(indices)

if asc.isabstract:
    
    # Let's define a primal 1D cochain on this abstract complex...
    k = 1
    Nk = asc[k].size
    c = Cochain(asc, k, [1 for i in range(Nk)])

    # ...and compute its exterior derivative
    dc = d(c)
    if dc.complex.isabstract:
        print('We computed the exterior derivative of a cochain defined on an abstrac complex.')

from dxtr.cochains import unit_cochain
from dxtr.operators import hodge_star

for k in range(3):
    cchn = unit_cochain(dim=k, manifold=True)
    mfld = cchn.complex
    
    cchn_2 = hodge_star(cchn)
    
    kvols, kcovols = mfld[k].volumes, mfld[k].covolumes
    
    s = -1 if k==2 else 1 # See note below.
    has_proper_values = (cchn_2.values == s * kcovols/kvols).all()
    has_proper_dim =  cchn_2.dim == mfld.dim - cchn.dim
    
    if cchn_2.isdual & has_proper_values & has_proper_dim:
        print(f'The Hodge star on {k}-cochains is properly computed')
    else:
        print(cchn_2.values)
        print( kcovols/kvols)

import numpy as np

for k in range(3):
    cchn = random_cochain(dim=k, manifold=True)
    cchn_bis = hodge_star(hodge_star(cchn))
    
    n = cchn.complex.dim
    s = (-1)**(k*(n-k))
    
    try:
        np.testing.assert_array_almost_equal(cchn.values, s*cchn_bis.values)
        print(f'it works for {k}-cochains.')
    except:
        print(f'it does not work {k}-cochains.')

from dxtr.cochains import normal_vector_field

normals = normal_vector_field('icosa')
nrl, mfld =  normals.values, normals.complex

ez = np.repeat([[0,0,1]], nrl.shape[0], axis=0)

vct = Cochain(mfld, dim=0, dual=True,
              values=np.cross(ez, nrl, axis=-1),
              name='e_theta')

if vct.isvectorvalued & vct.isdual & (vct.dim==0):
    print(f'{vct.name} is a discrete vector field on a {mfld.name}.')

from dxtr.operators import flat

flat_vct = flat(vct, name='Our first flatten vector field') 

if isinstance(flat_vct, Cochain) & (not flat_vct.isvectorvalued):
    dual = 'dual' if flat_vct.isdual else 'primal'
    k = flat_vct.dim

    print(f'{flat_vct.name} is a {dual}, {k}D, scalar-valued cochain.')

dual_flat_vct = flat(vct, dual=True, name='Our second flatten vector field')

if isinstance(dual_flat_vct, Cochain) & (not dual_flat_vct.isvectorvalued):
    dual = 'dual' if dual_flat_vct.isdual else 'primal'
    k = dual_flat_vct.dim

    print(f'This time, {dual_flat_vct.name} is a {dual}, {k}D, scalar-valued cochain.')

import pyvista as pv
from dxtr.visu import visualize

fig = pv.Plotter(shape=(1,3), border=False)

fig.subplot(0,0)
visualize(vct, fig=fig, display=False, 
          layout_parameters={'title':'Initial vector field'})

fig.subplot(0,1)
visualize(flat_vct, fig=fig, scaling_factor=10, display=False, 
          layout_parameters={'title':'Flatten as Primal 1-cochain'})

fig.subplot(0,2)
visualize(dual_flat_vct, fig=fig, scaling_factor=10, display=False, 
          layout_parameters={'title':'Flatten as Dual 1-cochain'})


fig.link_views()
fig.show()

from dxtr.operators import sharp

sharp_flat_vct = sharp(flat_vct, name='Our first sharpened cochain')

if isinstance(sharp_flat_vct, Cochain) & sharp_flat_vct.isdual & sharp_flat_vct.isvectorvalued & (sharp_flat_vct.dim == 0):
    print(f'{sharp_flat_vct.name} is indeed a proper discrete vector field.')

if sharp(dual_flat_vct) is None:
    print(f'WARNING: Sharpening a dual cochain is not possible yet.')

fig = pv.Plotter(shape=(1,3), border=False)

fig.subplot(0,0)
visualize(vct, fig=fig, display=False, 
          layout_parameters={'title':'Vector field'})

fig.subplot(0,1)
visualize(flat_vct, fig=fig, scaling_factor=10, display=False, 
          layout_parameters={'title':'Flat(vector field)'})

fig.subplot(0,2)
visualize(sharp_flat_vct, fig=fig, display=False, 
          layout_parameters={'title':'Sharp(flat(vector field))'})

fig.link_views()
fig.show()

import numpy.linalg as lng

err = [lng.norm(v2 - v1)/ lng.norm(v1) 
       for v1, v2 in zip(vct.values, sharp_flat_vct.values)]

print(f'The relative error between the amplitudes of the vectors and their sharp(flat()) counterparts is {np.mean(err):.2%} +/- {np.std(err):.2%}.')

from dxtr.operators import wedge
from dxtr.cochains import Cochain
from dxtr.complexes import icosahedron

mfld = icosahedron(manifold=True)

tol = 1e-15

N1 = mfld[1].size
N2 = mfld[2].size

cchn_1 = Cochain(mfld, 1, values= np.random.random(N1))
cchn_11 = wedge(cchn_1, cchn_1)

if (np.abs(cchn_11.values) < tol).all():
    print(f'The wedge of a {cchn_1.name} with it self yields a {cchn_11.name} of max value: {cchn_11.values.max():.2e}\n')

cchn_2 = Cochain(mfld, 1, values= 3*np.ones(N1))

cchn_12 = wedge(cchn_1, cchn_2)
cchn_21 = wedge(cchn_2, cchn_1)

if (np.abs(cchn_12.values + cchn_21.values) < tol).all():
    print('The wedge between two 1-cochains is indeed antisymmetric.')

N0 = mfld[0].size

cchn_0 = Cochain(mfld, 0, np.random.random(N0))
cchn_2 = Cochain(mfld, 2, np.random.random(N2))

cchn_02 = wedge(cchn_0, cchn_2)

error = []
for sidx, val in cchn_02.items():
    cfids = mfld.faces(2,0)[sidx]
    xpct_val = cchn_0.values[cfids].mean() * cchn_2.values[sidx]
    error.append(np.abs(val - xpct_val))

if (np.asarray(error) < tol).all():
    print('Wedge with a 0-cochain works as expected.')