Writing code to obtain TOC Calibration curves.

defs/Helper.py

@@ -15,6 +15,7 @@
 ### IMPORTS AND CONFIGURATIONS
 
 import pandas as pd
+import numpy as np
 
 ### ----------------------------------------------------------------------------------------------------------------------------------- ###
 ### Helper.py
@@ -64,4 +65,11 @@ class DATFile():
 
     def exclude_outliers():
             
-        pass
+        pass
+
+def compute_R2(y_true, y_est):
+
+    residuals = y_est - y_true
+    distances = y_true - np.mean(y_true)
+
+    return 1 - np.sum(np.power(residuals, 2)) / np.sum(np.power(distances, 2))

defs/TOC.py

@@ -16,39 +16,102 @@
 
 import numpy as np
 import pandas as pd
-import math
+from pathlib import Path, PurePath
+import datetime
 import warnings
+from scipy.optimize import curve_fit
+from defs.Helper import compute_R2
 
 ### ----------------------------------------------------------------------------------------------------------------------------------- ###
 ### TOC.py
 
-class TC_Calibration():
+class CalibrationCurve():
 
-    def init(self, cal_file):
+    def __init__(self):
+
+        self.exists = False
 
         self.slope = None
         self.intercept = None
         self.R2_value = None
+        self.excluded = None
+        self.lower_bound = None
+        self.upper_bound = None
+        self.computed_on = None
 
+    def __linear_function(self, x_values, a, b):
+
+        if not isinstance(x_values, np.ndarray):
+            x_values = np.array(x_values)
+
+        return a * x_values + b
+
+    def __optimise(self, x_data, y_data, p0):
+        return curve_fit(self.__linear_function, x_data, y_data, p0)
+
+    def __check_exist(self):
+
+        if self.exists == True:
+            return True
+        else:
+            return False
+
+    def __check(self):
+
+        if self.__check_exist():
+            warnings.warn("Calibration alreadyt exists, will be overwritten...")
+        else:
+            return
+
+    def calibrate(self, x_data, y_data):
+
+            self.__check()
+
+            parameters, _ = self.__optimise(x_data=x_data, y_data=y_data, p0=[10, 0.1])
+
+            self.exists = True
+            self.slope = parameters[0]
+            self.intercept = parameters[1]
+            self.R2_value = compute_R2(y_true = y_data, y_est = self.__linear_function(x_values=x_data, a = parameters[0], b=parameters[1]))
+            self.lower_bound = np.min(y_data)
+            self.upper_bound = np.max(y_data)
+            self.computed_on = pd.to_datetime(datetime.datetime.now())
+
+    def compute(self, sample_areas):
+        return self.__linear_function(x_values=sample_areas, a = self.slope, b = self.intercept)
+
+
+class CalibrationMetadata():
+
+    def __init__(self):
+
+        self.standard = None
+        self.unit = None
+        self.prepared = None
+        self.measured = None
+
+class Calibration():
+
+    def __init__(self, cal_file: Path):
+
+        self.metadata = CalibrationMetadata()
+        self.curve = CalibrationCurve()
         self.cal_file = cal_file
-        self.date = None
+        self.calibrant_data = None
+        self.__results_index = None
+        self.excluded_points = None
 
         self.__read_calibration_file()
 
-    def calibrate(self):
-        pass
-
     def __read_calibration_file(self):
 
-        concentrations = []
-        dilutions =      []
-        A1s =            []
-        A2s =            []
-        A3s =            []
+        conctr = []
+        diltns = []
+        area_1 = []
+        area_2 = []
+        area_3 = []
 
-        unit = None
-
-        calibration_results_flag = False
+        results_flag = False
 
         with open(self.cal_file, 'r', errors = 'ignore') as file:
 
@@ -56,43 +119,133 @@ class TC_Calibration():
 
             for index, line in enumerate(file.readlines()):
 
-                if index == 8:
-                    unit = line.split()[0]
+                # Reading in the metadata
+                if index == 2:
+                    self.metadata.standard = float(line.split()[1])
+                if index == 3:
+                    self.metadata.unit = line.split()[1]
+                if index == 4:
+                    self.metadata.prepared = line.split()[1]
+                if index == 5:
+                    self.metadata.measured = line.split()[1]
 
-                if index >= 9 and calibration_results_flag == False:
+                # Reading in the calibration data
+                if index >= 10 and results_flag == False:
 
                     if line.startswith('CALIBRATION RESULTS'):
 
-                        calibration_results_flag = True
+                        results_flag = True
                         continue
 
-                    else:
+                    elif not results_flag:
+
+                        if line == '\n':
+                            continue
 
                         items = line.split()
 
-                        concentrations.append(float(items[0]))
-                        dilutions.append(float(items[1]))
-                        A1s.append(float(items[2]))
-                        A2s.append(float(items[3]))
-                        A3s.append(float(items[4]))
+                        conctr.append(float(items[0]))
+                        diltns.append(float(items[1]))
+                        area_1.append(float(items[2]))
+                        area_2.append(float(items[3]))
+                        area_3.append(float(items[4]))
 
-                if calibration_results_flag == True:
+                if results_flag == True:
 
                     if line.startswith('-'):
 
                         if results_index is None:
-                            results_index = index + 1
+                            self.__results_index = index + 1
+                            continue
                         else:
                             continue
 
-                        continue
-
                     items = line.split()
 
-                    if index == results_index:
-                        self.slope = float(items[1])
-                    elif index == results_index + 1:
-                        self.intercept = float(items[2])
-                    elif index == results_index + 2:
-                        self.R2_value = float(items[3])
-                    
+                    if index == self.__results_index:
+                        self.curve.exists = True
+                        self.curve.slope = float(items[1])
+                    elif index == self.__results_index + 1:
+                        self.curve.intercept = float(items[1])
+                    elif index == self.__results_index + 2:
+                        self.curve.R2_value = float(items[1])
+                    elif index == self.__results_index + 3:
+                        if items[1] == 'None':
+                            continue
+                        else:
+                            del items[0]
+                            self.curve.excluded = [int(items[i].replace(',', '')) for i in range(len(items))]
+                    elif index == self.__results_index + 4:
+                        self.curve.lower_bound = float(items[1])
+                    elif index == self.__results_index + 5:
+                        self.curve.upper_bound = float(items[1])
+                    elif index == self.__results_index + 6:
+                        self.curve.computed_on = pd.to_datetime(items[1] + ' ' + items[2])
+
+        self.calibrant_data = pd.DataFrame({'concentration': conctr,
+                                            'dilution': diltns, 
+                                            'area_1': area_1, 
+                                            'area_2': area_2, 
+                                            'area_3': area_3,
+                                            'area_average': [(a1 + a2 + a3) / 3 for a1, a2, a3 in zip(area_1, area_2, area_3)]})
+    def produce_statistics():
+        # TODO
+        pass
+
+    def calibrate(self, excluded_points = [], write_to_new_file = False):
+
+        if excluded_points:
+            self.curve.excluded = excluded_points
+            mask = [False if i + 1 in self.curve.excluded else True for i in range(len(self.calibrant_data))]
+        else:
+            self.curve.excluded = None
+            mask = [True for i in range(len(self.calibrant_data))]
+
+        self.curve.calibrate(x_data=self.calibrant_data['area_average'].values[mask], y_data=self.calibrant_data['concentration'].values[mask])
+
+        with open(self.cal_file, 'r', errors = 'ignore') as file:
+            lines = file.readlines()
+
+        if self.__results_index is None:
+
+            lines.append('\nCALIBRATION RESULTS\n')
+            lines.append('----------------------------------------------------------------------\n')
+            lines.append('Slope\t\t' + f"{self.curve.slope:2.4f}\n")
+            lines.append('Intercept\t' + f"{self.curve.intercept:2.4f}\n")
+            lines.append('R2\t\t' + f"{self.curve.R2_value:1.4f}\n")
+            if self.curve.excluded is None:
+                lines.append('EXCLUDED\t' + f"None\n")
+            else:
+                lines.append('EXCLUDED\t' + ', '.join([str(n) for n in self.curve.excluded]) + '\n')
+            lines.append('FROM\t\t' + f"{self.curve.lower_bound:2.4f}\n")
+            lines.append('TO\t\t' + f"{self.curve.upper_bound:2.4f}\n")
+            lines.append('COMPUTED\t' + f"{self.curve.computed_on:%d.%m.%Y %H:%M:%S}")
+
+        else:
+            lines[self.__results_index] = 'Slope\t\t' + f"{self.curve.slope:2.4f}\n"
+            lines[self.__results_index + 1] ='Intercept\t' + f"{self.curve.intercept:2.4f}\n"
+            lines[self.__results_index + 2] ='R2\t\t' + f"{self.curve.R2_value:1.4f}\n"
+            if self.curve.excluded is None:
+                lines[self.__results_index + 3] = 'EXCLUDED\t' + f"None\n"
+            else:
+                lines[self.__results_index + 3] = 'EXCLUDED\t' + ', '.join([str(n) for n in self.curve.excluded]) + '\n'
+            lines[self.__results_index + 4] ='FROM\t\t' + f"{self.curve.lower_bound:2.4f}\n"
+            lines[self.__results_index + 5] ='TO\t\t' + f"{self.curve.upper_bound:2.4f}\n"
+            lines[self.__results_index + 6] ='COMPUTED\t' + f"{self.curve.computed_on:%d.%m.%Y %H:%M:%S}"
+
+        if write_to_new_file == False:
+
+            with open(self.cal_file, 'w', errors = 'ignore') as file:
+                file.writelines(lines)
+
+        elif write_to_new_file == True:
+
+            directory_path = self.cal_file.parents[0]
+            new_file = self.cal_file.stem + '_' + self.curve.computed_on.strftime('%d_%m_%Y_%H_%M_%S') + '.cal'
+            new_path = PurePath(directory_path, new_file)
+
+            with open(new_path, 'w', errors = 'ignore') as file:
+                            file.writelines(lines)
+
+    def compute(self, sample_areas):
+        return self.curve.compute(sample_areas=sample_areas)

main.py

@@ -1,22 +1,7 @@
 # Importing custom code
-from defs.SMPS import OFRAnalyser, OFRSample
-from defs.OzoneMonitor import OzoneAnalyser
-from defs.Helper import DATFile
-
-# Importing required libraries
-import numpy as np
+from defs.TOC import Calibration
 from pathlib import Path
-import pandas as pd
 
-# Running some test
-sample = OFRSample(source_file=Path(r"Z:\Data\SMPS\processed\2026\09_September\20260902\20260902_aPinene_SOA.txt"), sample_number=20)
-print(sample.Mass.integrate())
-
-# Specify the source file paths to for SMPS data for the different filters used
-smps_analyser = OFRAnalyser(data_source_file = Path(r"Z:\Data\SMPS\processed\2026\09_September\20260903\20260903_aPinene_SOA.txt"),
-                            flow_source_file = Path(r"Z:\People\David_Tschan\OFR_SMPS\flow_files\20260903_SET02_EXP14_flow_file.txt"),
-                            smpl_number_list = [20],
-                            on_time_interval = [pd.Timestamp('2026-09-03 15:20:00'), pd.Timestamp('2026-09-03 16:20:00')])
-
-smps_analyser.GetData()
-print(smps_analyser.data)
+calibration = Calibration(cal_file = Path(r"C:\Users\tschan_d\Documents\analysis\TOC_Analysis\20260909\IC_Calibration.cal"))
+calibration.calibrate(write_to_new_file=True)
+print(calibration.compute([1, 2, 3, 4, 20]))