Merge pull request #1080 from dirac-institute/Sanity-checks-from-orbits

Sanity checks from orbits and swapping --verbose with --log-level

docs/example_files/help_output.txt

@@ -31,5 +31,5 @@ Optional arguments:
                         Default False. (default: False)
   -s S, --survey S      Survey to simulate (default: rubin_sim)
   -t T, --stem T        Output file name stem. (default: SSPPOutput)
-  -v, --verbose         Verbosity. Default currently true; include to turn off
+  -l, --log-level         Verbosity. Default currently true; include to turn off
                         verbosity. (default: True)

src/sorcha/ephemeris/simulation_driver.py

@@ -102,7 +102,7 @@ def create_ephemeris(orbits_df, pointings_df, args, sconfigs):
     computed.  Details for those that actually do land within the field
     of view are passed along.
     """
-    verboselog = args.pplogger.info if args.verbose else lambda *a, **k: None
+    verboselog = args.pplogger.info if args.loglevel else lambda *a, **k: None
 
     ang_fov = sconfigs.simulation.ar_ang_fov
     buffer = sconfigs.simulation.ar_fov_buffer
@@ -351,7 +351,7 @@ def write_out_ephemeris_file(ephemeris_df, ephemeris_csv_filename, args, sconfig
     None.
     """
 
-    verboselog = args.pplogger.info if args.verbose else lambda *a, **k: None
+    verboselog = args.pplogger.info if args.loglevel else lambda *a, **k: None
 
     if sconfigs.input.eph_format == "csv":
         verboselog("Outputting ephemeris to CSV file...")

src/sorcha/modules/PPCommandLineParser.py

@@ -85,7 +85,7 @@ def PPCommandLineParser(args):
 
     cmd_args_dict["surveyname"] = args.s
     cmd_args_dict["outfilestem"] = args.t
-    cmd_args_dict["verbose"] = args.v
+    cmd_args_dict["loglevel"] = args.l
     cmd_args_dict["stats"] = args.st
 
     if cmd_args_dict["stats"] is not None:

src/sorcha/readers/OrbitAuxReader.py

@@ -92,6 +92,64 @@ def _process_and_validate_input_table(self, input_table, **kwargs):
             if not all(column in input_table.columns for column in keplerian_elements):
                 pplogger.error("ERROR: PPReadOrbitFile: Must provide all keplerian orbital elements.")
                 sys.exit("ERROR: PPReadOrbitFile: Must provide all keplerian orbital elements.")
+
+            filtered_df = input_table.loc[(input_table["inc"] < 0) | (input_table["inc"] > 180), ["ObjID"]]
+            if not filtered_df.empty:
+                pplogger.warning(
+                    f"WARNING: For Object(s) {', '.join(filtered_df['ObjID'].astype(str))}. Inclination (inc) is outside the valid range (0 180 degrees), which may cause the orbits to be mirrored affecting orbital calculations."
+                )
+            filtered_df = input_table.loc[
+                ~(
+                    (
+                        (input_table["argPeri"] >= 0)
+                        & (input_table["argPeri"] <= 360)
+                        & (input_table["node"] >= 0)
+                        & (input_table["node"] <= 360)
+                        & (input_table["ma"] >= 0)
+                        & (input_table["ma"] <= 360)
+                    )
+                    | (
+                        (input_table["argPeri"] >= -180)
+                        & (input_table["argPeri"] <= 180)
+                        & (input_table["node"] >= -180)
+                        & (input_table["node"] <= 180)
+                        & (input_table["ma"] >= -180)
+                        & (input_table["ma"] <= 180)
+                    )
+                ),
+                ["ObjID"],
+            ]
+            if not filtered_df.empty:
+                pplogger.warning(
+                    f"WARNING: For Object(s) {', '.join(filtered_df['ObjID'].astype(str))}. The argument of Perihelion (argPeri), the Longitude of the Ascending Node (node), and the Mean Anomaly (ma) are not in the same bounds (either [0 360] or [-180 180] degrees), which may lead to incorrect orbital calculations."
+                )
+            error_raised = False
+            # checks all objects before a system exit
+            filtered_df = input_table.loc[(input_table["e"] >= 1), ["ObjID"]]
+            if not filtered_df.empty:
+                error_raised = True
+                pplogger.error(
+                    f"ERROR: For Object(s) {', '.join(filtered_df['ObjID'].astype(str))}. Parabolic orbit (e == 1) and Hyperbolic orbit (e > 1) are not supported for Keplerian elements."
+                )
+            filtered_df = input_table.loc[(input_table["e"] < 0), ["ObjID"]]
+            if not filtered_df.empty:
+                error_raised = True
+                pplogger.error(
+                    f"ERROR: For Object(s) {', '.join(filtered_df['ObjID'].astype(str))}. Eccentricity (e) cannot be less than 0."
+                )
+            filtered_df = input_table.loc[(input_table["e"] < 1) & (input_table["a"] < 0), ["ObjID"]]
+            if not filtered_df.empty:
+                error_raised = True
+                pplogger.error(
+                    f"ERROR: For Object(s) {', '.join(filtered_df['ObjID'].astype(str))}. Bound orbit (e < 1) with negative semi-major axis (a < 0) is not physical."
+                )
+
+            if error_raised:
+                pplogger.error("ERROR: Invalid Keplerian elements detected for one or more objects.")
+                sys.exit(
+                    "ERROR: Invalid Keplerian elements detected for one or more objects (check log for information)"
+                )
+
         elif orb_format in ["COM", "BCOM"]:
             if not all(column in input_table.columns for column in cometary_elements):
                 pplogger.error("ERROR: PPReadOrbitFile: Must provide all cometary orbital elements.")
@@ -100,6 +158,38 @@ def _process_and_validate_input_table(self, input_table, **kwargs):
                 pplogger.warning(
                     "WARNING: At least one t_p_MJD_TDB is above 2400000.5 - make sure your t_p are MJD and not in JD"
                 )
+
+            error_raised = False
+            filtered_df = input_table.loc[(input_table["q"] == 0), ["ObjID"]]
+            if not filtered_df.empty:
+                pplogger.warning(
+                    f"WARNING: For Object(s) {', '.join(filtered_df['ObjID'].astype(str))}. q==0 is technically correct but suggests a collisional path with an object instead of an orbital path."
+                )
+
+            filtered_df = input_table.loc[(input_table["inc"] < 0) | (input_table["inc"] > 180), ["ObjID"]]
+            if not filtered_df.empty:
+                pplogger.warning(
+                    f"WARNING: For Object(s) {', '.join(filtered_df['ObjID'].astype(str))}. Inclination (inc) is outside the valid range (0 180 degrees), which may cause the orbits to be mirrored affecting orbital calculations."
+                )
+            filtered_df = input_table.loc[(input_table["q"] < 0), ["ObjID"]]
+            if not filtered_df.empty:
+                error_raised = True
+                pplogger.error(
+                    f"ERROR: For Object(s) {', '.join(filtered_df['ObjID'].astype(str))}. Perihelion distance (q) cannot be less than 0."
+                )
+            filtered_df = input_table.loc[(input_table["e"] < 0), ["ObjID"]]
+            if not filtered_df.empty:
+                error_raised = True
+                pplogger.error(
+                    f"ERROR: For Object(s) {', '.join(filtered_df['ObjID'].astype(str))}. Eccentricity (e) cannot be less than 0."
+                )
+
+            if error_raised:
+                pplogger.error("ERROR: Invalid cometary elements detected for one or more objects")
+                sys.exit(
+                    "ERROR: Invalid cometary elements detected for one or more objects (check log for information)"
+                )
+
         elif orb_format in ["CART", "BCART"]:
             if not all(column in input_table.columns for column in cartesian_elements):
                 pplogger.error("ERROR: PPReadOrbitFile: Must provide all cartesian coordinate values.")

src/sorcha/sorcha.py

@@ -111,7 +111,7 @@ def runLSSTSimulation(args, sconfigs):
 
     # if verbosity flagged, the verboselog function will log the message specified
     # if not, verboselog does absolutely nothing
-    verboselog = pplogger.info if args.verbose else lambda *a, **k: None
+    verboselog = pplogger.info if args.loglevel else lambda *a, **k: None
 
     sconfigs.filters.mainfilter, sconfigs.filters.othercolours = PPGetMainFilterAndColourOffsets(
         args.paramsinput, sconfigs.filters.observing_filters, sconfigs.input.aux_format
@@ -214,7 +214,7 @@ def runLSSTSimulation(args, sconfigs):
             sconfigs.filters.observing_filters,
             sconfigs.activity.comet_activity,
             lightcurve_choice=sconfigs.lightcurve.lc_model,
-            verbose=args.verbose,
+            verbose=args.loglevel,
         )
 
         if sconfigs.expert.trailing_losses_on:
@@ -239,15 +239,15 @@ def runLSSTSimulation(args, sconfigs):
         # Do NOT use trailedSourceMagTrue or PSFMagTrue, these are the unrandomised magnitudes.
         verboselog("Calculating astrometric and photometric uncertainties...")
         observations = PPAddUncertainties.addUncertainties(
-            observations, sconfigs, args._rngs, verbose=args.verbose
+            observations, sconfigs, args._rngs, verbose=args.loglevel
         )
 
         if sconfigs.expert.randomization_on:
             verboselog(
                 "Number of rows BEFORE randomizing astrometry and photometry: " + str(len(observations.index))
             )
             observations = PPRandomizeMeasurements.randomizeAstrometryAndPhotometry(
-                observations, sconfigs, args._rngs, verbose=args.verbose
+                observations, sconfigs, args._rngs, verbose=args.loglevel
             )
             verboselog(
                 "Number of rows AFTER randomizing astrometry and photometry: " + str(len(observations.index))
@@ -271,7 +271,7 @@ def runLSSTSimulation(args, sconfigs):
             verboselog("Applying field-of-view filters...")
             verboselog("Number of rows BEFORE applying FOV filters: " + str(len(observations.index)))
             observations = PPApplyFOVFilter(
-                observations, sconfigs, args._rngs, footprint=footprint, verbose=args.verbose
+                observations, sconfigs, args._rngs, footprint=footprint, verbose=args.loglevel
             )
             verboselog("Number of rows AFTER applying FOV filters: " + str(len(observations.index)))
 
@@ -299,7 +299,7 @@ def runLSSTSimulation(args, sconfigs):
                 sconfigs.fadingfunction.fading_function_peak_efficiency,
                 sconfigs.fadingfunction.fading_function_width,
                 args._rngs,
-                verbose=args.verbose,
+                verbose=args.loglevel,
             )
             verboselog("Number of rows AFTER applying fading function: " + str(len(observations.index)))
 
@@ -332,7 +332,7 @@ def runLSSTSimulation(args, sconfigs):
         if len(observations.index) > 0:
             pplogger.info("Post processing completed for this chunk")
             pplogger.info("Outputting results for this chunk")
-            PPWriteOutput(args, sconfigs, observations, verbose=args.verbose)
+            PPWriteOutput(args, sconfigs, observations, verbose=args.loglevel)
             if args.stats is not None:
                 stats(observations, args.stats, args.outpath, sconfigs)
         else:

src/sorcha/utilities/diffTestUtils.py

@@ -54,7 +54,7 @@ def compare_result_files(test_output, golden_output):
     "pointing_database": get_demo_filepath("baseline_v2.0_1yr.db"),
     "surveyname": "rubin_sim",
     "outfilestem": f"out_end2end",
-    "verbose": False,
+    "loglevel": False,
     "stats": None,
 }
 
@@ -66,7 +66,7 @@ def compare_result_files(test_output, golden_output):
     "output_ephemeris_file": "sorcha_ephemeris",
     "surveyname": "rubin_sim",
     "outfilestem": f"out_end2end_with_ephemeris_generation",
-    "verbose": False,
+    "loglevel": False,
     "stats": None,
 }
 
@@ -78,7 +78,7 @@ def compare_result_files(test_output, golden_output):
     "pointing_database": get_demo_filepath("baseline_v2.0_1yr.db"),
     "surveyname": "rubin_sim",
     "outfilestem": f"out_end2end_chunked",
-    "verbose": False,
+    "loglevel": False,
     "stats": None,
 }
 
@@ -90,7 +90,7 @@ def compare_result_files(test_output, golden_output):
     "pointing_database": get_demo_filepath("baseline_v2.0_1yr.db"),
     "surveyname": "rubin_sim",
     "outfilestem": f"out_end2end_unchunked",
-    "verbose": False,
+    "loglevel": False,
     "stats": None,
 }
 
@@ -103,7 +103,7 @@ def compare_result_files(test_output, golden_output):
     "output_ephemeris_file": "sorcha_ephemeris.csv",
     "surveyname": "rubin_sim",
     "outfilestem": f"verification_output",
-    "verbose": False,
+    "loglevel": False,
     "stats": None,
 }
 

src/sorcha/utilities/sorchaArguments.py

@@ -25,7 +25,7 @@ class sorchaArguments:
     outfilestem: str = ""
     """file system for output"""
 
-    verbose: bool = False
+    loglevel: bool = False
     """logger verbosity"""
 
     surveyname: str = ""
@@ -71,7 +71,7 @@ def read_from_dict(self, args):
         self.pointing_database = args["pointing_database"]
         self.output_ephemeris_file = args.get("output_ephemeris_file")
         self.ar_data_file_path = args.get("ar_data_path")
-        self.verbose = args["verbose"]
+        self.loglevel = args["loglevel"]
         self.stats = args["stats"]
 
         self.surveyname = args["surveyname"]

src/sorcha_cmdline/run.py

@@ -105,10 +105,10 @@ def main():
         "-t", "--stem", help="Output file name stem.", type=str, dest="t", default="SSPPOutput"
     )
     optional.add_argument(
-        "-v",
-        "--verbose",
+        "-l",
+        "--log-level",
         help="Print additional information to log while running",
-        dest="v",
+        dest="l",
         default=True,
         action="store_false",
     )

tests/data/orbit_test_files/orbit_sanity_check_com_e<0.csv

@@ -0,0 +1,3 @@
+ObjID,FORMAT,q,e,inc,node,argPeri,t_p_MJD_TDB,epochMJD_TDB
+S00000t,COM,0.952105479028,-0.504888475701,4.899098347472,148.881068605772,39.949789586436,54486.32292808,54466.0
+

tests/data/orbit_test_files/orbit_sanity_check_com_q<0.csv

@@ -0,0 +1,3 @@
+ObjID,FORMAT,q,e,inc,node,argPeri,t_p_MJD_TDB,epochMJD_TDB
+S00000t,COM,-0.952105479028,0.504888475701,4.899098347472,148.881068605772,39.949789586436,54486.32292808,54466.0
+

tests/data/orbit_test_files/orbit_sanity_check_kep_e<0.des

@@ -0,0 +1,2 @@
+ObjID		a		e		inc		node		argPeri		ma	  epochMJD_TDB	FORMAT
+2010_TU149	2.2103867	-1	1.9679		58.91911	92.60419	324.21145	60200	KEP

tests/data/orbit_test_files/orbit_sanity_check_kep_e<1_a<0.des

@@ -0,0 +1,2 @@
+ObjID		a		e		inc		node		argPeri		ma	  epochMJD_TDB	FORMAT
+2010_TU149	-2.2103867	0.5	1.9679		58.91911	92.60419	324.21145	60200	KEP

tests/data/orbit_test_files/orbit_sanity_check_kep_e=1.des

@@ -0,0 +1,2 @@
+ObjID		a		e		inc		node		argPeri		ma	  epochMJD_TDB	FORMAT
+2010_TU149	2.2103867	1	1.9679		58.91911	92.60419	324.21145	60200	KEP

tests/data/orbit_test_files/orbit_sanity_check_kep_e>1.des

@@ -0,0 +1,2 @@
+ObjID		a		e		inc		node		argPeri		ma	  epochMJD_TDB	FORMAT
+2010_TU149	2.2103867	2	1.9679		58.91911	92.60419	324.21145	60200	KEP

tests/ephemeris/test_ephemeris_generation.py

@@ -116,7 +116,7 @@ def test_ephemeris_end2end(single_synthetic_pointing, tmp_path):
         "outpath": tmp_path,
         "surveyname": "rubin_sim",
         "outfilestem": f"out_400k",
-        "verbose": False,
+        "loglevel": False,
         "stats": None,
     }
 
@@ -161,7 +161,7 @@ def test_ephemeris_writeread_csv(single_synthetic_ephemeris, tmp_path):
     params_in = get_demo_filepath("sspp_testset_colours.txt")
 
     class args(object):
-        verbose = False
+        loglevel = False
 
     cmd_args = args()
 
@@ -199,7 +199,7 @@ def test_ephemeris_writeread_whitespace(single_synthetic_ephemeris, tmp_path):
     params_in = get_demo_filepath("sspp_testset_colours.txt")
 
     class args(object):
-        verbose = False
+        loglevel = False
 
     cmd_args = args()
 
@@ -237,7 +237,7 @@ def test_ephemeris_writeread_hdf5(single_synthetic_ephemeris, tmp_path):
     params_in = get_demo_filepath("sspp_testset_colours.txt")
 
     class args(object):
-        verbose = False
+        loglevel = False
 
     cmd_args = args()
 

tests/ephemeris/test_pixdict.py

@@ -61,7 +61,7 @@ def test_pixeldict(tmp_path):
         "deleteTemporaryEphemerisDatabase": False,
         "surveyname": "rubin_sim",
         "outfilestem": f"out_400k",
-        "verbose": False,
+        "loglevel": False,
         "stats": None,
     }
 

tests/readers/test_OrbitAuxReader.py

@@ -306,3 +306,48 @@ def test_orbit_reader_wrong_delim():
 
     with pytest.raises(SystemExit) as err:
         _ = OrbitAuxReader(get_test_filepath("testorb.des"), "csv")
+
+def test_orbit_sanity_check_kep():
+    """If an orbit parameter is undefined, raise an exception"""
+
+    #for keplerian e==1
+    with pytest.raises(SystemExit) as err:
+        _ = OrbitAuxReader(get_test_filepath("orbit_test_files/orbit_sanity_check_kep_e=1.des"),"whitespace")
+        _.read_rows()
+    assert err.value.code == "ERROR: Invalid Keplerian elements detected for one or more objects (check log for information)"
+
+     #for keplerian e<0
+    with pytest.raises(SystemExit) as err:
+        _ = OrbitAuxReader(get_test_filepath("orbit_test_files/orbit_sanity_check_kep_e<0.des"),"whitespace")
+        _.read_rows()
+    assert err.value.code == "ERROR: Invalid Keplerian elements detected for one or more objects (check log for information)"
+
+     #for keplerian e>1
+    with pytest.raises(SystemExit) as err:
+        _ = OrbitAuxReader(get_test_filepath("orbit_test_files/orbit_sanity_check_kep_e>1.des"),"whitespace")
+        _.read_rows()
+    assert err.value.code == "ERROR: Invalid Keplerian elements detected for one or more objects (check log for information)"
+
+     #for keplerian e<1 a<0
+    with pytest.raises(SystemExit) as err:
+        _ = OrbitAuxReader(get_test_filepath("orbit_test_files/orbit_sanity_check_kep_e<1_a<0.des"),"whitespace")
+        _.read_rows()
+    assert err.value.code == "ERROR: Invalid Keplerian elements detected for one or more objects (check log for information)"
+
+
+
+def test_orbit_sanity_check_com():
+    """If an orbit parameter is undefined, raise an exception"""
+
+    #for comentary  q<0
+    with pytest.raises(SystemExit) as err:
+        _ = OrbitAuxReader(get_test_filepath("orbit_test_files/orbit_sanity_check_com_q<0.csv"),"csv")
+        _.read_rows()
+    assert err.value.code == "ERROR: Invalid cometary elements detected for one or more objects (check log for information)"
+
+    #for comentary  e<0
+    with pytest.raises(SystemExit) as err:
+        _ = OrbitAuxReader(get_test_filepath("orbit_test_files/orbit_sanity_check_com_e<0.csv"),"csv")
+        _.read_rows()
+    assert err.value.code == "ERROR: Invalid cometary elements detected for one or more objects (check log for information)"
+