Merge branch 'main' into issue/714

Conflicts:
	src/kbmod/search/stamp_creator.cpp

docs/source/user_manual/index.rst

@@ -69,7 +69,20 @@ You can run then run the tests to check that everything works:
    cd tests
    python -m unittest
 
-If tests fail or more than a few tests are skipped then it is possible that cmake was unable to find your GPU when compiling the code. Try running `cmake3 -B src/kbmod -S .` from the KBMOD directory. This will parse the `CMakeLists.txt` and produce more verbose output.
+If tests fail or more than a few tests are skipped then it is possible that cmake was unable to find your GPU when compiling the code. Depending on environment, you may need to install the CUDA libraries.
+
+.. code-block:: bash
+		
+   conda install -c conda-forge cudatoolkit-dev
+
+If you are seeing a warning saying "An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu." then you need to install the JAX cuda libraries ( `installation instructions here <https://jax.readthedocs.io/en/latest/installation.html>`_ ). For example you can install the libraries for an NVIDIA GPU on Linux with pip as:
+
+.. code-block:: bash
+		
+   pip install --upgrade "jax[cuda12]"
+
+If you still run into problems finding the GPU, try running `cmake3 -B src/kbmod -S .` from the KBMOD directory. This will parse the `CMakeLists.txt` and produce more verbose output.
+
 
 
 Running KBMOD

src/kbmod/fake_data/fake_data_creator.py

@@ -84,9 +84,9 @@ def add_fake_object(img, x, y, flux, psf=None):
     ----------
     img : `RawImage` or `LayeredImage`
         The image to modify.
-    x : `float`
+    x : `int` or `float`
         The x pixel location of the fake object.
-    y : `float`
+    y : `int` or `float`
         The y pixel location of the fake object.
     flux : `float`
         The flux value.
@@ -98,15 +98,24 @@ def add_fake_object(img, x, y, flux, psf=None):
     else:
         sci = img
 
+    # Explicitly cast to float because the indexing uses different order
+    # float integer and float.
+    x = int(x)
+    y = int(y)
+
     if psf is None:
-        sci.interpolated_add(x, y, flux)
+        if (x >= 0) and (y >= 0) and (x < sci.width) and (y < sci.height):
+            sci.set_pixel(y, x, flux + sci.get_pixel(y, x))
     else:
         dim = psf.get_dim()
-        initial_x = x - psf.get_radius()
-        initial_y = y - psf.get_radius()
+        initial_x = int(x - psf.get_radius())
+        initial_y = int(y - psf.get_radius())
         for i in range(dim):
             for j in range(dim):
-                sci.interpolated_add(float(initial_x + i), float(initial_y + j), flux * psf.get_value(i, j))
+                xp = initial_x + i
+                yp = initial_y + j
+                if (xp >= 0) and (yp >= 0) and (xp < sci.width) and (yp < sci.height):
+                    sci.set_pixel(yp, xp, flux * psf.get_value(i, j) + sci.get_pixel(yp, xp))
 
 
 def create_fake_times(num_times, t0=0.0, obs_per_day=1, intra_night_gap=0.01, inter_night_gap=1):

src/kbmod/fake_data/insert_fake_orbit.py

@@ -1,5 +1,6 @@
 from astropy.wcs import WCS
 
+from kbmod.fake_data.fake_data_creator import add_fake_object
 from kbmod.fake_data.pyoorb_helper import PyoorbOrbit
 from kbmod.search import ImageStack, LayeredImage, PSF, RawImage
 from kbmod.work_unit import WorkUnit
@@ -10,7 +11,7 @@ def safe_add_fake_detection(img, x, y, flux):
 
     Parameters
     ----------
-    img : `RawImage` or `LayeredImage`
+    img : `LayeredImage`
         The image to modify.
     x : `float`
         The x pixel location of the fake object.
@@ -35,16 +36,7 @@ def safe_add_fake_detection(img, x, y, flux):
     if img.get_mask().get_pixel(int(y), int(x)) != 0:
         return False
 
-    sci = img.get_science()
-    psf = img.get_psf()
-    dim = psf.get_dim()
-
-    initial_x = x - psf.get_radius()
-    initial_y = y - psf.get_radius()
-    for i in range(dim):
-        for j in range(dim):
-            sci.interpolated_add(float(initial_x + i), float(initial_y + j), flux * psf.get_value(i, j))
-
+    add_fake_object(img.get_science(), x, y, flux, img.get_psf())
     return True
 
 

src/kbmod/results.py

@@ -2,6 +2,7 @@
 and helper functions for filtering and maintaining consistency between different attributes in each row.
 """
 
+import copy
 import csv
 import logging
 import numpy as np
@@ -132,6 +133,10 @@ def get_num_times(self):
             return self.table["obs_valid"].shape[1]
         return 0
 
+    def copy(self):
+        """Return a deep copy of the current Results object."""
+        return copy.deepcopy(self)
+
     @classmethod
     def from_trajectories(cls, trajectories, track_filtered=False):
         """Extract data from a list of Trajectory objects.

src/kbmod/search/geom.h

@@ -157,56 +157,6 @@ inline Rectangle anchored_block(const Index& idx, const int r, const unsigned wi
     return {{top, left}, {anchor_top, anchor_left}, (unsigned)rangej, (unsigned)rangei};
 }
 
-// returns top-right-bot-left (clockwise) corners around an Index.
-inline auto manhattan_neighbors(const Index& idx, const unsigned width, const unsigned height) {
-    std::array<std::optional<Index>, 4> idxs;
-
-    // if conditions are not met, the element is not set. Because of optional
-    // type, casting a not-set element to bool returns false. Pybind11
-    // evaluates it as a None. See `interpolate` in RawImage for an example.
-    // top bot
-    if (idx.j >= 0 && idx.j < width) {
-        if (idx.i - 1 >= 0 && idx.i - 1 < height) idxs[0] = {idx.i - 1, idx.j};
-        if (idx.i + 1 >= 0 && idx.i + 1 < height) idxs[2] = {idx.i + 1, idx.j};
-    }
-
-    // right left
-    if (idx.i >= 0 && idx.i < height) {
-        if (idx.j + 1 >= 0 && idx.j + 1 < width) idxs[1] = {idx.i, idx.j + 1};
-        if (idx.j - 1 >= 0 && idx.j - 1 < width) idxs[3] = {idx.i, idx.j - 1};
-    }
-    return idxs;
-}
-
-// Note the distinct contextual distinction between manhattan neighborhood of
-// Index and Point. Point returns closes **pixel indices** that are neighbors.
-// This includes the pixel the Point is residing within. This is not the case
-// for Index, which will never return itself as a neighbor.
-inline auto manhattan_neighbors(const Point& p, const unsigned width, const unsigned height) {
-    std::array<std::optional<Index>, 4> idxs;
-
-    // The index in which the point resides.
-    // Almost always top-left corner, except when
-    // point is negative or (0, 0)
-    auto idx = p.to_index();
-
-    // top-left bot-left
-    if (idx.j >= 0 && idx.j < width) {
-        if (idx.i >= 0 && idx.i < height) idxs[0] = {idx.i, idx.j};
-        if (idx.i + 1 >= 0 && idx.i + 1 < height) idxs[3] = {idx.i + 1, idx.j};
-    }
-
-    // top-right
-    if (idx.i >= 0 && idx.i < height)
-        if (idx.j + 1 >= 0 && idx.j + 1 < width) idxs[1] = {idx.i, idx.j + 1};
-
-    // bot-right
-    if ((idx.i + 1 >= 0 && idx.i + 1 < width) && (idx.j + 1 >= 0 && idx.j + 1 < width))
-        idxs[2] = {idx.i + 1, idx.j + 1};
-
-    return idxs;
-}
-
 #ifdef Py_PYTHON_H
 static void index_bindings(py::module& m) {
     PYBIND11_NUMPY_DTYPE(Index, i, j);
@@ -312,23 +262,6 @@ static void geom_functions(py::module& m) {
                 return anchored_block({idx.first, idx.second}, r, shape.second, shape.first);
             },
             pydocs::DOC_anchored_block);
-
-    // Safe to cast to int as "ij" implies user is using indices, i.e. ints.
-    // CPP can be so odd, why not return a 1 or, you know... a bool? Mostly for
-    // testing purposes.
-    m.def(
-            "manhattan_neighbors",
-            [](const std::pair<float, float> coords, const std::pair<unsigned, unsigned> shape,
-               const std::string indexing) {
-                if (indexing.compare("ij") == 0)
-                    return manhattan_neighbors(Index{(int)coords.first, (int)coords.second}, shape.second,
-                                               shape.first);
-                else if (indexing.compare("xy") == 0)
-                    return manhattan_neighbors(Point{coords.first, coords.second}, shape.second, shape.first);
-                else
-                    throw std::domain_error("Expected 'ij' or 'xy' got " + indexing + " instead.");
-            },
-            pydocs::DOC_manhattan_neighbors);
 }
 #endif  // Py_PYTHON_H
 }  // namespace indexing

src/kbmod/search/pydocs/geom_docs.h

@@ -154,49 +154,6 @@ static const auto DOC_anchored_block = R"doc(
          [-1., 10., 11.]])
   )doc";
 
-static const auto DOC_manhattan_neighbors = R"doc(
-  Returns a list of nearest neighbors as determined by Manhattan distance of 1.
-
-  Indexing scheme ``ij`` handles the input as `Index`, i.e. the closest
-  neighbors are top, right, bot, and left indices when not on the edge of an
-  array. When ``xy``, the input is treated as a `Point`, i.e. real Cartesian
-  coordinates. In this case the closest neighbors are the closest pixels. Pixels
-  are array elements understood to span the whole integer range and which center
-  coordinates lie on a half-integer grid.
-  For example, origin of an image is the pixel with index ``(0, 0)``, it spans
-  the range ``[0, 1]`` and its center is the point ``(0.5, 0.5)``. So the
-  closest neighbor to a `Point(0.6, 0.6)` is `Index(0, 0)` and `Point(1, 1)` is
-  equidistant from indices ``[(0, 0), (0, 1), (1, 0), (1, 1)]``.
-
-  Parameters
-  ----------
-  coords : `tuple`
-      Center, around which the neighbors will be returned.
-  shape : `tuple`
-      Dimensions of the image/array.
-  indexing : `str`
-      Indexing scheme, ``ij`` or ``xy``.
-
-  Returns
-  -------
-  neighbors : `list[Index | None]`
-      List of indices in clockwise order starting from top or top-left (as
-      appropriate), or `None` when the returned `Index` would lie outside of the
-      array/
-
-  Raises
-  ------
-  ValueError - when indexing is not ``ij`` or ``xy``
-
-  Examples
-  --------
-  >>> shape = (10, 10)
-  >>> manhattan_neighbors((0, 0), shape, "ij")
-  [None,   (0, 1), (1, 0), None]
-  >>> manhattan_neighbors((1, 1), shape, "xy")
-  [(0, 0), (0, 1), (1, 1), (1, 0)]
-  )doc";
-
 }  // namespace pydocs
 
 #endif  // GEOM_DOCS

src/kbmod/search/pydocs/raw_image_docs.h

@@ -46,9 +46,9 @@ static const auto DOC_RawImage = R"doc(
   representing indices to a 2D matrix, usually expressed with the ``(i, j)``
   convention. Pixel accessing or setting methods of `RawImage`, such as
   `get_pixel`, use the indexing convention. This matches NumPy convention. Other
-  methods, such as `interpolate` or `add_fake_object`, however, use the `(x, y)`
-  convention which is the reversed NumPy convention. Refer to individual methods
-  signature and docstring to see which one they use.
+  methods, such as `add_fake_object`, however, use the `(x, y)` convention which
+  is the reversed NumPy convention. Refer to individual methods signature and docstring
+  to see which one they use.
 
   Examples
   --------
@@ -198,7 +198,7 @@ static const auto DOC_RawImage_compute_bounds = R"doc(
   bounds : `tuple`
       A ``(min, max)`` tuple.
   )doc";
-  
+
 static const auto DOC_RawImage_find_peak = R"doc(
   Returns the pixel coordinates of the maximum value.
 
@@ -264,52 +264,6 @@ static const auto DOC_RawImage_create_stamp = R"doc(
       The stamp.
   )doc";
 
-static const auto DOC_RawImage_interpolate = R"doc(
-  Get the interoplated value of a point.
-
-  Parameters
-  ----------
-  x : `float`
-      The x-coordinate, the abscissa.
-  y : `float`
-      The y-coordinate, the ordinate.
-
-  Returns
-  -------
-  value : `float`
-      Bilinearly interpolated value at that point.
-
-  )doc";
-
-static const auto DOC_RawImage_interpolated_add = R"doc(
-  Add the given value at a given point, to the image.
-
-  Addition is performed by determining the nearest Manhattan neighbors, weighing
-  the given value by the distance to these neighbors and then adding that value
-  at the index locations of the neighbors. Sort of like an inverse bilinear
-  interpolation.
-
-  Parameters
-  ----------
-  x : `float`
-      The x coordinate at which to add value.
-  y : `float`
-      Y coordinate.
-  value : `float`
-      Value to add.
-  )doc";
-
-static const auto DOC_RawImage_get_interp_neighbors_and_weights = R"doc(
-  Returns a tuple of Manhattan neighbors and interpolation weights.
-
-  Parameters
-  ----------
-  x : `float`
-      The x coordinate at which to add value.
-  y : `float`
-      Y coordinate.
-  )doc";
-
 static const auto DOC_RawImage_apply_mask = R"doc(
   Applies a mask to the RawImage by comparing the given bit vector with the
   values in the mask layer and marking pixels ``NO_DATA``.

src/kbmod/search/pydocs/stack_search_docs.h

@@ -169,7 +169,7 @@ static const auto DOC_StackSearch_get_number_total_results = R"doc(
   result : `int`
       The number of cached results.
   )doc";
-                                        
+
 static const auto DOC_StackSearch_get_results = R"doc(
   Get a batch of cached results.
 

src/kbmod/search/pydocs/trajectory_list_docs.h

@@ -143,7 +143,6 @@ static const auto DOC_TrajectoryList_sort_by_likelihood = R"doc(
   Raises a ``RuntimeError`` the data is on GPU.
   )doc";
 
-
 }  // namespace pydocs
 
 #endif /* TRAJECTORY_LIST_DOCS_ */