Mercurial > repos > kls286 > chap_test_20230328

comparison build/lib/CHAP/models/integration.py @ 0:cbbe42422d56 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/CHESSComputing/ChessAnalysisPipeline/tree/galaxy commit 1401a7e1ae007a6bda260d147f9b879e789b73e0-dirty
author kls286
date Tue, 28 Mar 2023 15:07:30 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:cbbe42422d56
1 import copy
2 from functools import cache, lru_cache
3 import json
4 import logging
5 import os
6 from time import time
7 from typing import Literal, Optional
8
9 from msnctools.general import input_menu
10 from multiprocessing.pool import ThreadPool
11 from nexusformat.nexus import (NXdata,
12 NXdetector,
13 NXfield,
14 NXprocess,
15 NXroot)
16 import numpy as np
17 from pydantic import (BaseModel,
18 validator,
19 constr,
20 conlist,
21 conint,
22 confloat,
23 FilePath)
24 import pyFAI, pyFAI.multi_geometry, pyFAI.units
25 from pyspec.file.tiff import TiffFile
26
27 from .map import MapConfig, SpecScans
28
29
30 class Detector(BaseModel):
31 """
32 Detector class to represent a single detector used in the experiment.
33
34 :param prefix: Prefix of the detector in the SPEC file.
35 :type prefix: str
36 :param poni_file: Path to the poni file.
37 :type poni_file: str
38 :param mask_file: Optional path to the mask file.
39 :type mask_file: str, optional
40 """
41 prefix: constr(strip_whitespace=True, min_length=1)
42 poni_file: FilePath
43 mask_file: Optional[FilePath]
44 @validator('poni_file', allow_reuse=True)
45 def validate_poni_file(cls, poni_file):
46 """
47 Validate the poni file by checking if it's a valid PONI file.
48
49 :param poni_file: Path to the poni file.
50 :type poni_file: str
51 :raises ValueError: If poni_file is not a valid PONI file.
52 :returns: Absolute path to the poni file.
53 :rtype: str
54 """
55 poni_file = os.path.abspath(poni_file)
56 try:
57 ai = azimuthal_integrator(poni_file)
58 except:
59 raise(ValueError(f'{poni_file} is not a valid PONI file'))
60 else:
61 return(poni_file)
62 @validator('mask_file', allow_reuse=True)
63 def validate_mask_file(cls, mask_file, values):
64 """
65 Validate the mask file. If a mask file is provided, it checks if it's a valid TIFF file.
66
67 :param mask_file: Path to the mask file.
68 :type mask_file: str or None
69 :param values: A dictionary of the Detector fields.
70 :type values: dict
71 :raises ValueError: If mask_file is provided and it's not a valid TIFF file.
72 :raises ValueError: If `'poni_file'` is not provided in `values`.
73 :returns: Absolute path to the mask file or None.
74 :rtype: str or None
75 """
76 if mask_file is None:
77 return(mask_file)
78 else:
79 mask_file = os.path.abspath(mask_file)
80 poni_file = values.get('poni_file')
81 if poni_file is None:
82 raise(ValueError('Cannot validate mask file without a PONI file.'))
83 else:
84 try:
85 mask_array = get_mask_array(mask_file, poni_file)
86 except BaseException as e:
87 raise(ValueError(f'Unable to open {mask_file} as a TIFF file'))
88 else:
89 return(mask_file)
90 @property
91 def azimuthal_integrator(self):
92 return(azimuthal_integrator(self.poni_file))
93 @property
94 def mask_array(self):
95 return(get_mask_array(self.mask_file, self.poni_file))
96
97 @cache
98 def azimuthal_integrator(poni_file:str):
99 if not isinstance(poni_file, str):
100 poni_file = str(poni_file)
101 return(pyFAI.load(poni_file))
102 @cache
103 def get_mask_array(mask_file:str, poni_file:str):
104 if mask_file is not None:
105 if not isinstance(mask_file, str):
106 mask_file = str(mask_file)
107 with TiffFile(mask_file) as tiff:
108 mask_array = tiff.asarray()
109 else:
110 mask_array = np.zeros(azimuthal_integrator(poni_file).detector.shape)
111 return(mask_array)
112
113 class IntegrationConfig(BaseModel):
114 """
115 Class representing the configuration for a raw detector data integration.
116
117 :ivar tool_type: type of integration tool; always set to "integration"
118 :type tool_type: str, optional
119 :ivar title: title of the integration
120 :type title: str
121 :ivar integration_type: type of integration, one of "azimuthal", "radial", or "cake"
122 :type integration_type: str
123 :ivar detectors: list of detectors used in the integration
124 :type detectors: List[Detector]
125 :ivar radial_units: radial units for the integration, defaults to `'q_A^-1'`
126 :type radial_units: str, optional
127 :ivar radial_min: minimum radial value for the integration range
128 :type radial_min: float, optional
129 :ivar radial_max: maximum radial value for the integration range
130 :type radial_max: float, optional
131 :ivar radial_npt: number of points in the radial range for the integration
132 :type radial_npt: int, optional
133 :ivar azimuthal_units: azimuthal units for the integration
134 :type azimuthal_units: str, optional
135 :ivar azimuthal_min: minimum azimuthal value for the integration range
136 :type azimuthal_min: float, optional
137 :ivar azimuthal_max: maximum azimuthal value for the integration range
138 :type azimuthal_max: float, optional
139 :ivar azimuthal_npt: number of points in the azimuthal range for the integration
140 :type azimuthal_npt: int, optional
141 :ivar error_model: error model for the integration, one of "poisson" or "azimuthal"
142 :type error_model: str, optional
143 """
144 tool_type: Literal['integration'] = 'integration'
145 title: constr(strip_whitespace=True, min_length=1)
146 integration_type: Literal['azimuthal', 'radial', 'cake']
147 detectors: conlist(item_type=Detector, min_items=1)
148 radial_units: str = 'q_A^-1'
149 radial_min: confloat(ge=0)
150 radial_max: confloat(gt=0)
151 radial_npt: conint(gt=0) = 1800
152 azimuthal_units: str = 'chi_deg'
153 azimuthal_min: confloat(ge=-180) = -180
154 azimuthal_max: confloat(le=360) = 180
155 azimuthal_npt: conint(gt=0) = 3600
156 error_model: Optional[Literal['poisson', 'azimuthal']]
157 sequence_index: Optional[conint(gt=0)]
158 @validator('radial_units', allow_reuse=True)
159 def validate_radial_units(cls, radial_units):
160 """
161 Validate the radial units for the integration.
162
163 :param radial_units: unvalidated radial units for the integration
164 :type radial_units: str
165 :raises ValueError: if radial units are not one of the recognized radial units
166 :return: validated radial units
167 :rtype: str
168 """
169 if radial_units in pyFAI.units.RADIAL_UNITS.keys():
170 return(radial_units)
171 else:
172 raise(ValueError(f'Invalid radial units: {radial_units}. Must be one of {", ".join(pyFAI.units.RADIAL_UNITS.keys())}'))
173 @validator('azimuthal_units', allow_reuse=True)
174 def validate_azimuthal_units(cls, azimuthal_units):
175 """
176 Validate that `azimuthal_units` is one of the keys in the
177 `pyFAI.units.AZIMUTHAL_UNITS` dictionary.
178
179 :param azimuthal_units: The string representing the unit to be validated.
180 :type azimuthal_units: str
181 :raises ValueError: If `azimuthal_units` is not one of the keys in `pyFAI.units.AZIMUTHAL_UNITS`
182 :return: The original supplied value, if is one of the keys in `pyFAI.units.AZIMUTHAL_UNITS`.
183 :rtype: str
184 """
185 if azimuthal_units in pyFAI.units.AZIMUTHAL_UNITS.keys():
186 return(azimuthal_units)
187 else:
188 raise(ValueError(f'Invalid azimuthal units: {azimuthal_units}. Must be one of {", ".join(pyFAI.units.AZIMUTHAL_UNITS.keys())}'))
189 def validate_range_max(range_name:str):
190 """Validate the maximum value of an integration range.
191
192 :param range_name: The name of the integration range (e.g. radial, azimuthal).
193 :type range_name: str
194 :return: The callable that performs the validation.
195 :rtype: callable
196 """
197 def _validate_range_max(cls, range_max, values):
198 """Check if the maximum value of the integration range is greater than its minimum value.
199
200 :param range_max: The maximum value of the integration range.
201 :type range_max: float
202 :param values: The values of the other fields being validated.
203 :type values: dict
204 :raises ValueError: If the maximum value of the integration range is not greater than its minimum value.
205 :return: The validated maximum range value
206 :rtype: float
207 """
208 range_min = values.get(f'{range_name}_min')
209 if range_min < range_max:
210 return(range_max)
211 else:
212 raise(ValueError(f'Maximum value of integration range must be greater than minimum value of integration range ({range_name}_min={range_min}).'))
213 return(_validate_range_max)
214 _validate_radial_max = validator('radial_max', allow_reuse=True)(validate_range_max('radial'))
215 _validate_azimuthal_max = validator('azimuthal_max', allow_reuse=True)(validate_range_max('azimuthal'))
216 def validate_for_map_config(self, map_config:MapConfig):
217 """
218 Validate the existence of the detector data file for all scan points in `map_config`.
219
220 :param map_config: The `MapConfig` instance to validate against.
221 :type map_config: MapConfig
222 :raises RuntimeError: If a detector data file could not be found for a scan point occurring in `map_config`.
223 :return: None
224 :rtype: None
225 """
226 for detector in self.detectors:
227 for scans in map_config.spec_scans:
228 for scan_number in scans.scan_numbers:
229 scanparser = scans.get_scanparser(scan_number)
230 for scan_step_index in range(scanparser.spec_scan_npts):
231 # Make sure the detector data file exists for all scan points
232 try:
233 detector_data_file = scanparser.get_detector_data_file(detector.prefix, scan_step_index)
234 except:
235 raise(RuntimeError(f'Could not find data file for detector prefix {detector.prefix} on scan number {scan_number} in spec file {scans.spec_file}'))
236 def get_azimuthal_adjustments(self):
237 """To enable a continuous range of integration in the azimuthal direction
238 for radial and cake integration, obtain adjusted values for this
239 `IntegrationConfig`'s `azimuthal_min` and `azimuthal_max` values, the
240 angle amount by which those values were adjusted, and the proper location
241 of the discontinuity in the azimuthal direction.
242
243 :return: Adjusted chi_min, adjusted chi_max, chi_offset, chi_discontinuity
244 :rtype: tuple[float,float,float,float]
245 """
246 return(get_azimuthal_adjustments(self.azimuthal_min, self.azimuthal_max))
247 def get_azimuthal_integrators(self):
248 """Get a list of `AzimuthalIntegrator`s that correspond to the detector
249 configurations in this instance of `IntegrationConfig`.
250
251 The returned `AzimuthalIntegrator`s are (if need be) artificially rotated
252 in the azimuthal direction to achieve a continuous range of integration
253 in the azimuthal direction.
254
255 :returns: A list of `AzimuthalIntegrator`s appropriate for use by this
256 `IntegrationConfig` tool
257 :rtype: list[pyFAI.azimuthalIntegrator.AzimuthalIntegrator]
258 """
259 chi_min, chi_max, chi_offset, chi_disc = self.get_azimuthal_adjustments()
260 return(get_azimuthal_integrators(tuple([detector.poni_file for detector in self.detectors]), chi_offset=chi_offset))
261 def get_multi_geometry_integrator(self):
262 """Get a `MultiGeometry` integrator suitable for use by this instance of
263 `IntegrationConfig`.
264
265 :return: A `MultiGeometry` integrator
266 :rtype: pyFAI.multi_geometry.MultiGeometry
267 """
268 poni_files = tuple([detector.poni_file for detector in self.detectors])
269 radial_range = (self.radial_min, self.radial_max)
270 azimuthal_range = (self.azimuthal_min, self.azimuthal_max)
271 return(get_multi_geometry_integrator(poni_files, self.radial_units, radial_range, azimuthal_range))
272 def get_azimuthally_integrated_data(self, spec_scans:SpecScans, scan_number:int, scan_step_index:int):
273 """Return azimuthally-integrated data for the scan step specified.
274
275 :param spec_scans: An instance of `SpecScans` containing the scan step requested.
276 :type spec_scans: SpecScans
277 :param scan_number: The number of the scan containing the scan step requested.
278 :type scan_number: int
279 :param scan_step_index: The index of the scan step requested.
280 :type scan_step_index: int
281 :return: A 1D array of azimuthally-integrated raw detector intensities.
282 :rtype: np.ndarray
283 """
284 detector_data = spec_scans.get_detector_data(self.detectors, scan_number, scan_step_index)
285 integrator = self.get_multi_geometry_integrator()
286 lst_mask = [detector.mask_array for detector in self.detectors]
287 result = integrator.integrate1d(detector_data, lst_mask=lst_mask, npt=self.radial_npt, error_model=self.error_model)
288 if result.sigma is None:
289 return(result.intensity)
290 else:
291 return(result.intensity, result.sigma)
292 def get_radially_integrated_data(self, spec_scans:SpecScans, scan_number:int, scan_step_index:int):
293 """Return radially-integrated data for the scan step specified.
294
295 :param spec_scans: An instance of `SpecScans` containing the scan step requested.
296 :type spec_scans: SpecScans
297 :param scan_number: The number of the scan containing the scan step requested.
298 :type scan_number: int
299 :param scan_step_index: The index of the scan step requested.
300 :type scan_step_index: int
301 :return: A 1D array of radially-integrated raw detector intensities.
302 :rtype: np.ndarray
303 """
304 # Handle idiosyncracies of azimuthal ranges in pyFAI
305 # Adjust chi ranges to get a continuous range of iintegrated data
306 chi_min, chi_max, chi_offset, chi_disc = self.get_azimuthal_adjustments()
307 # Perform radial integration on a detector-by-detector basis.
308 I_each_detector = []
309 variance_each_detector = []
310 integrators = self.get_azimuthal_integrators()
311 for i,(integrator,detector) in enumerate(zip(integrators,self.detectors)):
312 detector_data = spec_scans.get_detector_data([detector], scan_number, scan_step_index)[0]
313 result = integrator.integrate_radial(detector_data, self.azimuthal_npt,
314 unit=self.azimuthal_units, azimuth_range=(chi_min,chi_max),
315 radial_unit=self.radial_units, radial_range=(self.radial_min,self.radial_max),
316 mask=detector.mask_array) #, error_model=self.error_model)
317 I_each_detector.append(result.intensity)
318 if result.sigma is not None:
319 variance_each_detector.append(result.sigma**2)
320 # Add the individual detectors' integrated intensities together
321 I = np.nansum(I_each_detector, axis=0)
322 # Ignore data at values of chi for which there was no data
323 I = np.where(I==0, np.nan, I)
324 if len(I_each_detector) != len(variance_each_detector):
325 return(I)
326 else:
327 # Get the standard deviation of the summed detectors' intensities
328 sigma = np.sqrt(np.nansum(variance_each_detector, axis=0))
329 return(I, sigma)
330 def get_cake_integrated_data(self, spec_scans:SpecScans, scan_number:int, scan_step_index:int):
331 """Return cake-integrated data for the scan step specified.
332
333 :param spec_scans: An instance of `SpecScans` containing the scan step requested.
334 :type spec_scans: SpecScans
335 :param scan_number: The number of the scan containing the scan step requested.
336 :type scan_number: int
337 :param scan_step_index: The index of the scan step requested.
338 :type scan_step_index: int
339 :return: A 2D array of cake-integrated raw detector intensities.
340 :rtype: np.ndarray
341 """
342 detector_data = spec_scans.get_detector_data(self.detectors, scan_number, scan_step_index)
343 integrator = self.get_multi_geometry_integrator()
344 lst_mask = [detector.mask_array for detector in self.detectors]
345 result = integrator.integrate2d(detector_data, lst_mask=lst_mask,
346 npt_rad=self.radial_npt, npt_azim=self.azimuthal_npt,
347 method='bbox',
348 error_model=self.error_model)
349 if result.sigma is None:
350 return(result.intensity)
351 else:
352 return(result.intensity, result.sigma)
353 def get_integrated_data(self, spec_scans:SpecScans, scan_number:int, scan_step_index:int):
354 """Return integrated data for the scan step specified.
355
356 :param spec_scans: An instance of `SpecScans` containing the scan step requested.
357 :type spec_scans: SpecScans
358 :param scan_number: The number of the scan containing the scan step requested.
359 :type scan_number: int
360 :param scan_step_index: The index of the scan step requested.
361 :type scan_step_index: int
362 :return: An array of integrated raw detector intensities.
363 :rtype: np.ndarray
364 """
365 if self.integration_type == 'azimuthal':
366 return(self.get_azimuthally_integrated_data(spec_scans, scan_number, scan_step_index))
367 elif self.integration_type == 'radial':
368 return(self.get_radially_integrated_data(spec_scans, scan_number, scan_step_index))
369 elif self.integration_type == 'cake':
370 return(self.get_cake_integrated_data(spec_scans, scan_number, scan_step_index))
371
372 @property
373 def integrated_data_coordinates(self):
374 """
375 Return a dictionary of coordinate arrays for navigating the dimension(s)
376 of the integrated data produced by this instance of `IntegrationConfig`.
377
378 :return: A dictionary with either one or two keys: 'azimuthal' and/or
379 'radial', each of which points to a 1-D `numpy` array of coordinate
380 values.
381 :rtype: dict[str,np.ndarray]
382 """
383 if self.integration_type == 'azimuthal':
384 return(get_integrated_data_coordinates(radial_range=(self.radial_min,self.radial_max),
385 radial_npt=self.radial_npt))
386 elif self.integration_type == 'radial':
387 return(get_integrated_data_coordinates(azimuthal_range=(self.azimuthal_min,self.azimuthal_max),
388 azimuthal_npt=self.azimuthal_npt))
389 elif self.integration_type == 'cake':
390 return(get_integrated_data_coordinates(radial_range=(self.radial_min,self.radial_max),
391 radial_npt=self.radial_npt,
392 azimuthal_range=(self.azimuthal_min,self.azimuthal_max),
393 azimuthal_npt=self.azimuthal_npt))
394 @property
395 def integrated_data_dims(self):
396 """Return a tuple of the coordinate labels for the integrated data
397 produced by this instance of `IntegrationConfig`.
398 """
399 directions = list(self.integrated_data_coordinates.keys())
400 dim_names = [getattr(self, f'{direction}_units') for direction in directions]
401 return(dim_names)
402 @property
403 def integrated_data_shape(self):
404 """Return a tuple representing the shape of the integrated data
405 produced by this instance of `IntegrationConfig` for a single scan step.
406 """
407 return(tuple([len(coordinate_values) for coordinate_name,coordinate_values in self.integrated_data_coordinates.items()]))
408
409 @cache
410 def get_azimuthal_adjustments(chi_min:float, chi_max:float):
411 """
412 Fix chi discontinuity at 180 degrees and return the adjusted chi range,
413 offset, and discontinuty.
414
415 If the discontinuity is crossed, obtain the offset to artificially rotate
416 detectors to achieve a continuous azimuthal integration range.
417
418 :param chi_min: The minimum value of the azimuthal range.
419 :type chi_min: float
420 :param chi_max: The maximum value of the azimuthal range.
421 :type chi_max: float
422 :return: The following four values: the adjusted minimum value of the
423 azimuthal range, the adjusted maximum value of the azimuthal range, the
424 value by which the chi angle was adjusted, the position of the chi
425 discontinuity.
426 """
427 # Fix chi discontinuity at 180 degrees for now.
428 chi_disc = 180
429 # If the discontinuity is crossed, artificially rotate the detectors to
430 # achieve a continuous azimuthal integration range
431 if chi_min < chi_disc and chi_max > chi_disc:
432 chi_offset = chi_max - chi_disc
433 else:
434 chi_offset = 0
435 return(chi_min-chi_offset, chi_max-chi_offset, chi_offset, chi_disc)
436 @cache
437 def get_azimuthal_integrators(poni_files:tuple, chi_offset=0):
438 """
439 Return a list of `AzimuthalIntegrator` objects generated from PONI files.
440
441 :param poni_files: Tuple of strings, each string being a path to a PONI file. : tuple
442 :type poni_files: tuple
443 :param chi_offset: The angle in degrees by which the `AzimuthalIntegrator` objects will be rotated, defaults to 0.
444 :type chi_offset: float, optional
445 :return: List of `AzimuthalIntegrator` objects
446 :rtype: list[pyFAI.azimuthalIntegrator.AzimuthalIntegrator]
447 """
448 ais = []
449 for poni_file in poni_files:
450 ai = copy.deepcopy(azimuthal_integrator(poni_file))
451 ai.rot3 += chi_offset * np.pi/180
452 ais.append(ai)
453 return(ais)
454 @cache
455 def get_multi_geometry_integrator(poni_files:tuple, radial_unit:str, radial_range:tuple, azimuthal_range:tuple):
456 """Return a `MultiGeometry` instance that can be used for azimuthal or cake
457 integration.
458
459 :param poni_files: Tuple of PONI files that describe the detectors to be
460 integrated.
461 :type poni_files: tuple
462 :param radial_unit: Unit to use for radial integration range.
463 :type radial_unit: str
464 :param radial_range: Tuple describing the range for radial integration.
465 :type radial_range: tuple[float,float]
466 :param azimuthal_range:Tuple describing the range for azimuthal integration.
467 :type azimuthal_range: tuple[float,float]
468 :return: `MultiGeometry` instance that can be used for azimuthal or cake
469 integration.
470 :rtype: pyFAI.multi_geometry.MultiGeometry
471 """
472 chi_min, chi_max, chi_offset, chi_disc = get_azimuthal_adjustments(*azimuthal_range)
473 ais = copy.deepcopy(get_azimuthal_integrators(poni_files, chi_offset=chi_offset))
474 multi_geometry = pyFAI.multi_geometry.MultiGeometry(ais,
475 unit=radial_unit,
476 radial_range=radial_range,
477 azimuth_range=(chi_min,chi_max),
478 wavelength=sum([ai.wavelength for ai in ais])/len(ais),
479 chi_disc=chi_disc)
480 return(multi_geometry)
481 @cache
482 def get_integrated_data_coordinates(azimuthal_range:tuple=None, azimuthal_npt:int=None, radial_range:tuple=None, radial_npt:int=None):
483 """
484 Return a dictionary of coordinate arrays for the specified radial and/or
485 azimuthal integration ranges.
486
487 :param azimuthal_range: Tuple specifying the range of azimuthal angles over
488 which to generate coordinates, in the format (min, max), defaults to
489 None.
490 :type azimuthal_range: tuple[float,float], optional
491 :param azimuthal_npt: Number of azimuthal coordinate points to generate,
492 defaults to None.
493 :type azimuthal_npt: int, optional
494 :param radial_range: Tuple specifying the range of radial distances over
495 which to generate coordinates, in the format (min, max), defaults to
496 None.
497 :type radial_range: tuple[float,float], optional
498 :param radial_npt: Number of radial coordinate points to generate, defaults
499 to None.
500 :type radial_npt: int, optional
501 :return: A dictionary with either one or two keys: 'azimuthal' and/or
502 'radial', each of which points to a 1-D `numpy` array of coordinate
503 values.
504 :rtype: dict[str,np.ndarray]
505 """
506 integrated_data_coordinates = {}
507 if azimuthal_range is not None and azimuthal_npt is not None:
508 integrated_data_coordinates['azimuthal'] = np.linspace(*azimuthal_range, azimuthal_npt)
509 if radial_range is not None and radial_npt is not None:
510 integrated_data_coordinates['radial'] = np.linspace(*radial_range, radial_npt)
511 return(integrated_data_coordinates)