# Author: David Tschan # Institute: Paul-Scherrer Institute, Laboratory for Atmospheric Chemistry, Multiphase Chemistry Group, # Forschungsstrasse 111, 5232 Villigen PSI, Switzerland # Role: Master Student # License: n/a # Date: 2026 ### ----------------------------------------------------------------------------------------------------------------------------------- ### ### SUBMODULE DOCSTRING """This submodule contains functionality to calibrate TOC measurement data and analyse them. """ ### ----------------------------------------------------------------------------------------------------------------------------------- ### ### IMPORTS AND CONFIGURATIONS import numpy as np import pandas as pd from pathlib import Path, PurePath import datetime import warnings from scipy.optimize import curve_fit from defs.Helper import compute_R2 import re ### ----------------------------------------------------------------------------------------------------------------------------------- ### ### TOC.py class CalibrationCurve(): r"""An object representing a calibration curve for TOC, IC or NPOC Measurements. Description ----------- This object instanciates a class that relates the area $A$ measured in the TOC Analyser to the concentration $c$ used in the calibration solutions. The calibration curve is affine, $c=aA+b$ where $a$ and $b$ denote slope and intercept, respectively, of the affine map. This object can be used to represent linear calibration curves for TC (total carbon)-, IC (inorganic carbon)- and NPOC (non-purgeable organic carbon) measurements. Attributes ---------- exists: Bool This value is set to `True` if a calibration curve is already available, else it is set to `False`. slope: float Value of the slope for the calibration of concentration (dependent variable) with measured area (independent variable). intercept: float Value of the intercep for the calibration of concentration (dependent variable) with measured area (independent variable). R2_value: float The coefficient of determination of the calibration curve. For a good calibration curve, this should be above 0.995.\ excluded: list A list of integers denoting the samples in the calibrant data that are excluded. lower_bound: float The lower bound value for which this calibration curve can be used. upper_bound: float The upper bound value for which this calibration curve can be used. computed_on: pd.Timestamp A timestamp indicating the time at which the calibration curve was computed. """ def __init__(self): # initialises all member variables 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): r"""Computes the value of the calibration curve for the passed parameters. Description ----------- Computes the affine map, $c=aA+b$, given data $A$ and parameters $a$ (slope) and $b$ (intercept). Parameters ---------- x_values: array-like An array-like (array or list) object of length $n$ containing the data for the independent variable that are to be mapped. a: float The slope of the calibration curve. b: float The intercept of the calibration curve. Returns ------- y_values: np.ndarray An array (of size $n$) containing the image under the affine map of `x_values`. """ # check if the input data is an instance of a np.ndarray. If not, cast it to np.ndarray if not isinstance(x_values, np.ndarray): x_values = np.array(x_values) # evaluate the image of the input under the affine map and return the results return a * x_values + b def __optimise(self, x_data, y_data, p0): r"""Finds the optimal parameters of the calibration curve Description ----------- Uses `scipy.optimize.curve_fit` to find optimal parameters for the calibrant data. Parameters ---------- x_data: array-like The data representing the independent variable. y_data: array-like The data representing the dependent variable. p0: list A list containing the initial guesses of the parameters. Returns ------- p: np.ndarray An array containing the optimal parameter values. pcov: The covariance matrix associated with the optimal parameters. """ # evaluate the optimisation problem and return optimal parameters and their covariance matrix. return curve_fit(self.__linear_function, x_data, y_data, p0) def __check_exist(self): r"""Checks if attribute `exists` is set Description ----------- Checks of the attribute `exists` is set to `None` (its instanciation value) and returns `False` in this case, or returns `True` if attribute value is not `None` (i.e. if it has been set). Returns ------- t: Bool `True` if attribute has a value, else returns `None`. """ # this line is equivalent to "if None is True" if attribute is not set, which evaluates to 'False'. Else evaluates to 'True'. if self.exists == True: return True else: return False def __check(self): r"""Gives a warning message when overwriting an existing calibration curve. Description ----------- Calls the private `check_exists` method to see if a calibration curve has already been computed. If yes, this will throw a warning message, but not an error. """ if self.__check_exist(): warnings.warn("Calibration alreadyt exists, will be overwritten...") else: return def calibrate(self, x_data, y_data): """Compute a calibration curve based on the available calibrant data. Description ----------- Bla Bla Bla Parameters ---------- x_data: np.ndarray Bla bla bla y_data: np.ndarray Bla Bla Bla """ 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.calibrant_data = None self.__results_index = None self.excluded_points = None self.__read_calibration_file() self.mask = None def __read_calibration_file(self): conctr = [] diltns = [] area_1 = [] area_2 = [] area_3 = [] results_flag = False with open(self.cal_file, 'r', errors = 'ignore') as file: results_index = None for index, line in enumerate(file.readlines()): # 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] # Reading in the calibration data if index >= 10 and results_flag == False: if line.startswith('CALIBRATION RESULTS'): results_flag = True continue elif not results_flag: if line == '\n': continue items = line.split() 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 results_flag == True: if line.startswith('-'): if results_index is None: self.__results_index = index + 1 continue else: continue items = line.split() 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 self.mask = [False if i + 1 in self.curve.excluded else True for i in range(len(self.calibrant_data))] else: self.curve.excluded = None self.mask = [True for i in range(len(self.calibrant_data))] self.curve.calibrate(x_data=self.calibrant_data['area_average'].values[self.mask], y_data=self.calibrant_data['concentration'].values[self.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) def TXTFileTransformer(file: Path, headline = 11, paketsize = 3, to_clipboard = True): methods = [] # column 1 sample_names = [] # column 2 comments = [] # column 6 vials = [] # column 21 date_times_1 = [] # column 7 date_times_2 = [] # column 7 date_times_3 = [] # column 7 areas_1 = [] # column 11 areas_2 = [] # column 11 areas_3 = [] # column 11 concentrations_1 = [] # column 13 concentrations_2 = [] # column 13 concentrations_3 = [] # column 13 with open(file = file, errors='ignore', mode = 'r') as f: for index, line in enumerate(f.readlines()): if index < headline: continue else: items = re.split(r'\t', line)[:-1] # first sample if (index - headline) % paketsize == 0: methods.append(str(items[1])) sample_names.append(str(items[2])) comments.append(str(items[6])) vials.append(int(items[21])) date_times_1.append(pd.to_datetime(str(items[7]), format=r'%m/%d/%Y %I:%M:%S %p')) areas_1.append(float(items[11])) concentrations_1.append(float(items[13])) # second sample elif (index - headline) % paketsize == 1: date_times_2.append(pd.to_datetime(str(items[7]), format=r'%m/%d/%Y %I:%M:%S %p')) areas_2.append(float(items[11])) concentrations_2.append(float(items[13])) # third sample elif (index - headline) % paketsize == 2: date_times_3.append(pd.to_datetime(str(items[7]), format=r'%m/%d/%Y %I:%M:%S %p')) areas_3.append(float(items[11])) concentrations_3.append(float(items[13])) out_frame = pd.DataFrame({'method': methods, 'sample': sample_names, 'comment': comments, 'vial': vials, 'datetime_1': date_times_1, 'datetime_2': date_times_2, 'datetime_3': date_times_3, 'area_1': areas_1, 'area_2': areas_2, 'area_3': areas_3, 'concentration_1': concentrations_1, 'concentration_2': concentrations_2, 'concentration_3': concentrations_3}) if to_clipboard: out_frame.to_clipboard(excel=True, index=False) else: return out_frame