435 lines
17 KiB
Python
435 lines
17 KiB
Python
# 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 |