Coverage for src/nmrlineshapeanalyser/core.py: 81%

1import nmrglue as ng 

2import numpy as np 

3import scipy 

4from scipy.optimize import curve_fit 

5import matplotlib.pyplot as plt 

6import matplotlib as mpl 

7from typing import List, Tuple, Dict, Optional, Union 

8import warnings 

9import pandas as pd 

10import sys 

11import glob 

12import os 

13 

14print(f"nmrglue: {ng.__version__}") 

15print(f"numpy: {np.__version__}") 

16print(f"scipy: {scipy.__version__}") 

17print(f"matplotlib: {mpl.__version__}") 

18print(f"pandas: {pd.__version__}") 

19 

20class NMRProcessor: 

21 """ 

22 A comprehensive class for processing and analyzing NMR data. 

23 Combines data loading, processing, peak fitting, and visualization. 

24 """ 

25 

26 def __init__(self): 

27 """Initialize the NMR processor with default plot style.""" 

28 self.data = None 

29 self.number = None 

30 self.nucleus = None 

31 self.carrier_freq = None 

32 self.uc = None 

33 self.ppm = None 

34 self.ppm_limits = None 

35 self.fixed_params = None 

36 self.set_plot_style() 

37 

38 @staticmethod 

39 def set_plot_style() -> None: 

40 """Set up the matplotlib plotting style.""" 

41 mpl.rcParams['font.family'] = "sans-serif" 

42 plt.rcParams['font.sans-serif'] = ['Arial'] 

43 plt.rcParams['font.size'] = 14 

44 plt.rcParams['axes.linewidth'] = 2 

45 mpl.rcParams['xtick.major.size'] = mpl.rcParams['ytick.major.size'] = 8 

46 mpl.rcParams['xtick.major.width'] = mpl.rcParams['ytick.major.width'] = 1 

47 mpl.rcParams['xtick.direction'] = mpl.rcParams['ytick.direction'] = 'out' 

48 mpl.rcParams['xtick.major.top'] = mpl.rcParams['ytick.major.right'] = False 

49 mpl.rcParams['xtick.minor.size'] = mpl.rcParams['ytick.minor.size'] = 5 

50 mpl.rcParams['xtick.minor.width'] = mpl.rcParams['ytick.minor.width'] = 1 

51 mpl.rcParams['xtick.top'] = mpl.rcParams['ytick.right'] = True 

52 

53 def load_data(self, filepath: str) -> None: 

54 """ 

55 Load and process Bruker NMR data from the specified filepath. 

56  

57 Args: 

58 filepath (str): Path to the Bruker data directory 

59 """ 

60 # Read the Bruker data 

61 dic, self.data = ng.bruker.read_pdata(filepath) 

62 

63 # Set the spectral parameters 

64 udic = ng.bruker.guess_udic(dic, self.data) 

65 

66 nuclei = udic[0]['label'] 

67 carrier_freq = udic[0]['obs'] 

68 self.carrier_freq = carrier_freq 

69 

70 # Extract number and nucleus symbols 

71 self.number = ''.join(filter(str.isdigit, nuclei)) 

72 self.nucleus = ''.join(filter(str.isalpha, nuclei)) 

73 

74 # Create converter and get scales 

75 self.uc = ng.fileiobase.uc_from_udic(udic, dim=0) 

76 self.ppm = self.uc.ppm_scale() 

77 self.ppm_limits = self.uc.ppm_limits() 

78 

79 def load_csv(self, filepath: str, atomic_no: str, nucleus: str, larmor_freq: float) -> None: 

80 """ 

81 Load CSV data exported from MNova (columns: 'ppm', 'Intensity'). 

82  

83 Args: 

84 filepath (str): Directory path containing the CSV file 

85 atomic_no (str): Atomic number of the nucleus 

86 nucleus (str): Nuclear symbol 

87 larmor_freq (float): Larmor frequency in MHz 

88  

89 Raises: 

90 FileNotFoundError: If no CSV files are found in the specified directory 

91 """ 

92 csv_files = glob.glob(os.path.join(filepath, '*.csv')) 

93 

94 if not csv_files: 

95 raise FileNotFoundError(f"No CSV files found in {filepath}") 

96 

97 data = pd.read_csv(csv_files[0], sep='[,\t]', engine='python') 

98 

99 required_columns = ['ppm', 'Intensity'] 

100 if not all(col in data.columns for col in required_columns): 

101 raise ValueError(f"CSV must contain columns: {required_columns}") 

102 

103 x_data = data['ppm'].values 

104 y_data = data['Intensity'].values 

105 

106 # Create pseudo-complex data to match the Bruker data format 

107 y_data = y_data + 1j * np.zeros_like(y_data) 

108 

109 self.ppm = x_data 

110 self.data = y_data 

111 self.number = str(atomic_no) 

112 self.nucleus = nucleus 

113 self.carrier_freq = float(larmor_freq) 

114 

115 def select_region(self, ppm_start: float, ppm_end: float) -> Tuple[np.ndarray, np.ndarray]: 

116 """Select a specific region of the NMR spectrum for analysis.""" 

117 if self.data is None: 

118 raise ValueError("No data loaded. Call load_data first.") 

119 

120 region_mask = (self.ppm >= ppm_start) & (self.ppm <= ppm_end) 

121 x_region = self.ppm[region_mask] 

122 y_real = self.data.real 

123 y_region = y_real[region_mask] 

124 

125 if x_region.size == 0: 

126 raise ValueError(f"No data found in region {ppm_start} to {ppm_end} ppm.") 

127 

128 return x_region, y_region 

129 

130 def normalize_data(self, x_data: np.ndarray, y_data: np.ndarray) -> Tuple[np.ndarray, np.ndarray, float, float]: 

131 """Normalize the data for processing to 0-1 range.""" 

132 y_ground = np.min(y_data) 

133 y_normalized = y_data - y_ground 

134 y_amp = np.max(y_normalized) 

135 y_normalized = y_normalized / y_amp if y_amp != 0 else np.zeros_like(y_normalized, dtype=float) 

136 return x_data, y_normalized, y_amp, y_ground 

137 

138 @staticmethod 

139 def pseudo_voigt(x: np.ndarray, x0: float, amp: float, width: float, eta: float) -> np.ndarray: 

140 """Calculate the Pseudo-Voigt function.""" 

141 sigma = width / (2 * np.sqrt(2 * np.log(2))) 

142 gamma = width / 2 

143 lorentzian = amp * (gamma**2 / ((x - x0)**2 + gamma**2)) 

144 gaussian = amp * np.exp(-0.5 * ((x - x0) / sigma)**2) 

145 return eta * lorentzian + (1 - eta) * gaussian 

146 

147 def pseudo_voigt_multiple(self, x: np.ndarray, *params) -> np.ndarray: 

148 """ 

149 Calculate multiple Pseudo-Voigt peaks summed on top of a single shared offset. 

150 The shared offset is always the last entry in params. 

151 """ 

152 n_peaks = len(self.fixed_params) 

153 

154 expected_len = sum( 

155 (fixed_x0 is None) + (fixed_amp is None) + (fixed_width is None) + 1 

156 for fixed_x0, fixed_amp, fixed_width, _, _ in self.fixed_params 

157 ) + 1 # +1 for the shared offset 

158 if len(params) != expected_len: 

159 raise ValueError( 

160 f"Expected {expected_len} parameters for {n_peaks} peak(s), got {len(params)}" 

161 ) 

162 

163 param_idx = 0 

164 y = np.zeros_like(x) 

165 

166 for i in range(n_peaks): 

167 fixed_x0, fixed_amp, fixed_width, _, _ = self.fixed_params[i] 

168 

169 if fixed_x0 is not None: 

170 x0 = fixed_x0 

171 else: 

172 x0 = params[param_idx] 

173 param_idx += 1 

174 

175 if fixed_amp is not None: 

176 amp = fixed_amp 

177 else: 

178 amp = params[param_idx] 

179 param_idx += 1 

180 

181 if fixed_width is not None: 

182 width = fixed_width 

183 else: 

184 width = params[param_idx] 

185 param_idx += 1 

186 

187 eta = params[param_idx] 

188 param_idx += 1 

189 

190 y += self.pseudo_voigt(x, x0, amp, width, eta) 

191 

192 # Single shared offset applied once to the summed peaks 

193 offset = params[param_idx] 

194 y += offset 

195 

196 return y 

197 

198 def fit_peaks(self, x_data: np.ndarray, y_data: np.ndarray, 

199 initial_params: List[float], fixed_x0: Optional[List[bool]] = None, 

200 fixed_amp: Optional[List[bool]] = None, fixed_width: Optional[List[bool]] = None, 

201 y_scale: float = 1.0, y_offset: float = 0.0) -> Tuple[np.ndarray, List[Dict], np.ndarray]: 

202 """ 

203 Fit multiple Pseudo-Voigt peaks to the data, sharing a single baseline offset. 

204  

205 Args: 

206 fixed_x0: per-peak flags to fix the peak position at its initial value 

207 fixed_amp: per-peak flags to fix the peak amplitude at its initial value 

208 fixed_width: per-peak flags to fix the peak width at its initial value 

209  

210 Note: all peaks share one fitted offset (baseline), seeded from the mean of 

211 the per-peak offsets given in initial_params. The returned popt/peak_metrics 

212 report that same shared offset value for every peak. 

213  

214 IMPORTANT: Initial offsets should be in normalized scale (0-1), e.g., use 0.0. 

215 """ 

216 if len(initial_params) % 5 != 0: 

217 raise ValueError("Number of initial parameters must be divisible by 5") 

218 

219 n_peaks = len(initial_params) // 5 

220 

221 if fixed_x0 is None: 

222 fixed_x0 = [False] * n_peaks 

223 if fixed_amp is None: 

224 fixed_amp = [False] * n_peaks 

225 if fixed_width is None: 

226 fixed_width = [False] * n_peaks 

227 

228 self.fixed_params = [] 

229 fit_params = [] 

230 lower_bounds = [] 

231 upper_bounds = [] 

232 

233 # Process each peak's shape parameters (x0, amp, width, eta). 

234 # The offset is handled separately below as a single shared parameter. 

235 for i in range(n_peaks): 

236 x0, amp, width, eta, offset = initial_params[5*i:5*(i+1)] 

237 

238 self.fixed_params.append(( 

239 x0 if fixed_x0[i] else None, 

240 amp if fixed_amp[i] else None, 

241 width if fixed_width[i] else None, 

242 None, 

243 None, 

244 )) 

245 

246 if not fixed_x0[i]: 

247 fit_params.append(x0) 

248 lower_bounds.append(x0 - width/2) 

249 upper_bounds.append(x0 + width/2) 

250 

251 if not fixed_amp[i]: 

252 fit_params.append(amp) 

253 lower_bounds.append(0) 

254 upper_bounds.append(np.inf) 

255 

256 if not fixed_width[i]: 

257 fit_params.append(width) 

258 lower_bounds.append(1) 

259 upper_bounds.append(np.inf) 

260 

261 # eta is always free 

262 fit_params.append(eta) 

263 lower_bounds.append(0) 

264 upper_bounds.append(1) 

265 

266 # Single shared offset, appended once at the end, seeded from the 

267 # mean of the per-peak offsets in initial_params 

268 shared_offset_init = float(np.mean(initial_params[4::5])) 

269 fit_params.append(shared_offset_init) 

270 lower_bounds.append(-1.0) 

271 upper_bounds.append(2.0) 

272 

273 # Perform the fit 

274 with warnings.catch_warnings(): 

275 warnings.filterwarnings('ignore', category=RuntimeWarning) 

276 popt, pcov = curve_fit(self.pseudo_voigt_multiple, x_data, y_data, 

277 p0=fit_params, bounds=(lower_bounds, upper_bounds), 

278 maxfev=10000, method='trf') 

279 

280 # Process results 

281 full_popt = self._process_fit_results(popt, initial_params, fixed_x0, fixed_amp, fixed_width) 

282 peak_metrics = self.calculate_peak_metrics(full_popt, pcov, fixed_x0, fixed_amp, fixed_width) 

283 fitted_data = self.pseudo_voigt_multiple(x_data, *popt) 

284 

285 return full_popt, peak_metrics, fitted_data 

286 

287 def _process_fit_results(self, popt: np.ndarray, initial_params: List[float], 

288 fixed_x0: List[bool], fixed_amp: List[bool], 

289 fixed_width: List[bool]) -> np.ndarray: 

290 """Process and organize fitting results, broadcasting the shared offset to every peak.""" 

291 n_peaks = len(initial_params) // 5 

292 param_idx = 0 

293 peak_shape_params = [] 

294 

295 for i in range(n_peaks): 

296 x0_init, amp_init, width_init, eta_init, offset_init = initial_params[5*i:5*(i+1)] 

297 

298 if fixed_x0[i]: 

299 x0 = x0_init 

300 else: 

301 x0 = popt[param_idx] 

302 param_idx += 1 

303 

304 if fixed_amp[i]: 

305 amp = amp_init 

306 else: 

307 amp = popt[param_idx] 

308 param_idx += 1 

309 

310 if fixed_width[i]: 

311 width = width_init 

312 else: 

313 width = popt[param_idx] 

314 param_idx += 1 

315 

316 eta = popt[param_idx] 

317 param_idx += 1 

318 

319 peak_shape_params.append((x0, amp, width, eta)) 

320 

321 # Last remaining parameter is the single shared offset 

322 shared_offset = popt[param_idx] 

323 

324 full_popt = [] 

325 for x0, amp, width, eta in peak_shape_params: 

326 full_popt.extend([x0, amp, width, eta, shared_offset]) 

327 

328 return np.array(full_popt) 

329 

330 def calculate_peak_metrics(self, popt: np.ndarray, pcov: np.ndarray, 

331 fixed_x0: List[bool], fixed_amp: List[bool], 

332 fixed_width: List[bool]) -> List[Dict]: 

333 """Calculate metrics for each fitted peak, using the shared offset's error for all peaks.""" 

334 n_peaks = len(popt) // 5 

335 errors = np.sqrt(np.diag(pcov)) if pcov.size else np.zeros(len(popt)) 

336 error_idx = 0 

337 

338 peak_shape_errors = [] 

339 for i in range(n_peaks): 

340 if fixed_x0[i]: 

341 x0_err = 0 

342 else: 

343 x0_err = errors[error_idx] 

344 error_idx += 1 

345 

346 if fixed_amp[i]: 

347 amp_err = 0 

348 else: 

349 amp_err = errors[error_idx] 

350 error_idx += 1 

351 

352 if fixed_width[i]: 

353 width_err = 0 

354 else: 

355 width_err = errors[error_idx] 

356 error_idx += 1 

357 

358 eta_err = errors[error_idx] 

359 error_idx += 1 

360 

361 peak_shape_errors.append((x0_err, amp_err, width_err, eta_err)) 

362 

363 # Last remaining error entry belongs to the single shared offset 

364 offset_err = errors[error_idx] if error_idx < len(errors) else 0 

365 

366 peak_results = [] 

367 for i in range(n_peaks): 

368 x0, amp, width, eta, offset = popt[5*i:5*(i+1)] 

369 x0_err, amp_err, width_err, eta_err = peak_shape_errors[i] 

370 

371 sigma = width / (2 * np.sqrt(2 * np.log(2))) 

372 gamma = width / 2 

373 

374 gauss_area = (1 - eta) * amp * sigma * np.sqrt(2 * np.pi) 

375 lorentz_area = eta * amp * np.pi * gamma 

376 total_area = gauss_area + lorentz_area 

377 

378 gauss_area_err = np.sqrt( 

379 ((1 - eta) * sigma * np.sqrt(2 * np.pi) * amp_err) ** 2 + 

380 (amp * sigma * np.sqrt(2 * np.pi) * eta_err) ** 2 + 

381 ((1 - eta) * amp * np.sqrt(2 * np.pi) * (width_err / (2 * np.sqrt(2 * np.log(2))))) ** 2 

382 ) 

383 

384 lorentz_area_err = np.sqrt( 

385 (eta * np.pi * gamma * amp_err) ** 2 + 

386 (amp * np.pi * gamma * eta_err) ** 2 + 

387 (eta * amp * np.pi * (width_err / 2)) ** 2 

388 ) 

389 

390 total_area_err = np.sqrt(gauss_area_err ** 2 + lorentz_area_err ** 2) 

391 

392 peak_results.append({ 

393 'x0': (x0, x0_err), 

394 'amplitude': (amp, amp_err), 

395 'width': (width, width_err), 

396 'eta': (eta, eta_err), 

397 'offset': (offset, offset_err), 

398 'gaussian_area': (gauss_area, gauss_area_err), 

399 'lorentzian_area': (lorentz_area, lorentz_area_err), 

400 'total_area': (total_area, total_area_err) 

401 }) 

402 

403 return peak_results 

404 

405 def plot_results(self, x_data: np.ndarray, y_data: np.ndarray, 

406 fitted_data: np.ndarray, 

407 popt: np.ndarray) -> Tuple[plt.Figure, plt.Axes, List[np.ndarray]]: 

408 """Plot the fitting results with components (all in normalized 0-1 scale).""" 

409 fig, ax1 = plt.subplots(1, 1, figsize=(12, 6)) 

410 

411 ax1.plot(x_data, y_data, 'ok', ms=1, label='Data') 

412 ax1.plot(x_data, fitted_data, '-r', lw=2, label='Fit') 

413 residuals = y_data - fitted_data 

414 ax1.plot(x_data, residuals-0.05, '-g', lw=2, label='Residuals', alpha=0.5) 

415 

416 n_peaks = len(popt) // 5 

417 components = [] 

418 

419 for i in range(n_peaks): 

420 x0, amp, width, eta, offset = popt[5*i:5*(i+1)] 

421 component = self.pseudo_voigt(x_data, x0, amp, width, eta) + offset 

422 components.append(component) 

423 

424 ax1.fill(x_data, component, alpha=0.5, label=f'Component {i+1}') 

425 peak_height = self.pseudo_voigt(np.array([x0]), x0, amp, width, eta)[0] + offset 

426 ax1.plot(x0, peak_height, 'ob', markersize=8, label='Peak Position' if i == 0 else '') 

427 

428 ax1.invert_xaxis() 

429 ax1.legend(ncol=2, fontsize=10) 

430 ax1.set_title('NMR Peak Fit (Normalized Scale 0-1)') 

431 ax1.set_xlabel(f'$^{{{self.number}}} \\ {self.nucleus}$ chemical shift (ppm)') 

432 ax1.set_ylabel('Intensity (normalized 0-1)') 

433 ax1.hlines(0, x_data[0], x_data[-1], colors='blue', linestyles='dashed', alpha=0.5) 

434 ax1.set_ylim(-0.15, 1.15) 

435 

436 return fig, ax1, components 

437 

438 def _print_detailed_results(self, peak_metrics: List[Dict]) -> None: 

439 """Print detailed fitting results and statistics.""" 

440 print("\nPeak Fitting Results:") 

441 print("===================") 

442 

443 area_of_peaks = [] 

444 for i, metrics in enumerate(peak_metrics, 1): 

445 print(f"\nPeak {i} (Position: {metrics['x0'][0]:.2f} ± {metrics['x0'][1]:.2f}):") 

446 print(f" Amplitude: {metrics['amplitude'][0]:.3f} ± {metrics['amplitude'][1]:.3f}") 

447 print(f" Width (FWHM): {metrics['width'][0]:.2f} ± {metrics['width'][1]:.2f} ppm") 

448 print(f" Width (FWHM): {metrics['width'][0]*self.carrier_freq:.2f} ± {metrics['width'][1]*self.carrier_freq:.2f} Hz") 

449 print(f" Carrier Frequency: {self.carrier_freq} MHz") 

450 print(f" Eta: {metrics['eta'][0]:.2f} ± {metrics['eta'][1]:.2f}") 

451 print(f" Offset: {metrics['offset'][0]:.4f} ± {metrics['offset'][1]:.4f}") 

452 print(f" Gaussian Area: {metrics['gaussian_area'][0]:.2f} ± {metrics['gaussian_area'][1]:.2f}") 

453 print(f" Lorentzian Area: {metrics['lorentzian_area'][0]:.2f} ± {metrics['lorentzian_area'][1]:.2f}") 

454 print(f" Total Area: {metrics['total_area'][0]:.2f} ± {metrics['total_area'][1]:.2f}") 

455 area_of_peaks.append(metrics['total_area']) 

456 

457 self._calculate_and_print_percentages(area_of_peaks) 

458 

459 def _calculate_and_print_percentages(self, area_of_peaks: List[Tuple[float, float]]) -> None: 

460 """Calculate and print percentage contributions of each peak.""" 

461 total_area_sum = sum(area[0] for area in area_of_peaks) 

462 total_area_sum_err = np.sqrt(sum(area[1]**2 for area in area_of_peaks)) 

463 

464 for i, (area, area_err) in enumerate(area_of_peaks, 1): 

465 percentage = (area / total_area_sum) * 100 

466 percentage_err = percentage * np.sqrt((area_err / area) ** 2 + 

467 (total_area_sum_err / total_area_sum) ** 2) 

468 print(f'Peak {i} Percentage: {percentage:.2f}% ± {percentage_err:.2f}%') 

469 

470 overall_percentage = sum((area[0] / total_area_sum) * 100 for area in area_of_peaks) 

471 print(f'Total: {overall_percentage:.2f}%') 

472 

473 def save_results(self, filepath: str, x_data: np.ndarray, y_data: np.ndarray, 

474 fitted_data: np.ndarray, peak_metrics: List[Dict], 

475 popt: np.ndarray, components: List[np.ndarray]) -> None: 

476 """Save all results to files.""" 

477 self._save_peak_data(filepath, x_data, y_data, fitted_data, components) 

478 self._save_metrics(filepath, peak_metrics) 

479 self._save_plot(filepath, x_data, y_data, fitted_data, popt) 

480 self._print_detailed_results(peak_metrics) 

481 

482 def _save_peak_data(self, filepath: str, x_data: np.ndarray, y_data: np.ndarray, 

483 fitted_data: np.ndarray, components: List[np.ndarray]) -> None: 

484 """Save peak data to CSV file.""" 

485 df = pd.DataFrame({'x_data': x_data, 'y_data': y_data, 'y_fit': fitted_data}) 

486 for i, component in enumerate(components): 

487 df[f'component_{i+1}'] = component 

488 df.to_csv(filepath + 'peak_data.csv', index=False) 

489 

490 def _save_metrics(self, filepath: str, peak_metrics: List[Dict]) -> None: 

491 """Save peak metrics to text file.""" 

492 with open(filepath + 'pseudoVoigtPeak_metrics.txt', 'w') as file: 

493 area_of_peaks = [] 

494 for i, metrics in enumerate(peak_metrics, 1): 

495 file.write(f"\nPeak {i} (Position: {metrics['x0'][0]:.2f} ± {metrics['x0'][1]:.2f}):\n") 

496 file.write(f"Amplitude: {metrics['amplitude'][0]:.3f} ± {metrics['amplitude'][1]:.3f}\n") 

497 file.write(f"Width (FWHM): {metrics['width'][0]:.2f} ± {metrics['width'][1]:.2f} ppm\n") 

498 file.write(f"Width (FWHM): {metrics['width'][0]*self.carrier_freq:.2f} ± {metrics['width'][1]*self.carrier_freq:.2f} Hz\n") 

499 file.write(f"Eta: {metrics['eta'][0]:.2f} ± {metrics['eta'][1]:.2f}\n") 

500 file.write(f"Offset: {metrics['offset'][0]:.4f} ± {metrics['offset'][1]:.4f}\n") 

501 file.write(f"Total Area: {metrics['total_area'][0]:.2f} ± {metrics['total_area'][1]:.2f}\n\n") 

502 area_of_peaks.append(metrics['total_area']) 

503 

504 total_area_sum = sum(area[0] for area in area_of_peaks) 

505 for i, (area, area_err) in enumerate(area_of_peaks, 1): 

506 percentage = (area / total_area_sum) * 100 

507 file.write(f'Peak {i} Percentage: {percentage:.2f}%\n') 

508 

509 def _save_plot(self, filepath: str, x_data: np.ndarray, y_data: np.ndarray, 

510 fitted_data: np.ndarray, popt: np.ndarray) -> None: 

511 """Save the plot to a file.""" 

512 fig, _, _ = self.plot_results(x_data, y_data, fitted_data, popt) 

513 fig.savefig(filepath + 'pseudoVoigtPeakFit.png', bbox_inches='tight') 

514 plt.close(fig)