Skills › Scientific Agent Skills › nmrglue › process_1d.py
process_1d.py
190 lines · skills/nmrglue/scripts/process_1d.py at 5f2bc03
1#!/usr/bin/env python3
2"""Process a calibrated complex 1D FID; no acquisition settings are inferred."""
3from __future__ import annotations
5import argparse
6import csv
7import hashlib
8import json
9from pathlib import Path
11import nmrglue as ng
12import numpy as np
13import scipy
14from scipy.signal import find_peaks
17def finite_number(value, name, positive=False):
18 if isinstance(value, bool) or not isinstance(value, (int, float)):
19 raise ValueError(f"{name} must be a number")
20 number = float(value)
21 if not np.isfinite(number) or (positive and number <= 0):
22 raise ValueError(f"{name} must be finite" + (" and positive" if positive else ""))
23 return number
26def regions_in_axis(regions, ppm):
27 result = []
28 for region in regions:
29 if len(region) != 2:
30 raise ValueError("Each region needs two ppm bounds")
31 low, high = sorted(finite_number(v, "region bound") for v in region)
32 if low == high or low < ppm.min() or high > ppm.max():
33 raise ValueError("Region must have width and lie wholly inside the ppm axis")
34 result.append((low, high))
35 return result
38def process(fid, settings):
39 fid = np.asarray(fid)
40 if fid.ndim != 1 or len(fid) < 16 or not np.iscomplexobj(fid):
41 raise ValueError("fid must be a one-dimensional complex array with at least 16 points")
42 if not np.isfinite(fid).all():
43 raise ValueError("fid contains nonfinite values")
44 allowed = {"spectral_width_hz", "observation_mhz", "carrier_ppm", "nucleus", "fid_sign",
45 "line_broadening_hz", "zero_fill_points", "first_point_scale", "phase0_deg",
46 "phase1_deg", "baseline", "baseline_regions_ppm", "integration_regions_ppm",
47 "peak_prominence_fraction"}
48 if set(settings) - allowed:
49 raise ValueError(f"Unknown settings: {sorted(set(settings) - allowed)}")
50 required = allowed - {"baseline_regions_ppm"}
51 if required - set(settings):
52 raise ValueError(f"Missing explicit settings: {sorted(required - set(settings))}")
53 sw = finite_number(settings["spectral_width_hz"], "spectral_width_hz", True)
54 obs = finite_number(settings["observation_mhz"], "observation_mhz", True)
55 car = finite_number(settings["carrier_ppm"], "carrier_ppm")
56 if not isinstance(settings["nucleus"], str) or not settings["nucleus"].strip():
57 raise ValueError("nucleus must identify the observed nucleus")
58 if settings["fid_sign"] not in ("-i", "+i"):
59 raise ValueError("fid_sign must be '-i' or '+i'")
60 lb = finite_number(settings["line_broadening_hz"], "line_broadening_hz")
61 if lb < 0:
62 raise ValueError("Negative line broadening is unsupported")
63 size = settings["zero_fill_points"]
64 if isinstance(size, bool) or not isinstance(size, int) or size < len(fid) or size % 2:
65 raise ValueError("zero_fill_points must be an even integer at least the FID length")
66 first = finite_number(settings["first_point_scale"], "first_point_scale")
67 if not 0 <= first <= 1:
68 raise ValueError("first_point_scale must be between zero and one")
69 p0 = finite_number(settings["phase0_deg"], "phase0_deg")
70 p1 = finite_number(settings["phase1_deg"], "phase1_deg")
71 fraction = finite_number(settings["peak_prominence_fraction"], "peak_prominence_fraction", True)
72 if fraction > 1:
73 raise ValueError("peak_prominence_fraction must be at most one")
74 data = fid.astype(np.complex128, copy=True)
75 if settings["fid_sign"] == "+i":
76 data = data.conj()
77 data[0] *= first
78 data = ng.proc_base.em(data, lb=lb / sw)
79 spectrum = ng.proc_base.ps(ng.proc_base.fft(ng.proc_base.zf_size(data, size)), p0=p0, p1=p1)
80 ppm = ng.fileiobase.unit_conversion(size, False, sw, obs, car * obs).ppm_scale()
81 raw_real = spectrum.real.copy()
82 baseline = np.zeros(size)
83 regions = regions_in_axis(settings.get("baseline_regions_ppm", []), ppm)
84 if settings["baseline"] == "linear":
85 mask = np.zeros(size, dtype=bool)
86 for low, high in regions:
87 mask |= (ppm >= low) & (ppm <= high)
88 if mask.sum() < 4 or len(regions) < 2:
89 raise ValueError("Linear baseline needs at least two signal-free regions and four points")
90 baseline = np.polyval(np.polyfit(ppm[mask], raw_real[mask], 1), ppm)
91 elif settings["baseline"] != "none":
92 raise ValueError("baseline must be 'none' or 'linear'")
93 elif regions:
94 raise ValueError("baseline_regions_ppm requires baseline='linear'")
95 real = raw_real - baseline
96 prominence = fraction * max(float(real.max()), 0)
97 indices, properties = find_peaks(real, height=0, prominence=max(prominence, np.finfo(float).tiny))
98 peaks = [{"ppm": float(ppm[i]), "height": float(real[i]), "prominence": float(p)}
99 for i, p in zip(indices, properties["prominences"])]
100 integrals = []
101 for low, high in regions_in_axis(settings["integration_regions_ppm"], ppm):
102 axis, values = ppm[::-1], real[::-1]
103 inside = (axis > low) & (axis < high)
104 x = np.concatenate(([low], axis[inside], [high]))
105 y = np.concatenate(([np.interp(low, axis, values)], values[inside], [np.interp(high, axis, values)]))
106 integrals.append({"low_ppm": low, "high_ppm": high, "area_signal_ppm": float(np.trapezoid(y, x))})
107 return {"ppm": ppm, "real": real, "imaginary": spectrum.imag, "baseline": baseline,
108 "peaks": peaks, "integrals": integrals, "acquired_points": len(fid),
109 "acquisition_time_s": len(fid) / sw, "digital_spacing_hz": sw / size}
112def load_fid(input_path, settings, input_format):
113 if input_format == "npz":
114 with np.load(input_path, allow_pickle=False) as archive:
115 if set(archive.files) != {"fid"}:
116 raise ValueError("NPZ archive must contain exactly the complex array 'fid'")
117 return archive["fid"], {"format": "npz"}
118 if input_format != "nmrpipe":
119 raise ValueError("input_format must be npz or nmrpipe")
120 header, fid = ng.pipe.read(str(input_path))
121 if header["FDDIMCOUNT"] != 1 or header["FDDIMORDER"][0] != 2 or fid.ndim != 1:
122 raise ValueError("Only canonical 1D NMRPipe direct-dimension FDF2 files are supported")
123 if header["FDF2FTFLAG"] != 0 or header["FDF2QUADFLAG"] != 0 or not np.iscomplexobj(fid):
124 raise ValueError("NMRPipe input must be complex and time-domain (FTFLAG=0, QUADFLAG=0)")
125 if len(fid) < 16:
126 raise ValueError("NMRPipe FID must contain at least 16 complex points")
127 if header["FDSIZE"] != len(fid):
128 raise ValueError("NMRPipe header size disagrees with the decoded FID")
129 if header["FDF2TDSIZE"] != len(fid):
130 raise ValueError("NMRPipe FDF2TDSIZE must equal the stored FID length; prior zero filling or truncation is unsupported")
131 for field, key in (("FDF2SW", "spectral_width_hz"), ("FDF2OBS", "observation_mhz"), ("FDF2CAR", "carrier_ppm")):
132 value = finite_number(settings[key], key)
133 if not np.isclose(header[field], value, rtol=1e-6, atol=1e-6):
134 raise ValueError(f"NMRPipe {field} disagrees with explicit {key}")
135 # Match the centered, unextracted axis constructed by pipe.create_dic.
136 # CAR alone is insufficient: make_uc derives its carrier from ORIG instead.
137 center = len(fid) // 2 + 1
138 origin = header["FDF2CAR"] * header["FDF2OBS"] - header["FDF2SW"] * (len(fid) - center) / len(fid)
139 origin_atol = 1e-6 * max(abs(header["FDF2CAR"] * header["FDF2OBS"]), abs(header["FDF2SW"]), 1.0)
140 if header["FDF2CENTER"] != center or not np.isclose(header["FDF2ORIG"], origin, rtol=1e-6, atol=origin_atol):
141 raise ValueError("NMRPipe FDF2CENTER/FDF2ORIG must describe a canonical centered axis")
142 recorded = {key: header[key] for key in ("FDDIMCOUNT", "FDSIZE", "FDF2SW", "FDF2OBS", "FDF2CAR",
143 "FDF2LABEL", "FDF2FTFLAG", "FDF2QUADFLAG", "FDF2P0", "FDF2P1",
144 "FDF2TDSIZE", "FDF2CENTER", "FDF2ORIG")}
145 return fid, {"format": "nmrpipe", "header": recorded,
146 "reader_validation": "synthetic 1D round trip and upstream NMRPipe-generated fixture; no experimental vendor validation"}
149def run(input_path, settings_path, output_dir, input_format="npz"):
150 input_path, settings_path, output_dir = map(Path, (input_path, settings_path, output_dir))
151 if output_dir.exists():
152 raise ValueError("Output directory already exists; choose a new output directory")
153 settings = json.loads(settings_path.read_text())
154 fid, input_metadata = load_fid(input_path, settings, input_format)
155 result = process(fid, settings)
156 report = {"schema_version": 1, "input": input_metadata, "input_sha256": hashlib.sha256(input_path.read_bytes()).hexdigest(),
157 "settings_sha256": hashlib.sha256(settings_path.read_bytes()).hexdigest(),
158 "versions": {"nmrglue": ng.__version__, "numpy": np.__version__, "scipy": scipy.__version__},
159 "settings": settings, "acquired_points": result["acquired_points"],
160 "acquisition_time_s": result["acquisition_time_s"],
161 "digital_spacing_hz": result["digital_spacing_hz"],
162 "ppm_order": "descending", "fft_normalization": "unnormalized",
163 "phase_law_deg": "phase0_deg + phase1_deg * index / zero_fill_points",
164 "peaks": result["peaks"], "integrals": result["integrals"],
165 "limits": ["Peak candidates are not assignments", "Integrals are signed arbitrary signal times ppm; no concentration calibration"]}
166 output_dir.mkdir(parents=True)
167 with (output_dir / "spectrum.csv").open("w", newline="") as handle:
168 writer = csv.writer(handle)
169 writer.writerow(["ppm", "real", "imaginary", "baseline"])
170 writer.writerows(zip(result["ppm"], result["real"], result["imaginary"], result["baseline"]))
171 (output_dir / "report.json").write_text(json.dumps(report, indent=2, allow_nan=False) + "\n")
172 return report
175def main():
176 parser = argparse.ArgumentParser(description=__doc__)
177 parser.add_argument("input", type=Path, help="Complex FID: NPZ or canonical 1D NMRPipe")
178 parser.add_argument("--input-format", choices=("npz", "nmrpipe"), default="npz")
179 parser.add_argument("settings", type=Path, help="Explicit acquisition and processing JSON")
180 parser.add_argument("output", type=Path, help="New output directory")
181 args = parser.parse_args()
182 try:
183 report = run(args.input, args.settings, args.output, args.input_format)
184 except (ValueError, OSError, KeyError, TypeError) as error:
185 parser.error(str(error))
186 print(json.dumps({"output": str(args.output), "peak_candidates": len(report["peaks"])}))
189if __name__ == "__main__":
190 main()