some code quarry for homology computation

Homology.cs

using System;
using System.Collections.Generic;

using System.Collections.ObjectModel;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

using MathNet.Numerics;
using System.Numerics;
using MathNet.Numerics.Providers.LinearAlgebra.Mkl;
using LibmatCore;
using Rhino.Geometry;

using System.Diagnostics;

namespace LibmatGeometry
{

    // https://www.ayasdi.com/blog/topology/measuring-shape/

    public class Homology
    {
        // holds a bunch of stuff

        // a tri mesh version
        public Mesh _M { get; private set; }


        // a custom version of the edges
        MeshTopology _MT;

        public List<Tuple<int, int>> myEdgeList { get; private set; }

        // the most important part: the simplicial complex' boundary operators
        public MathNet.Numerics.LinearAlgebra.Double.DenseMatrix d1 { get; private set; }

        // much faster than sparse matrix, for some reason
        public MathNet.Numerics.LinearAlgebra.Double.DenseMatrix d2 { get; private set; }


        // the Betti numbers
        public int[] betti { get; private set; }

        public String _laprovider;

        public MathNet.Numerics.LinearAlgebra.Vector<double>[] _harmonicOneChains;
        //   public DataFrame _U;
// #if experimental
        public double[] conformalCoordinate;
// #endif
        // exists on Mac only ..... no comment.
        public static Vector3d CrossProduct(Vector3d a, Vector3d b)
        {
            return new Vector3d(
                a.Y * b.Z - a.Z * b.Y,
                a.Z * b.X - a.X * b.Z,
                a.X * b.Y - a.Y * b.X
                );
        }


        // sets up the matrices computing the simplicial complex

        public Homology(Mesh M)
        {
            // sanitize mesh
            Control.UseNativeMKL();
            _laprovider = Control.LinearAlgebraProvider.ToString();

            M.Weld(0.1);
            M.Faces.ConvertQuadsToTriangles();
            M.UnifyNormals();
            M.Normals.ComputeNormals();

            this._M = M;


            // Rhino.Geometry.Collections.MeshTopologyEdgeList edgeList = M.TopologyEdges;
            _MT = new MeshTopology(this._M, true);
            myEdgeList = _MT._list;
            // System.Diagnostics.Debug.Assert(myEdgeList.Count == v1); // doesn't hold for big meshes

            // set up the boundary operator https://en.wikipedia.org/wiki/Chain_(algebraic_topology)
            // https://en.wikipedia.org/wiki/Simplicial_homology
            // the boundary of a one-chain is a zero-chain
            d1 = new MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(this._M.Vertices.Count, myEdgeList.Count);
            d2 = new MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(myEdgeList.Count, this._M.Faces.Count);

            for (int edgeIndex = 0; edgeIndex < myEdgeList.Count; edgeIndex++)
            {

                //Rhino.IndexPair ip = edgeList.GetTopologyVertices(edgeIndex);
                var myip = myEdgeList[edgeIndex];
                int firstIndex = myip.Item1;
                int secondIndex = myip.Item2;

                d1[firstIndex, edgeIndex] = 1;
                d1[secondIndex, edgeIndex] = -1;

#if extraDebug
                System.Diagnostics.Debug.WriteLine("In the first loop, edge index " + edgeIndex.ToString() + " has " +
                        "contributed a plus one to vertex " + firstIndex.ToString() + " " +
                        "with coordinates (" +
                        M.Vertices[firstIndex].X.ToString() + ", " +
                        M.Vertices[firstIndex].Y.ToString() + ", " +
                        M.Vertices[firstIndex].Z.ToString() + ") and a minus one to vertex " + secondIndex.ToString() +
                        " with coordinates (" +
                        M.Vertices[secondIndex].X.ToString() + ", " +
                        M.Vertices[secondIndex].Y.ToString() + ", " +
                        M.Vertices[secondIndex].Z.ToString() + ")."
                        );
#endif
            }

            // ok, now we can set up the matrix d2
            Rhino.Geometry.Collections.MeshFaceList m = this._M.Faces;

            // we need to have a fixed order, so we have to convert the face list to an actual list
            List<MeshFace> faces = m.ToList();

            for (int faceIndex = 0; faceIndex < this._M.Faces.Count; faceIndex++)
            {
                MeshFace f = faces[faceIndex];
#if extraDebug
                Debug.WriteLine("looking at face " + faceIndex.ToString() +
                    " with points" +
                    "(" + _M.Vertices[f.A].X.ToString() + " " + _M.Vertices[f.A].Y.ToString() + " " + _M.Vertices[f.A].Z.ToString() + ") " +
                    "(" + _M.Vertices[f.B].X.ToString() + " " + _M.Vertices[f.B].Y.ToString() + " " + _M.Vertices[f.B].Z.ToString() + ") " +
                    "(" + _M.Vertices[f.C].X.ToString() + " " + _M.Vertices[f.C].Y.ToString() + " " + _M.Vertices[f.C].Z.ToString() + ")." +
                    "(" + _M.Vertices[f.D].X.ToString() + " " + _M.Vertices[f.D].Y.ToString() + " " + _M.Vertices[f.D].Z.ToString() + ").")
                    ;
#endif
                // int[] edges = edgeList.GetEdgesForFace(faceIndex); // this one doesn't give correct edges. Bug in Rhino!!
                // instead, the following returns the correct code.
                var edges = new int[3];
                edges[0] = _MT.GetEdgeIndex(f.A, f.B);
                edges[1] = _MT.GetEdgeIndex(f.B, f.C);
                edges[2] = _MT.GetEdgeIndex(f.A, f.C);

                // have to assert all edges have been found
                System.Diagnostics.Debug.Assert(edges[0] != -1 && edges[1] != -1 && edges[2] != -1);

                foreach (int edgeIndex in edges)
                {
                    // have to keep track of the parity of the iteration of this loop
                    //    edgeCounter++;
                    //    bool edgeCounterOdd = (edgeCounter % 2 == 1);

                    // each edge generates an entry in the matrix representing the simplicial homology boundary operator.
                    // however, it is clumsy to find its sign ... plur or minus one.
                    // We have to decide whether the orientation of the face induces the same
                    // orientation on the edge as the one given by its associated indexPair.
                    // It does induce the right one if the scalar product between the face's normal with the
                    // crossed product of the vector from the first to the second index with the vector from the first
                    // to the center of the triangle is positive. Otherwise, it doesn't. 

                    var myip = myEdgeList[edgeIndex];
                    Point3d firstPointOnEdge = this._M.Vertices[myip.Item1];
                    Point3d secondPointOnEdge = this._M.Vertices[myip.Item2];
                    Vector3d firstToSecond = secondPointOnEdge - firstPointOnEdge;

                    // one could also use the face's third vertex here
                    Point3d faceCenter = new Point3d(
                        (_M.Vertices[f.A].X + _M.Vertices[f.B].X + _M.Vertices[f.C].X) / 3,
                        (_M.Vertices[f.A].Y + _M.Vertices[f.B].Y + _M.Vertices[f.C].Y) / 3,
                        (_M.Vertices[f.A].Z + _M.Vertices[f.B].Z + _M.Vertices[f.C].Z) / 3
                    );

                    Vector3d firstToCenter = faceCenter - firstPointOnEdge;

                    Vector3d crossProduct = CrossProduct(firstToSecond, firstToCenter);
                    Vector3d DirectionOfFaceNormal = CrossProduct(
                        ((Point3d)_M.Vertices[f.B] - (Point3d)_M.Vertices[f.A]),
                        ((Point3d)_M.Vertices[f.C] - (Point3d)_M.Vertices[f.A])
                        );

                    double scalarProduct = crossProduct * DirectionOfFaceNormal;

                    //   if the scalar product is zero, they are perpendicular which must not happen, by the geometry of the situation.
                    System.Diagnostics.Debug.Assert(scalarProduct != 0);
                    bool OrientationsCoincide = (scalarProduct > 0);

                    // finally, we can decide about the sign of the matrix entry.
                    int sign = (OrientationsCoincide) ? 1 : -1;

                    // done! We can finally set the corresponding matrix entry.                    
                    d2[edgeIndex, faceIndex] += sign;
#if extraDebug
                    System.Diagnostics.Debug.WriteLine("In the second loop, face index " +
                    faceIndex.ToString() + " with center (" +
                    faceCenter.X.ToString() + ", " +
                    faceCenter.Y.ToString() + ", " +
                    faceCenter.Z.ToString() + ") has " +
                    "contributed a sign " + sign.ToString() + " to edge index " + edgeIndex.ToString() + " " +
                    "which leads from (" +
                    firstPointOnEdge.X.ToString() + ", " +
                    firstPointOnEdge.Y.ToString() + ", " +
                    firstPointOnEdge.Z.ToString() + ") to (" +
                    secondPointOnEdge.X.ToString() + ", " +
                    secondPointOnEdge.Y.ToString() + ", " +
                    secondPointOnEdge.Z.ToString() + ").");
#endif
                }
            }
            // ok, now we have the de-Rham complex! That's already been the bulk of the work.


            // the number of connected components of the mesh must be equal to v0 minus the rank of d1,
            // because both are equal to the dimension of the zero-th simplicial homology group with real
            // coefficients.
            //            System.Diagnostics.Debug.Assert(M.DisjointMeshCount == v0 - rank);



            // now starts the computation-intensive part
            //IRHostSession _session = RHostSession.Create("one");

            //TU(d1, _session).Wait();
            int d1rank = d1.Rank();

            // int d1rank = qrr.Diagonal().Find(g => g < 0.00001).Item1;
            int d2rank = d2.Rank();

            //     int d1nullity = d1.Nullity();
            int d1nullity = this.myEdgeList.Count - d1rank;
            betti = new int[3];
            betti[0] = d1.RowCount - d1rank;

            // the dimension of the kernel of d2 is the dimension of Z2 minus the rank
            betti[2] = d2.ColumnCount - d2rank;

            // dimension of kernel of d1 minus dimension of image of d2
            betti[1] = d1nullity - d2rank;

            var dSquare = d1 * d2;

            double test = dSquare.FrobeniusNorm(); // has to be zero
            System.Diagnostics.Debug.Assert(test < 0.01);


            //  GetHomologyBasis


            //var projImaged1T = d1.Transpose() * (d1 * d1.Transpose()).PseudoInverse() * d1;
            // var projKerneld1 = MathNet.Numerics.LinearAlgebra.Double.SparseMatrix.CreateDiagonal(this.myEdgeList.Count, this.myEdgeList.Count, 1) - projImaged1T;
            //int kerneld1 = projKerneld1.Rank();
            //var homBasis = projKerneld1.Range();

            //int homBasisCount = homBasis.Length;
            //int a = 2;
            // we now need an ONB for the null space of d1.
            // since ker(A) = orthogonal complement of (image of t(A)), we first build the latter image

            // var s = d1.Svd();
            // we need to get the d1nullity last columns of s.VT
            // var kernelOfd1 = MathNet.Numerics.LinearAlgebra.

            // strategy: find a matrix whose column space is the kernel of d1
            // then, project each column onto the space orthogonal to the image of d2
            // this will not lead outside the kernel of d1
            // then, compute the .Range() of the matrix thus obtained. This is a homology basis.


            // we can then find a harmonic one-form

            //            var imaged1T = d1.Transpose().Range();
            // need to find its orthogonal complement

            var TopologicalEdgeLaplacian = d2 * d2.Transpose() + d1.Transpose() * d1; // formula (3) of http://www3.cs.stonybrook.edu/~gu/publications/pdf/2003/sgp_global_conformal_param.pdf

            //            var projectionKernelLaplacian = MathNet.Numerics.LinearAlgebra.Double.SparseMatrix.CreateDiagonal(this.myEdgeList.Count, this.myEdgeList.Count, 1) - TopologicalEdgeLaplacian.PseudoInverse() * TopologicalEdgeLaplacian;
            //            var hom = projectionKernelLaplacian.Range();
            //            int kkk = projectionKernelLaplacian.Rank();
            //            int hum = TopologicalEdgeLaplacian.Nullity();
            _harmonicOneChains = TopologicalEdgeLaplacian.Kernel();

            // out of interest: is it really harmonic, i.e., is it closed and co-closed?

            // d1 * _harmonicOneChains;
            // d2.Transpose() * _harmonicOneChains;


            // we can now set up the linear system of equations for the harmonic one-form.
            // thre requirements are: the cotangent-weighted laplacian has to be zero

            // then, it is harmonic itself, meaning that

            // 1) Closedness
            // secondly, it is closed meaning that for every face, its integral around the face is zero.

            // 2) Duality
            // its pairing with the harmonic one-chains has to be zero (duality)

            // 3) Harmonicity


            // var cotWeightedLaplacian = new MathNet.Numerics.LinearAlgebra.Double.SparseMatrix(this.myEdgeList.Count);
            // this gives a single requirement row
            // this gives a requirement for every face
            int twog = _harmonicOneChains.Length;

            //            var Requirement = new MathNet.Numerics.LinearAlgebra.Double.SparseMatrix(_M.Faces.Count + twog + 1, myEdgeList.Count);
            // to clad it with zeros
            var Requirement = new MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(_M.Faces.Count + twog + _M.Vertices.Count, myEdgeList.Count);



            // closedness
            // gives #F requirements
            int meshCounter = 0;

            foreach (MeshFace f in _M.Faces)
            {
                Requirement[meshCounter, _MT.GetEdgeIndex(f.A, f.B)] = 1;
                Requirement[meshCounter, _MT.GetEdgeIndex(f.B, f.C)] = 1;
                Requirement[meshCounter, _MT.GetEdgeIndex(f.A, f.C)] = 1;
                meshCounter++;
            }

            

            // 2) Duality
            // gives 2g requirements ... one for each cohomology class of degree 1
            for (int j = 0; j < twog; j++)
            {
                for (int i = 0; i < myEdgeList.Count; i++)
                {
                    Requirement[j + _M.Faces.Count, i] = _harmonicOneChains[j][i];
                }
            }

            // 3) Harmonicity in the sense of being in the kernel of the contangent-weighted Laplacian
            // gives a requirement for any vertex
            // this part of the requirement matrix corresponds to 
            MathNet.Numerics.LinearAlgebra.Double.DenseMatrix cotangents = _MT.EdgeCotangentWeight();

            for (int j = 0; j < _MT._list.Count; j++)
            {
                double cotgt = cotangents[_MT._list[j].Item1, _MT._list[j].Item2];
                Requirement[_MT._list[j].Item1 + _M.Faces.Count + twog, j] =  -cotgt;
                Requirement[_MT._list[j].Item2 + _M.Faces.Count + twog, j] =  -cotgt;
            }

            // in total, we get #F + 2g + #V = #E + 2 requirements (because of #V - #E + #F = 2 - 2g, the Euler characteristic), so these are just two requirements too many, which is not too bad.

            // int a = Requirement.Nullity();
            

            // we wanna do var x = A.Solve(b); for this, we need b
            var b = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(_M.Faces.Count + twog + _M.Vertices.Count);

            // try something dual to the first basis vector
            b[_M.Faces.Count] = 1;
            //            var firstConformalOneForm = Requirement.Solve(b);
            // var n = Requirement.Nullity();
            var someSolution = Requirement.QR().Solve(b);

            // #if experimentalvar someSolution = Requirement.QR().Solve(b);

//            var someSolution = Requirement.PseudoInverse() * b;

            
            conformalCoordinate = new double[_M.Vertices.Count];
            bool[] visited = new bool[_M.Vertices.Count];

            for (int i = 0; i < _M.Vertices.Count; i++) visited[i] = false;

            conformalCoordinate[0] = 0;
            visited[0] = true;

            // while there is any unvisited vertex left
            // works only if there is only one connected component ....
            while (visited.ToList().Any(g => !g)) {
                foreach (var e in _MT._enum)
                {
                    if (visited[e.Item1] & !visited[e.Item2])
                    {
                        conformalCoordinate[e.Item2] = conformalCoordinate[e.Item1] + someSolution[_MT.GetEdgeIndex(e.Item1, e.Item2)];
                        visited[e.Item2] = true;
                    }
                    if (visited[e.Item2] & !visited[e.Item1])
                    {
                        conformalCoordinate[e.Item1] = conformalCoordinate[e.Item2] - someSolution[_MT.GetEdgeIndex(e.Item1, e.Item2)];
                        visited[e.Item1] = true;
                    }
                }
            }
// #endif
        }
    }
}