Merge pull request #77 from tberlok/actions

Actions

.gitignore

@@ -13,3 +13,4 @@ dist/
 *.rst
 *.png
 *.mp4
+*.log

paicos/__init__.py

@@ -44,6 +44,7 @@
 from .image_creators.nested_projector import NestedProjector
 from .image_creators.tree_projector import TreeProjector
 from .image_creators.slicer import Slicer
+from .image_creators.image_creator_actions import Actions
 
 # Histograms
 from .histograms.histogram import Histogram

paicos/cython/sph_projectors.pyx

@@ -51,9 +51,10 @@ def project_image(real_t[:] xvec, real_t[:] yvec, real_t[:] variable,
     cdef int Np = xvec.shape[0]
 
     # Shape of projection array
-    cdef int ny = <int> (sidelength_y/sidelength_x * nx)
+    cdef int ny = <int> round(sidelength_y/sidelength_x * nx)
     msg = '(sidelength_y/sidelength_x * nx) needs to be an integer'
-    assert (sidelength_y/sidelength_x * nx) == <float> ny, msg
+
+    assert abs((sidelength_y/sidelength_x * nx) - <float> ny) < 1e-8, msg
 
     # Create projection array
     cdef real_t[:, :] projection = np.zeros((nx, ny), dtype=np.float64)
@@ -67,7 +68,7 @@ def project_image(real_t[:] xvec, real_t[:] yvec, real_t[:] variable,
     cdef real_t x, y, norm
     cdef real_t x0, y0
 
-    assert sidelength_x/nx == sidelength_y/ny
+    # assert sidelength_x/nx == sidelength_y/ny
 
     # Lower left corner of image in Arepo coordinates
     x0 = xc - sidelength_x / 2.0
@@ -159,9 +160,10 @@ def project_image_omp(real_t[:] xvec, real_t[:] yvec, real_t[:] variable,
     cdef int Np = xvec.shape[0]
 
     # Shape of projection array
-    cdef int ny = <int> (sidelength_y/sidelength_x * nx)
+    cdef int ny = <int> round(sidelength_y/sidelength_x * nx)
     msg = '(sidelength_y/sidelength_x * nx) needs to be an integer'
-    assert (sidelength_y/sidelength_x * nx) == <float> ny, msg
+
+    assert abs((sidelength_y/sidelength_x * nx) - <float> ny) < 1e-8, msg
 
     # Loop integers and other variables
     cdef int ip, ix, iy, ih, threadnum
@@ -172,7 +174,7 @@ def project_image_omp(real_t[:] xvec, real_t[:] yvec, real_t[:] variable,
     cdef real_t x, y, norm
     cdef real_t x0, y0
 
-    assert sidelength_x/nx == sidelength_y/ny
+    # assert sidelength_x/nx == sidelength_y/ny
 
     # Lower left corner of image in Arepo coordinates
     x0 = xc - sidelength_x/2.0
@@ -282,9 +284,10 @@ def project_oriented_image(real_t[:] xvec, real_t[:] yvec, real_t[:] zvec, real_
     cdef int Np = xvec.shape[0]
 
     # Shape of projection array
-    cdef int ny = <int> (sidelength_y/sidelength_x * nx)
+    cdef int ny = <int> round(sidelength_y/sidelength_x * nx)
     msg = '(sidelength_y/sidelength_x * nx) needs to be an integer'
-    assert (sidelength_y/sidelength_x * nx) == <float> ny, msg
+
+    assert abs((sidelength_y/sidelength_x * nx) - <float> ny) < 1e-8, msg
 
     # Create projection array
     cdef real_t[:, :] projection = np.zeros((nx, ny), dtype=np.float64)
@@ -298,7 +301,7 @@ def project_oriented_image(real_t[:] xvec, real_t[:] yvec, real_t[:] zvec, real_
     cdef real_t x, y, norm
     cdef real_t cen_x, cen_y, cen_z
 
-    assert sidelength_x/nx == sidelength_y/ny
+    # assert sidelength_x/nx == sidelength_y/ny
 
     for ip in range(Np):
         # Centered coordinates
@@ -396,9 +399,10 @@ def project_oriented_image_omp(real_t[:] xvec, real_t[:] yvec, real_t[:] zvec, r
     cdef int Np = xvec.shape[0]
 
     # Shape of projection array
-    cdef int ny = <int> (sidelength_y/sidelength_x * nx)
+    cdef int ny = <int> round(sidelength_y/sidelength_x * nx)
     msg = '(sidelength_y/sidelength_x * nx) needs to be an integer'
-    assert (sidelength_y/sidelength_x * nx) == <float> ny, msg
+
+    assert abs((sidelength_y/sidelength_x * nx) - <float> ny) < 1e-8, msg
 
     # Loop integers and other variables
     cdef int ip, ix, iy, ih, threadnum
@@ -409,7 +413,7 @@ def project_oriented_image_omp(real_t[:] xvec, real_t[:] yvec, real_t[:] zvec, r
     cdef real_t x, y, norm
     cdef real_t cen_x, cen_y, cen_z
 
-    assert sidelength_x/nx == sidelength_y/ny
+    #assert sidelength_x/nx == sidelength_y/ny
 
     # Create projection arrays
     cdef real_t[:, :, :] tmp_var = np.zeros((nx, ny, numthreads), dtype=np.float64)

paicos/image_creators/gpu_ray_projector.py

@@ -257,17 +257,22 @@ def _send_small_data_to_gpu(self):
             self.gpu_variables['widths'] = cp.array(self.widths)
             self.gpu_variables['center'] = cp.array(self.center)
 
-    def _gpu_project(self, variable_str):
+    def _gpu_project(self, variable_str, additive):
         """
         Private method for projecting using cuda code
         """
-        rotation_matrix = self.gpu_variables['rotation_matrix']
-        widths = self.gpu_variables['widths']
-        center = self.gpu_variables['center']
-        hsml = self.gpu_variables['hsml']
+        gpu_vars = self.gpu_variables
+        rotation_matrix = gpu_vars['rotation_matrix']
+        widths = gpu_vars['widths']
+        center = gpu_vars['center']
+        hsml = gpu_vars['hsml']
         nx = self.npix_width
         ny = self.npix_height
-        variable = self.gpu_variables[variable_str]
+
+        if additive:
+            variable = gpu_vars[variable_str] / gpu_vars[f'{self.parttype}_Volume']
+        else:
+            variable = gpu_vars[variable_str]
 
         tree_scale_factor = self.tree.conversion_factor
         tree_offsets = self.tree.off_sets
@@ -286,29 +291,7 @@ def _gpu_project(self, variable_str):
         projection = cp.asnumpy(image)
         return projection
 
-    def project_variable(self, variable, additive=False):
-        """
-        projects a given variable onto a 2D plane.
-
-        Parameters
-        ----------
-        variable : str, function, numpy array
-            variable, it can be passed as string or an array
-
-        Returns
-        -------
-        numpy array
-            The image of the projected variable
-        """
-
-        # This calls _do_region_selection if resolution, Orientation,
-        # widths or center changed
-        self._check_if_properties_changed()
-
-        if additive:
-            err_msg = "GPU ray tracer does not yet support additive=True"
-            raise RuntimeError(err_msg)
-
+    def _send_variable_to_gpu(self, variable, gpu_key='projection_variable'):
         if isinstance(variable, str):
             variable_str = str(variable)
             err_msg = 'projector uses a different parttype'
@@ -323,9 +306,10 @@ def project_variable(self, variable, additive=False):
 
         # Select same part of array that the projector has selected
         if self.do_pre_selection:
+            # TODO: Check that this is not applied more than once...
             variable = variable[self.index]
 
-        if variable_str in self.gpu_variables and variable_str != 'projection_variable':
+        if variable_str in self.gpu_variables and variable_str != gpu_key:
             pass
         else:
             # Send variable to gpu
@@ -338,25 +322,69 @@ def project_variable(self, variable, additive=False):
             self.gpu_variables[variable_str] = self.gpu_variables[variable_str][
                 self.tree.sort_index]
 
+        if isinstance(variable, units.PaicosQuantity):
+            unit_quantity = variable.unit_quantity
+        else:
+            unit_quantity = None
+
+        return variable_str, unit_quantity
+
+    def project_variable(self, variable, additive=False):
+        """
+        projects a given variable onto a 2D plane.
+
+        Parameters
+        ----------
+        variable : str, function, numpy array
+            variable, it can be passed as string or an array
+
+        Returns
+        -------
+        numpy array
+            The image of the projected variable
+        """
+
+        # This calls _do_region_selection if resolution, Orientation,
+        # widths or center changed
+        self._check_if_properties_changed()
+
+        variable_str, unit_quantity = self._send_variable_to_gpu(variable)
+        if additive:
+            _, vol_unit_quantity = self._send_variable_to_gpu(f'{self.parttype}_Volume')
+
         # Do the projection
-        projection = self._gpu_project(variable_str)
+        projection = self._gpu_project(variable_str, additive)
 
         # Transpose
         projection = projection.T
 
         assert projection.shape[0] == self.npix_height
         assert projection.shape[1] == self.npix_width
 
-        if isinstance(variable, units.PaicosQuantity):
+        if unit_quantity is not None:
             unit_length = self.snap['0_Coordinates'].uq
-            projection = projection * variable.unit_quantity * unit_length
+            projection = projection * unit_quantity * unit_length
+            if additive:
+                projection = projection / vol_unit_quantity
 
-        return projection / self.depth
+        if additive:
+            return projection
+        else:
+            return projection / self.depth
 
     def __del__(self):
         """
         Clean up like this? Not sure it is needed...
         """
-        del self.gpu_variables
-        del self.tree
+        self.release_gpu_memory()
+
+    def release_gpu_memory(self):
+        if hasattr(self, 'gpu_variables'):
+            for key in list(self.gpu_variables):
+                del self.gpu_variables[key]
+            del self.gpu_variables
+        if hasattr(self, 'tree'):
+            self.tree.release_gpu_memory()
+            del self.tree
+
         cp._default_memory_pool.free_all_blocks()

paicos/image_creators/gpu_sph_projector.py

@@ -439,7 +439,12 @@ def __del__(self):
         """
         Clean up like this? Not sure it is needed...
         """
-        del self.gpu_variables
-        # for key in list(self.gpu_variables.keys):
-        #     del self.gpu_variables[key]
+        self.release_gpu_memory()
+
+    def release_gpu_memory(self):
+        if hasattr(self, 'gpu_variables'):
+            for key in list(self.gpu_variables):
+                del self.gpu_variables[key]
+            del self.gpu_variables
+
         cp._default_memory_pool.free_all_blocks()

paicos/image_creators/image_creator.py

@@ -190,6 +190,14 @@ def center(self, value):
             self._center = np.array(value)
         self._properties_changed = True
 
+    def reset_center(self):
+        """
+        Resets the center of the image creator to the
+        value it had at initialization.
+        """
+        self.center = self.center - self._diff_center
+        self._diff_center = np.zeros_like(self._diff_center)
+
     @property
     def npix(self):
         """
@@ -365,41 +373,41 @@ def width(self, value):
         if settings.use_units:
             assert hasattr(value, 'unit')
         if self.direction == 'x':
-            self._widths[1] = value.copy
+            self._widths[1] = np.copy(value)
 
         elif self.direction == 'y':
-            self._widths[2] = value.copy
+            self._widths[2] = np.copy(value)
 
         elif self.direction == 'z' or self.direction == 'orientation':
-            self._widths[0] = value.copy
+            self._widths[0] = np.copy(value)
         self._properties_changed = True
 
     @height.setter
     def height(self, value):
         if settings.use_units:
             assert hasattr(value, 'unit')
         if self.direction == 'x':
-            self._widths[2] = value.copy
+            self._widths[2] = np.copy(value)
 
         elif self.direction == 'y':
-            self._widths[0] = value.copy
+            self._widths[0] = np.copy(value)
 
         elif self.direction == 'z' or self.direction == 'orientation':
-            self._widths[1] = value.copy
+            self._widths[1] = np.copy(value)
         self._properties_changed = True
 
     @depth.setter
     def depth(self, value):
         if settings.use_units:
             assert hasattr(value, 'unit')
         if self.direction == 'x':
-            self._widths[0] = value.copy
+            self._widths[0] = np.copy(value)
 
         elif self.direction == 'y':
-            self._widths[1] = value.copy
+            self._widths[1] = np.copy(value)
 
         elif self.direction == 'z' or self.direction == 'orientation':
-            self._widths[2] = value.copy
+            self._widths[2] = np.copy(value)
         self._properties_changed = True
 
     def double_resolution(self):
@@ -485,9 +493,9 @@ def npix_height(self):
         The number of pixels of an image in the vertical direction.
         """
         if settings.use_units:
-            return int((self.height / self.width).value * self.npix_width)
+            return int(round((self.height / self.width).value * self.npix_width))
         else:
-            return int((self.height / self.width) * self.npix_width)
+            return int(round((self.height / self.width) * self.npix_width))
 
     @property
     def area(self):