Merge pull request #298 from saraedum/saf

Expose SAF invariant on flow components

doc/news/saf.rst

@@ -0,0 +1,7 @@
+**Added:**
+
+* Added `FlowComponent::safInvariant()` to expose the SAF invariant of the underlying Interval Exchange Transformation.
+
+**Fixed:**
+
+* Fixed access to intervalxt objects through pyflatsurf, e.g., `FlowComoponent::dynamicalComponent`.

libflatsurf/flatsurf/flow_component.hpp

@@ -1,7 +1,7 @@
 /**********************************************************************
  *  This file is part of flatsurf.
  *
- *        Copyright (C) 2019-2020 Julian Rüth
+ *        Copyright (C) 2019-2022 Julian Rüth
  *
  *  Flatsurf is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -27,6 +27,7 @@
 #include <iosfwd>
 #include <list>
 #include <vector>
+#include <gmpxx.h>
 
 #include "copyable.hpp"
 
@@ -81,6 +82,8 @@ class FlowComponent : boost::equality_comparable<FlowComponent<Surface>> {
 
   FlowTriangulation<Surface> triangulation() const;
 
+  std::vector<mpq_class> safInvariant() const;
+
   const intervalxt::Component& dynamicalComponent() const;
 
   [[deprecated("This method returns a inconsistent version of the IET from which this component was created originally. Use dynamicalComponent().iet() instead.")]] const IntervalExchangeTransformation<FlatTriangulationCollapsed<T>>& intervalExchangeTransformation() const;

libflatsurf/src/flow_component.cc

@@ -1,7 +1,7 @@
 /**********************************************************************
  *  This file is part of flatsurf.
  *
- *        Copyright (C) 2019-2021 Julian Rüth
+ *        Copyright (C) 2019-2022 Julian Rüth
  *
  *  Flatsurf is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -26,6 +26,7 @@
 #include <intervalxt/decomposition_step.hpp>
 #include <intervalxt/fmt.hpp>
 #include <intervalxt/label.hpp>
+#include <intervalxt/interval_exchange_transformation.hpp>
 #include <memory>
 #include <ostream>
 #include <unordered_map>
@@ -414,6 +415,11 @@ FlowDecomposition<Surface> FlowComponent<Surface>::decomposition() { return Impl
 template <typename Surface>
 const FlowDecomposition<Surface> FlowComponent<Surface>::decomposition() const { return ImplementationOf<FlowDecomposition<Surface>>::make(self->state); }
 
+template <typename Surface>
+std::vector<mpq_class> FlowComponent<Surface>::safInvariant() const {
+  return dynamicalComponent().iet().safInvariant();
+}
+
 template <typename Surface>
 FlowComponent<Surface> ImplementationOf<FlowComponent<Surface>>::make(std::shared_ptr<FlowDecompositionState<Surface>> state, FlowComponentState<Surface>* component) {
   return FlowComponent<Surface>(PrivateConstructor{}, state, component);

libflatsurf/test/flow_decomposition.test.cc

@@ -85,11 +85,50 @@ TEMPLATE_TEST_CASE("Flow Decomposition", "[FlowDecomposition]", (long long), (mp
     CAPTURE(flowDecomposition);
     REQUIRE(area(flowDecomposition.components()) == surface->area());
 
+    // Return the sum of the SAF invariants of the components that make up
+    // a flow decomposition.
+    // Note that the length of the SAF invariants can differ because some
+    // components can be rational and therefore their SAF invariant is
+    // trivial.
+    const auto safSum = [](const auto& components) {
+      std::vector<mpq_class> sum;
+      for (const auto& component : components) {
+        const auto saf = component.safInvariant();
+        if (saf.size() == 0)
+          continue;
+
+        if (sum.size() == 0) {
+          sum = saf;
+          continue;
+        }
+
+        REQUIRE(sum.size() == saf.size());
+
+        for (int i = 0; i < sum.size(); i++)
+          sum[i] += saf[i];
+      }
+
+      return sum;
+    };
+
+    auto saf = safSum(flowDecomposition.components());
+
     REQUIRE(flowDecomposition.decompose());
 
     CAPTURE(flowDecomposition);
     REQUIRE(area(flowDecomposition.components()) == surface->area());
 
+    SECTION("The SAF Invariants are Consistent") {
+      auto saf2 = safSum(flowDecomposition.components());
+
+      if (saf.size() == 0)
+        saf = std::vector<mpq_class>(saf2.size());
+      if (saf2.size() == 0)
+        saf2 = std::vector<mpq_class>(saf.size());
+
+      REQUIRE_THAT(saf2, Catch::Matchers::Equals(saf));
+    }
+
     SECTION("The Area of the Cylinders is Consistent with Their Width & Height (modulo some scaling)") {
       for (const auto& component : flowDecomposition.components()) {
         if (component.cylinder()) {

pyflatsurf/configure.ac

@@ -32,6 +32,7 @@ AS_IF([test "x$with_pytest" != "xno" && test "x$have_python" = "xyes"],
        with_pytest=yes
        AS_IF([$PYTHON -c 'import pytest'], , AC_MSG_ERROR([pytest for make check not found; run --without-pytest to disable Python tests in make check]))
        AS_IF([$PYTHON -c 'import pyexactreal'], , AC_MSG_ERROR([pyexactreal for make check not found; run --without-pytest to disable Python tests in make check]))
+       AS_IF([$PYTHON -c 'import pyintervalxt'], , AC_MSG_ERROR([pyintervalxt for make check not found; run --without-pytest to disable Python tests in make check]))
        AS_IF([$PYTHON -c 'import cppyy'], , AC_MSG_ERROR([cppyy for make check not found; run --without-pytest to disable Python tests in make check]))
        AS_IF([$PYTHON -c 'import cppyythonizations'], , AC_MSG_ERROR([cppyythonizations for make check not found; run --without-pytest to disable Python tests in make check]))
       ], [])

pyflatsurf/environment.yml

@@ -14,6 +14,7 @@ dependencies:
   - ipywidgets
   - libtool
   - make
+  - pyintervalxt>=3.2.0,<4
   - pyexactreal>=3.0.1,<5
   - pytest
   - python

pyflatsurf/recipe/meta.yaml

@@ -24,13 +24,15 @@ requirements:
     - python
     - setuptools
     - gmpxxyy >=1.0.1,<2
+    - pyintervalxt >=3.2.0,<4
     - pyexactreal >=3.0.1,<4
     - cppyythonizations
     - libflatsurf
   run:
     - python
     - cppyy >=2.1.0,<3
     - cppyythonizations
+    - pyintervalxt >=3.2.0,<4
     - pyexactreal >=3.0.1,<4
     - boost-cpp
 

pyflatsurf/src/pyflatsurf/__init__.py

@@ -27,7 +27,7 @@
 #*********************************************************************
 #  This file is part of flatsurf.
 #
-#        Copyright (C) 2019-2020 Julian Rüth
+#        Copyright (C) 2019-2022 Julian Rüth
 #
 #  Flatsurf is free software: you can redistribute it and/or modify
 #  it under the terms of the GNU General Public License as published by
@@ -43,6 +43,11 @@
 #  along with flatsurf. If not, see <https://www.gnu.org/licenses/>.
 #*********************************************************************
 
+# The API of flatsurf exposes intervalxt objects such as
+# intervalxt::Component. We therefore need the intervalxt headers to
+# be loaded so cppyy can make these parts of the API availalbe.
+import pyintervalxt
+
 from .cppyy_flatsurf import flatsurf
 
 from .factory import make_surface as Surface

pyflatsurf/src/pyflatsurf/cppyy_flatsurf.py

@@ -109,6 +109,9 @@ def unwrapTracked(self):
 # We have to work around issues with complex std::function parameters in cppyy
 cppyy.py.add_pythonization(filtered(re.compile("FlatTriangulation<.*>"))(add_method("isomorphism")(lambda self, other, kind=None, filter_matrix=lambda a,b,c,d: a == 1 and b == 0 and c == 0 and d == 1, filter_map=lambda a, b: True: cppyy.gbl.flatsurf.isomorphism[type(self).Coordinate.__cpp_name__](self, other, kind or cppyy.gbl.flatsurf.ISOMORPHISM.FACES, lambda m: filter_matrix(m.a, m.b, m.c, m.d), filter_map))), "flatsurf")
 
+# Return the SAF invariant as a Python list
+cppyy.py.add_pythonization(filtered(re.compile("FlowComponent<.*>"))(wrap_method("safInvariant")(lambda self, cpp: list(cpp()))), "flatsurf")
+
 for path in os.environ.get('PYFLATSURF_INCLUDE','').split(':'):
     if path: cppyy.add_include_path(path)
 

pyflatsurf/test/flow_decomposition.py

@@ -4,7 +4,7 @@
 ######################################################################
 # This file is part of flatsurf.
 #
-#       Copyright (C) 2020 Julian Rüth
+#       Copyright (C) 2020-2022 Julian Rüth
 #
 # flatsurf is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -67,6 +67,8 @@ def test_hexagon_eantic():
                 assert area1 == 2 * area2
                 assert area3 == area2 * (vx*vx + vy*vy)
 
+            assert component.safInvariant() == [0]
+
         assert all(component.cylinder() for component in decomposition.cylinders())
         assert len(decomposition.cylinders()) == n_cylinders
 
@@ -93,16 +95,21 @@ def test_D33():
         assert str(decomposition).startswith("FlowDecomposition")
         decomposition.decompose(-1)
         n_cylinders, n_without_periodic_trajectory, n_indeterminate = 0, 0, 0
-        for dec in decomposition.components():
-            if dec.cylinder() == True:
+        for component in decomposition.components():
+            if component.cylinder() == True:
                 n_cylinders += 1
-                h = dec.circumferenceHolonomy()
+                h = component.circumferenceHolonomy()
                 assert h.x() or h.y()
                 assert h.x() * v.y() == h.y() * v.x()
-            elif dec.cylinder() == False:
+
+                assert all(saf == 0 for saf in component.safInvariant()), component.safInvariant()
+            elif component.cylinder() == False:
                 n_without_periodic_trajectory += 1
+
+                assert any(saf != 0 for saf in component.safInvariant()), component.safInvariant()
             else:
                 n_indeterminate = 0
+
         assert n_indeterminate == 0
         decompositions[(n_cylinders, n_without_periodic_trajectory)] += 1