Merge pull request #4 from thomgrand/cleanup

Cleanup branch

.coveragerc

@@ -22,10 +22,6 @@ exclude_lines =
     if 0:
     if __name__ == .__main__.:
 
-    #Ignore numba njit code
-    \@njit
-    \@vectorize
-
 ignore_errors = True
 
 [html]

.coveragerc_cpu

@@ -24,10 +24,6 @@ exclude_lines =
     if cupy_enabled:
     if not cupy_available:
 
-    #Ignore numba njit code
-    \@njit
-    \@vectorize
-
 ignore_errors = True
 
 [html]

README.md

@@ -37,14 +37,14 @@ pip install fim-python #CPU version
 
 # Usage
 
-The main interface to create a solver object to use is [`FIMPY.create_fim_solver`](https://fim-python.readthedocs.io/en/latest/interface.html#fimpy.solver.FIMPY.create_fim_solver)
+The main interface to create a solver object to use is [`create_fim_solver`](https://fim-python.readthedocs.io/en/latest/interface.html#fimpy.solver.create_fim_solver)
 
 ```python
-from fimpy.solver import FIMPY
+from fimpy.solver import create_fim_solver
 
 #Create a FIM solver, by default the GPU solver will be called with the active list
 #Set device='cpu' to run on cpu and use_active_list=False to use Jacobi method
-fim = FIMPY.create_fim_solver(points, elems, D)
+fim = create_fim_solver(points, elems, D)
 ```
 
 Example
@@ -55,7 +55,7 @@ The following code reproduces the [above example](#details)
 ```python
 import numpy as np
 import cupy as cp
-from fimpy.solver import FIMPY
+from fimpy.solver import create_fim_solver
 from scipy.spatial import Delaunay
 import matplotlib.pyplot as plt
 
@@ -76,7 +76,7 @@ x0 = np.array([np.argmin(np.linalg.norm(points, axis=-1), axis=0)])
 x0_vals = np.array([0.])
 
 #Create a FIM solver, by default the GPU solver will be called with the active list
-fim = FIMPY.create_fim_solver(points, elems, D)
+fim = create_fim_solver(points, elems, D)
 phi = fim.comp_fim(x0, x0_vals)
 
 #Plot the data of all points to the given x0 at the center of the domain
@@ -99,7 +99,7 @@ On the CPU, `use_active_list=True` outperforms the Jacobi approach for almost al
 
 # Citation
 
-If you find this work useful in your research, please consider citing the paper in the [Journal of Open Source Software](https://joss.theoj.org/)
+If you find this work useful in your research, please consider citing the [paper](https://doi.org/10.21105/joss.03641) in the [Journal of Open Source Software](https://joss.theoj.org/)
 ```bibtex
 @article{grandits_fast_2021,
   doi = {10.21105/joss.03641},

docs/requirements_docs.txt

@@ -1,5 +1,4 @@
 numpy
-numba
 cython
 sphinx
 pydata_sphinx_theme

fimpy/__init__.py

@@ -1,5 +1,5 @@
 #TODO: Import all symbols and methods here
-from .solver import FIMPY
+from .solver import create_fim_solver
 
-__version__ = "1.2"
+__version__ = "1.2.1"
 __author__ = "Thomas Grandits"

fimpy/cupy_kernels.py

@@ -101,4 +101,4 @@
         }
       }
   }
-}''') #: CUDA kernel to compute a mask of all element permutations containing at least one active index. New, more inefficient version using shared memory.
+}''') #: CUDA kernel to compute a mask of all element permutations containing at least one active index. New, more efficient version using shared memory.

fimpy/fim_base.py

@@ -6,7 +6,6 @@
 from itertools import permutations
 from abc import abstractmethod
 from .utils.tsitsiklis import norm_map
-#import pyximport; pyximport.install()
 from .fim_cutils import compute_point_elem_map_c, compute_neighborhood_map_c
 
 class FIMBase():
@@ -71,7 +70,6 @@ def __init__(self, points, elems, metrics=None, precision=np.float32, comp_conne
             self.point_elem_map = self.compute_point_elem_map()
 
         self.precision = precision
-        self.norm_map = norm_map
         self.dims = self.points.shape[-1]
         self.choose_update_alg()
 
@@ -107,7 +105,7 @@ def check_metrics_argument(self, metrics):
         assert(metrics.shape[0] == self.nr_elems) #One constant metric for each element
         assert(metrics.ndim == 3)
         assert(metrics.shape[-1] == metrics.shape[-2] and metrics.shape[-1] == self.points.shape[-1])
-        assert(np.allclose(metrics - np.transpose(metrics, axes=(0, 2, 1)), 0.)) #Symmetric
+        assert(np.allclose(metrics - np.transpose(metrics, axes=(0, 2, 1)), 0., atol=1e-4)) #Symmetric
         assert(np.all(np.linalg.eigh(metrics)[0] > 1e-4)) #Positive definite
 
 
@@ -186,7 +184,7 @@ def tsitsiklis_update_line(self, x1, x2, D, u1, lib=np):
         ndarray (precision)
             An [..., N] array that holds :math:`||\\mathbf{x}_2 - \\mathbf{x}_1||_M`
         """
-        norm_f = self.norm_map[lib][D.shape[-1]][0]
+        norm_f = norm_map[lib][D.shape[-1]][0]
         a1 = x2 - x1
         return u1 + norm_f(D, a1, a1)
 
@@ -199,7 +197,7 @@ def tsitsiklis_update_point_sol(self, x1, x2, x3, D, u1, u2, lib=np):
         in a broadcasted way.
         For more information on the type and shape of the parameters and return value, see :meth:`tsitsiklis_update_line`.
         """
-        norm_f = self.norm_map[lib][D.shape[-1]][0]
+        norm_f = norm_map[lib][D.shape[-1]][0]
 
         a1 = x3 - x1
         a2 = x3 - x2
@@ -220,7 +218,7 @@ def tsitsiklis_update_triang(self, x1, x2, x3, D, u1, u2, lib=np):
         z2 = x2 - x3
         z1 = x1 - x2
 
-        norm_f, norm_sqr_f = self.norm_map[lib][D.shape[-1]]
+        norm_f, norm_sqr_f = norm_map[lib][D.shape[-1]]
 
         p11 = norm_sqr_f(D, x1=z1, x2=z1)
         p12 = norm_sqr_f(D, x1=z1, x2=z2)
@@ -251,7 +249,7 @@ def tsitsiklis_update_triang(self, x1, x2, x3, D, u1, u2, lib=np):
     def tsitsiklis_update_tetra_quadr(self, D, k, z1, z2, lib=np):
         """Computes the quadratic equation for the tetrahedra update in :meth:`calculate_tet_update`.
         """
-        norm_f, norm_sqr_f = self.norm_map[lib][D.shape[-1]]
+        norm_f, norm_sqr_f = norm_map[lib][D.shape[-1]]
         p11 = norm_sqr_f(D, z1, z1)
         p12 = norm_sqr_f(D, z1, z2)
         p22 = norm_sqr_f(D, z2, z2)
@@ -294,7 +292,7 @@ def calculate_tet_update(self, x1, x2, x3, x4, D, u1, u2, u3, lib=np):
         """
         xs = lib.stack([x1, x2, x3, x4], axis=-1)
         us = lib.stack([u1, u2, u3], axis=-1)
-        norm_f, norm_sqr_f = self.norm_map[lib][D.shape[-1]]
+        norm_f, norm_sqr_f = norm_map[lib][D.shape[-1]]
 
         y3 = x4 - x3
         y1 = x3 - x1
@@ -357,10 +355,7 @@ def calculate_specific_triang_updates(self, elems_perm, xs_perm, D, us, lib=np):
                                                   D, us_perm[..., 0], us_perm[..., 1], lib=lib)
 
         #Now we need to take the minimum result of old and all new
-        #if lib == np:
         lib.minimum.at(us_new, elems_perm[..., -1], us_result)
-        #elif lib == cp:
-        #    cpx.scatter_min(us_new, elems_perm[..., -1], us_result)
 
         return us_new
 
@@ -483,9 +478,10 @@ def comp_fim(self, x0, x0_vals, metrics=None, max_iterations=int(1e10)):
         """
         #Suppress warnings of the computations, since they should be handled internally
         with np.errstate(divide='ignore', invalid='ignore', over='ignore'):
-            assert(metrics is not None or self.metrics is not None) #Metrics were neither in the constructor, nor here given
+            #TODO: Maybe use identity in this case?
+            assert metrics is not None or self.metrics is not None, f"Metrics (D) need to be provided in comp_fim, or at construction in __init__"
             if metrics is not None:
                 self.check_metrics_argument(metrics)
                 metrics = np.linalg.inv(metrics).astype(self.precision) #The inverse metric is used in the FIM algorithm
 
-            return self._comp_fim(x0, x0_vals, metrics)
+            return self._comp_fim(x0, x0_vals, metrics, max_iterations=max_iterations)

fimpy/fim_cupy.py

@@ -139,9 +139,8 @@ def __init__(self, points, elems, metrics=None, precision=np.float32):
     self.streams = [cp.cuda.Stream(non_blocking=False) for i in range(4)]
     self.mempool = cp.get_default_memory_pool()
 
-  def __del__(self):
-    #self.mempool.free_all_blocks()
-    pass
+  def free_gpu_mem(self):
+    self.mempool.free_all_blocks()
 
   def calculate_all_line_updates(self, elems_perm, xs_perm, D, us, lib=np):
     us_new = us.copy() 
@@ -150,7 +149,7 @@ def calculate_all_line_updates(self, elems_perm, xs_perm, D, us, lib=np):
 
     us_result = self.tsitsiklis_update_line(xs_perm[..., 0, :], xs_perm[..., 1, :],
                             D_broadcasted, us_perm[..., 0], lib=cp)
-    cpx.scatter_min(us_new, elems_perm[..., -1], us_result)
+    cp.minimum.at(us_new, elems_perm[..., -1], us_result)
 
     return us_new
 
@@ -162,7 +161,7 @@ def calculate_all_triang_updates(self, elems_perm, xs_perm, D, us, lib=np):
     us_result = self.tsitsiklis_update_triang(xs_perm[..., 0, :], xs_perm[..., 1, :], xs_perm[..., 2, :],
                                               D_broadcasted, us_perm[..., 0], us_perm[..., 1], lib=lib)
 
-    cpx.scatter_min(us_new, elems_perm[..., -1], us_result)
+    cp.minimum.at(us_new, elems_perm[..., -1], us_result)
     return us_new
 
   def tsitsiklis_update_tetra(self, x1, x2, x3, x4, D, u1, u2, u3, lib=np):
@@ -191,7 +190,7 @@ def calculate_all_tetra_updates(self, elems_perm, xs_perm, D, us, lib=np):
 
     us_result = self.tsitsiklis_update_tetra(xs_perm[..., 0, :], xs_perm[..., 1, :], xs_perm[..., 2, :], xs_perm[..., 3, :],
                             D_broadcasted, us_perm[..., 0], us_perm[..., 1], us_perm[..., 2], lib=cp)
-    cpx.scatter_min(us_new, elems_perm[..., -1], us_result)
+    cp.minimum.at(us_new, elems_perm[..., -1], us_result)
 
     return us_new
 
@@ -252,8 +251,6 @@ def __init__(self, points, elems, metrics=None, precision=np.float32):
 
     self.active_list = cp.array(self.active_list)
     self.mempool = cp.get_default_memory_pool()
-    #self.mempool.set_limit(size=self.nr_elems * 4 * 8 * 250)
-    self.norm = self.norm_map[cp][self.dims]
 
     if self.dims == 2:
       self.u3_comp_cupy = u3_comp_cupy_2D
@@ -295,27 +292,20 @@ def comp_unique_map(self, inds, elem_map=True):
         An [?] array holding the indices of the unique element/point indices.
     """
     if elem_map:
-      cpx.scatter_add(self.elem_unique_map, inds, self.elem_ones[inds])
-      #self.elem_unique_map[inds] = 1
+      cp.add.at(self.elem_unique_map, inds, self.elem_ones[inds])
       unique_inds = self.elem_unique_map.nonzero()[0]
-      #unique_mask = self.elem_unique_map != 0
       self.elem_unique_map[:] = 0 #Reset for next run
     else:
-      cpx.scatter_add(self.points_unique_map, inds, self.points_ones[inds])
-      #self.points_unique_map[inds] = 1
+      cp.add.at(self.points_unique_map, inds, self.points_ones[inds])
       unique_inds = self.points_unique_map.nonzero()[0]
-      #unique_mask = self.points_unique_map != 0
       self.points_unique_map[:] = 0 #Reset for next run
 
-    #return unique_mask
     return unique_inds
 
-  def __del__(self):
+  def free_gpu_mem(self):
     self.mempool.free_all_blocks()
-    #pass
 
   def tsitsiklis_update_triang(self, x1, x2, x3, D, u1, u2, lib=np, use_streams=None):
-        #norm_f, norm_sqr_f = self.norm_map[lib][D.shape[-1]]
 
         #Called for non tetra meshes -> Streams are available
         if self.elem_dims == 3:
@@ -349,7 +339,7 @@ def calculate_specific_triang_updates(self, elems_perm, xs_perm, D, us, lib=np):
                                               D, us_perm[..., 0], us_perm[..., 1], lib=lib)
 
     #Now we need to take the minimum result of old and all new
-    cpx.scatter_min(us_new, elems_perm[..., -1], us_result)
+    cp.minimum.at(us_new, elems_perm[..., -1], us_result)
 
     return us_new
 
@@ -358,7 +348,7 @@ def calculate_specific_line_updates(self, elems_perm, xs_perm, D, us, lib=np):
     us_perm = us[elems_perm]
     us_result = self.tsitsiklis_update_line(xs_perm[..., 0, :], xs_perm[..., 1, :],
                             D, us_perm[..., 0], lib=lib)
-    cpx.scatter_min(us_new, elems_perm[..., -1], us_result)
+    cp.minimum.at(us_new, elems_perm[..., -1], us_result)
 
     return us_new
 
@@ -369,11 +359,8 @@ def tsitsiklis_update_tetra(self, x1, x2, x3, x4, D, u1, u2, u3, lib=np):
       u_tet = lib.where(cp.isnan(u_tet), lib.inf, u_tet)
 
     #Face calculations (Includes possible line calculations)
-    #with self.streams[1]:
     triang1 = self.tsitsiklis_update_triang(x1, x2, x4, D, u1, u2, cp, (self.streams[1], self.streams[2]))
-    #with self.streams[2]:
     triang2 = self.tsitsiklis_update_triang(x1, x3, x4, D, u1, u3, cp, (self.streams[3], self.streams[4]))
-    #with self.streams[3]:
     triang3 = self.tsitsiklis_update_triang(x2, x3, x4, D, u2, u3, cp, (self.streams[5], self.streams[6]))
 
     u_triang = lib.minimum(triang1, triang2)
@@ -387,16 +374,12 @@ def tsitsiklis_update_tetra_quadr(self, D, k, z1, z2, lib=np):
     else:
       return super().tsitsiklis_update_tetra_quadr(D, k, z1, z2, cp)
 
-  #def tsitsiklis_update_tetra(self, x1, x2, x3, x4, D, u1, u2, u3, lib=np):
-
-
-
   def calculate_specific_tetra_updates(self, elems_perm, xs_perm, D, us, lib=np):
     us_new = us.copy() 
     us_perm = us[elems_perm]
     us_result = self.tsitsiklis_update_tetra(xs_perm[..., 0, :], xs_perm[..., 1, :], xs_perm[..., 2, :], xs_perm[..., 3, :],
                             D, us_perm[..., 0], us_perm[..., 1], us_perm[..., 2], lib=cp)
-    cpx.scatter_min(us_new, elems_perm[..., -1], us_result)
+    cp.minimum.at(us_new, elems_perm[..., -1], us_result)
 
     return us_new