Revisions to score.py

Author: joseph-hellerstein

.gitignore

@@ -8,6 +8,8 @@ __pycache__/
 
 # C extensions
 *.so
+*.py.sav
+err_*
 
 # Distribution / packaging
 .Python

scripts/calculate_linear_prediction_scores.py

@@ -146,7 +146,7 @@ def processModels(first_model_num: int, last_model_num: int,
                     threshold=threshold)
             discovery.fit()
             prediction_df = discovery.predict()
-            score.add(timecourse.timecourse_df, prediction_df, label=model_name)
+            score.add(timecourse.timecourse_df, prediction_df, system_id=model_name)
         except Exception as e:
             print(f"Error occurred while processing model {model_name}: {e}")
             continue

src/constants.py

@@ -39,6 +39,21 @@
 COL_NUM_PERTURBATION = "num_perturbation"
 COL_NUM_TIMEPOINT = "num_timepoint"
 COL_NAMES = [COL_MODEL_NAME, COL_MAXCV, COL_ENDTIME, COL_ENDTIME_SOURCE]
+COL_SYSTEM_ID = "system_id"  # Unique identifier for the system, e.g. model name or species name.
+COL_MEAN = "mean" # Mean value of the valid values used in calculating the statistics.
+COL_MIN = "min"  # Minimum value of the valid values used in calculating the statistics.
+COL_MAX = "max"  # Maximum value of the valid values used in calculating the statistics.
+COL_COUNT = "count"  # Number of valid values used in calculating the statistics.
+COL_INVALID_COUNT = "invalid_count"  # Number of invalid (sentinel -1) values excluded from aggregation.
+COL_LABEL = "label"  # Unique identifier for the row of statistics, e.g. model name or species name.
+COL_P25 = "p25"  # 25th percentile of the valid values used in calculating the statistics.
+COL_P30 = "p30"  # 30th percentile of the valid values used in calculating the statistics.
+COL_P50 = "p50"  # 50th percentile of the valid values used in calculating the statistics.
+COL_P80 = "p80"  # 80th percentile of the valid values used in calculating the statistics.
+COL_P95 = "p95"  # 95th percentile of the valid values used in calculating the statistics.
+COL_P99 = "p99"  # 99th percentile of the valid values used
+COL_PERCENTILES = [COL_P25, COL_P30, COL_P50, COL_P80, COL_P95, COL_P99]
+STATISTICS = [COL_MEAN, COL_MIN, COL_MAX, COL_COUNT, COL_INVALID_COUNT] + COL_PERCENTILES
 
 # Symbolic values
 ENDTIME_SOURCE_RECIROCAL_MIN_EIGENVALUE = "reciprocal_min_eigenvalue"

src/score.py

@@ -7,106 +7,15 @@
     MAPE = max(0, 1 - (abs(prediction - true) / true))
 """
 
-from src.dataframe_serializer import DataframeSerializer  # type: ignore
 from src.plot_options import PlotOptions  # type: ignore
 import src.constants as cn  # type: ignore
+from src.statistic_calculator import StatisticCalculator # type: ignore
 
 import matplotlib.pyplot as plt  # type: ignore
 import numpy as np  # type: ignore
 import pandas as pd  # type: ignore
 
-MEAN = "mean"
-MIN = "min"
-MAX = "max"
-COUNT = "count"
 
-#########################################
-class _StatisticAccumulator(object):
-    STATISTICS = ["mean", "min", "max", "count", "invalid_count", "p25", "p30",
-            "p50", "p80", "p95", "p99"]
-    # count is the number of valid values used in calculating the statistics.
-    # invalid_count is the number of invalid (sentinel -1) values excluded from aggregation.
-
-    """A container for storing statistics accumulated in a dictionary"""
-    def __init__(self) -> None:
-        self.statistic_dct: dict = {n: [] for n in self.STATISTICS}
-        self.statistic_dct[cn.AGGREGATION_TYPE] = []
-        self.statistic_dct[cn.DESCRIPTION] = []
-
-    def add(self,
-            value_arr: np.ndarray = np.array([]),
-            label: str = "",
-            aggregation_type: str = "",  # "model" or species name
-    ) -> None:
-        """Computes statistics and adds to cumulative values.
-
-        Parameters
-        ----------
-        value_arr : np.ndarray
-            1D Array of MAPE metric values. Negative sentinel values (-1)
-            indicate invalid measurements and are excluded from aggregation.
-        label : str
-            Descriptive label for this aggregation.
-        aggregation_type : str
-            Type of aggregation, either "model" or a species name.
-
-        Returns
-        -------
-        None
-            Modifies internal state (self.statistic_dct) only.
-        """
-        value_arr = value_arr.flatten()  # Flatten to 1D for aggregation
-        LARGE_VAL = 1e6
-        # Filter out negative sentinel values (invalid/undefined true values).
-        valid_mask = value_arr >= 0
-        invalid_count = int(np.sum(~valid_mask))
-        valid_arr = value_arr[valid_mask].copy()
-
-        count = int(len(valid_arr))
-
-        # If no input values at all, still record a row so the CSV reflects that data was collected.
-        if len(value_arr) == 0:
-            self.statistic_dct[cn.DESCRIPTION].append(label)
-            self.statistic_dct[cn.AGGREGATION_TYPE].append(aggregation_type)
-            self.statistic_dct[MEAN].append(0.0)
-            self.statistic_dct[MIN].append(0.0)
-            self.statistic_dct[MAX].append(0.0)
-            self.statistic_dct[COUNT].append(0)
-            self.statistic_dct["invalid_count"].append(invalid_count)
-            for p in [p for p in self.STATISTICS if p.startswith("p")]:
-                self.statistic_dct[p].append(0.0)
-            return
-
-        # If all values are invalid, still record a row so the CSV reflects that data was collected.
-        if count == 0:
-            self.statistic_dct[cn.DESCRIPTION].append(label)
-            self.statistic_dct[cn.AGGREGATION_TYPE].append(aggregation_type)
-            self.statistic_dct[MEAN].append(0.0)
-            self.statistic_dct[MIN].append(0.0)
-            self.statistic_dct[MAX].append(0.0)
-            self.statistic_dct[COUNT].append(0)
-            self.statistic_dct["invalid_count"].append(invalid_count)
-            for p in [p for p in self.STATISTICS if p.startswith("p")]:
-                self.statistic_dct[p].append(0.0)
-            return
-
-        # Replace remaining NaN/inf/large values with LARGE_VAL for aggregation.
-        sel = np.isnan(valid_arr) | np.isinf(valid_arr) | (valid_arr > LARGE_VAL)
-        valid_arr[sel] = LARGE_VAL
-        # Update dictionary
-        self.statistic_dct[cn.DESCRIPTION].append(label)
-        self.statistic_dct[cn.AGGREGATION_TYPE].append(aggregation_type)
-        self.statistic_dct[MEAN].append(float(np.nanmean(valid_arr)))
-        self.statistic_dct[MIN].append(float(np.nanmin(valid_arr)))
-        self.statistic_dct[MAX].append(float(np.nanmax(valid_arr)))
-        self.statistic_dct[COUNT].append(count)
-        self.statistic_dct["invalid_count"].append(invalid_count)
-        percentiles = [p for p in self.STATISTICS if p.startswith("p")]
-        for p in percentiles:
-            self.statistic_dct[p].append(float(np.nanpercentile(valid_arr, int(p[1:]))))
-
-
-#########################################
 class Score:
     """Scores prediction timecourses against true timecourses using zero-floor MAPE.
 
@@ -129,20 +38,17 @@ def __init__(self, serialization_path: str = "",
         is_persist : bool
             Whether to persist the DataFrame to the CSV file.
         """
+        self._is_persist = is_persist
         if len(serialization_path) == 0:
             serialization_path = self.SERIALIZATION_PATH
-        self._serializer = DataframeSerializer(serialization_path,
-                is_initialize=is_initialize, is_persist=is_persist)
+        self._serialization_path = serialization_path
         #
-        self.statistic_accumulator = _StatisticAccumulator()
+        self.score_df = pd.DataFrame()
+        self.statistic_calculator = StatisticCalculator()
 
     @property
     def serialization_path(self) -> str:
-        return self._serializer.serialization_path
-
-    @property
-    def score_df(self) -> pd.DataFrame:
-        return self._serializer.dataframe
+        return self._serialization_path
 
     @staticmethod 
     def calculateMAPE(true_df: pd.DataFrame,
@@ -180,14 +86,14 @@ def calculateMAPE(true_df: pd.DataFrame,
     def add(self,
             true_timecourse_df: pd.DataFrame,
             prediction_timecourse_df: pd.DataFrame,
-            label: str = "",
+            system_id: str = "",
             ) -> pd.DataFrame:
-        """Computes MAPE scores and accumulates statistics for model-level and per-species aggregations.
+        """Computes scores and accumulates statistics for model-level and per-species aggregations.
 
         Parameters
         ----------
-        label : str
-            Descriptive label stored in each aggregation row.
+        system_id : str
+            Descriptive ID stored in each aggregation row.
         true_timecourse_df : pd.DataFrame
             True timecourse with timepoints as index and species as columns.
         prediction_timecourse_df : pd.DataFrame
@@ -198,38 +104,18 @@ def add(self,
         pd.DataFrame
             The full score DataFrame after this addition.
         """
-        score_df: pd.DataFrame = self.calculateMAPE(true_timecourse_df, prediction_timecourse_df)
-
-        # Record how many rows were in the accumulator BEFORE this call.
-        dct = self.statistic_accumulator.statistic_dct
-        start_count = len(dct[cn.DESCRIPTION])
-
-        # Model level aggregation (all species and timepoints combined)
-        model_arr = np.asarray(score_df.values, dtype=float)
-        self.statistic_accumulator.add(model_arr,
-                aggregation_type=cn.AGGREGATION_TYPE_MODEL,
-                label=label)
-
+        mape_df: pd.DataFrame = self.calculateMAPE(true_timecourse_df, prediction_timecourse_df)
+        self.statistic_calculator.add(cn.AGGREGATION_TYPE_MODEL, mape_df.values.flatten())
         # Species level aggregations (one per species column, across all timepoints)
-        species_names = list(score_df.columns)
+        species_names = list(mape_df.columns)
         for species_name in species_names:
-            species_arr = np.asarray(score_df[species_name].values, dtype=float)
-            self.statistic_accumulator.add(species_arr,
-                    aggregation_type=species_name,
-                    label=label)
-
-        # Build list of dicts from only the NEW rows added in this call.
-        end_count = len(dct[cn.DESCRIPTION])
-        new_rows: list[dict] = []
-        for i in range(start_count, end_count):
-            row: dict = {}
-            for key in dct:
-                if isinstance(dct[key], list):
-                    row[key] = dct[key][i]
-            new_rows.append(row)
-
-        # Persist the new rows to the DataFrame/CSV.
-        self._serializer.serializeDct(new_rows)
+            self.statistic_calculator.add(species_name, mape_df[species_name].to_numpy())
+        # Add the system ID
+        self.score_df = self.statistic_calculator.dataframe.copy()
+        self.score_df[cn.COL_SYSTEM_ID] = system_id
+        # Serialize the accumulated statistics
+        if self._is_persist:
+            self.score_df.to_csv(self.serialization_path, index=False)
         return self.score_df
 
     def plotCDF(self, 

src/statistic_calculator.py

@@ -0,0 +1,104 @@
+'''Calculates univariate descriptive statistics'''
+
+# Statistics that begin with 'p' are percentiles, e.g. p25 is the 25th percentile
+# of the valid values used in calculating the statistics.
+# Assumes that all statistics
+"""
+Usage:
+    statistics = StatisticCalculator()
+    statistics.add(label="model1", value_arr=np.array([1, 2, 3, -1, 4]))
+    statistics.add(label="model2", value_arr=np.array([5, 6, 7, 8, 9]))
+    df = statistics.dataframe
+"""
+
+
+import src.constants as cn  # type: ignore
+
+import numpy as np  # type: ignore
+import pandas as pd  # type: ignore
+
+LARGE_VAL = 1e6
+
+
+#########################################
+class StatisticCalculator(object):
+    # count is the number of valid values used in calculating the statistics.
+    # invalid_count is the number of invalid (sentinel -1) values excluded from aggregation.
+
+    """A container for storing statistics accumulated in a dictionary"""
+    def __init__(self) -> None:
+        self.statistic_dct: dict = {n: [] for n in cn.STATISTICS}
+        self.statistic_dct[cn.COL_LABEL] = []
+
+    @property
+    def dataframe(self) -> pd.DataFrame:
+        """Returns a DataFrame of the accumulated statistics.
+
+        Returns
+        -------
+        pd.DataFrame
+            DataFrame containing the accumulated statistics.
+        """
+        return pd.DataFrame(self.statistic_dct)
+    
+    def _is_percentile(self, stat_name: str) -> bool:
+        """Returns True if the statistic name is a percentile (e.g. p25, p50, etc.)"""
+        return stat_name.startswith("p") and stat_name[1:].isdigit() and 0 <= int(stat_name[1:]) <= 100
+
+    def add(self,
+            label: str,
+            value_arr: np.ndarray = np.array([]),
+            is_non_negative: bool = True,
+    ) -> None:
+        """Computes statistics and accumulates them in the internal dictionary (self.statistic_dct).
+
+        Parameters
+        ----------
+        label : str
+            Descriptive label for this aggregation.
+        value_arr : np.ndarray
+            1D Array of values. Negative sentinel values (-1)
+            indicate invalid measurements and are excluded from aggregation.
+        is_non_negative : bool
+            If True, only non-negative values are considered valid for aggregation.
+
+        Returns
+        -------
+        None
+            Modifies internal state (self.statistic_dct) only.
+        """
+        value_arr = value_arr.flatten()  # Flatten to 1D for aggregation
+        # Filter out NaN, inf, and optionally negative values.
+        if is_non_negative:
+            valid_mask = [not (np.isnan(v) or np.isinf(v)) and v >= 0 for v in value_arr]
+        else:
+            valid_mask = [not (np.isnan(v) or np.isinf(v)) for v in value_arr]
+        count = int(np.sum(valid_mask))  # Count of valid values used in aggregation
+        invalid_count = len(value_arr) - int(np.sum(valid_mask))
+
+        # If no input values at all, still record a row so the CSV reflects that data was collected.
+        if (len(value_arr) == 0) or (count == 0):
+            self.statistic_dct[cn.COL_LABEL].append(label)
+            self.statistic_dct[cn.COL_MEAN].append(np.nan)
+            self.statistic_dct[cn.COL_MIN].append(np.nan)
+            self.statistic_dct[cn.COL_MAX].append(np.nan)
+            self.statistic_dct[cn.COL_COUNT].append(0)
+            self.statistic_dct[cn.COL_INVALID_COUNT].append(invalid_count)
+            for p in [p for p in cn.STATISTICS if self._is_percentile(p)]:
+                self.statistic_dct[p].append(np.nan)
+            return
+
+        # Replace large values with LARGE_VAL for aggregation.
+        valid_arr = value_arr[valid_mask].copy()
+        sel = valid_arr > LARGE_VAL
+        valid_arr[sel] = LARGE_VAL
+        # Update dictionary
+        self.statistic_dct[cn.COL_LABEL].append(label)
+        self.statistic_dct[cn.COL_MEAN].append(float(np.nanmean(valid_arr)))
+        self.statistic_dct[cn.COL_MIN].append(float(np.nanmin(valid_arr)))
+        self.statistic_dct[cn.COL_MAX].append(float(np.nanmax(valid_arr)))
+        self.statistic_dct[cn.COL_COUNT].append(count)
+        self.statistic_dct[cn.COL_INVALID_COUNT].append(invalid_count)
+        percentiles = [p for p in cn.STATISTICS if self._is_percentile(p)]
+        for p in percentiles:
+            self.statistic_dct[p].append(float(np.nanpercentile(valid_arr, int(p[1:]))))