{{ message }}
-
Notifications
You must be signed in to change notification settings - Fork 5
Expand file tree
/
Copy pathinputs.py
More file actions
603 lines (490 loc) · 23.3 KB
/
Copy pathinputs.py
File metadata and controls
603 lines (490 loc) · 23.3 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
"""Converts python models to the necessary inputs for the compiled RAT code."""
import importlib
import os
import pathlib
import warnings
from collections.abc import Callable
import numpy as np
import ratapi
import ratapi.wrappers
from ratapi.rat_core import Checks, Control, NameStore, ProblemDefinition
from ratapi.utils.enums import Calculations, Languages, LayerModels, Procedures, TypeOptions
parameter_field = {
"parameters": "params",
"bulk_in": "bulkIns",
"bulk_out": "bulkOuts",
"scalefactors": "scalefactors",
"domain_ratios": "domainRatios",
"background_parameters": "backgroundParams",
"resolution_parameters": "resolutionParams",
}
def get_python_handle(file_name: str, function_name: str, path: str | pathlib.Path = "") -> Callable:
"""Get the function handle from a function defined in a python module located anywhere within the filesystem.
Parameters
----------
file_name : str
The name of the file containing the function of interest.
function_name : str
The name of the function we wish to obtain the handle for within the module.
path : str
The path to the file containing the function (default is "", which represent the working directory).
Returns
-------
handle : Callable
The handle of the function defined in the python module file.
"""
spec = importlib.util.spec_from_file_location(pathlib.Path(file_name).stem, os.path.join(path, file_name))
custom_module = importlib.util.module_from_spec(spec)
spec.loader.exec_module(custom_module)
handle = getattr(custom_module, function_name)
return handle
class FileHandles:
"""Class to defer creation of custom file handles.
Parameters
----------
files : ClassList[CustomFile]
A list of custom file models.
"""
def __init__(self, files=None):
self.index = 0
self.files = [] if files is None else [file.model_dump() for file in files]
def __iter__(self):
self.index = 0
return self
def get_handle(self, index: int):
"""Return file handle for a given custom file.
Parameters
----------
index : int
The index of the custom file.
"""
custom_file = self.files[index]
full_path = os.path.join(custom_file["path"], custom_file["filename"])
if not os.path.isfile(full_path):
raise FileNotFoundError(f"The custom file ({custom_file['name']}) does not have a valid path.")
if not custom_file["function_name"] and custom_file["language"] != Languages.Matlab:
raise ValueError(f"The custom file ({custom_file['name']}) does not have a valid function name.")
if custom_file["language"] == Languages.Python:
file_handle = get_python_handle(custom_file["filename"], custom_file["function_name"], custom_file["path"])
elif custom_file["language"] == Languages.Matlab:
file_handle = ratapi.wrappers.MatlabWrapper(full_path).getHandle()
elif custom_file["language"] == Languages.Cpp:
file_handle = ratapi.wrappers.DylibWrapper(full_path, custom_file["function_name"]).getHandle()
return file_handle
def copy(self) -> "FileHandles":
"""Create a copy of the FileHandles object.
Returns
-------
FileHandles
The copy of this FileHandles object.
"""
handles = FileHandles()
handles.files = [file.copy() for file in self.files]
return handles
def __next__(self):
if self.index < len(self.files):
custom_file = self.get_handle(self.index)
self.index += 1
return custom_file
else:
raise StopIteration
def __len__(self):
return len(self.files)
def make_input(project: ratapi.Project, controls: ratapi.Controls) -> tuple[ProblemDefinition, Control]:
"""Construct the inputs required for the compiled RAT code using the data defined in the input project and controls.
Parameters
----------
project : RAT.Project
The project model, which defines the physical system under study.
controls : RAT.Controls
The controls model, which defines algorithmic properties.
Returns
-------
problem : RAT.rat_core.ProblemDefinition
The problem input used in the compiled RAT code.
cpp_controls : RAT.rat_core.Control
The controls object used in the compiled RAT code.
"""
problem = make_problem(project, controls.procedure != Procedures.Calculate)
cpp_controls = make_controls(controls)
return problem, cpp_controls
def make_problem(project: ratapi.Project, validate_range: bool = False) -> ProblemDefinition:
"""Construct the problem input required for the compiled RAT code.
Parameters
----------
project : RAT.Project
The project model, which defines the physical system under study.
validate_range : bool, default True
Whether parameter range should be validated.
Returns
-------
problem : RAT.rat_core.ProblemDefinition
The problem input used in the compiled RAT code.
"""
prior_id = {"uniform": 1, "gaussian": 2, "jeffreys": 3}
# Ensure all contrast fields are properly defined
for contrast in project.contrasts:
contrast_fields = ["data", "background", "bulk_in", "bulk_out", "scalefactor", "resolution"]
if project.calculation == Calculations.Domains:
contrast_fields.append("domain_ratio")
for field in contrast_fields:
if getattr(contrast, field) == "":
raise ValueError(
f'In the input project, the "{field}" field of contrast "{contrast.name}" does not have a '
f"value defined. A value must be supplied before running the project."
)
# Ensure backgrounds and resolutions have a source defined
background = project.backgrounds[contrast.background]
resolution = project.resolutions[contrast.resolution]
if background.source == "":
raise ValueError(
f"All backgrounds must have a source defined. For a {background.type} type background, "
f"the source must be defined in "
f'"{ratapi.project.values_defined_in[f"backgrounds.{background.type}.source"]}"'
)
if resolution.source == "" and resolution.type != TypeOptions.Data:
raise ValueError(
f"Constant resolutions must have a source defined. The source must be defined in "
f'"{ratapi.project.values_defined_in[f"resolutions.{resolution.type}.source"]}"'
)
# Set contrast parameters according to model type
if project.model == LayerModels.StandardLayers:
if project.calculation == Calculations.Domains:
contrast_models = [
[project.domain_contrasts.index(domain_contrast, True) for domain_contrast in contrast.model]
for contrast in project.contrasts
]
else:
contrast_models = [
[project.layers.index(layer, True) for layer in contrast.model] for contrast in project.contrasts
]
else:
contrast_models = [[]] * len(project.contrasts)
# Set contrast parameters according to model type
if project.model == LayerModels.StandardLayers:
contrast_custom_files = [float("NaN")] * len(project.contrasts)
else:
contrast_custom_files = [project.custom_files.index(contrast.model[0], True) for contrast in project.contrasts]
# Get details of defined layers
layer_details = get_layer_details(project)
contrast_background_params = []
contrast_background_types = []
all_data = []
data_limits = []
simulation_limits = []
contrast_resolution_params = []
contrast_resolution_types = []
# Set data, background and resolution for each contrast
for contrast in project.contrasts:
# Set data
data_index = project.data.index(contrast.data)
data = project.data[data_index].data
data_range = project.data[data_index].data_range
simulation_range = project.data[data_index].simulation_range
if data_range:
data_limits.append(data_range)
else:
data_limits.append([0.0, 0.0])
if simulation_range:
simulation_limits.append(simulation_range)
else:
simulation_limits.append([0.0, 0.0])
# Set background parameters
background = project.backgrounds[contrast.background]
contrast_background_types.append(background.type)
contrast_background_param = []
if background.type == TypeOptions.Data:
# Background data is appended to contrast data so empty index
if background.value_1 != "":
contrast_background_param.append(project.background_parameters.index(background.value_1))
# If we are using a data background, we add the background data to the contrast data
data = append_data_background(data, project.data[background.source].data)
elif background.type == TypeOptions.Function:
contrast_background_param.append(project.custom_files.index(background.source, True))
contrast_background_param.extend(
[
project.background_parameters.index(value, True)
for value in [
background.value_1,
background.value_2,
background.value_3,
background.value_4,
background.value_5,
]
if value != ""
]
)
else:
contrast_background_param.append(project.background_parameters.index(background.source, True))
contrast_background_params.append(contrast_background_param)
# Set resolution parameters
resolution = project.resolutions[contrast.resolution]
contrast_resolution_types.append(resolution.type)
contrast_resolution_param = []
if resolution.type == TypeOptions.Function:
contrast_resolution_param.append(project.custom_files.index(resolution.source, True))
contrast_resolution_param.extend(
[
project.resolution_parameters.index(value, True)
for value in [
resolution.value_1,
resolution.value_2,
resolution.value_3,
resolution.value_4,
resolution.value_5,
]
if value != ""
]
)
elif resolution.type == TypeOptions.Constant:
contrast_resolution_param.append(project.resolution_parameters.index(resolution.source, True))
contrast_resolution_params.append(contrast_resolution_param)
# Contrast data has exactly six columns to include background data if relevant
all_data.append(np.column_stack((data, np.zeros((data.shape[0], 6 - data.shape[1])))))
problem = ProblemDefinition()
problem.TF = project.calculation
problem.resample = make_resample(project)
problem.data = all_data
problem.dataPresent = make_data_present(project)
problem.dataLimits = data_limits
problem.simulationLimits = simulation_limits
problem.numberOfContrasts = len(project.contrasts)
problem.geometry = project.geometry
problem.useImaginary = project.absorption
problem.repeatLayers = [contrast.repeat_layers for contrast in project.contrasts]
problem.contrastBackgroundParams = contrast_background_params
problem.contrastBackgroundTypes = contrast_background_types
problem.contrastBackgroundActions = [contrast.background_action for contrast in project.contrasts]
problem.contrastScalefactors = [
project.scalefactors.index(contrast.scalefactor, True) for contrast in project.contrasts
]
problem.contrastBulkIns = [project.bulk_in.index(contrast.bulk_in, True) for contrast in project.contrasts]
problem.contrastBulkOuts = [project.bulk_out.index(contrast.bulk_out, True) for contrast in project.contrasts]
problem.contrastResolutionParams = contrast_resolution_params
problem.contrastResolutionTypes = contrast_resolution_types
problem.backgroundParams = [param.value for param in project.background_parameters]
problem.scalefactors = [param.value for param in project.scalefactors]
problem.bulkIns = [param.value for param in project.bulk_in]
problem.bulkOuts = [param.value for param in project.bulk_out]
problem.resolutionParams = [param.value for param in project.resolution_parameters]
problem.params = [param.value for param in project.parameters]
problem.numberOfLayers = len(project.layers)
problem.contrastLayers = [contrast_model if contrast_model else [] for contrast_model in contrast_models]
problem.layersDetails = layer_details if project.model == LayerModels.StandardLayers else []
problem.customFiles = FileHandles(project.custom_files)
problem.modelType = project.model
problem.contrastCustomFiles = contrast_custom_files
problem.contrastDomainRatios = [
project.domain_ratios.index(contrast.domain_ratio, True) if hasattr(contrast, "domain_ratio") else 0
for contrast in project.contrasts
]
problem.domainRatios = [param.value for param in project.domain_ratios]
problem.numberOfDomainContrasts = len(project.domain_contrasts)
domain_contrast_models = [
[project.layers.index(layer, True) for layer in domain_contrast.model]
for domain_contrast in project.domain_contrasts
]
problem.domainContrastLayers = [
domain_contrast_model if domain_contrast_model else [] for domain_contrast_model in domain_contrast_models
]
fit_params = []
fit_limits = []
prior_names = []
prior_values = []
problem.checks = Checks()
for class_list in ratapi.project.parameter_class_lists:
field = parameter_field[class_list]
check_list = []
for param in getattr(project, class_list):
prior_names.append(param.name)
prior_values.append([prior_id[param.prior_type], param.mu, param.sigma])
check_list.append(int(param.fit))
if param.fit:
min_range = 1e-6 if param.value == 0 else abs(param.value) * 1e-6
if validate_range and (param.max - param.min) < min_range:
warnings.warn(
f'{class_list.replace("_", " ").title()} "{param.name}" was removed from the '
f"fit because its range is too small (< {min_range:g}).",
stacklevel=2,
)
check_list[-1] = 0
else:
fit_params.append(param.value)
fit_limits.append([param.min, param.max])
setattr(problem.checks, field, check_list)
problem.fitParams = fit_params
problem.fitLimits = fit_limits
problem.priorNames = prior_names
problem.priorValues = prior_values
# Names
problem.names = NameStore()
for class_list in ratapi.project.parameter_class_lists:
setattr(problem.names, parameter_field[class_list], [param.name for param in getattr(project, class_list)])
problem.names.contrasts = [contrast.name for contrast in project.contrasts]
check_indices(problem)
return problem
def get_layer_details(project: ratapi.Project) -> list[int]:
"""Get parameter indices for all layers defined in the project."""
hydrate_id = {"bulk in": 0, "bulk out": 1}
layer_details = []
# Get the thickness, SLD, roughness fields from the appropriate model
if project.absorption:
layer_fields = list(ratapi.models.AbsorptionLayer.model_fields.keys())[1:-2]
else:
layer_fields = list(ratapi.models.Layer.model_fields.keys())[1:-2]
for layer in project.layers:
for field in layer_fields:
if getattr(layer, field) == "":
raise ValueError(
f'In the input project, the "{field}" field of layer {layer.name} does not have a value '
f"defined. A value must be supplied before running the project."
)
layer_params = [project.parameters.index(getattr(layer, attribute), True) for attribute in list(layer_fields)]
layer_params.append(project.parameters.index(layer.hydration, True) if layer.hydration else float("NaN"))
layer_params.append(hydrate_id[layer.hydrate_with])
layer_details.append(layer_params)
return layer_details
def make_resample(project: ratapi.Project) -> list[int]:
"""Construct the "resample" field of the problem input required for the compiled RAT code.
Parameters
----------
project : RAT.Project
The project model, which defines the physical system under study.
Returns
-------
list[int]
The "resample" field of the problem input used in the compiled RAT code.
"""
return [contrast.resample for contrast in project.contrasts]
def make_data_present(project: ratapi.Project) -> list[int]:
"""Construct the "dataPresent" field of the problem input required for the compiled RAT code.
Parameters
----------
project : RAT.Project
The project model, which defines the physical system under study.
Returns
-------
list[int]
The "dataPresent" field of the problem input used in the compiled RAT code.
"""
return [1 if project.data[contrast.data].data.size != 0 else 0 for contrast in project.contrasts]
def check_indices(problem: ProblemDefinition) -> None:
"""Check the indices given in a problem's contrasts lie within the range of the corresponding parameter lists.
Parameters
----------
problem : RAT.rat_core.ProblemDefinition
The problem input used in the compiled RAT code.
"""
index_list = {
"scalefactors": "contrastScalefactors",
"bulkIns": "contrastBulkIns",
"bulkOuts": "contrastBulkOuts",
"domainRatios": "contrastDomainRatios",
}
# Check the indices -- note we have switched to 1-based indexing at this point
for params in index_list:
param_list = getattr(problem, params)
if len(param_list) > 0:
elements = [
element
for element in getattr(problem, index_list[params])
if (element != -1) and not (0 < element <= len(param_list))
]
if elements:
raise IndexError(
f'The problem field "{index_list[params]}" contains: {", ".join(str(i) for i in elements)}'
f', which lie{"s" * (len(elements) == 1)} outside of the range of "{params}"',
)
# backgroundParams has a different structure, so is handled separately:
# it is of type list[list[int]], where each list[int] is the indices for
# source, value_1, value_2, value_3, value_4, value_5 where they are defined
# e.g. for a data background with offset it is [source value_1], for a function
# with 3 values it is [source value_1 value_2 value_3], etc.
source_param_lists = {
"constant": "backgroundParams",
"data": "backgroundParams",
"function": "customFiles",
}
for i, background_data in enumerate(problem.contrastBackgroundParams):
background_type = problem.contrastBackgroundTypes[i]
# check source param is in range of the relevant parameter list
param_list = getattr(problem, source_param_lists[background_type])
# For data, background data is appended to the contrast data so first index could
# either be empty or be the index of an offset.
if background_data and not 0 < background_data[0] <= len(param_list):
first_entry_name = "data offset" if background_type == "data" else "source"
raise IndexError(
f'Entry {i} of contrastBackgroundParams has type "{background_type}" '
f"and {first_entry_name} index {background_data[0]}, "
f'which is outside the range of "{source_param_lists[background_type]}".'
)
# check value params are in range for background params
if len(background_data) > 1:
elements = [element for element in background_data[1:] if not 0 < element <= len(problem.backgroundParams)]
if elements:
raise IndexError(
f"Entry {i} of contrastBackgroundParams contains: {', '.join(str(i) for i in elements)}"
f', which lie{"s" * (len(elements) == 1)} outside of the range of "backgroundParams"',
)
def append_data_background(data: np.array, background: np.array) -> np.array:
"""Append background data to contrast data.
Parameters
----------
data : np.array
The contrast data to which we are appending a background.
background : np.array
The background data to append to the contrast.
Returns
-------
np.array
The contrast data with background data appended as two additional columns.
"""
if not np.allclose(data[:, 0], background[:, 0]):
raise ValueError("The q-values of the data and background must be equal.")
return np.hstack((data, np.zeros((data.shape[0], 4 - data.shape[1])), background[:, 1:]))
def make_controls(input_controls: ratapi.Controls) -> Control:
"""Convert the controls object to the format required by the compiled RAT code.
Parameters
----------
input_controls : RAT.Controls
The controls model, which defines algorithmic properties.
Returns
-------
controls : RAT.rat_core.Control
The controls object used in the compiled RAT code.
"""
controls = Control()
controls.procedure = input_controls.procedure
controls.parallel = input_controls.parallel
controls.numSimulationPoints = input_controls.numSimulationPoints
controls.resampleMinAngle = input_controls.resampleMinAngle
controls.resampleNPoints = input_controls.resampleNPoints
controls.display = input_controls.display
# Simplex
controls.xTolerance = input_controls.xTolerance
controls.funcTolerance = input_controls.funcTolerance
controls.maxFuncEvals = input_controls.maxFuncEvals
controls.maxIterations = input_controls.maxIterations
controls.updateFreq = input_controls.updateFreq
controls.updatePlotFreq = input_controls.updatePlotFreq
# DE
controls.populationSize = input_controls.populationSize
controls.fWeight = input_controls.fWeight
controls.crossoverProbability = input_controls.crossoverProbability
controls.strategy = int(input_controls.strategy) # RAT core expects strategy as an integer
controls.targetValue = input_controls.targetValue
controls.numGenerations = input_controls.numGenerations
# NS
controls.nLive = input_controls.nLive
controls.nMCMC = input_controls.nMCMC
controls.propScale = input_controls.propScale
controls.nsTolerance = input_controls.nsTolerance
# Dream
controls.nSamples = input_controls.nSamples
controls.nChains = input_controls.nChains
controls.jumpProbability = input_controls.jumpProbability
controls.pUnitGamma = input_controls.pUnitGamma
controls.boundHandling = input_controls.boundHandling
controls.adaptPCR = input_controls.adaptPCR
controls.IPCFilePath = input_controls._IPCFilePath
return controls
You can’t perform that action at this time.