musical

edge_nsimplex_couples(complex)

Couples together edge ids with their highest-dimensional cofaces ids.

Parameters:

Name	Type	Description	Default
complex	SimplicialComplex	The SimplicialComplex to work on.	required

Returns:

Type	Description
ndarray	A Nc * 2 arrays each row corresponds to a couple (eid, nsid).

Note

Each edge id (eid) can appear multiple times, depending on its n-valence (i.e. number of n-cofaces). The number of row corresponds to the sum of the edge n-valences for all edges.

Source code in src/dxtr/operators/musical.py

def edge_nsimplex_couples(complex:SimplicialComplex)->np.ndarray[int]:
    """
    Couples together edge ids with their highest-dimensional cofaces ids.

    Parameters
    ----------
    complex : SimplicialComplex
        The SimplicialComplex to work on.

    Returns
    -------
    np.ndarray
        A Nc * 2 arrays each row corresponds to a couple (eid, nsid). 

    Note
    ----
    Each edge id (eid) can appear multiple times, depending on 
    its n-valence (i.e. number of n-cofaces). The number of row
    corresponds to the sum of the edge n-valences for all edges.
    """
    n = complex.dim

    eid_nsid_couples = [np.array([[eid]*len(sids), sids]).T
                        for eid, sids in enumerate(complex.cofaces(1, n))]

    return np.vstack(eid_nsid_couples)

edge_nsimplex_weight(on, for_operator='dual_to_primal_flat')

Computes the edge_nsimplex_weight for each edge-n-Simplex couple.

Parameters:

Name	Type	Description	Default
on	SimplicialManifold	The simplicial manifold we work on.	required
for_operator	str	A flag to change the normalization factor depending on the use case. Should either be 'dual_to_primal_flat', '_dual_to_dual_flat', 'sharp' or 'div'. Default is 'dual_to_primal_flat'.	'dual_to_primal_flat'

Returns:

Type	Description
csr_matrix	A (N1 x Nn)-sparse matrix that contains the sought edge_nsimplex_weight.

Notes

All edges (i.e. shared by multiple n-cofaces) have one weight for each of their n-cofaces. * The values of these weights depends on the operator to compute, c.f. See also section for more details. * In the sharp and flat cases: For each edge, the sum of its edge_nsimplex_weight must egals one. * This condition seems validated only if the n-simplices are well-centered. * With this definition, the divergence theorem for vector fields defined as vector-values dual 0-forms is verified.

See also

Hirani's PhD Thesis for more details: * For the sharp weight formula: Proposition 5.5.1, p.51. * For the flat weight formula: Definition 5.8.1, p.58. * For the div weight formula: Lemma 6.1.2 p.60.

Source code in src/dxtr/operators/musical.py

def edge_nsimplex_weight(on: SimplicialManifold, 
                         for_operator: str = 'dual_to_primal_flat'
                         ) -> sp.csr_matrix[float]:
    """
    Computes the edge_nsimplex_weight for each edge-n-Simplex couple.

    Parameters
    ----------
    on : SimplicialManifold
        The simplicial manifold we work on.
    for_operator : str, optional
        A flag to change the normalization factor depending on the use case.
        Should either be 'dual_to_primal_flat', '_dual_to_dual_flat', 'sharp' 
        or 'div'. Default is 'dual_to_primal_flat'.

    Returns
    -------
    sp.csr_matrix
        A (N1 x Nn)-sparse matrix that contains the sought edge_nsimplex_weight.

    Notes
    -----
    All edges (i.e. shared by multiple n-cofaces) have one weight for 
        each of their n-cofaces. 
      * The values of these weights depends on the operator to compute, 
        c.f. See also section for more details.
      * In the sharp and flat cases: For each edge, the sum of its
        edge_nsimplex_weight must egals one.
      * This condition seems validated only if the n-simplices 
        are well-centered.
      * With this definition, the divergence theorem for vector fields defined 
        as vector-values dual 0-forms is verified.

    See also
    --------
    Hirani's PhD Thesis for more details:
      * For the sharp weight formula:  Proposition 5.5.1, p.51. 
      * For the flat weight formula: Definition 5.8.1, p.58.
      * For the div weight formula: Lemma 6.1.2 p.60.
    """
    mfld, operator = on,  for_operator

    n = mfld.dim
    N1 = mfld.shape[1]
    Nn = mfld.shape[-1]
    chn_cplx = mfld._chain_complex

    if operator == 'sharp':
        splx_ids_chains = complex_simplicial_partition(chn_cplx)[0]
        ccenters = [mfld[i].circumcenters for i in range(n+1)]
        top_simplex_volumes = mfld[-1].volumes

    else:
        k = n-1 if operator=='dual_to_dual_flat' else 1

        splx_ids_chains = complex_simplicial_partition(chn_cplx,lowest_dim=k)[0]
        ccenters = [mfld[i].circumcenters for i in range(k, n+1)] 
        edge_covolumes = mfld[k].covolumes

    weight_matrix = sp.lil_matrix((N1, Nn), dtype=float)
    for (eid, nsid) in edge_nsimplex_couples(mfld):
        if operator in ['dual_to_primal_flat', 'dual_to_dual_flat', 'div']:

            ids = np.where((splx_ids_chains[:, 0] == eid) & 
                           (splx_ids_chains[:,-1] == nsid))

            covolume_intersection = volume_polytope(ccenters, 
                                                    splx_ids_chains[ids])

            weight_matrix[eid, nsid] = covolume_intersection

            if operator != 'div':
                weight_matrix[eid, nsid] /= edge_covolumes[eid]

        elif operator == 'sharp':
            ids = np.where((splx_ids_chains[:, 1] == eid) & 
                           (splx_ids_chains[:,-1] == nsid))

            support_volume_intersection = volume_polytope(ccenters, 
                                                        splx_ids_chains[ids])

            weight_matrix[eid, nsid] = support_volume_intersection \
                / top_simplex_volumes[nsid]

        else:
            logger.warning(
                'Wrong operator name. Should be either sharp, flat or div.')
            return None

    return weight_matrix.tocsr()

flat(vector_field, manifold=None, dual=False, name='1-cochain, flat from vector field')

Transforms a vector field into a 1-Cochain.

Parameters:

Name	Type	Description	Default
vector_field	ndarray or Cochain	The vector field to convert.	required
manifold	SimplicialManifold	The Simplicial Manifold where the vector field is defined. If the provided vector field is a vector-valued Cochain, this argument is not needed. Default is None.	None
dual	bool	If True, a dual 1-Cochain is returned; otherwise, it is a primal one. Default is False.	False
name	str	The name to give to the flattened cochain. Default is '1-cochain, flat from vector field'.	'1-cochain, flat from vector field'

Returns:

Type	Description
Cochain or None	The resulting 1-cochain defined on the prescribed manifold.

Notes

The input vector field must be defined on the circumcenters of the top simplices of the considered SimplicialManifold.

This implementation of the discrete flat operator is derived from the definitions provided by Anil Hirani in his PhD manuscript.

We chose to implement two types of discrete flat: one returns a primal 1-cochain, the other a dual 1-cochain. The dual argument enables the selection of the desired one.

TODO: Add the ability to flatten a vector field defined on the 0-simplices.

See Also

dual_to_primal_flat : Transforms a vector field into a primal 1-cochain. dual_to_dual_flat : Transforms a vector field into a dual 1-cochain.

Source code in src/dxtr/operators/musical.py

def flat(vector_field: np.ndarray[float] | Cochain, 
         manifold: Optional[SimplicialManifold] = None,
         dual: bool = False,
         name: Optional[str] = '1-cochain, flat from vector field'
         ) -> Optional[Cochain]:
    """
    Transforms a vector field into a 1-`Cochain`.

    Parameters
    ----------
    vector_field : np.ndarray or Cochain
        The vector field to convert.
    manifold : SimplicialManifold, optional
        The Simplicial Manifold where the vector field is defined. If the 
        provided vector field is a vector-valued `Cochain`, this argument 
        is not needed. Default is None.
    dual : bool, optional
        If True, a dual 1-`Cochain` is returned; otherwise, it is a primal one. 
        Default is False.
    name : str, optional
        The name to give to the flattened cochain. Default is '1-cochain, flat from vector field'.

    Returns
    -------
    Cochain or None
        The resulting 1-cochain defined on the prescribed manifold.

    Notes
    -----
    The input vector field must be defined on the circumcenters of the 
    top simplices of the considered `SimplicialManifold`.

    This implementation of the discrete flat operator is derived from
    the definitions provided by Anil Hirani in his PhD manuscript.

    We chose to implement two types of discrete flat: one returns a primal 
    1-cochain, the other a dual 1-cochain. The `dual` argument enables the 
    selection of the desired one.

    TODO: Add the ability to flatten a vector field defined on 
    the 0-simplices.

    See Also
    --------
    dual_to_primal_flat : Transforms a vector field into a primal 1-cochain.
    dual_to_dual_flat : Transforms a vector field into a dual 1-cochain.
    """
    if isinstance(vector_field, Cochain):
        manifold = vector_field.complex
        vector_field = vector_field.values

    if not _isvalid_for_flat(vector_field, manifold):
        return None
    elif dual:
        return _dual_to_dual_flat(vector_field, manifold, name)
    else:
        return _dual_to_primal_flat(vector_field, manifold, name)        

interpolate_cochain_on_nsimplices(cochain)

Interpolates a k-cochain on the circumcenters of the top-simplices.

Parameters:

Name	Type	Description	Default
cochain	Cochain	The cochain to interpolate.	required

Returns:

Type	Description
ndarray	(Nn, Nk, D)-array containing the corresponding k-vector field. Nn: number of n-simplices, Nk: number of k-simplices, D: embedding dimension.

Notes

In the computation we added a factor 2 in order to get the proper results, but we don't know how or why it should be here. We just noticed that without it the sharp(flat(vector_field)) ended up being scaled by a factor 2 for regular complexes.

TODO: understand where this factor 2 comes from, maybe: look into the divergence theorem, see chapter 6 of Hirani's thesis.

Source code in src/dxtr/operators/musical.py

def interpolate_cochain_on_nsimplices(cochain: Cochain) -> np.ndarray[float]:
    """
    Interpolates a k-cochain on the circumcenters of the top-simplices.

    Parameters
    ----------
    cochain : Cochain
        The cochain to interpolate.

    Returns
    -------
    np.ndarray
        (Nn, Nk, D)-array containing the corresponding k-vector field.
        Nn: number of n-simplices, Nk: number of k-simplices, 
        D: embedding dimension.

    Notes
    -----
    In the computation we added a factor 2 in order to get the proper 
    results, but we don't know how or why it should be here. 
    We just noticed that without it the sharp(flat(vector_field)) 
    ended up being scaled by a factor 2 for regular complexes.

    TODO: understand where this factor 2 comes from, maybe:
    look into the divergence theorem, see chapter 6 of Hirani's thesis.
    """
    k = cochain.dim
    coef = cochain.values
    cplx = cochain.complex

    D = cplx.emb_dim
    Nn = cplx.shape[-1]
    Nk = cplx.shape[k]

    w_map = WhitneyMap.of(cplx, k, normalized=True)
    weights = edge_nsimplex_weight(cplx, 'sharp')

    kvct_shape= (D,D) if cochain.isvectorvalued else (1,D)
    kvector_field = np.zeros((Nn, Nk, *kvct_shape)) 
    for (i, j), base_vector in zip(w_map.indices, w_map.base):
        kvector_field[i,j,:] += 2 * np.outer(coef[j], 
                                             weights.T[i,j] * base_vector)

    if not cochain.isvectorvalued: kvector_field = kvector_field[:,:,0,:]

    return kvector_field

project(vector_field, on, of)

Projects a vector field onto the edges of a simplicial/cellular complex.

Parameters:

Name	Type	Description	Default
vector_field	ndarray	The vector field to project.	required
on	str	Should either be 'primal' or 'dual'. Specifies on which complex the vector field should be projected: Either on the edges formed by the 1-simplices or the edges formed by the dual cells of the (n-1)-simplices.	required
of	SimplicialManifold	The simplicial manifold where to project the vector field.	required

Returns:

Type	Description
ndarray	An array of shape (Nk,); Nk = manifold[n-1].size containing the projected values.

Notes

The provided vector field must be defined as a vector-valued dual 0-cochain. Meaning that it must contain 1 vector per top-simplex of the considered SimplicialManifold.

Source code in src/dxtr/operators/musical.py

def project(vector_field: np.array[float], on: str, of: SimplicialManifold
            ) -> np.ndarray[float]:
    """
    Projects a vector field onto the edges of a simplicial/cellular complex.

    Parameters
    ----------
    vector_field : np.ndarray
        The vector field to project.
    on : str
        Should either be 'primal' or 'dual'. Specifies on which complex
        the vector field should be projected: Either on the edges formed
        by the 1-simplices or the edges formed by the dual cells of the 
        (n-1)-simplices.
    of : SimplicialManifold
        The simplicial manifold where to project the vector field.

    Returns
    -------
    np.ndarray
        An array of shape (Nk,); Nk = manifold[n-1].size 
        containing the projected values.

    Notes
    -----
    The provided vector field must be defined as a vector-valued dual 
    0-cochain. Meaning that it must contain 1 vector per top-simplex of 
    the considered `SimplicialManifold`.
    """
    manifold = of
    n = manifold.dim

    if on == 'primal':
        edges = edge_vectors(manifold, normalized=True) 
        weights = edge_nsimplex_weight(manifold)

        return np.einsum('ij,ij->i', edges, weights @ vector_field)

    elif on == 'dual':
        vectors = [vector_field[cfids] for cfids in manifold.cofaces(n-1)]
        half_edges = _dual_half_edges_vectors(manifold, normalized=True)
        weights = _dual_half_edge_nsimplex_weights(manifold)
        signs = _dual_edges_orientation(manifold, half_edges)

        coefs = [sign * np.einsum('ij, i, ij -> i', hedgs, wghts, vcts) 
                for sign, hedgs, wghts, vcts 
                in zip(signs, half_edges, weights, vectors)]

        return np.asarray([coef.sum() for coef in coefs])

sharp(cochain, name='Discrete vector field')

Transforms a primal 1-Cochain into a discrete vector field.

Parameters:

Name	Type	Description	Default
cochain	Cochain	The primal 1-Cochain to transform into a vector field.	required
name	str	The name to give to the output Cochain instance. Default is 'Discrete vector field'.	'Discrete vector field'

Returns:

Type	Description
Cochain or None	A dual vector-valued 0-Cochain corresponding to the sought vector field.

Notes

In the present version, we can only transform a primal 1-Cochain into a dual vector-valued 0-Cochain.

Source code in src/dxtr/operators/musical.py

def sharp(cochain: Cochain, 
          name: str = 'Discrete vector field') -> Optional[Cochain]:
    """
    Transforms a primal 1-`Cochain` into a discrete vector field.

    Parameters
    ----------
    cochain : Cochain
        The primal 1-`Cochain` to transform into a vector field.
    name : str, optional
        The name to give to the output `Cochain` instance. Default is 'Discrete vector field'.

    Returns
    -------
    Cochain or None
        A dual vector-valued 0-`Cochain` corresponding to the sought vector field.

    Notes
    -----
    In the present version, we can only transform a primal 1-`Cochain`
    into a dual vector-valued 0-`Cochain`.
    """
    if not _isvalid_for_sharp(cochain): return None

    vector_field = interpolate_cochain_on_nsimplices(cochain).sum(axis=1)

    return Cochain(cochain.complex, values=vector_field, 
                   dim=0, dual=True, name=name)