From 9ae41e7baf18490b969c3319db842571235ab4c9 Mon Sep 17 00:00:00 2001 From: Ryan Lyttle Date: Tue, 1 Apr 2025 17:14:51 +0100 Subject: [PATCH 1/7] Added config file linkingfilter.ssp_number_observations into PPMiniDifi Added linkingfilter.ssp_number_observations into MiniDifi function hasTracklet. Now the function is no longer hard coded for 2 observations to make a tracklet and depends on config variable ssp_number_observations. Updated test_PPLinkingFilter unit test and ran notebook demo_miniDifiValidation to check these changes with the brute force method to ensure the linking filter still works the same way. --- docs/notebooks/demo_miniDifiValidation.ipynb | 10 +++--- src/sorcha/modules/PPLinkingFilter.py | 1 + src/sorcha/modules/PPMiniDifi.py | 34 +++++++++++++++----- tests/sorcha/test_PPLinkingFilter.py | 2 ++ 4 files changed, 34 insertions(+), 13 deletions(-) diff --git a/docs/notebooks/demo_miniDifiValidation.ipynb b/docs/notebooks/demo_miniDifiValidation.ipynb index 56bc55fe..45147b23 100644 --- a/docs/notebooks/demo_miniDifiValidation.ipynb +++ b/docs/notebooks/demo_miniDifiValidation.ipynb @@ -223,8 +223,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 s, sys: 50.8 ms, total: 11.1 s\n", + "Wall time: 11.1 s\n" ] }, { @@ -276,8 +276,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 746 ms, sys: 151 ms, total: 897 ms\n", + "Wall time: 1.22 s\n" ] } ], @@ -297,7 +297,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)" ] }, { diff --git a/src/sorcha/modules/PPLinkingFilter.py b/src/sorcha/modules/PPLinkingFilter.py index 87489381..3f842713 100755 --- a/src/sorcha/modules/PPLinkingFilter.py +++ b/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, diff --git a/src/sorcha/modules/PPMiniDifi.py b/src/sorcha/modules/PPMiniDifi.py index 926e314c..cd6f9003 100755 --- a/src/sorcha/modules/PPMiniDifi.py +++ b/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,11 +81,14 @@ 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 >= 1 and min_observations == 1: # edge case when 1 observation is needed for a tracklet + return True if nobs < 2: return False maxdt = maxdt_minutes / (60 * 24) minlen = minlen_arcsec / 3600 + detection_pair_count = 0 # counting number of detection pairs for i in range(nobs): for j in range(nobs): @@ -90,13 +96,15 @@ def hasTracklet(mjd, ra, dec, maxdt_minutes, minlen_arcsec): if diff > 0 and diff < maxdt: sep = haversine_np(ra[i], dec[i], ra[j], dec[j]) if sep > minlen: - return True - - return False + detection_pair_count += 1 + if detection_pair_count >= (min_observations - 1): + return True + else: + 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 +129,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 +150,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 +243,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 +275,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 +329,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) diff --git a/tests/sorcha/test_PPLinkingFilter.py b/tests/sorcha/test_PPLinkingFilter.py index 5858b7f8..a5f88bdb 100755 --- a/tests/sorcha/test_PPLinkingFilter.py +++ b/tests/sorcha/test_PPLinkingFilter.py @@ -111,6 +111,7 @@ def test_PPLinkingFilter_discoveryChances(): objectId="ssObjectId", maxdt_minutes=max_time_separation * 24 * 60, minlen_arcsec=min_angular_separation, + min_observations = min_observations, window=min_tracklet_window, nlink=min_tracklets, p=detection_efficiency, @@ -158,6 +159,7 @@ def test_PPLinkingFilter_nlink1(): objectId="ssObjectId", maxdt_minutes=max_time_separation * 24 * 60, minlen_arcsec=min_angular_separation, + min_observations = min_observations, window=min_tracklet_window, nlink=min_tracklets, p=detection_efficiency, From 83f36c9e0f402b8df8e405aba380ea1dcbb8889c Mon Sep 17 00:00:00 2001 From: Ryan Lyttle Date: Wed, 2 Apr 2025 17:20:56 +0100 Subject: [PATCH 2/7] Changed logic in hasTracklet function in MiniDifi and created unit tests to check this works --- src/sorcha/modules/PPMiniDifi.py | 10 ++- tests/sorcha/test_PPLinkingFilter.py | 114 +++++++++++++++++++++++++++ 2 files changed, 121 insertions(+), 3 deletions(-) diff --git a/src/sorcha/modules/PPMiniDifi.py b/src/sorcha/modules/PPMiniDifi.py index cd6f9003..78915556 100755 --- a/src/sorcha/modules/PPMiniDifi.py +++ b/src/sorcha/modules/PPMiniDifi.py @@ -88,7 +88,7 @@ def hasTracklet(mjd, ra, dec, maxdt_minutes, minlen_arcsec, min_observations): maxdt = maxdt_minutes / (60 * 24) minlen = minlen_arcsec / 3600 - detection_pair_count = 0 # counting number of detection pairs + detection_indexes = [] # index of objects linked to eachother in pairs for i in range(nobs): for j in range(nobs): @@ -96,8 +96,12 @@ def hasTracklet(mjd, ra, dec, maxdt_minutes, minlen_arcsec, min_observations): if diff > 0 and diff < maxdt: sep = haversine_np(ra[i], dec[i], ra[j], dec[j]) if sep > minlen: - detection_pair_count += 1 - if detection_pair_count >= (min_observations - 1): + detection_indexes.append(i) + detection_indexes.append(j) # store indexes of tracklet pairs that have pairs + + if ( + len(set(detection_indexes)) >= min_observations + ): # removes duplicates and checks minimum number of observations return True else: return False diff --git a/tests/sorcha/test_PPLinkingFilter.py b/tests/sorcha/test_PPLinkingFilter.py index a5f88bdb..0d9d92af 100755 --- a/tests/sorcha/test_PPLinkingFilter.py +++ b/tests/sorcha/test_PPLinkingFilter.py @@ -1,6 +1,7 @@ import pandas as pd import numpy as np + from sorcha.utilities.dataUtilitiesForTests import get_test_filepath @@ -395,3 +396,116 @@ def test_PPLinkingFilter_nodrop(): assert all(~linked_observations[linked_observations["ObjID"] == "unlinked_object"]["object_linked"]) assert len(linked_observations[linked_observations["ObjID"] == "linked_object"]) == 6 assert len(linked_observations[linked_observations["ObjID"] == "unlinked_object"]) == 6 + + +def test_PPLinkingFilter_minobs(): + """ + Testing min_observations for one night to consider a tracklet + """ + from sorcha.modules.PPLinkingFilter import PPLinkingFilter + + # min obs 5 and 5 detections on 1 night. should return all 5 observations + min_observations = 5 + min_angular_separation = 0.5 + max_time_separation = 0.0625 + min_tracklets = 1 + min_tracklet_window = 15 + detection_efficiency = 1 + night_start_utc = 17.0 + + obj_id = ["pretend_object"] * 5 + field_id = np.arange(1, 6) + t0 = 60000.0 + times = np.asarray([0.03, 0.06,0.07,0.08,0.09]) + t0 + ra = [142, 142.1,142.2,142.3,142.4] + dec = [8, 8.1,8.2,8.3,8.4] + + observations = pd.DataFrame( + {"ObjID": obj_id, "FieldID": field_id, "fieldMJD_TAI": times, "RA_deg": ra, "Dec_deg": dec} + ) + + linked_observations = PPLinkingFilter( + observations, + detection_efficiency, + min_observations, + min_tracklets, + min_tracklet_window, + min_angular_separation, + max_time_separation, + night_start_utc, + ) + + assert len(linked_observations) == 5 + + # min obs 6 and 5 detections on 1 night. Should return empty dataframe + + + min_observations = 6 + min_angular_separation = 0.5 + max_time_separation = 0.0625 + min_tracklets = 1 + min_tracklet_window = 15 + detection_efficiency = 1 + night_start_utc = 17.0 + + # create object that should link 6 observations on one night to tracklet but only give it 5 observations for that night. + obj_id = ["pretend_object"] * 5 + field_id = np.arange(1, 6) + t0 = 60000.0 + times = np.asarray([0.03, 0.06,0.07,0.08,0.09]) + t0 + ra = [142, 142.1,142.2,142.3,142.4] + dec = [8, 8.1,8.2,8.3,8.4] + + observations = pd.DataFrame( + {"ObjID": obj_id, "FieldID": field_id, "fieldMJD_TAI": times, "RA_deg": ra, "Dec_deg": dec} + ) + + linked_observations = PPLinkingFilter( + observations, + detection_efficiency, + min_observations, + min_tracklets, + min_tracklet_window, + min_angular_separation, + max_time_separation, + night_start_utc, + ) + + assert len(linked_observations) == 0 + + # min obs 1 and 2 detections on two different nights + + min_observations = 1 + min_angular_separation = 0 + max_time_separation = 0 + min_tracklets = 1 + min_tracklet_window = 15 + detection_efficiency = 1 + night_start_utc = 17.0 + + # create object that should link 1 observation on one night to tracklet. + obj_id = ["pretend_object"] * 2 + field_id = np.arange(1,3) + t0 = 60000.0 + times = np.asarray([0.03,1000]) + t0 + ra = [142,152] + dec = [8,9] + + observations = pd.DataFrame( + {"ObjID": obj_id, "FieldID": field_id, "fieldMJD_TAI": times, "RA_deg": ra, "Dec_deg": dec} + ) + + linked_observations = PPLinkingFilter( + observations, + detection_efficiency, + min_observations, + min_tracklets, + min_tracklet_window, + min_angular_separation, + max_time_separation, + night_start_utc, + ) + print(linked_observations) + assert len(linked_observations) == 2 + + \ No newline at end of file From b2676ae29724fc257dbace590a2c92092d833192 Mon Sep 17 00:00:00 2001 From: Ryan Lyttle Date: Thu, 3 Apr 2025 16:02:05 +0100 Subject: [PATCH 3/7] Update PPMiniDifi.py --- src/sorcha/modules/PPMiniDifi.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/sorcha/modules/PPMiniDifi.py b/src/sorcha/modules/PPMiniDifi.py index 78915556..f7a4c184 100755 --- a/src/sorcha/modules/PPMiniDifi.py +++ b/src/sorcha/modules/PPMiniDifi.py @@ -81,11 +81,11 @@ def hasTracklet(mjd, ra, dec, maxdt_minutes, minlen_arcsec, min_observations): ## tracklets by taking all observations in a night and computing ## all of theirs pairwise distances, then selecting on that. nobs = len(ra) - if nobs >= 1 and min_observations == 1: # edge case when 1 observation is needed for a tracklet - return True - 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 detection_indexes = [] # index of objects linked to eachother in pairs @@ -97,11 +97,9 @@ def hasTracklet(mjd, ra, dec, maxdt_minutes, minlen_arcsec, min_observations): sep = haversine_np(ra[i], dec[i], ra[j], dec[j]) if sep > minlen: detection_indexes.append(i) - detection_indexes.append(j) # store indexes of tracklet pairs that have pairs + detection_indexes.append(j) # store indexes of observations that are linked. - if ( - len(set(detection_indexes)) >= min_observations - ): # removes duplicates and checks minimum number of observations + if len((detection_indexes)) >= min_observations: return True else: return False From 073bafce25c19e66350bf0ccbd68882c71c5ee2c Mon Sep 17 00:00:00 2001 From: Ryan Lyttle Date: Thu, 3 Apr 2025 16:22:24 +0100 Subject: [PATCH 4/7] Changed linkingfilter config class checks Changed linkingfilter config class checks to so that ssp_separation_threshold can be zero when ssp_number_observations =1 --- src/sorcha/utilities/sorchaConfigs.py | 11 +++++++---- tests/sorcha/test_sorchaConfigs.py | 13 ++++++++++++- 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/src/sorcha/utilities/sorchaConfigs.py b/src/sorcha/utilities/sorchaConfigs.py index f1a7ea78..ad4ef9f0 100644 --- a/src/sorcha/utilities/sorchaConfigs.py +++ b/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.") diff --git a/tests/sorcha/test_sorchaConfigs.py b/tests/sorcha/test_sorchaConfigs.py index 989f8f4c..0ff73c60 100644 --- a/tests/sorcha/test_sorchaConfigs.py +++ b/tests/sorcha/test_sorchaConfigs.py @@ -725,12 +725,23 @@ def test_linkingfilter_bounds(key_name): test_configs = linkingfilterConfigs(**linkingfilter_configs) assert error_text.value.code == f"ERROR: {key_name} is negative." - elif key_name == "ssp_separation_threshold" or key_name == "ssp_number_observations": + elif key_name == "ssp_number_observations": linkingfilter_configs[key_name] = -5 with pytest.raises(SystemExit) as error_text: test_configs = linkingfilterConfigs(**linkingfilter_configs) assert error_text.value.code == f"ERROR: {key_name} is zero or negative." + elif key_name == "ssp_separation_threshold": + linkingfilter_configs[key_name] = -5 + with pytest.raises(SystemExit) as error_text: + test_configs = linkingfilterConfigs(**linkingfilter_configs) + assert error_text.value.code == f"ERROR: {key_name} is negative." + + linkingfilter_configs[key_name] = 0 + with pytest.raises(SystemExit) as error_text: + test_configs = linkingfilterConfigs(**linkingfilter_configs) + assert error_text.value.code == f"ERROR: {key_name} is zero." + elif key_name == "ssp_number_tracklets": linkingfilter_configs[key_name] = -5 with pytest.raises(SystemExit) as error_text: From 0a299009da27caab55732ea60f5c99a2f76a9057 Mon Sep 17 00:00:00 2001 From: Ryan Lyttle Date: Tue, 8 Apr 2025 12:01:34 +0100 Subject: [PATCH 5/7] Fixed if statement in miniDifi to not count duplicate detections Fixed if statement in miniDifi to not count duplicate detections using set function. Added another unit test to check case where not enough observations are linked to count as a tracklet. --- src/sorcha/modules/PPMiniDifi.py | 2 +- tests/sorcha/test_PPLinkingFilter.py | 36 ++++++++++++++++++++++++++++ 2 files changed, 37 insertions(+), 1 deletion(-) diff --git a/src/sorcha/modules/PPMiniDifi.py b/src/sorcha/modules/PPMiniDifi.py index f7a4c184..68ce9b7d 100755 --- a/src/sorcha/modules/PPMiniDifi.py +++ b/src/sorcha/modules/PPMiniDifi.py @@ -99,7 +99,7 @@ def hasTracklet(mjd, ra, dec, maxdt_minutes, minlen_arcsec, min_observations): detection_indexes.append(i) detection_indexes.append(j) # store indexes of observations that are linked. - if len((detection_indexes)) >= min_observations: + if len(set(detection_indexes)) >= min_observations: return True else: return False diff --git a/tests/sorcha/test_PPLinkingFilter.py b/tests/sorcha/test_PPLinkingFilter.py index 0d9d92af..42e65b52 100755 --- a/tests/sorcha/test_PPLinkingFilter.py +++ b/tests/sorcha/test_PPLinkingFilter.py @@ -437,6 +437,42 @@ def test_PPLinkingFilter_minobs(): assert len(linked_observations) == 5 + + # min obs 5 and 5 detections on 1 night, but only 4 observations can link into a tracklet. + # should return 0 + min_observations = 5 + min_angular_separation = 0.5 + max_time_separation = 0.0625 + min_tracklets = 1 + min_tracklet_window = 15 + detection_efficiency = 1 + night_start_utc = 17.0 + + obj_id = ["pretend_object"] * 5 + field_id = np.arange(1, 6) + t0 = 60000.0 + times = np.asarray([0.03, 0.06,0.07,0.08,0.15]) + t0 #last time has difference > max_time_separation so unlinkable + ra = [142, 142.1,142.2,142.3,142.15] + dec = [8, 8.1,8.2,8.3,8.4] + + observations = pd.DataFrame( + {"ObjID": obj_id, "FieldID": field_id, "fieldMJD_TAI": times, "RA_deg": ra, "Dec_deg": dec} + ) + + linked_observations = PPLinkingFilter( + observations, + detection_efficiency, + min_observations, + min_tracklets, + min_tracklet_window, + min_angular_separation, + max_time_separation, + night_start_utc, + ) + + assert len(linked_observations) == 0 + + # min obs 6 and 5 detections on 1 night. Should return empty dataframe From d692cb5d63ee7d86cb621b46c9cd7df478595569 Mon Sep 17 00:00:00 2001 From: Ryan Lyttle Date: Tue, 8 Apr 2025 14:50:42 +0100 Subject: [PATCH 6/7] Updated demo_miniDifiValidation to allow different minimum observations for a tracklet Updated demo_miniDifiValidation to allow different minimum observations for a tracklet in the brute force linking test. Changed miniDifi.hasTracklet function to be more efficient. --- docs/notebooks/demo_miniDifiValidation.ipynb | 34 ++++++++++---------- src/sorcha/modules/PPMiniDifi.py | 8 ++--- 2 files changed, 21 insertions(+), 21 deletions(-) diff --git a/docs/notebooks/demo_miniDifiValidation.ipynb b/docs/notebooks/demo_miniDifiValidation.ipynb index 45147b23..aea80c2a 100644 --- a/docs/notebooks/demo_miniDifiValidation.ipynb +++ b/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 11 s, sys: 50.8 ms, total: 11.1 s\n", - "Wall time: 11.1 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 746 ms, sys: 151 ms, total: 897 ms\n", - "Wall time: 1.22 s\n" + "CPU times: user 1.1 s, sys: 190 ms, total: 1.29 s\n", + "Wall time: 1.47 s\n" ] } ], @@ -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": { diff --git a/src/sorcha/modules/PPMiniDifi.py b/src/sorcha/modules/PPMiniDifi.py index 68ce9b7d..a2416677 100755 --- a/src/sorcha/modules/PPMiniDifi.py +++ b/src/sorcha/modules/PPMiniDifi.py @@ -99,10 +99,10 @@ def hasTracklet(mjd, ra, dec, maxdt_minutes, minlen_arcsec, min_observations): detection_indexes.append(i) detection_indexes.append(j) # store indexes of observations that are linked. - if len(set(detection_indexes)) >= min_observations: - return True - else: - return False + if len(set(detection_indexes)) >= min_observations: + return True + + return False @njit(cache=True) From d8888300825dcda04fd537f4153622170a9b48a5 Mon Sep 17 00:00:00 2001 From: Ryan Lyttle Date: Wed, 30 Apr 2025 13:00:37 -0700 Subject: [PATCH 7/7] Updated miniDifi changes discussed with Mario --- src/sorcha/modules/PPMiniDifi.py | 8 +--- tests/sorcha/test_PPLinkingFilter.py | 70 +++++++++++++++++++++++++++- 2 files changed, 70 insertions(+), 8 deletions(-) diff --git a/src/sorcha/modules/PPMiniDifi.py b/src/sorcha/modules/PPMiniDifi.py index a2416677..6bf71fea 100755 --- a/src/sorcha/modules/PPMiniDifi.py +++ b/src/sorcha/modules/PPMiniDifi.py @@ -83,23 +83,17 @@ def hasTracklet(mjd, ra, dec, maxdt_minutes, minlen_arcsec, min_observations): nobs = len(ra) 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 - detection_indexes = [] # index of objects linked to eachother in pairs - 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: - detection_indexes.append(i) - detection_indexes.append(j) # store indexes of observations that are linked. - - if len(set(detection_indexes)) >= min_observations: + if (i + 1) - (j) >= min_observations: return True return False diff --git a/tests/sorcha/test_PPLinkingFilter.py b/tests/sorcha/test_PPLinkingFilter.py index 42e65b52..dba2fdf8 100755 --- a/tests/sorcha/test_PPLinkingFilter.py +++ b/tests/sorcha/test_PPLinkingFilter.py @@ -544,4 +544,72 @@ def test_PPLinkingFilter_minobs(): print(linked_observations) assert len(linked_observations) == 2 - \ No newline at end of file + + + # min obs 4 and 4 detections on 1 night, but only 2 observations can link into a tracklet either side. + # should return 0 + min_observations = 4 + min_angular_separation = 0.5 + max_time_separation = 0.0625 + min_tracklets = 1 + min_tracklet_window = 15 + detection_efficiency = 1 + night_start_utc = 17.0 + + obj_id = ["pretend_object"] * 4 + field_id = np.arange(1, 5) + t0 = 60000.0 + times = np.asarray([0.03, 0.04,0.15,0.16]) + t0 # obs split in two with time difference > max_time_separation so unlinkable + ra = [142, 142.1,142.5,142.6] + dec = [8, 8.1,8.5,8.6] + + observations = pd.DataFrame( + {"ObjID": obj_id, "FieldID": field_id, "fieldMJD_TAI": times, "RA_deg": ra, "Dec_deg": dec} + ) + + linked_observations = PPLinkingFilter( + observations, + detection_efficiency, + min_observations, + min_tracklets, + min_tracklet_window, + min_angular_separation, + max_time_separation, + night_start_utc, + ) + + assert len(linked_observations) == 0 + + # min obs 4 and 6 detections on 1 night, where only 4 observations can link into a tracklet. + # should return 6 + min_observations = 4 + min_angular_separation = 0.5 + max_time_separation = 0.0625 + min_tracklets = 1 + min_tracklet_window = 15 + detection_efficiency = 1 + night_start_utc = 17.0 + + obj_id = ["pretend_object"] * 6 + field_id = np.arange(1, 7) + t0 = 60000.0 + times = np.asarray([0.03, 0.04,0.15,0.16,0.17,0.18]) + t0 # obs split in two with time difference > max_time_separation 2 on left and 4 on right + ra = [142, 142.1,142.5,142.6,142.7,142.8] + dec = [8, 8.1,8.5,8.6,8.7,8.8] + + observations = pd.DataFrame( + {"ObjID": obj_id, "FieldID": field_id, "fieldMJD_TAI": times, "RA_deg": ra, "Dec_deg": dec} + ) + + linked_observations = PPLinkingFilter( + observations, + detection_efficiency, + min_observations, + min_tracklets, + min_tracklet_window, + min_angular_separation, + max_time_separation, + night_start_utc, + ) + + assert len(linked_observations) == 6 \ No newline at end of file