symplectic basis for triangulated 3 manifolds (#117)

cython/SnapPy.pxi

@@ -705,6 +705,9 @@ cdef extern from "SnapPea.h":
     extern void two_bridge(c_Triangulation *manifold, Boolean *is_two_bridge, long int *p, long int *q) except *
     extern Real volume(c_Triangulation *manifold, int *precision) except *
     extern Boolean mark_fake_cusps(c_Triangulation   *manifold) except *
+
+    extern int** get_symplectic_basis(c_Triangulation *manifold, int *, int *, int) except *
+    extern void free_symplectic_basis(int **, int) except *
     extern void register_callbacks(void (*begin_callback)(),
                                    void (*middle_callback)(),
                                    void (*end_callback)())

cython/core/triangulation.pyx

@@ -3030,3 +3030,90 @@ cdef class Triangulation():
                 ignore_curve_orientations = ignore_curve_orientations,
                 ignore_orientation = ignore_orientation)
         return self._cache.save(result, 'triangulation_isosig', *args)
+
+    def symplectic_basis(self, verify=False):
+        """
+        Extend the Neumann-Zagier Matrix to one which is symplectic (up to factors of 2)
+        using oscillating curves. Verify parameter explicitly tests if the resulting matrix is symplectic.
+        Only accepts triangulations with 1 cusp.
+
+        >>> M = Manifold("4_1")
+        >>> M.symplectic_basis()
+        [-1  0 -1 -1]
+        [ 2  0 -2  0]
+        [-2 -1 -2 -1]
+        [ 0 -1 -2 -1]
+
+        <https://arxiv.org/abs/2208.06969>
+        """
+        def is_symplectic(M):
+            """
+            Test if the matrix M is symplectic
+            :param M: square matrix
+            :return: true or false
+            """
+            n = len(M)
+
+            for i in range(n):
+                for j in range(i, n):
+                    omega = abs(symplectic_form(M[i], M[j]))
+
+                    if i % 2 == 0 and j % 2 == 1 and j == i + 1:
+                        if omega != 2:
+                            return False
+                    elif omega:
+                        return False
+
+            return True
+
+        def symplectic_form(u, v):
+            return sum([u[2 * i] * v[2 * i + 1] - u[2 * i + 1] * v[2 * i] for i in range(len(u) // 2)])
+
+        cdef int **c_eqns;
+        cdef int **g_eqns;
+        cdef int num_rows, num_cols, dual_rows;
+        cdef int* eqn;
+
+        if self.c_triangulation is NULL:
+            raise ValueError('The Triangulation is empty.')
+
+        # current get_symplectic_eqns() implementation requires 1 cusp
+        if self.num_cusps() > 1:
+            raise ValueError('Triangulation contains {} cusps, only accepts triangulations with 1 cusp'.format(self.num_cusps()))
+
+        eqns = []
+
+        peripheral_curves(self.c_triangulation)
+
+        # Cusp Equations
+        for i in range(self.num_cusps()):
+            cusp_info = self.cusp_info(i)
+            if cusp_info.is_complete:
+                to_do = [(1,0), (0,1)]
+            else:
+                to_do = [cusp_info.filling]
+            for (m, l) in to_do:
+                eqn = get_cusp_equation(self.c_triangulation,
+                                        i, int(m), int(l), &num_rows)
+                eqns.append([eqn[j] for j in range(num_rows)])
+                free_cusp_equation(eqn)
+
+        # Dual Curve Equations
+        g_eqns = get_symplectic_basis(self.c_triangulation, &dual_rows, &num_cols, 0)
+
+        for i in range(dual_rows):
+            eqns.append([g_eqns[i][j] for j in range(num_cols)])
+
+        free_symplectic_basis(g_eqns, dual_rows)
+
+        # Convert to Neumann Zagier Matrix
+        rows = len(eqns)
+        retval = [[eqns[i][3 * (j // 2) + j % 2] - eqns[i][3 * (j // 2) + 2] for j in range(rows)] for i in range(rows)]
+
+        if verify:
+            if is_symplectic(retval):
+                print("Result is symplectic (up to factors of 2)")
+            else:
+                print("Warning: Result is not symplectic")
+
+        return matrix(retval)

dev/symplectic_basis/CMakeLists.txt

@@ -0,0 +1,63 @@
+# To compile project 
+# - run 'cmake -S . -B build/' in the dev/symplectic_basis dir 
+# - run 'cmake --build build/'
+
+# For clangd copy 'dev/symplectic_basis/build/compile_commands.json' to 'build/'
+
+cmake_minimum_required(VERSION 3.12)
+set(CMAKE_C_COMPILER "gcc")
+set(CMAKE_CXX_COMPILER "g++")
+
+project(SnapPea)
+set(CMAKE_EXPORT_COMPILE_COMMANDS 1)
+
+# Set the name of your program
+set(YOUR_PROGRAM symplectic_basis)
+
+# Set the location of the SnapPea kernel code
+set(SNAPPEA_KERNEL ${PROJECT_SOURCE_DIR}/../../kernel)
+
+# Compiler options
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -g -Wall -fanalyzer")
+
+# Add include directories
+include_directories(${SNAPPEA_KERNEL}/addl_code ${SNAPPEA_KERNEL}/headers ${SNAPPEA_KERNEL}/real_type ${SNAPPEA_KERNEL}/unix_kit ${SNAPPEA_KERNEL}/)
+
+# Add source files
+file(GLOB KERNEL_SOURCES ${SNAPPEA_KERNEL}/*/*.c)
+file(GLOB HEADER_FILES ${SNAPPEA_KERNEL}/*/*.h)
+file(COPY ${CMAKE_SOURCE_DIR}/CuspedCensusData DESTINATION ${CMAKE_SOURCE_DIR}/build)      # Copy census data to be used by program
+
+# Create executable target
+add_executable(${YOUR_PROGRAM} ${YOUR_PROGRAM}_main.c ${KERNEL_SOURCES} ${HEADER_FILES} ${CUSPED_CENSUS})
+
+# Link math library
+target_link_libraries(${YOUR_PROGRAM} m)
+
+# Define custom target to create BuildDate file
+#add_custom_target(BuildDate ALL
+#        COMMAND date > ${CMAKE_SOURCE_DIR}/BuildDate
+#        DEPENDS ${KERNEL_SOURCES} ${HEADER_FILES}
+#        WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+#        )
+
+# testing binary
+add_executable(${YOUR_PROGRAM}_test ${YOUR_PROGRAM}_test.c ${KERNEL_SOURCES} ${HEADER_FILES} ${CUSPED_CENSUS})
+target_link_libraries(${YOUR_PROGRAM}_test m)
+
+# enable testing functionality
+enable_testing()
+
+# define tests
+add_test(
+        NAME symplectic_test
+        COMMAND $<TARGET_FILE:${YOUR_PROGRAM}_test>
+)
+
+# Clean target
+add_custom_target(clean-all
+        COMMAND ${CMAKE_COMMAND} -E remove ${YOUR_PROGRAM} ${YOUR_PROGRAM}.o
+        COMMAND ${CMAKE_COMMAND} -E remove_directory KernelObjects
+        COMMAND ${CMAKE_COMMAND} -E remove *.pyc
+        )
+

dev/symplectic_basis/file_formats/GeneratorsFileFormat

@@ -0,0 +1,111 @@
+				SnapPea Matrix Generators File Format
+
+A generators file must begin with the line
+
+	% Generators
+
+to distinguish it from a triangulation file, which begins with
+"% Triangulation", or a link projection file which begins with
+"% Link Projection".  Next comes an integer telling how many matrices
+are present.  The matrices may by in either O(3,1) or SL(2,C).
+Orientation-reversing generators are allowed in O(3,1) but not in
+PSL(2,C).  If the matrices are in SL(2,C), read_generators() will
+convert them to O(3,1).  read_generators() can tell which format
+you are using by comparing the total number of matrix entries to
+the total number of matrices.  (SL(2,C) matrices contain 8 real entries
+each, while O(3,1) matrices contain 16 real entries each.)
+
+In PSL(2,C), the entries of a matrix
+
+                a   b
+                c   d
+
+should be written as
+
+        Re(a)    Im(a)
+        Re(b)    Im(b)
+        Re(c)    Im(c)
+        Re(d)    Im(d).
+
+Actually, the arrangement of the white space (blanks, tabs and returns)
+is irrelevant, so if you prefer you may write a PSL(2,C) matrix as, say,
+
+        Re(a)    Im(a)        Re(b)    Im(b)
+        Re(c)    Im(c)        Re(d)    Im(d).
+
+In O(3,1) the entries of each matrix should be written as
+
+        m00  m01  m02  m03
+        m10  m11  m12  m13
+        m20  m21  m22  m23
+        m30  m31  m32  m33
+
+where the 0-th coordinate is the timelike one.  Again, the arrangement
+of the white space is irrelevant.
+
+Here are two sample files.
+
+Sample #1.   PSL(2,C) generators for the Borromean rings complement.
+
+% Generators
+6
+   0.000000000000000    0.000000000000000
+   0.000000000000000   -1.000000000000000
+   0.000000000000000   -1.000000000000000
+   2.000000000000000    0.000000000000000
+
+   0.000000000000000    0.000000000000000
+   0.000000000000000    1.000000000000000
+   0.000000000000000    1.000000000000000
+   2.000000000000000    0.000000000000000
+
+   1.000000000000000   -1.000000000000000
+   0.000000000000000   -1.000000000000000
+   0.000000000000000    1.000000000000000
+   1.000000000000000    1.000000000000000
+
+   1.000000000000000   -1.000000000000000
+   0.000000000000000    1.000000000000000
+   0.000000000000000   -1.000000000000000
+   1.000000000000000    1.000000000000000
+
+   1.000000000000000    0.000000000000000
+  -2.000000000000000    0.000000000000000
+   0.000000000000000    0.000000000000000
+   1.000000000000000    0.000000000000000
+
+   1.000000000000000    0.000000000000000
+   0.000000000000000    0.000000000000000
+  -2.000000000000000    0.000000000000000
+   1.000000000000000    0.000000000000000
+
+Sample #2.  O(3,1) generators for a mirrored regular ideal tetrahedron.
+
+% Generators
+
+4
+
+ 1.25        -0.433012    -0.433012    -0.433012
+ 0.433012     0.25        -0.75        -0.75    
+ 0.433012    -0.75         0.25        -0.75    
+ 0.433012    -0.75        -0.75         0.25    
+
+ 1.25        -0.433012    +0.433012    +0.433012
+ 0.433012     0.25        +0.75        +0.75    
+-0.433012    +0.75         0.25        -0.75    
+-0.433012    +0.75        -0.75         0.25    
+
+ 1.25        +0.433012    -0.433012    +0.433012
+-0.433012     0.25        +0.75        -0.75    
+ 0.433012    +0.75         0.25        +0.75    
+-0.433012    -0.75        +0.75         0.25    
+
+ 1.25        +0.433012    +0.433012    -0.433012
+-0.433012     0.25        -0.75        +0.75    
+-0.433012    -0.75         0.25        +0.75    
+ 0.433012    +0.75        +0.75         0.25    
+
+(Note:  I truncated sqrt(3)/4 = 0.433012701892219323... to 0.433012
+to fit the above matrices within the width of this window.
+If you want to try out this example, please restore the
+high-precision value.)