Added config file linkingfilter.ssp_number_observations into PPMiniDifi

docs/notebooks/demo_miniDifiValidation.ipynb

@@ -12,7 +12,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "id": "878bd071-d694-4bdf-91d2-16e29e090df8",
    "metadata": {},
    "outputs": [],
@@ -83,7 +83,7 @@
     "    tshifted = df[\"midPointTai\"] - night_start_utc_days\n",
     "    return tshifted.astype(int)\n",
     "\n",
-    "def has_tracklet(df, minlen_arcsec=1, maxdt_minutes=90):\n",
+    "def has_tracklet(df, minlen_arcsec=1, maxdt_minutes=90, min_observations=2):\n",
     "    # this computes whether a given night has /any/ observations close enough to be\n",
     "    # made into a tracklet.\n",
     "    #\n",
@@ -94,7 +94,12 @@
     "\n",
     "    maxdt = maxdt_minutes / 60. / 24.\n",
     "    minlen = minlen_arcsec / 3600.\n",
+    "    detections = []\n",
     "\n",
+    "    if len(df)< min_observations:\n",
+    "        return False\n",
+    "    if min_observations == 1:  # for edge case where one observation is needed for a tracklet\n",
+    "        return True\n",
     "    for i in range(len(df)):\n",
     "        for j in range(i+1, len(df)):\n",
     "            t0, t1 = df[\"midPointTai\"].iloc[i], df[\"midPointTai\"].iloc[j]\n",
@@ -118,8 +123,11 @@
     "\n",
     "            # this combo will make a tracklet; return true!\n",
     "            ##print(f'dt = {dt*24*60:.2f} min, sep = {dist*3600:.2f}\"')\n",
-    "            return True\n",
-    "    \n",
+    "            detections.append(i)\n",
+    "            detections.append(j)\n",
+    "            if len(set(detections)) >= min_observations:\n",
+    "                return True\n",
+    "\n",
     "    return False\n",
     "\n",
     "def nights_with_tracklets(df):\n",
@@ -129,7 +137,7 @@
     "        assert len(df[\"ssObjectId\"].unique()) == 1, \"This function must be run on a single object\"\n",
     "\n",
     "    # Filter the DataFrame to only keep rows where the 'night' value appears 2 or more times\n",
-    "    df_filtered = df[df.groupby('night')['night'].transform('count') >= 2].reset_index()\n",
+    "    df_filtered = df[df.groupby('night')['night'].transform('count') >= 1].reset_index()\n",
     "\n",
     "    # Now, for each night compute if it has a tracklet by calling has_tracklet()\n",
     "    nightHasTracklet = df_filtered.groupby(\"night\").apply(has_tracklet, include_groups=False)\n",
@@ -223,8 +231,8 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "CPU times: user 20.1 s, sys: 48.9 ms, total: 20.2 s\n",
-      "Wall time: 20.3 s\n"
+      "CPU times: user 11.4 s, sys: 36.9 ms, total: 11.5 s\n",
+      "Wall time: 11.5 s\n"
      ]
     },
     {
@@ -276,8 +284,8 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "CPU times: user 4.81 s, sys: 358 ms, total: 5.16 s\n",
-      "Wall time: 10.5 s\n"
+      "CPU times: user 1.1 s, sys: 190 ms, total: 1.29 s\n",
+      "Wall time: 1.47 s\n"
      ]
     }
    ],
@@ -297,7 +305,7 @@
     "    )\n",
     "\n",
     "# run minidifi\n",
-    "minidifi = md.linkObservations(obsv, seed=42, p=1)"
+    "minidifi = md.linkObservations(obsv, seed=42, p=1, objectId=\"ssObjectId\", maxdt_minutes=90, minlen_arcsec=1, min_observations=2, window=14, nlink=3, night_start_utc_days=17.0 / 24.0)"
    ]
   },
   {
@@ -390,14 +398,6 @@
     "print(\"Differences between miniDifi and brute force (note: empty is good!):\")\n",
     "res[~res[\"equal\"]]"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "4f785000",
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {

src/sorcha/modules/PPLinkingFilter.py

@@ -76,6 +76,7 @@ def PPLinkingFilter(
         objectId="ssObjectId",
         maxdt_minutes=maximum_time * 24 * 60,
         minlen_arcsec=minimum_separation,
+        min_observations=min_observations,
         window=tracklet_interval,
         nlink=min_tracklets,
         p=detection_efficiency,

src/sorcha/modules/PPMiniDifi.py

@@ -45,7 +45,7 @@ def haversine_np(lon1, lat1, lon2, lat2):
 
 
 @njit(cache=True)
-def hasTracklet(mjd, ra, dec, maxdt_minutes, minlen_arcsec):
+def hasTracklet(mjd, ra, dec, maxdt_minutes, minlen_arcsec, min_observations):
     """
     Given a set of observations in one night, calculate it has
     at least onedetectable tracklet.
@@ -67,6 +67,9 @@ def hasTracklet(mjd, ra, dec, maxdt_minutes, minlen_arcsec):
     minlen_arcsec : float
         Minimum allowable distance separation between observations [Units: arcsec]
 
+    min_observations (int):
+        the minimum number of observations in a night required to form a tracklet.
+
     Returns
     --------
         : boolean
@@ -78,25 +81,26 @@ def hasTracklet(mjd, ra, dec, maxdt_minutes, minlen_arcsec):
     ## tracklets by taking all observations in a night and computing
     ## all of theirs pairwise distances, then selecting on that.
     nobs = len(ra)
-    if nobs < 2:
+    if nobs < min_observations:
         return False
-
+    if min_observations == 1:  # for edge case where one observation is needed for a tracklet
+        return True
     maxdt = maxdt_minutes / (60 * 24)
     minlen = minlen_arcsec / 3600
-
     for i in range(nobs):
         for j in range(nobs):
             diff = mjd[i] - mjd[j]
             if diff > 0 and diff < maxdt:
                 sep = haversine_np(ra[i], dec[i], ra[j], dec[j])
                 if sep > minlen:
-                    return True
+                    if (i + 1) - (j) >= min_observations:
+                        return True
 
     return False
 
 
 @njit(cache=True)
-def trackletsInNights(night, mjd, ra, dec, maxdt_minutes, minlen_arcsec):
+def trackletsInNights(night, mjd, ra, dec, maxdt_minutes, minlen_arcsec, min_observations):
     """
     Calculate, for a given set of observations sorted by observation time,
     whether or not it has at least one discoverable tracklet in each night.
@@ -121,6 +125,9 @@ def trackletsInNights(night, mjd, ra, dec, maxdt_minutes, minlen_arcsec):
     minlen_arcsec : float
         Minimum allowable distance separation between observations [Units: arcsec]
 
+    min_observations (int):
+        the minimum number of observations in a night required to form a tracklet.
+
     Returns
     --------
     nights : float or array of floats
@@ -139,7 +146,7 @@ def trackletsInNights(night, mjd, ra, dec, maxdt_minutes, minlen_arcsec):
     # for each night, test if it has a tracklet
     b = 0
     for k, e in enumerate(i):
-        hasTrk[k] = hasTracklet(mjd[b:e], ra[b:e], dec[b:e], maxdt_minutes, minlen_arcsec)
+        hasTrk[k] = hasTracklet(mjd[b:e], ra[b:e], dec[b:e], maxdt_minutes, minlen_arcsec, min_observations)
         b = e
 
     return nights, hasTrk
@@ -232,7 +239,9 @@ def discoveryOpportunities(nights, nightHasTracklets, window, nlink, p, rng):
     return discIdx, disc
 
 
-def linkObject(obsv, seed, maxdt_minutes, minlen_arcsec, window, nlink, p, night_start_utc_days):
+def linkObject(
+    obsv, seed, maxdt_minutes, minlen_arcsec, min_observations, window, nlink, p, night_start_utc_days
+):
     """
     For a set of observations of a single object, calculate if there are any tracklets,
     if there are enough tracklets to form a discovery window, and then report back all of
@@ -262,6 +271,9 @@ def linkObject(obsv, seed, maxdt_minutes, minlen_arcsec, window, nlink, p, night
     minlen_arcsec : float
         Minimum allowable distance separation between observations [Units: arcsec]
 
+    min_observations (int):
+        the minimum number of observations in a night required to form a tracklet.
+
     window : float
         Number of tracklets required with <= this window to complete a detection
 
@@ -313,7 +325,9 @@ def linkObject(obsv, seed, maxdt_minutes, minlen_arcsec, window, nlink, p, night
         seed %= 0xFFFF_FFFF
         rng = np.random.default_rng(seed)
 
-        nights, hasTrk = trackletsInNights(night, mjd, ra, dec, maxdt_minutes, minlen_arcsec)
+        nights, hasTrk = trackletsInNights(
+            night, mjd, ra, dec, maxdt_minutes, minlen_arcsec, min_observations
+        )
         discIdx, discNights = discoveryOpportunities(nights, hasTrk, window, nlink, p, rng)
         if discIdx != -1:
             discoveryChances = len(discNights)

src/sorcha/utilities/sorchaConfigs.py

@@ -532,6 +532,13 @@ def _validate_linkingfilter_configs(self):
                 logging.error("ERROR: ssp_number_observations is zero or negative.")
                 sys.exit("ERROR: ssp_number_observations is zero or negative.")
 
+            if self.ssp_number_observations > 1 and self.ssp_separation_threshold == 0.0:
+                logging.error("ERROR: ssp_separation_threshold is zero.")
+                sys.exit("ERROR: ssp_separation_threshold is zero.")
+            if self.ssp_separation_threshold < 0.0:
+                logging.error("ERROR: ssp_separation_threshold is negative.")
+                sys.exit("ERROR: ssp_separation_threshold is negative.")
+
             if self.ssp_number_tracklets < 1:
                 logging.error("ERROR: ssp_number_tracklets is zero or less.")
                 sys.exit("ERROR: ssp_number_tracklets is zero or less.")
@@ -544,10 +551,6 @@ def _validate_linkingfilter_configs(self):
                 logging.error("ERROR: ssp_detection_efficiency out of bounds (should be between 0 and 1).")
                 sys.exit("ERROR: ssp_detection_efficiency out of bounds (should be between 0 and 1).")
 
-            if self.ssp_separation_threshold <= 0.0:
-                logging.error("ERROR: ssp_separation_threshold is zero or negative.")
-                sys.exit("ERROR: ssp_separation_threshold is zero or negative.")
-
             if self.ssp_maximum_time < 0:
                 logging.error("ERROR: ssp_maximum_time is negative.")
                 sys.exit("ERROR: ssp_maximum_time is negative.")