Merge branch 'feature/CO2_fitting_refinement' into 'master'

CO2 fitting algorithm refinement

See merge request caimira/caimira!503

Author: lrdossan

caimira/apps/calculator/__init__.py

@@ -33,6 +33,7 @@
 from . import markdown_tools
 from . import model_generator, co2_model_generator
 from .report_generator import ReportGenerator, calculate_report_data
+from .co2_report_generator import CO2ReportGenerator
 from .user import AuthenticatedUser, AnonymousUser
 
 # The calculator version is based on a combination of the model version and the
@@ -404,7 +405,10 @@ async def post(self, endpoint: str) -> None:
 
         requested_model_config = tornado.escape.json_decode(self.request.body)
         try:
-            form = co2_model_generator.CO2FormData.from_dict(requested_model_config, data_registry)
+            form: co2_model_generator.CO2FormData = co2_model_generator.CO2FormData.from_dict(
+                requested_model_config, 
+                data_registry
+            )
         except Exception as err:
             if self.settings.get("debug", False):
                 import traceback
@@ -414,29 +418,21 @@ async def post(self, endpoint: str) -> None:
             self.finish(json.dumps(response_json))
             return
 
+        CO2_report_generator: CO2ReportGenerator = CO2ReportGenerator()
         if endpoint.rstrip('/') == 'plot':
-            transition_times = co2_model_generator.CO2FormData.find_change_points_with_pelt(form.CO2_data)
-            self.finish({'CO2_plot': co2_model_generator.CO2FormData.generate_ventilation_plot(form.CO2_data, transition_times),
-                        'transition_times': [round(el, 2) for el in transition_times]})
+            report = CO2_report_generator.build_initial_plot(form)
+            self.finish(report)
         else:
             executor = loky.get_reusable_executor(
                 max_workers=self.settings['handler_worker_pool_size'],
                 timeout=300,
             )
             report_task = executor.submit(
-                co2_model_generator.CO2FormData.build_model, form,
+                CO2_report_generator.build_fitting_results, form,
             )
+            
             report = await asyncio.wrap_future(report_task)
-
-            result = dict(report.CO2_fit_params())
-            ventilation_transition_times = report.ventilation_transition_times
-
-            result['fitting_ventilation_type'] = form.fitting_ventilation_type
-            result['transition_times'] = ventilation_transition_times
-            result['CO2_plot'] = co2_model_generator.CO2FormData.generate_ventilation_plot(CO2_data=form.CO2_data,
-                                                                transition_times=ventilation_transition_times[:-1],
-                                                                predictive_CO2=result['predictive_CO2'])
-            self.finish(result)
+            self.finish(report)
 
 
 def get_url(app_root: str, relative_path: str = '/'):

caimira/apps/calculator/co2_model_generator.py

@@ -2,8 +2,9 @@
 import logging
 import typing
 import numpy as np
-import ruptures as rpt
 import matplotlib.pyplot as plt
+from scipy.signal import find_peaks
+import pandas as pd
 import re
 
 from caimira import models
@@ -21,22 +22,20 @@
 class CO2FormData(FormData):
     CO2_data: dict
     fitting_ventilation_states: list
-    fitting_ventilation_type: str
     room_capacity: typing.Optional[int]
 
     #: The default values for undefined fields. Note that the defaults here
     #: and the defaults in the html form must not be contradictory.
     _DEFAULTS: typing.ClassVar[typing.Dict[str, typing.Any]] = {
         'CO2_data': '{}',
+        'fitting_ventilation_states': '[]',
         'exposed_coffee_break_option': 'coffee_break_0',
         'exposed_coffee_duration': 5,
         'exposed_finish': '17:30',
         'exposed_lunch_finish': '13:30',
         'exposed_lunch_option': True,
         'exposed_lunch_start': '12:30',
         'exposed_start': '08:30',
-        'fitting_ventilation_states': '[]',
-        'fitting_ventilation_type': 'fitting_natural_ventilation',
         'infected_coffee_break_option': 'coffee_break_0',
         'infected_coffee_duration': 5,
         'infected_dont_have_breaks_with_exposed': False,
@@ -97,55 +96,75 @@ def validate(self):
                             raise TypeError(f'Unable to fetch "finish_time" key. Got "{dict_keys[1]}".')
                         for time in input_break.values():
                             if not re.compile("^(2[0-3]|[01]?[0-9]):([0-5]?[0-9])$").match(time):
-                                raise TypeError(f'Wrong time format - "HH:MM". Got "{time}".')
+                                raise TypeError(f'Wrong time format - "HH:MM". Got "{time}".')        
 
-    @classmethod
-    def find_change_points_with_pelt(self, CO2_data: dict):
+    def find_change_points(self) -> list:
         """
-        Perform change point detection using Pelt algorithm from ruptures library with pen=15.
-        Returns a list of tuples containing (index, X-axis value) for the detected significant changes.
+        Perform change point detection using scipy library (find_peaks method) with rolling average of data.
+        Incorporate existing state change candidates and adjust the result accordingly.
+        Returns a list of the detected ventilation state changes, discarding any occupancy state change.
         """
-
-        times: list = CO2_data['times']
-        CO2_values: list = CO2_data['CO2']
+        times: list = self.CO2_data['times']
+        CO2_values: list = self.CO2_data['CO2']
 
         if len(times) != len(CO2_values):
             raise ValueError("times and CO2 values must have the same length.")
 
-        # Convert the input list to a numpy array for use with the ruptures library
-        CO2_np = np.array(CO2_values)
+        # Time difference between two consecutive time data entries, in seconds
+        diff = (times[1] - times[0]) * 3600 # Initial data points in absolute hours, e.g. 14.78
+
+        # Calculate minimum interval for smoothing technique
+        smooth_min_interval_in_minutes = 1 # Minimum time difference for smooth technique
+        window_size = max(int((smooth_min_interval_in_minutes * 60) // diff), 1)
 
-        # Define the model for change point detection (Radial Basis Function kernel)
-        model = "rbf"
+        # Applying a rolling average to smooth the initial data
+        smoothed_co2 = pd.Series(CO2_values).rolling(window=window_size, center=True).mean()
 
-        # Fit the Pelt algorithm to the data with the specified model
-        algo = rpt.Pelt(model=model).fit(CO2_np)
+        # Calculate minimum interval for peaks and valleys detection
+        peak_valley_min_interval_in_minutes = 15 # Minimum time difference between two peaks or two valleys
+        min_distance_points = max(int((peak_valley_min_interval_in_minutes * 60) // diff), 1)
 
-        # Predict change points using the Pelt algorithm with a penalty value of 15
-        result = algo.predict(pen=15)
+        # Calculate minimum width of datapoints for valley detection
+        width_min_interval_in_minutes = 20 # Minimum duration of a valley
+        min_valley_width = max(int((width_min_interval_in_minutes * 60) // diff), 1)
 
-        # Find local minima and maxima
-        segments = np.split(np.arange(len(CO2_values)), result)
-        merged_segments = [np.hstack((segments[i], segments[i + 1])) for i in range(len(segments) - 1)]
-        result_set = set()
-        for segment in merged_segments[:-2]:
-            result_set.add(times[CO2_values.index(min(CO2_np[segment]))])
-            result_set.add(times[CO2_values.index(max(CO2_np[segment]))])
-        return list(result_set)
+        # Find peaks (maxima) in the smoothed data applying the distance factor
+        peaks, _ = find_peaks(smoothed_co2.values, prominence=100, distance=min_distance_points)
+        
+        # Find valleys (minima) by inverting the smoothed data and applying the width and distance factors
+        valleys, _ = find_peaks(-smoothed_co2.values, prominence=50, width=min_valley_width, distance=min_distance_points)
 
-    @classmethod
-    def generate_ventilation_plot(self, CO2_data: dict,
-                                  transition_times: typing.Optional[list] = None,
-                                  predictive_CO2: typing.Optional[list] = None):
-            times_values = CO2_data['times']
-            CO2_values = CO2_data['CO2']
+        # Extract peak and valley timestamps
+        timestamps = np.array(times)
+        peak_timestamps = timestamps[peaks]
+        valley_timestamps = timestamps[valleys]
+
+        return sorted(np.concatenate((peak_timestamps, valley_timestamps)))
+
+    def generate_ventilation_plot(self,
+                                  ventilation_transition_times: typing.Optional[list] = None,
+                                  occupancy_transition_times: typing.Optional[list] = None,
+                                  predictive_CO2: typing.Optional[list] = None) -> str:
+            
+            # Plot data (x-axis: times; y-axis: CO2 concentrations)
+            times_values: list = self.CO2_data['times']
+            CO2_values: list = self.CO2_data['CO2']
 
             fig = plt.figure(figsize=(7, 4), dpi=110)
             plt.plot(times_values, CO2_values, label='Input CO₂')
 
-            if (transition_times):
-                for time in transition_times:
-                    plt.axvline(x = time, color = 'grey', linewidth=0.5, linestyle='--')
+            # Add occupancy state changes:
+            if (occupancy_transition_times):
+                for i, time in enumerate(occupancy_transition_times):
+                    plt.axvline(x = time, color = 'grey', linewidth=0.5, linestyle='--', label='Occupancy change (from input)' if i == 0 else None)
+            # Add ventilation state changes:
+            if (ventilation_transition_times):
+                for i, time in enumerate(ventilation_transition_times):
+                    if i == 0:
+                        label = 'Ventilation change (detected)' if occupancy_transition_times else 'Ventilation state changes'
+                    else: label = None
+                    plt.axvline(x = time, color = 'red', linewidth=0.5, linestyle='--', label=label)
+
             if (predictive_CO2):
                 plt.plot(times_values, predictive_CO2, label='Predictive CO₂')
             plt.xlabel('Time of day')
@@ -158,14 +177,18 @@ def population_present_changes(self, infected_presence: models.Interval, exposed
         state_change_times.update(exposed_presence.transition_times())
         return sorted(state_change_times)
 
-    def ventilation_transition_times(self) -> typing.Tuple[float, ...]:
-        # Check what type of ventilation is considered for the fitting
-        if self.fitting_ventilation_type == 'fitting_natural_ventilation':
-            vent_states = self.fitting_ventilation_states
-            vent_states.append(self.CO2_data['times'][-1])
-            return tuple(vent_states)
-        else:
-            return tuple((self.CO2_data['times'][0], self.CO2_data['times'][-1]))
+    def ventilation_transition_times(self) -> typing.Tuple[float]:
+        '''
+        Check if the last time from the input data is
+        included in the ventilation ventilations state.
+        Given that the last time is a required state change, 
+        if not included, this method adds it.
+        '''
+        vent_states = self.fitting_ventilation_states
+        last_time_from_input = self.CO2_data['times'][-1]
+        if (vent_states and last_time_from_input != vent_states[-1]): # The last time value is always needed for the last ACH interval.
+            vent_states.append(last_time_from_input)
+        return tuple(vent_states)
 
     def build_model(self, size=None) -> models.CO2DataModel: # type: ignore
         size = size or self.data_registry.monte_carlo['sample_size']
@@ -184,7 +207,7 @@ def build_model(self, size=None) -> models.CO2DataModel: # type: ignore
             activity=None, # type: ignore
         )
 
-        all_state_changes=self.population_present_changes(infected_presence, exposed_presence)
+        all_state_changes = self.population_present_changes(infected_presence, exposed_presence)
         total_people = [infected_population.people_present(stop) + exposed_population.people_present(stop)
                         for _, stop in zip(all_state_changes[:-1], all_state_changes[1:])]
 

caimira/apps/calculator/co2_report_generator.py

@@ -0,0 +1,68 @@
+import dataclasses
+import typing
+
+from caimira.models import CO2DataModel, Interval, IntPiecewiseConstant
+from .co2_model_generator import CO2FormData
+
+
+@dataclasses.dataclass
+class CO2ReportGenerator:
+    
+    def build_initial_plot(
+            self,
+            form: CO2FormData,
+    ) -> dict:
+        '''
+        Initial plot with the suggested ventilation state changes. 
+        This method receives the form input and returns the CO2
+        plot with the respective transition times.
+        '''
+        CO2model: CO2DataModel = form.build_model()
+
+        occupancy_transition_times = list(CO2model.occupancy.transition_times)
+        
+        ventilation_transition_times: list = form.find_change_points()
+        # The entire ventilation changes consider the initial and final occupancy state change
+        all_vent_transition_times: list = sorted(
+            [occupancy_transition_times[0]] + 
+            [occupancy_transition_times[-1]] + 
+            ventilation_transition_times)
+
+        ventilation_plot: str = form.generate_ventilation_plot(
+            ventilation_transition_times=all_vent_transition_times,
+            occupancy_transition_times=occupancy_transition_times
+        )
+
+        context = {
+            'CO2_plot': ventilation_plot,
+            'transition_times': [round(el, 2) for el in all_vent_transition_times],
+        }
+
+        return context
+    
+    def build_fitting_results(
+            self,
+            form: CO2FormData,
+    ) -> dict:
+        '''
+        Final fitting results with the respective predictive CO2. 
+        This method receives the form input and returns the fitting results
+        along with the CO2 plot with the predictive CO2.
+        '''
+        CO2model: CO2DataModel = form.build_model()
+
+        # Ventilation times after user manipulation from the suggested ventilation state change times.
+        ventilation_transition_times = list(CO2model.ventilation_transition_times)
+
+        # The result of the following method is a dict with the results of the fitting
+        # algorithm, namely the breathing rate and ACH values. It also returns the
+        # predictive CO2 result based on the fitting results.
+        context: typing.Dict = dict(CO2model.CO2_fit_params())
+
+        # Add the transition times and CO2 plot to the results.
+        context['transition_times'] = ventilation_transition_times
+        context['CO2_plot'] = form.generate_ventilation_plot(ventilation_transition_times=ventilation_transition_times[:-1], 
+                                                       predictive_CO2=context['predictive_CO2'])
+
+        return context
+    

caimira/apps/calculator/model_generator.py

@@ -330,13 +330,10 @@ def ventilation(self) -> models._VentilationBase:
                                                     min(self.infected_start, self.exposed_start)/60)
         if self.ventilation_type == 'from_fitting':
             ventilations = []
-            if self.CO2_fitting_result['fitting_ventilation_type'] == 'fitting_natural_ventilation':
-                transition_times = self.CO2_fitting_result['transition_times']
-                for index, (start, stop) in enumerate(zip(transition_times[:-1], transition_times[1:])):
-                    ventilations.append(models.AirChange(active=models.SpecificInterval(present_times=((start, stop), )),
-                                                         air_exch=self.CO2_fitting_result['ventilation_values'][index]))
-            else:
-                ventilations.append(models.AirChange(active=always_on, air_exch=self.CO2_fitting_result['ventilation_values'][0]))
+            transition_times = self.CO2_fitting_result['transition_times']
+            for index, (start, stop) in enumerate(zip(transition_times[:-1], transition_times[1:])):
+                ventilations.append(models.AirChange(active=models.SpecificInterval(present_times=((start, stop), )),
+                                                        air_exch=self.CO2_fitting_result['ventilation_values'][index]))
             return models.MultipleVentilation(tuple(ventilations))
 
         # Initializes a ventilation instance as a window if 'natural_ventilation' is selected, or as a HEPA-filter otherwise

caimira/apps/calculator/report_generator.py

@@ -368,9 +368,9 @@ def readable_minutes(minutes: int) -> str:
 
 def hour_format(hour: float) -> str:
     # Convert float hour to HH:MM format
-    hours = int(hour)
-    minutes = int(hour % 1 * 60)
-    return f"{hours}:{minutes if minutes != 0 else '00'}"
+    hours = f"{int(hour):02}"
+    minutes = f"{int(hour % 1 * 60):02}"
+    return f"{hours}:{minutes}"
 
 
 def percentage(absolute: float) -> float: