Merge pull request #17 from esa/pickle-support

Pickle support

.github/workflows/build-and-test.yml

@@ -0,0 +1,65 @@
+name: Build & Test
+on:
+  # Triggers the workflow on push or pull request events but only for the master branch
+  push:
+    branches: [ main ]
+  pull_request:
+    branches: [ main ]
+  # Allows you to run this workflow manually from the Actions tab
+  workflow_dispatch:
+
+# A workflow run is made up of one or more jobs that can run sequentially or in parallel
+jobs:
+  # Build & Test on Linux
+  build-and-test-linux:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v3
+      - name: Build
+        run: |
+          mkdir build
+          cd build
+          cmake -DCMAKE_BUILD_TYPE=Release ..
+          cmake --build .
+      - name: Test
+        run: cd build/test && ./polyhedralGravity_test
+
+  # Install the python interface by building from source and run pytest
+  pytest-linux:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v3
+      - name: Install Ninja Build
+        run: |
+          sudo apt-get install ninja-build -y
+      - name: Install Conda environment from environment.yml
+        uses: mamba-org/provision-with-micromamba@main
+        with:
+          environment-file: environment.yml
+          cache-downloads: true
+          cache-env: true
+      - name: Install & Test polyhedral-gravity
+        shell: bash -l {0}
+        run: |
+          pip install . -vv --no-build-isolation
+          pytest -n 3
+          
+
+  # Builds the polyhedral gravity python interface the interface
+  # No testing on Windows since it takes too long
+  # Enables the long paths feature on Windows (newer thrust versions require this)
+  build-windows:
+    runs-on: windows-latest
+    steps:
+      - uses: actions/checkout@v3
+      - uses: ilammy/msvc-dev-cmd@v1
+      - name: Build
+        run: |
+          New-ItemProperty -Path "HKLM:\SYSTEM\CurrentControlSet\Control\FileSystem" -Name "LongPathsEnabled" -Value 1 -PropertyType DWORD -Force
+          git config --system core.longpaths true
+          mkdir build
+          cd build
+          cmake -DCMAKE_BUILD_TYPE=Release -DBUILD_POLYHEDRAL_PYTHON_INTERFACE=ON ..
+          cmake --build . --config Release
+
+      

.github/workflows/wheels.yml

@@ -7,9 +7,11 @@ on:
     # See https://docs.github.com/en/actions/reference/events-that-trigger-workflows#release
     # for more information on the release event
     prerelease: false
+  # Allows you to run this workflow manually from the Actions tab
+  workflow_dispatch:
 
 jobs:
-  # 1.a Buidl the wheels on a matrix of Windows, MacOS, and Linux plattforms usinbg cibuildwheel
+  # 1.a Buidl the wheels on a matrix of Windows, MacOS, and Linux platforms using cibuildwheel
   build_wheels:
     name: Build wheels on ${{ matrix.os }}
     runs-on: ${{ matrix.os }}
@@ -44,7 +46,7 @@ jobs:
           package-dir: .
           output-dir: dist
         if: matrix.os == 'macos-latest'
-      # Set up the Visual Studio environment on Windows (required, so that CNaje finds the compiler)
+      # Set up the Visual Studio environment on Windows (required, so that CMake finds the compiler)
       - uses: ilammy/msvc-dev-cmd@v1
         if: matrix.os == 'windows-latest'
       - name: Build wheels (Windows)
@@ -93,7 +95,7 @@ jobs:
           repository-url: https://test.pypi.org/legacy/
 
   # 3. Check if the package can be installed from testpypi
-  #    Notice, that this is more of an install test since the
+  #    Notice, that this is more of an installation test since the
   #    import check has already been done in the build_wheels job
   check_testpypi:
     needs: [upload_testpypi]
@@ -126,7 +128,7 @@ jobs:
       - name: Check import
         run: python -c "import polyhedral_gravity"
 
-  # 4. Upload the wheels to the actualy Python Package Index
+  # 4. Upload the wheels to the actually Python Package Index
   #    using trusted publishing
   upload_pypi:
     needs: [build_wheels, make_sdist, check_testpypi]

environment.yml

@@ -0,0 +1,7 @@
+name: polyhedral-gravity-env
+channels:
+  - conda-forge
+dependencies:
+  - numpy==1.23.4
+  - pytest
+  - pytest-xdist

pytest.ini

@@ -0,0 +1,4 @@
+[pytest]
+minversion = 6.0
+testpaths =
+    test

script/time.py

@@ -4,6 +4,7 @@
 import matplotlib.pyplot as plt
 import mesh_utility
 import timeit
+import pickle
 
 vertices, faces = mesh_utility.read_pk_file("/Users/schuhmaj/Programming/polyhedral-gravity-model/script/mesh/Eros.pk")
 vertices, faces = np.array(vertices), np.array(faces)
@@ -67,3 +68,24 @@
 print(f"--> Time taken: {delta:.3f} microseconds per point")
 
 
+########################################################################################################################
+# PICKLE SUPPORT
+########################################################################################################################
+
+
+with open("evaluable.pk", "wb") as f:
+    pickle.dump(evaluable, f, pickle.HIGHEST_PROTOCOL)
+
+with open("evaluable.pk", "rb") as f:
+    evaluable2 = pickle.load(f)
+
+
+start_time = timeit.default_timer()
+evaluable2(computation_points)
+end_time = timeit.default_timer()
+
+delta = (end_time - start_time) / N * 1e6
+# Print the time in milliseconds
+print("######## Pickle ########")
+print(f"--> 1 time {N} points with GravityEvaluable")
+print(f"--> Time taken: {delta:.3f} microseconds per point")

setup.py

@@ -147,7 +147,7 @@ def build_extension(self, ext):
 # --------------------------------------------------------------------------------
 setup(
     name="polyhedral_gravity",
-    version="2.0.3",
+    version="2.1",
     author="Jonas Schuhmacher",
     author_email="[email protected]",
     description="Package to compute full gravity tensor of a given constant density polyhedron for arbitrary points",

src/polyhedralGravity/calculation/GravityEvaluable.cpp

@@ -52,10 +52,10 @@ namespace polyhedralGravity {
 
         if constexpr (Parallelization) {
             result = thrust::transform_reduce(thrust::device, zip1, zip2, &GravityEvaluable::evaluateFace, result,
-                                              util::operator+<double, Array3, Array6>);
+                                              util::operator+ < double, Array3, Array6 > );
         } else {
             result = thrust::transform_reduce(thrust::host, zip1, zip2, &GravityEvaluable::evaluateFace, result,
-                                              util::operator+<double, Array3, Array6>);
+                                              util::operator+ < double, Array3, Array6 > );
         }
 
         SPDLOG_LOGGER_DEBUG(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(),
@@ -203,7 +203,11 @@ namespace polyhedralGravity {
             // indicators for numerical magnitudes appearing during the calculation:
             // planeDistance gets very big when far away, sum2 remains independently very small
             SPDLOG_LOGGER_WARN(PolyhedralGravityLogger::DEFAULT_LOGGER.getLogger(),
-                               "The results of point [{}, {}, {}] may be wrong due to floating point arithmetic");
+                               "While evaluating the plane with coordinates v1 = [{}, {}, {}], v2 = [{}, {}, {}], "
+                               "v3 = [{}, {}, {}] (with computation point re-located at the origin) a "
+                               "significant difference of magnitudes occurred during the evaluation. "
+                               "This may lead to numerically unstable results!", face[0][0], face[0][1], face[0][2],
+                               face[1][0], face[1][1], face[1][2], face[2][0], face[2][1], face[2][2]);
         }
 
         //6. Step: Sum for tensor
@@ -229,10 +233,14 @@ namespace polyhedralGravity {
 
     std::string GravityEvaluable::toString() const {
         std::stringstream ss;
-        ss  << "<polyhedral_gravity.GravityEvaluable, density=" << _density
-            << ", vertices= " << _polyhedron.countVertices()
-            << ", faces= " << _polyhedron.countFaces() << ">";
+        ss << "<polyhedral_gravity.GravityEvaluable, density=" << _density << ", vertices= "
+           << _polyhedron.countVertices() << ", faces= " << _polyhedron.countFaces() << ">";
         return ss.str();
     }
 
+    std::tuple<Polyhedron, double, std::vector<Array3Triplet>, std::vector<Array3>, std::vector<Array3Triplet>>
+    GravityEvaluable::getState() const {
+        return std::make_tuple(_polyhedron, _density, _segmentVectors, _planeUnitNormals, _segmentUnitNormals);
+    }
+
 }

src/polyhedralGravity/calculation/GravityEvaluable.h

@@ -59,16 +59,32 @@ namespace polyhedralGravity {
          * @param polyhedralSource - the vertices & faces of the polyhedron as tuple or the filenames of the files
          * @param density - the constant density in [kg/m^3]
          */
-        GravityEvaluable(
-                const PolyhedralSource &polyhedralSource,
-                double density) : _polyhedron{
+        GravityEvaluable(const PolyhedralSource &polyhedralSource, double density) : _polyhedron{
                 std::holds_alternative<std::tuple<std::vector<Array3>, std::vector<IndexArray3>>>(polyhedralSource)
                 ? Polyhedron{std::get<std::tuple<std::vector<Array3>, std::vector<IndexArray3>>>(polyhedralSource)}
-                : TetgenAdapter(std::get<std::vector<std::string>>(polyhedralSource)).getPolyhedron()
-        }, _density{density} {
+                : TetgenAdapter(std::get<std::vector<std::string>>(polyhedralSource)).getPolyhedron()},
+                                                                                     _density{density} {
             this->prepare();
         }
 
+        /**
+         * Instantiates a GravityEvaluable with a given constant density polyhedron and caches.
+         * This is for restoring a GravityEvaluable from a previous state.
+         * @param polyhedron - the polyhedron
+         * @param density - the constant density in [kg/m^3]
+         * @param segmentVectors - the segment vectors
+         * @param planeUnitNormals - the plane unit normals
+         * @param segmentUnitNormals - the segment unit normals
+         */
+        GravityEvaluable(const Polyhedron &polyhedron, double density, const std::vector<Array3Triplet> &segmentVectors,
+                         const std::vector<Array3> &planeUnitNormals,
+                         const std::vector<Array3Triplet> &segmentUnitNormals) : _polyhedron{polyhedron},
+                                                                                 _density{density},
+                                                                                 _segmentVectors{segmentVectors},
+                                                                                 _planeUnitNormals{planeUnitNormals},
+                                                                                 _segmentUnitNormals{
+                                                                                         segmentUnitNormals} {}
+
         /**
          * Evaluates the polyhedrale gravity model for a given constant density polyhedron at computation
          * point P. Wrapper for evaluate<parallelization>.
@@ -100,6 +116,14 @@ namespace polyhedralGravity {
          */
         std::string toString() const;
 
+        /**
+         * Returns the polyhedron, the density and the internal caches.
+         *
+         * @return tuple of polyhedron, density, segmentVectors, planeUnitNormals and segmentUnitNormals
+         */
+        std::tuple<Polyhedron, double, std::vector<Array3Triplet>, std::vector<Array3>, std::vector<Array3Triplet>>
+        getState() const;
+
     private:
 
         /**
@@ -120,6 +144,7 @@ namespace polyhedralGravity {
         template<bool Parallelization = true>
         GravityModelResult evaluate(const Array3 &computationPoint) const;
 
+
         /**
          * Evaluates the polyhedrale gravity model for a given constant density polyhedron at computation
          * at multiple computation points.
@@ -130,7 +155,6 @@ namespace polyhedralGravity {
         template<bool Parallelization = true>
         std::vector<GravityModelResult> evaluate(const std::vector<Array3> &computationPoints) const;
 
-
         /**
          * Evaluates the polyhedrale gravity model for a given constant density polyhedron at computation a certain face.
          * @param tuple - consisting of face, segmentVectors, planeUnitNormal and segmentUnitNormals

src/polyhedralGravityPython/PolyhedralGravityPython.cpp

@@ -83,7 +83,29 @@ PYBIND11_MODULE(polyhedral_gravity, m) {
                     Returns:
                         Either a tuple of potential, acceleration and second derivatives at the computation points or
                         if multiple computation points are given, the result is a list of tuples
-                )mydelimiter", py::arg("computation_points"), py::arg("parallel") = true);
+                )mydelimiter", py::arg("computation_points"), py::arg("parallel") = true)
+            .def(py::pickle(
+                    [](const GravityEvaluable &evaluable) {
+                        const auto&[polyhedron, density, segmentVectors, planeUnitNormals, segmentUnitNormals] =
+                                evaluable.getState();
+                        return py::make_tuple(polyhedron.getVertices(), polyhedron.getFaces(), density, segmentVectors,
+                                              planeUnitNormals, segmentUnitNormals);
+                    },
+                    [](const py::tuple &tuple) {
+                        constexpr size_t tupleSize = 6;
+                        if (tuple.size() != tupleSize) {
+                            throw std::runtime_error("Invalid state!");
+                        }
+                        Polyhedron polyhedron {
+                            tuple[0].cast<std::vector<Array3>>(), tuple[1].cast<std::vector<IndexArray3>>()
+                        };
+                        GravityEvaluable evaluable{
+                            polyhedron, tuple[2].cast<double>(), tuple[3].cast<std::vector<Array3Triplet>>(),
+                            tuple[4].cast<std::vector<Array3>>(), tuple[5].cast<std::vector<Array3Triplet>>()
+                        };
+                        return evaluable;
+                    }
+            ));
 
     py::module_ utility = m.def_submodule("utility",
                                           "This submodule contains useful utility functions like parsing meshes "