From 9c362429a63c9b19bbb65a7a033a88627cddfc1e Mon Sep 17 00:00:00 2001 From: Barbara Frosik Date: Wed, 24 Apr 2024 10:57:47 -0500 Subject: [PATCH] Dev (#24) * Strain from CoupledRec (#16) * Additional documentation for multipeak phasing * Update setup.py * Update meta.yaml * Refactored get_centered to reduce number of numpy function calls. * Added a bunch of functions to cplib.py * Get strain from final reconstruction. This commit also contains commented code that would implement the partial fourier constraint with free-varying pixels according to a preset mask. * Added several new methods to the cohlib signature: - diff - gradient - argmin - take_along_axis - moveaxis - lstsq - zeros - indices - concatenate These methods have been implemented in cplib and nplib, but will raise NotImplementedError if called from torchlib or aflib. * Removed some experimental code that accidentally got included in this version. --------- Co-authored-by: Barbara Frosik * fixed interface issue with cohere-ui after merge * fixed interface issue after merge * Delete reconstruction_GA_HPC.py * moved cleaning memory to lib * Delete reconstruction_HPC_GA.py * fixed problrm with parameter * fixed imoprt name * fixed docs for developers * supressed warnings in verifier when parsing device * fixed crop issue in ui * Autodata (#18) * changes for autodata option * synch with ui --------- Co-authored-by: bfrosik * added debug option * moved utils from cohere-ui to cohere core and cleane up * cleaned up utils * support for auto data for separate scan ranges * aligned with cohere-ui * align with cohere-ui * relaxed verifier for auto data * relaxed verifier * aligning with ui * added missing save_dir in reconstruction scripts * added for hpc * sync with cohere-ui module * removed warning * added cluster capability * removed print statements * cleaned up, bug fixes * added description for cluster, fixed host naming issue * fixed typo * Update .readthedocs.yaml * Update .readthedocs.yaml * Update .readthedocs.yaml * Update .readthedocs.yaml * Update .readthedocs.yaml * Update README.md * Update README.md * Update README.md * fixed verifier, aligned cohere-ui * aligned with cohere-ui module * formatted text, fixed verifier for multipeak * blocked auto binning for mp and aligned ui * return correlation error in addition to aligned array * aligned ui * aligned with cohere-ui * added description in verifier, aligned with ui * Updated multipeak capabilities (Dev) (#22) * Additional documentation for multipeak phasing * Update setup.py * Update meta.yaml * Refactored get_centered to reduce number of numpy function calls. * Added a bunch of functions to cplib.py * Get strain from final reconstruction. This commit also contains commented code that would implement the partial fourier constraint with free-varying pixels according to a preset mask. * Added several new methods to the cohlib signature: - diff - gradient - argmin - take_along_axis - moveaxis - lstsq - zeros - indices - concatenate These methods have been implemented in cplib and nplib, but will raise NotImplementedError if called from torchlib or aflib. * Removed some experimental code that accidentally got included in this version. * assigned release tag * fixed doc for developers * fixed crop issue in ui * Added the following functions to cplib.py: amin(), affine_transform(), pad(), histogram2d(), calc_nmi(), calc_ehd(). These have only been implemented in cplib, other libraries are only implemented as stubs. * Significant changes that should have been committed earlier: - Resampling is now included in the phasing process. - Added several new error metrics, including normalized mutual information (NMI) and expected histogram deviation (EHD) - Apply the support constraint to the full object during ER iterations (to prevent buildup of HIO feedback) - Fixed incorrect normalization when projecting to each Bragg peak. - Added a "control_peak" option to exclude one peak from the phasing process to use for unbiased error calculations. - Added a "calc_strain" option which can skip the slow strain calculation. - Added a "fast_resample" option to toggle whether the original or resampled data is used. --------- Co-authored-by: Barbara Frosik * aligned with ui * updated torch lib * fixed torchlib * removed unused function, aligned ui * synch ui * revised fftconvolve function * added comments to op flow, and simplified * modified get_ratio * removed defaults in cohere-core, added checks for misconfiguration * added error checks and prints of error messages if fail * removed code to activate lpf at last iteration * aligned with cohere-ui * added tqdm package in setup.py * added example of different sw type in features * moved checks based on type * change definition of algorithms dict to include to_direct and to_reciprocal * refactored cluster handling * changed API to receive hostfile instead of arbitrary naming it * changed api * Removed Support class and moved support array to Rec * renamed phase modulus to phase constrain * change description, slightly modified code * updated documentation * corrected docs * modified for new release * modified to do beta buid * refactored to set defaults through get, corrected decs * modified doc conf.py * modified doc conf.py * updated docs * cleaned up docs * aligning with cohere-ui --------- Co-authored-by: Nick Porter <34489968+jacione@users.noreply.github.com> Co-authored-by: Barbara Frosik Co-authored-by: pfrosik Co-authored-by: cxduser --- .readthedocs.yaml | 30 +- README.md | 4 +- build.sh | 2 +- cohere_core/controller/AI_guess.py | 87 +-- cohere_core/controller/__init__.py | 2 + cohere_core/controller/features.py | 271 ++++--- cohere_core/controller/mpi_cmd.py | 45 ++ cohere_core/controller/op_flow.py | 338 +++++---- cohere_core/controller/phasing.py | 676 +++++++++++------ cohere_core/controller/reconstruction_GA.py | 672 +++++----------- .../controller/reconstruction_coupled.py | 24 +- .../controller/reconstruction_multi.py | 239 ++---- .../controller/reconstruction_populous_GA.py | 384 ++++++++++ .../controller/reconstruction_single.py | 22 +- cohere_core/data/alien_tools.py | 41 +- cohere_core/data/auto_prep.py | 0 cohere_core/data/standard_preprocess.py | 87 ++- cohere_core/lib/aflib.py | 49 ++ cohere_core/lib/cohlib.py | 122 ++- cohere_core/lib/cplib.py | 119 ++- cohere_core/lib/nplib.py | 76 +- cohere_core/lib/torchlib.py | 137 +++- cohere_core/utilities/config_errors_dict.py | 12 +- cohere_core/utilities/config_verifier.py | 114 +-- cohere_core/utilities/dvc_utils.py | 264 ++++++- cohere_core/utilities/ga_utils.py | 183 +++++ cohere_core/utilities/host_utils.py | 36 + cohere_core/utilities/utils.py | 717 ++++++++++++++---- docs/source/about.rst | 4 +- docs/source/conf.py | 13 +- docs/source/config_data.rst | 30 +- docs/source/config_disp.rst | 20 +- docs/source/config_instr.rst | 9 +- docs/source/config_main.rst | 29 +- docs/source/config_mp.rst | 29 +- docs/source/config_prep.rst | 19 +- docs/source/config_rec.rst | 136 ++-- docs/source/define_alg_seq.rst | 58 +- docs/source/experiment_space.rst | 78 +- docs/source/for_developers.rst | 34 +- docs/source/how_to_use.rst | 4 +- docs/source/index.rst | 2 +- docs/source/installation.rst | 30 +- meta.yaml | 27 - requirements.txt | 2 +- setup.py | 7 +- 46 files changed, 3451 insertions(+), 1833 deletions(-) create mode 100644 cohere_core/controller/mpi_cmd.py create mode 100644 cohere_core/controller/reconstruction_populous_GA.py delete mode 100644 cohere_core/data/auto_prep.py create mode 100644 cohere_core/utilities/ga_utils.py create mode 100644 cohere_core/utilities/host_utils.py delete mode 100644 meta.yaml diff --git a/.readthedocs.yaml b/.readthedocs.yaml index f924969..3ecaf61 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -5,24 +5,28 @@ # Required version: 2 -# Set the version of Python and other tools you might need -#build: -# python: "3.8" +# Set the OS, Python version and other tools you might need +build: + os: "ubuntu-22.04" + tools: + python: "3.12" # You can also specify other tool versions: - # nodejs: "16" - # rust: "1.55" - # golang: "1.17" + # nodejs: "19" + # rust: "1.64" + # golang: "1.19" -# Build documentation in the docs/ directory with Sphinx +# Build documentation in the "docs/" directory with Sphinx sphinx: configuration: docs/source/conf.py -# If using Sphinx, optionally build your docs in additional formats such as PDF +# Optionally build your docs in additional formats such as PDF and ePub # formats: # - pdf +# - epub -# Optionally declare the Python requirements required to build your docs -python: - version: 3.8 - install: - - requirements: requirements.txt +# Optional but recommended, declare the Python requirements required +# to build your documentation +# See https://docs.readthedocs.io/en/stable/guides/reproducible-builds.html +# python: +# install: +# - requirements: docs/requirements.txt diff --git a/README.md b/README.md index 8a5450f..3471bc6 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ -[![Documentation Status](https://readthedocs.org/projects/cohere/badge/?version=latest)](http://cohere.readthedocs.io/en/latest/?badge=latest) +[![Documentation Status](https://readthedocs.org/projects/cohere-dev/badge/?version=latest)](http://cohere-dev.readthedocs.io/en/latest/?badge=latest) -Project home page: https://cohere.readthedocs.io/ +Project home page: [https://cohere-dev.readthedocs.io/](https://cohere-dev.readthedocs.io/) The cohere package provides tools for reconstruction of image of a nanoscale structures from data obtained using Bragg Coherent Diffraction Imaging technique. diff --git a/build.sh b/build.sh index 40b0193..762d43b 100644 --- a/build.sh +++ b/build.sh @@ -1 +1 @@ -python setup.py install +python setup.py install diff --git a/cohere_core/controller/AI_guess.py b/cohere_core/controller/AI_guess.py index 8107952..324dfdc 100644 --- a/cohere_core/controller/AI_guess.py +++ b/cohere_core/controller/AI_guess.py @@ -51,84 +51,6 @@ def __init__(self, model_file): outputs=model.get_layer('phi').output) -def threshold_by_edge(fp: np.ndarray) -> np.ndarray: - # threshold by left edge value - mask = np.ones_like(fp, dtype=bool) - mask[tuple([slice(1, None)] * fp.ndim)] = 0 - zero = 1e-6 - cut = np.max(fp[mask]) - binary = np.zeros_like(fp) - binary[(np.abs(fp) > zero) & (fp > cut)] = 1 - return binary - - -def select_central_object(fp: np.ndarray) -> np.ndarray: - import scipy.ndimage as ndimage - zero = 1e-6 - binary = np.abs(fp) - binary[binary > zero] = 1 - binary[binary <= zero] = 0 - - # cluster by connectivity - struct = ndimage.morphology.generate_binary_structure(fp.ndim, - 1).astype("uint8") - label, nlabel = ndimage.label(binary, structure=struct) - - # select largest cluster - select = np.argmax(np.bincount(np.ravel(label))[1:]) + 1 - - binary[label != select] = 0 - - fp[binary == 0] = 0 - return fp - - -def get_central_object_extent(fp: np.ndarray) -> list: - fp_cut = threshold_by_edge(np.abs(fp)) - need = select_central_object(fp_cut) - - # get extend of cluster - extent = [np.max(s) + 1 - np.min(s) for s in np.nonzero(need)] - return extent - - -def get_oversample_ratio(fp: np.ndarray) -> np.ndarray: - """ get oversample ratio - fp = diffraction pattern - """ - # autocorrelation - acp = np.fft.fftshift(np.fft.ifftn(np.abs(fp)**2.)) - aacp = np.abs(acp) - - # get extent - blob = get_central_object_extent(aacp) - - # correct for underestimation due to thresholding - correction = [0.025, 0.025, 0.0729][:fp.ndim] - - extent = [ - min(m, s + int(round(f * aacp.shape[i], 1))) - for i, (s, f, m) in enumerate(zip(blob, correction, aacp.shape)) - ] - - # oversample ratio - oversample = [ - 2. * s / (e + (1 - s % 2)) for s, e in zip(aacp.shape, extent) - ] - return np.round(oversample, 3) - - -def Resize(IN, dim): - ft = np.fft.fftshift(np.fft.fftn(IN)) / np.prod(IN.shape) - - pad_value = np.array(dim) // 2 - np.array(ft.shape) // 2 - pad = [[pad_value[0], pad_value[0]], [pad_value[1], pad_value[1]], - [pad_value[2], pad_value[2]]] - ft_resize = ut.adjust_dimensions(ft, pad) - output = np.fft.ifftn(np.fft.ifftshift(ft_resize)) * np.prod(dim) - return output - - def match_oversample_diff( diff: np.ndarray, fr: Union[list, np.ndarray, None] = None, @@ -144,7 +66,6 @@ def match_oversample_diff( # adjustment needed to match oversample ratio change = [np.round(f / t).astype('int32') for f, t in zip(fr, to)] change = [np.max([1, c]) for c in change] - diff = ut.binning(diff, change) # crop diff to match output shape shape_arr = np.array(shape) @@ -293,7 +214,7 @@ def run_AI(data, model_file, dir): # prepare data to make the oversampling ratio ~3 wos = 3.0 - orig_os = get_oversample_ratio(data) + orig_os = ut.get_oversample_ratio(data) # match oversampling to wos wanted_os = [wos, wos, wos] # match diff os @@ -314,7 +235,7 @@ def run_AI(data, model_file, dir): pred_obj = preds_amp * np.exp(1j * preds_phi) # match object size with the input data - pred_obj = Resize(pred_obj, inshape) + pred_obj = ut.Resize(pred_obj, inshape) pad_value = np.array(data.shape) // 2 - np.array(pred_obj.shape) // 2 pad = [[pad_value[0], pad_value[0]], [pad_value[1], pad_value[1]], @@ -352,7 +273,7 @@ def start_AI(pars, datafile, dir): print ('no AI_trained_model in config') return None if not os.path.isfile(pars['AI_trained_model']): - print('there is no file', pars['AI_trained_model']) + print(f'there is no file {pars["AI_trained_model"]}') return None if datafile.endswith('tif') or datafile.endswith('tiff'): @@ -372,7 +293,7 @@ def start_AI(pars, datafile, dir): return None # The results will be stored in the directory /AI_guess - ai_dir = dir + '/results_AI' + ai_dir = ut.join(dir, 'results_AI') if os.path.exists(ai_dir): pass else: diff --git a/cohere_core/controller/__init__.py b/cohere_core/controller/__init__.py index ca8b1bd..c35f00e 100644 --- a/cohere_core/controller/__init__.py +++ b/cohere_core/controller/__init__.py @@ -1,3 +1,5 @@ +from .mpi_cmd import * +from .reconstruction_populous_GA import * from .reconstruction_GA import * from .reconstruction_single import * from .reconstruction_multi import * diff --git a/cohere_core/controller/features.py b/cohere_core/controller/features.py index e4c14c9..5a6033d 100644 --- a/cohere_core/controller/features.py +++ b/cohere_core/controller/features.py @@ -1,3 +1,4 @@ +import cohere_core.utilities.utils as ut import cohere_core.utilities.dvc_utils as dvut from abc import ABC, abstractmethod @@ -5,9 +6,9 @@ __copyright__ = "Copyright (c) 2016, UChicago Argonne, LLC." __docformat__ = 'restructuredtext en' __all__ = ['Pcdi', - 'Feature', + 'TriggeredOp', 'ShrinkWrap', - 'PhaseMod', + 'PhaseConstrain', 'LowPassFilter'] @@ -19,39 +20,35 @@ def set_lib(dlib): class Pcdi: def __init__(self, params, data, dir=None): - if 'pc_type' in params: - self.type = params['pc_type'] - else: - self.type = "LUCY" + self.type = params.get('pc_type', 'LUCY') if 'pc_LUCY_iterations' in params: self.iterations = params['pc_LUCY_iterations'] else: - self.iterations = 20 - if 'pc_normalize' in params: - self.normalize = params['pc_normalize'] - else: - self.normalize = True + print('pc_LUCY_iterations parameter not defined') + raise + self.normalize = params.get('pc_normalize', True) if 'pc_LUCY_kernel' in params: self.kernel_area = params['pc_LUCY_kernel'] else: - self.kernel_area = (16, 16, 16) + print('pc_LUCY_kernel parameter not defined') + raise if dir is None: self.kernel = None else: try: - self.kernel = devlib.load(dir + '/coherence.npy') + self.kernel = devlib.load(ut.join(dir, 'coherence.npy')) except: self.kernel = None self.dims = devlib.dims(data) - self.roi_data = dvut.crop_center(devlib.fftshift(data), self.kernel_area) + self.roi_data = ut.crop_center(devlib.fftshift(data), self.kernel_area) if self.normalize: self.sum_roi_data = devlib.sum(devlib.square(self.roi_data)) if self.kernel is None: self.kernel = devlib.full(self.kernel_area, 0.5, dtype=devlib.dtype(data)) def set_previous(self, abs_amplitudes): - self.roi_amplitudes_prev = dvut.crop_center(devlib.fftshift(abs_amplitudes), self.kernel_area) + self.roi_amplitudes_prev = ut.crop_center(devlib.fftshift(abs_amplitudes), self.kernel_area) def apply_partial_coherence(self, abs_amplitudes): abs_amplitudes_2 = devlib.square(abs_amplitudes) @@ -60,7 +57,7 @@ def apply_partial_coherence(self, abs_amplitudes): return converged def update_partial_coherence(self, abs_amplitudes): - roi_amplitudes = dvut.crop_center(devlib.fftshift(abs_amplitudes), self.kernel_area) + roi_amplitudes = ut.crop_center(devlib.fftshift(abs_amplitudes), self.kernel_area) roi_combined_amp = 2 * roi_amplitudes - self.roi_amplitudes_prev if self.normalize: amplitudes_2 = devlib.square(roi_combined_amp) @@ -86,12 +83,12 @@ def lucy_deconvolution(self, amplitudes, data, iterations): self.kernel = devlib.absolute(self.kernel) / coh_sum -class Feature(ABC): +class TriggeredOp(ABC): """ - Base class for features. This class creates feature objects and manages the trigger or subtrigger depending on + Base class for Triggered operation. This class creates feature objects and manages the trigger or subtrigger depending on configuration. """ - def __init__(self, sub_rows_trigs, params): + def __init__(self, trig_op_name): # the self.objs will hold either one concrete feature object if configuration is for general trigger, # i.e. spanning across all iterations, # or it will hold a list of objects, each being a concrete sub-feature, appearing in order that will be @@ -102,12 +99,8 @@ def __init__(self, sub_rows_trigs, params): # or to self.apply_trigger_seq if sub-triggers are configured self.f = None - # this function sets self.objs and self.f and creates all objects - self.create_objs(sub_rows_trigs, params) + self.trig_op_name = trig_op_name # ex: 'shrink_wrap_trigger' - self.key = None # ex: 'SW' - - self.trig_name = None # ex: 'shrink_wrap_trigger' def apply_trigger_obj(self, *args): """ @@ -149,7 +142,7 @@ def create_obj(self, *args): """ pass - def create_objs(self, sub_rows_trigs, params): + def create_objs(self, params, sub_rows_trigs): """ The params map contains value of the feature trigger. If the trigger is configured as a general one, i.e. a single trigger, the self.objs is set to the created instance of the feature, and the self.f @@ -162,14 +155,17 @@ def create_objs(self, sub_rows_trigs, params): :param params: configuration parameters :return: """ - if self.key in sub_rows_trigs.keys(): - row = sub_rows_trigs[self.key][0] - sub_trigs = sub_rows_trigs[self.key][1] + trigger_name = f'{self.trig_op_name}_trigger' + + if trigger_name in sub_rows_trigs.keys(): + row = sub_rows_trigs[trigger_name][0] + sub_trigs = sub_rows_trigs[trigger_name][1] sub_objs = {} for sub_t in sub_trigs: (beg, end, idx) = sub_t index = int(idx) - sub_objs[index] = self.create_obj(params, index=index, beg=beg, end=end) + if index not in sub_objs.keys(): + sub_objs[index] = self.create_obj(params, index=index, beg=beg, end=end) trigs = [i-1 for i in row.tolist() if i != 0] # the operation of creating object might fail # if conditions are not met, ex: the sw_type is not supported @@ -177,16 +173,14 @@ def create_objs(self, sub_rows_trigs, params): self.objs = [sub_objs[idx] for idx in trigs] self.f = self.apply_trigger_seq else: - if self.trig_name in params: + if trigger_name in params: self.objs = self.create_obj(params) self.f = self.apply_trigger_obj -class ShrinkWrap(Feature): - def __init__(self, params, sub_rows_trigs): - self.key = 'SW' - self.trig_name = 'shrink_wrap_trigger' - super().__init__(sub_rows_trigs, params) +class ShrinkWrap(TriggeredOp): + def __init__(self, trig_op): + super().__init__(trig_op) class GaussSW: def __init__(self, gauss_sigma, threshold): @@ -194,88 +188,131 @@ def __init__(self, gauss_sigma, threshold): self.threshold = threshold def apply_trigger(self, *args): -# print('shrink wrap, sigma', self.gauss_sigma) ds_image = args[0] return dvut.shrink_wrap(ds_image, self.threshold, self.gauss_sigma) + class Gauss1SW: + def __init__(self, gauss_sigma, threshold): + self.gauss_sigma = gauss_sigma + self.threshold = threshold + + def apply_trigger(self, *args): + ds_image = args[0] + return dvut.shrink_wrap(ds_image, self.threshold + .01, self.gauss_sigma) + + def create_obj(self, params, index=None, beg=None, end=None): + def check_Gauss_type(): + # for now cohere supports only Gauss type, so the following parameters are mandatory + if 'shrink_wrap_gauss_sigma' not in params: + print('shrink_wrap_gauss_sigma parameter not defined') + raise + if 'shrink_wrap_threshold' not in params: + print('shrink_wrap_threshold parameter not defined') + raise + + def check_Gauss1_type(): + # for now cohere supports only Gauss type, so the following parameters are mandatory + if 'shrink_wrap_gauss_sigma' not in params: + print('shrink_wrap_gauss_sigma parameter not defined') + raise + if 'shrink_wrap_threshold' not in params: + print('shrink_wrap_threshold parameter not defined') + raise + + if 'shrink_wrap_type' not in params: + print('shrink_wrap_type parameter not defined') + raise + if index is None: - if 'shrink_wrap_type' not in params: - sw_type = 'GAUSS' - else: - sw_type = params['shrink_wrap_type'] + sw_type = params['shrink_wrap_type'] if sw_type == 'GAUSS': - if 'shrink_wrap_gauss_sigma' not in params: - sigma = 1.0 - else: - sigma = params['shrink_wrap_gauss_sigma'] - if 'shrink_wrap_threshold' not in params: - threshold = 0.1 - else: - threshold = params['shrink_wrap_threshold'] - else: - if 'shrink_wrap_type' not in params or len(params['shrink_wrap_type']) - 1 < index: - sw_type = 'GAUSS' + check_Gauss_type() + sigma = params['shrink_wrap_gauss_sigma'] + threshold = params['shrink_wrap_threshold'] + return self.GaussSW(sigma, threshold) + elif sw_type == 'GAUSS1': + check_Gauss1_type() + sigma = params['shrink_wrap_gauss_sigma'] + threshold = params['shrink_wrap_threshold'] + return self.Gauss1SW(sigma, threshold) else: - sw_type = params['shrink_wrap_type'][index] + print(f'{sw_type} shrink wrap type is not supported') + raise + else: + if len(params['shrink_wrap_type']) - 1 < index: + print(f'shrink_wrap_type not defined for sub-trigger {index}') + raise + sw_type = params['shrink_wrap_type'][index] if sw_type == 'GAUSS': - if 'shrink_wrap_gauss_sigma' not in params or len(params['shrink_wrap_gauss_sigma']) - 1 < index: - sigma = 1.0 - else: - sigma = params['shrink_wrap_gauss_sigma'][index] - if 'shrink_wrap_threshold' not in params or len(params['shrink_wrap_threshold']) - 1 < index: - threshold = 0.1 - else: - threshold = params['shrink_wrap_threshold'][index] - return self.GaussSW(sigma, threshold) - - -class PhaseMod(Feature): - def __init__(self, params, sub_rows_trigs): - self.key = 'PHM' - self.trig_name = 'phm_trigger' - super().__init__(sub_rows_trigs, params) - - class PhasePHM: - def __init__(self, phm_phase_min, phm_phase_max): - self.phm_phase_min = phm_phase_min - self.phm_phase_max = phm_phase_max + check_Gauss_type() + if len(params['shrink_wrap_gauss_sigma']) - 1 < index: + print(f'shrink_wrap_gauss_sigma not defined for sub-trigger {index}') + raise + sigma = params['shrink_wrap_gauss_sigma'][index] + if len(params['shrink_wrap_threshold']) - 1 < index: + print(f'shrink_wrap_threshold not defined for sub-trigger {index}') + raise + threshold = params['shrink_wrap_threshold'][index] + return self.GaussSW(sigma, threshold) + elif sw_type == 'GAUSS1': + check_Gauss1_type() + if len(params['shrink_wrap_gauss_sigma']) - 1 < index: + print(f'shrink_wrap_gauss_sigma not defined for sub-trigger {index}') + raise + sigma = params['shrink_wrap_gauss_sigma'][index] + if len(params['shrink_wrap_threshold']) - 1 < index: + print(f'shrink_wrap_threshold not defined for sub-trigger {index}') + raise + threshold = params['shrink_wrap_threshold'][index] + return self.Gauss1SW(sigma, threshold) + else: + print(f'{sw_type} shrink wrap type is not supported') + raise + + + +class PhaseConstrain(TriggeredOp): + def __init__(self, trig_op): + super().__init__(trig_op) + + class PhasePHC: + def __init__(self, phc_phase_min, phc_phase_max): + self.phc_phase_min = phc_phase_min + self.phc_phase_max = phc_phase_max def apply_trigger(self, *args): -# print('phm trig, phase min', self.phm_phase_min) ds_image = args[0] phase = devlib.angle(ds_image) - return (phase > self.phm_phase_min) & (phase < self.phm_phase_max) + return (phase > self.phc_phase_min) & (phase < self.phc_phase_max) def create_obj(self, params, index=None, beg=None, end=None): + if 'phc_phase_min' not in params: + print('phc_phase_min parameter not defined') + raise + if 'phc_phase_max' not in params: + print('phc_phase_max parameter not defined') + raise if index is None: - if 'phm_phase_min' not in params: - phase_min = -1.57 - else: - phase_min = params['phm_phase_min'] - if 'phm_phase_max' not in params: - phase_max = 1.57 - else: - phase_max = params['phm_phase_max'] + phase_min = params['phc_phase_min'] + phase_max = params['phc_phase_max'] else: - if 'phm_phase_min' not in params or len(params['phm_phase_min']) - 1 < index: - phase_min = -1.57 - else: - phase_min = params['phm_phase_min'][index] - if 'phm_phase_max' not in params or len(params['phm_phase_max']) - 1 < index: - phase_max = 1.57 - else: - phase_max = params['phm_phase_max'][index] - return self.PhasePHM(phase_min, phase_max) + if len(params['phc_phase_min']) - 1 < index: + print(f'phc_phase_min not defined for sub-trigger {index}') + raise + phase_min = params['phc_phase_min'][index] + if len(params['phc_phase_max']) - 1 < index: + print(f'phc_phase_max not defined for sub-trigger {index}') + raise + phase_max = params['phc_phase_max'][index] + return self.PhasePHC(phase_min, phase_max) -class LowPassFilter(Feature): - def __init__(self, params, sub_rows_trigs): - self.key = 'LPF' - self.trig_name = 'lowpass_filter_trigger' - super().__init__(sub_rows_trigs, params) +class LowPassFilter(TriggeredOp): + def __init__(self, trig_op): + super().__init__(trig_op) class Filter: def __init__(self, filter_range, iter_range, threshold): @@ -292,30 +329,52 @@ def apply_trigger(self, *args): iter = args[1] ds_image = args[2] filter_sigma = self.filter_sigmas[iter - self.iter_offset] - # print('in apply_low_filter, filter_sigma', filter_sigma) return (devlib.gaussian_filter(data, filter_sigma), dvut.shrink_wrap(ds_image, self.threshold, 1/filter_sigma)) def create_obj(self, params, index=None, beg=None, end=None): + if 'lowpass_filter_sw_threshold' not in params: + print('lowpass_filter_sw_threshold parameter not defined') + raise + if 'lowpass_filter_range' not in params: + print('lowpass_filter_range parameter not defined') + raise # in case of a all_iteration configuration the lowpass filter # iteration can be read from lowpass_filter_trigger if index is None: iter_range = (params['lowpass_filter_trigger'][0], params['lowpass_filter_trigger'][2]) filter_range = params['lowpass_filter_range'] - if 'lowpass_filter_threshold' in params: - threshold = params['lowpass_filter_threshold'] - else: - threshold = .1 + threshold = params['lowpass_filter_sw_threshold'] # otherwise the iterations to apply this instance are given as arguments else: iter_range = (beg, end) + if len(params['lowpass_filter_range']) - 1 < index: + print(f'lowpass_filter_range not defined for sub-trigger {index}') + raise filter_range = params['lowpass_filter_range'][index] - if 'lowpass_filter_threshold' not in params or len(params['lowpass_filter_threshold']) - 1 < index: - threshold = .1 - else: - threshold = params['lowpass_filter_threshold'][index] + if len(params['lowpass_filter_sw_threshold']) - 1 < index: + print(f'lowpass_filter_sw_threshold not defined for sub-trigger {index}') + raise + threshold = params['lowpass_filter_sw_threshold'][index] if len(filter_range) == 1: filter_range.append(1.0) return self.Filter(filter_range, iter_range, threshold) +def create(trig_op, params, trig_op_info): + if trig_op == 'shrink_wrap': + to = ShrinkWrap(trig_op) + if trig_op == 'phc': + to = PhaseConstrain(trig_op) + if trig_op == 'lowpass_filter': + to = LowPassFilter(trig_op) + + # this function sets self.objs and self.f and creates all objects + try: + to.create_objs(params, trig_op_info) + return to + except: + return None + + + diff --git a/cohere_core/controller/mpi_cmd.py b/cohere_core/controller/mpi_cmd.py new file mode 100644 index 0000000..b4e7475 --- /dev/null +++ b/cohere_core/controller/mpi_cmd.py @@ -0,0 +1,45 @@ +import time +import os +import subprocess +import argparse + + +def run_with_mpi(ga_method, lib, conf_file, datafile, dir, devices, hostfile=None): + start_time = time.time() + if ga_method is None: + script = '/reconstruction_multi.py' + else: + script = '/reconstruction_GA.py' + + script = os.path.realpath(os.path.dirname(__file__)).replace(os.sep, '/') + script + if hostfile is None: + command = ['mpiexec', '-n', str(len(devices)), 'python', script, lib, conf_file, datafile, dir, str(devices)] + else: + command = ['mpiexec', '-n', str(len(devices)), '--hostfile', hostfile, 'python', script, lib, conf_file, datafile, dir, str(devices)] + + subprocess.run(command, capture_output=True) + + run_time = time.time() - start_time + if ga_method is None: + print(f'multiple reconstructions took {run_time} seconds') + else: # fast GA + print(f'GA reconstructions for directory {dir} took {run_time} seconds') + + +def main(): + import ast + + parser = argparse.ArgumentParser() + parser.add_argument("lib", help="lib") + parser.add_argument("conf_file", help="conf_file") + parser.add_argument("datafile", help="datafile") + parser.add_argument('dir', help='dir') + parser.add_argument('dev', help='dev') + + args = parser.parse_args() + dev = ast.literal_eval(args.dev) + run_with_mpi(args.lib, args.conf_file, args.datafile, args.dir, dev) + + +if __name__ == "__main__": + exit(main()) diff --git a/cohere_core/controller/op_flow.py b/cohere_core/controller/op_flow.py index 34e4926..8dbc3bb 100644 --- a/cohere_core/controller/op_flow.py +++ b/cohere_core/controller/op_flow.py @@ -1,47 +1,61 @@ import numpy as np import re -# -algs = {'ER': ('er', 'modulus'), - 'HIO': ('hio', 'modulus'), - 'ERpc': ('er', 'pc_modulus'), - 'HIOpc': ('hio', 'pc_modulus'), - 'SF' : ('new_alg', 'pc_modulus'), +# This dict maps the mnemonic used when defining algorithm sequence parameter to the four steps of +# phase retrieval algorithm functions. + +algs = {'ER': ('to_reciprocal_space', 'modulus', 'to_direct_space', 'er'), + 'HIO': ('to_reciprocal_space', 'modulus', 'to_direct_space', 'hio'), + 'ERpc': ('to_reciprocal_space', 'pc_modulus', 'to_direct_space', 'er'), + 'HIOpc': ('to_reciprocal_space', 'pc_modulus', 'to_direct_space', 'hio'), } -# this map keeps the names of triggers that can be configured as sub-trigger, i.e. be a trigger for the iteration span +# This map keeps the names of triggers that can be configured as sub-trigger, i.e. be a trigger for the iteration span # defined by preceding algorithm. The key is the trigger name and value is the mnemonic. The mnemonic is used in the # configuration. -sub_triggers = {'shrink_wrap_trigger': 'SW', - 'phm_trigger': 'PHM', - 'lowpass_filter_trigger': 'LPF'} - -def get_algs(): - return algs +sub_triggers = {'SW' : 'shrink_wrap_trigger', + 'PHC' : 'phc_trigger', + 'LPF' : 'lowpass_filter_trigger'} +# This list contains triggers that will be active at the last iteration defined by trigger, despite +# not being a trigger calculated by the step formula. +# It applies to sub-triggers, setting the last iteration to that of sub-trigger. +last_iter_op_triggers = ['progress_trigger', + 'switch_peaks_trigger', + 'switch_resampling_trigger'] def get_alg_rows(s, pc_conf_start): + """ + Parses algorithm sequence string into structures being: algorithm rows, and sub-trigger operations info. + + :param s: str + algorithm sequence + :param pc_conf_start: boolean or None + if None, no partial coherence is scheduled + if True, the configured partial coherence will be scheduled + if False, the partial coherence started ago (in GA case) and will continue here + :return: tuple + rows : ndarry + ndarray that depicts algorithms (modulus, pc_modulus, hio, er) operations + sub_rows : dict + dictionary with entries of k : v, where + k is the trigger name that is being configured as sub-triggers + v is a list of sub-trigger operations + iter_no : int + number of iterations + pc_start : None or int + starting iteration of partial coherence if any + """ seq = [] accum_iter = 0 def parse_entry(ent, accum_iter): + # parses elementary part of the algorithm sequence r_e = ent.split('*') seq.append([int(r_e[0]), r_e[1], accum_iter]) accum_iter += int(r_e[0]) return accum_iter - if pc_conf_start is None: # no pc in this - # this is kind of hackish, but will be replaced by sequence for each generation - if s =='ERpc': - s = 'ER' - if s == 'HIOpc': - s = 'HIO' - elif not pc_conf_start: # GA case, the coherence will start at first iteration - if s == 'ER': - s = 'ERpc' - if s == 'HIO': - s = 'HIOpc' - s = s.replace(' ', '') entries = s.split('+') i = 0 @@ -67,39 +81,66 @@ def parse_entry(ent, accum_iter): accum_iter = parse_entry(entry, accum_iter) i += 1 iter_no = sum([e[0] for e in seq]) - rows = {} + alg_rows = {} sub_rows = {} row = np.zeros(iter_no, dtype=int) fs = set([i for sub in algs.values() for i in sub]) for f in fs: - rows[f] = row.copy() - # for each possible subtrigger add entry - for f in sub_triggers.values(): - sub_rows[f] = [] + alg_rows[f] = row.copy() i = 0 pc_start = None for entry in seq: repeat = entry[0] funs = entry[1].split('.') if funs[0] not in algs: - return 'undefined algorithm ' + funs[0] + print(f'algorithm {funs[0]} is not defined in op_flow.py file, algs dict.') + raise + # the pc will not be executed if pc_conf_start is None + # this code will be revised after each generation has separate config + if pc_conf_start is None: + if funs[0].endswith('pc'): + funs[0] = funs[0][:-2] + elif not pc_conf_start: + if not funs[0].endswith('pc'): + funs[0] = funs[0] + 'pc' + row_keys = algs[funs[0]] for row_key in row_keys: - rows[row_key][i:i+repeat] = 1 - if 'pc' in row_key and pc_start is None: - pc_start = i + alg_rows[row_key][i:i+repeat] = 1 + if 'pc' in row_key and pc_start == None: + if pc_conf_start == True: + pc_start = i + else: + pc_start = 1 # find sub-triggers for row_key in funs[1:]: match = re.match(r"([A-Z]+)([0-9]+)", row_key, re.I) if match: - (feature, idx) = match.groups(0) - # sub_rows[trigs[feature]].append((entry[2], entry[0] + entry[2], idx)) - sub_rows[feature].append((entry[2], entry[0] + entry[2], idx)) + (trig_op, idx) = match.groups(0) + sub_t = sub_triggers[trig_op] + if trig_op not in sub_triggers.keys(): + print(f'the sub-trigger {trig_op} must be defined in op_flow.py file, sub_triggers dict.') + raise + if sub_t not in sub_rows: + sub_rows[sub_t] = [] + sub_rows[sub_t].append((entry[2], entry[0] + entry[2], idx)) i += repeat - return rows, sub_rows, iter_no, pc_start + + return alg_rows, sub_rows, iter_no, pc_start -def trigger_row(trig, iter_no, row=None): +def fill_trigger_row(trig, iter_no, last_trig, row=None): + """ + This functions creates ndarray that depicts triggered operations for a given trigger. + + :param trig: list + a list with 1, 2, or 3 elements defining trigger + :param iter_no: int + total number of iterations + :param row: ndarray + if given, the row will be used to fill the trigger + :return: + """ if row is None: row = np.zeros(iter_no, dtype=int) if len(trig) ==1: @@ -120,11 +161,56 @@ def trigger_row(trig, iter_no, row=None): trig_stop = iter_no for it in range(trig_start, trig_stop, step): row[it] = 1 + if last_trig: + row[trig_stop - 1] = 1 return row -def get_flow_arr(params, flow_items_list, curr_gen=None, first_run=False): - # get information about GA/pc from config_map +def fill_sub_trigger_row(sub_iters, sub_trigs, iter_no, last_trig): + """ + Based on iterations allocated to sub-triggers and sub-triggers definitions this functions + creates ndarray that depicts triggered operations. + + :param sub_iters: list + contains entry for each sub-trigger + the entry consisting of starting and ending iterations where the sub-trigger + is active and index value (+1) that specifies sub-trigger. + :param sub_trigs: list + list of sub-trigger, defined in configuration + :param iter_no: int + total number of iterations + :return: ndarray + array of int, value of zero meaning no trigger operation in this iteration + value greater than zero meaning the sub-trigger operation related to the value + will be triggered in this iteration + """ + # create array indicating triggered operation (1) or no action (0) along iterations + sub_trig_row = np.zeros(iter_no, dtype=int) + # create array indicating with index which sub-triggered operation may happen in the iterations + sub_trig_idx_row = np.zeros(iter_no, dtype=int) + # for each defined sub iteration chunk apply corresponding sub-trigger + for (b, e, idx) in sub_iters: + index = int(idx) + sub_trig_idx_row[b:e] = index + 1 + if len(sub_trigs) - 1 < index: + print('not enough entries in sub-trigger') + raise + trigger = sub_trigs[index].copy() + trigger[0] += b + if len(trigger) == 2: + trigger.append(e) + elif len(trigger) == 3: + trigger[2] = min(e, trigger[0] + trigger[2]) + sub_trig_row = fill_trigger_row(trigger, iter_no, last_trig, sub_trig_row) + + return sub_trig_row * sub_trig_idx_row + + +def get_flow_arr(params, flow_items_list, curr_gen=None): + # get information about GA and partial coherence from config_map + # pc_conf_start is None if partial coherence is inactive in this reconstruction + # it is True if the partial coherence starts in this generation + # and False if partial coherence started in previous generation and is continued if 'pc_interval' in params: if curr_gen is None: pc_conf_start = True @@ -137,96 +223,92 @@ def get_flow_arr(params, flow_items_list, curr_gen=None, first_run=False): pc_conf_start = False else: pc_conf_start = None + if pc_conf_start is None: + params.pop('pc_interval', None) - parsed_seq = get_alg_rows(params['algorithm_sequence'], pc_conf_start) - if type(parsed_seq) == str: - return None, None - alg_rows, sub_trig_rows, iter_no, pc_start = get_alg_rows(params['algorithm_sequence'], pc_conf_start) - flow_arr = np.zeros((len(flow_items_list), iter_no), dtype=int) + # parse algorithm sequence to get the algorithm rows and sub-triggers rows, number iterations, + # and partial coherence starting iteration + try: + (alg_rows, sub_iters, iter_no, pc_start) = get_alg_rows(params['algorithm_sequence'], pc_conf_start) + except: + return False, None, None - sub_feats = {} - is_res = False - # special case; during lowpass filter a lpf shrink wrap is applied, and shrink wrap suppressed + # do some checks to find if the sequence and configuration are runnable + # and special cases + + # this array is a mask blocking shrink wrap when low pass filter is active apply_sw_row = np.ones(iter_no, dtype=int) - # last iteration when lowpass filter is applied - last_lpf = 0 - for i in range(len(flow_items_list)): - flow_item = flow_items_list[i] - if flow_item == 'next' or flow_item == 'to_reciprocal_space' or flow_item == 'to_direct_space': + last_lpf = None + if 'lowpass_filter_trigger' in params: + if 'lowpass_filter_trigger' in sub_iters: + for b, e, i in sub_iters['lowpass_filter_trigger']: + apply_sw_row[b:e] = 0 + last_lpf = e + elif len(params['lowpass_filter_trigger']) < 2: + print('Low pass trigger misconfiguration error. This trigger should have upper bound.') + raise + elif params['lowpass_filter_trigger'][2] >= iter_no: + print('Low pass trigger misconfiguration error. The upper bound should be less than total iterations.') + raise + else: + apply_sw_row[params['lowpass_filter_trigger'][0]:params['lowpass_filter_trigger'][2]] = 0 + last_lpf = params['lowpass_filter_trigger'][2] + + if pc_start is not None: + if pc_start == 0: + print('partial coherence is configured in first iteration, allow several ER before') + raise + else: + pc_interval = params['pc_interval'] + params['pc_trigger'] = [pc_start, pc_interval] + + # initialize + sub_trig_op = {} + + # create empty array with the size of number of all functions by number of all iterations + flow_arr = np.zeros((len(flow_items_list), iter_no), dtype=int) + + # fill the flow array with ones if function should execute in iteration + for i, flow_item in enumerate(flow_items_list): + if flow_item == 'next': + # these functions are executed in each iteration flow_arr[i, :] = 1 - elif flow_item == 'lowpass_filter_trigger': - if first_run \ - and 'lowpass_filter_trigger' in params \ - and len(params['lowpass_filter_trigger']) == 3 \ - and type(params['lowpass_filter_trigger'][0]) == int: - flow_arr[i] = trigger_row(params['lowpass_filter_trigger'], iter_no) - apply_sw_row[params['lowpass_filter_trigger'][0]:params['lowpass_filter_trigger'][2]] = 0 - is_res = True - elif flow_item == 'reset_resolution': - if is_res: - flow_arr[i] = trigger_row([params['lowpass_filter_trigger'][-1],], iter_no) - if last_lpf > 0: - flow_arr[i][last_lpf] = 1 - elif flow_item == 'shrink_wrap_trigger': - if 'shrink_wrap_trigger' in params and type(params['shrink_wrap_trigger'][0]) == int: - flow_arr[i] = trigger_row(params['shrink_wrap_trigger'], iter_no) * apply_sw_row - elif flow_item == 'phm_trigger': - if first_run and 'phm_trigger' in params and type(params['phm_trigger'][0]) == int: - flow_arr[i] = trigger_row(params['phm_trigger'], iter_no) - elif flow_item == 'new_func_trigger': - if 'new_func_trigger' in [algs]: - flow_arr[i] = trigger_row(params['new_func_trigger'], iter_no) - elif flow_item == 'pc_trigger': - if pc_start is not None: - pc_interval = params['pc_interval'] - pc_trigger = [pc_start, pc_interval] - flow_arr[i] = trigger_row(pc_trigger, iter_no) - pc_row = i - elif flow_item == 'set_prev_pc_trigger': - if pc_start is not None: - flow_arr[i, : -1] = flow_arr[pc_row, 1:] elif flow_item in alg_rows.keys(): + # fill out the algorithm rows flow_arr[i] = alg_rows[flow_item] - elif flow_item == 'twin_trigger': - if first_run and 'twin_trigger' in params: - flow_arr[i] = trigger_row(params['twin_trigger'], iter_no) - elif flow_item == 'average_trigger': - if 'average_trigger' in params and curr_gen is not None and curr_gen == params['ga_generations'] -1: - flow_arr[i] = trigger_row(params['average_trigger'], iter_no) - elif flow_item == 'progress_trigger': - if 'progress_trigger' in params: - flow_arr[i] = trigger_row(params['progress_trigger'], iter_no) - flow_arr[i][-1] = 1 - elif flow_item == 'switch_peaks': - if 'switch_peak_trigger' in params: - flow_arr[i] = trigger_row(params['switch_peak_trigger'], iter_no) - flow_arr[i][-1] = 1 - - # Determine features based on sub-triggers - if flow_item in sub_triggers.keys(): - mne = sub_triggers[flow_item] - if len(sub_trig_rows[sub_triggers[flow_item]]): - sub_trig_row = np.zeros(iter_no, dtype=int) - sub_feats_row = np.zeros(iter_no, dtype=int) - for (b, e, idx) in sub_trig_rows[mne]: - index = int(idx) - sub_feats_row[b:e] = index + 1 - # special case for lowpass filter feature that suppresses shrink wrap - if flow_item == 'lowpass_filter_trigger': - apply_sw_row[b:e] = 0 - last_lpf = e - trigger = params[flow_item][index].copy() - trigger[0] += b - if len(trigger) == 2: - trigger.append(e) - elif len(trigger) == 3: - trigger[2] = min(e, trigger[0] + trigger[2]) - sub_trig_row = trigger_row(trigger, iter_no, sub_trig_row) - flow_arr[i] = sub_trig_row - sub_feats_row = sub_feats_row * sub_trig_row - # suppress the shrink wrap during lowpass filtering - if flow_item == 'shrink_wrap_trigger': - sub_feats_row *= apply_sw_row - flow_arr[i] *= apply_sw_row - sub_feats[mne] = (sub_feats_row, sub_trig_rows[mne]) - return pc_start is not None, flow_arr, sub_feats + elif flow_item.endswith('operation'): + # fill out trigger/sub-trigger operations rows + # The function name and associated trigger differ in prefix. + # the function name ends with 'operation', and trigger ends with 'trigger' + trigger_name = flow_item.replace('operation', 'trigger') + if trigger_name in params: + # set the switch last_trig if the trigger should end with operation + last_trig = trigger_name in last_iter_op_triggers + + # determined in algorithm sequence parsing if the triggered operation is configured + # with sub-triggers or trigger + if trigger_name in sub_iters.keys(): + try: + flow_arr[i] = fill_sub_trigger_row(sub_iters[trigger_name], params[trigger_name], iter_no, last_trig) + except: + flow_arr = None + break + # add entry to sub trigger operation dict with key of the trigger mnemonic + # and the value of a list with the row and sub triggers iterations chunks + sub_trig_op[trigger_name] = (flow_arr[i], sub_iters[trigger_name]) + else: + flow_arr[i] = fill_trigger_row(params[trigger_name], iter_no, last_trig) + + # handle special cases + if flow_item == 'shrink_wrap_operation': + # shrink wrap is blocked by low pass filter + flow_arr[i] = flow_arr[i] * apply_sw_row + elif flow_item == 'set_prev_pc' and pc_start is not None: + # set_prev_pc is executed one iteration before pc_trigger + pc_row = flow_items_list.index('pc_operation') + flow_arr[i, : -1] = flow_arr[pc_row, 1:] + elif flow_item == 'reset_resolution' and last_lpf is not None: + # reset low pass filter (i.e. data set to original) after the last LPF operation + flow_arr[i][last_lpf] = 1 + + return pc_start is not None, flow_arr, sub_trig_op diff --git a/cohere_core/controller/phasing.py b/cohere_core/controller/phasing.py index a7eb57a..dcfd9d7 100644 --- a/cohere_core/controller/phasing.py +++ b/cohere_core/controller/phasing.py @@ -14,11 +14,18 @@ """ +from pathlib import Path import time import os from math import pi import random import importlib +import tqdm + +import numpy as np +from matplotlib import pyplot as plt +import tifffile as tf + import cohere_core.utilities.dvc_utils as dvut import cohere_core.utilities.utils as ut import cohere_core.utilities.config_verifier as ver @@ -28,7 +35,10 @@ __author__ = "Barbara Frosik" __copyright__ = "Copyright (c) 2016, UChicago Argonne, LLC." __docformat__ = 'restructuredtext en' -__all__ = ['reconstruction', +__all__ = ['set_lib', + 'get_norm', + 'reconstruction', + 'Breeder', 'Rec'] @@ -43,31 +53,20 @@ def get_norm(arr): return devlib.sum(devlib.square(devlib.absolute(arr))) -class Support: - """ - Support class is a state holder for the support array. It initializes the support at creation time. - Features that modify support: shrink wrap, phase modifier, lowpass filter. - """ - def __init__(self, params, dims, dir=None): - # create or get initial support - if dir is None or not os.path.isfile(dir + '/support.npy'): - initial_support_area = params['initial_support_area'] - init_support = [] - for i in range(len(initial_support_area)): - if type(initial_support_area[0]) == int: - init_support.append(initial_support_area[i]) - else: - init_support.append(int(initial_support_area[i] * dims[i])) - center = devlib.full(init_support, 1) - self.support = dvut.pad_around(center, dims, 0) - else: - self.support = devlib.load(dir + '/support.npy') +class Breeder: + def __init__(self, alpha_dir): + self.alpha_dir = alpha_dir - def get_support(self): - return self.support + def breed(self, phaser, breed_mode, threshold, sigma): + if breed_mode != 'none': + try: + phaser.ds_image = dvut.breed(breed_mode, self.alpha_dir, phaser.ds_image) + phaser.support = dvut.shrink_wrap(phaser.ds_image, threshold, sigma) + except: + print('exception during breeding') + return -1 - def flip(self): - self.support = devlib.flip(self.support) + return 0 class Rec: @@ -85,39 +84,35 @@ class Rec: __all__ = [] def __init__(self, params, data_file): self.iter_functions = [self.next, - self.lowpass_filter_trigger, - self.reset_resolution, - self.shrink_wrap_trigger, - self.phm_trigger, - self.to_reciprocal_space, - self.new_func_trigger, - self.pc_trigger, - self.pc_modulus, - self.modulus, - self.set_prev_pc_trigger, - self.to_direct_space, - self.er, - self.hio, - self.new_alg, - self.twin_trigger, - self.average_trigger, - self.progress_trigger] - - if 'init_guess' not in params: - params['init_guess'] = 'random' - elif params['init_guess'] == 'AI_guess': + self.lowpass_filter_operation, + self.reset_resolution, + self.shrink_wrap_operation, + self.phc_operation, + self.to_reciprocal_space, + self.new_func_operation, + self.pc_operation, + self.pc_modulus, + self.modulus, + self.set_prev_pc, + self.to_direct_space, + self.er, + self.hio, + self.new_alg, + self.twin_operation, + self.average_operation, + self.progress_operation] + + params['init_guess'] = params.get('init_guess', 'random') + if params['init_guess'] == 'AI_guess': if 'AI_threshold' not in params: params['AI_threshold'] = params['shrink_wrap_threshold'] if 'AI_sigma' not in params: params['AI_sigma'] = params['shrink_wrap_gauss_sigma'] - if 'reconstructions' not in params: - params['reconstructions'] = 1 - if 'hio_beta' not in params: - params['hio_beta'] = 0.9 - if 'initial_support_area' not in params: - params['initial_support_area'] = (.5, .5, .5) - if 'twin_trigger' in params and 'twin_halves' not in params: - params['twin_halves'] = (0, 0) + params['reconstructions'] = params.get('reconstructions', 1) + params['hio_beta'] = params.get('hio_beta', 0.9) + params['initial_support_area'] = params.get('initial_support_area', (.5, .5, .5)) + if 'twin_trigger' in params: + params['twin_halves'] = params.get('twin_halves', (0, 0)) if 'pc_interval' in params and 'pc' in params['algorithm_sequence']: self.is_pcdi = True else: @@ -128,16 +123,21 @@ def __init__(self, params, data_file): self.ds_image = None self.need_save_data = False self.saved_data = None + self.er_iter = False # Indicates whether the last iteration done was ER, used in CoupledRec + def init_dev(self, device_id): os.environ["CUDA_DEVICE_ORDER"] = "PCI_BUS_ID" - self.dev = device_id - try: - devlib.set_device(device_id) - except Exception as e: - print(e) - print('may need to restart GUI') - return -1 + if device_id != -1: + self.dev = device_id + if device_id != -1: + try: + devlib.set_device(device_id) + except Exception as e: + print(e) + print('may need to restart GUI') + return -1 + if self.data_file.endswith('tif') or self.data_file.endswith('tiff'): try: data_np = ut.read_tif(self.data_file) @@ -165,64 +165,46 @@ def init_dev(self, device_id): return 0 - def fast_ga(self, worker_qin, worker_qout): - if self.params['low_resolution_generations'] > 0: - self.need_save_data = True - - functions_dict = {} - functions_dict[self.init_dev.__name__] = self.init_dev - functions_dict[self.init.__name__] = self.init - functions_dict[self.breed.__name__] = self.breed - functions_dict[self.iterate.__name__] = self.iterate - functions_dict[self.get_metric.__name__] = self.get_metric - functions_dict[self.save_res.__name__] = self.save_res - - done = False - while not done: - # cmd is a tuple containing name of the function, followed by arguments - cmd = worker_qin.get() - if cmd == 'done': - done = True - else: - if type(cmd) == str: - # the function does not have any arguments - try: - ret = functions_dict[cmd]() - except: - print('cought exception') - ret = -3 - else: - try: - # the function name is the first element in the cmd tuple, and the other elements are arguments - ret = functions_dict[cmd[0]](*cmd[1:]) - except: - print('cought exception') - ret = -3 - - worker_qout.put(ret) - def init(self, dir=None, alpha_dir=None, gen=None): - def create_feat_objects(params, feats): + def create_feat_objects(params, trig_op_info): if 'shrink_wrap_trigger' in params: - self.shrink_wrap_obj = ft.ShrinkWrap(params, feats) - if 'phm_trigger' in params: - self.phm_obj = ft.PhaseMod(params, feats) + self.shrink_wrap_obj = ft.create('shrink_wrap', params, trig_op_info) + if self.shrink_wrap_obj is None: + print('failed to create shrink wrap object') + return False + if 'phc_trigger' in params: + self.phc_obj = ft.create('phc', params, trig_op_info) + if self.phc_obj is None: + print('failed to create phase constrain object') + return False if 'lowpass_filter_trigger' in params: - self.lowpass_filter_obj = ft.LowPassFilter(params, feats) + self.lowpass_filter_obj = ft.create('lowpass_filter', params, trig_op_info) + if self.lowpass_filter_obj is None: + print('failed to create lowpass filter object') + return False + return True if self.ds_image is not None: first_run = False - elif dir is None or not os.path.isfile(dir + '/image.npy'): + elif dir is None or not os.path.isfile(ut.join(dir, 'image.npy')): self.ds_image = devlib.random(self.dims, dtype=self.data.dtype) first_run = True else: - self.ds_image = devlib.load(dir + '/image.npy') + self.ds_image = devlib.load(ut.join(dir, 'image.npy')) first_run = False + # When running GA the lowpass filter, phc, and twin triggers should be active only during first + # generation. The code below inactivates the triggers in subsequent generations. + # This will be removed in the future when each generation will have own configuration. + if not first_run: + self.params.pop('lowpass_filter_trigger', None) + self.params.pop('phc_trigger', None) + self.params.pop('twin_trigger', None) + self.flow_items_list = [f.__name__ for f in self.iter_functions] - self.is_pc, flow, feats = of.get_flow_arr(self.params, self.flow_items_list, gen, first_run) + self.is_pc, flow, feats = of.get_flow_arr(self.params, self.flow_items_list, gen) if flow is None: return -1 @@ -230,7 +212,7 @@ def create_feat_objects(params, feats): (op_no, self.iter_no) = flow.shape for i in range(self.iter_no): for j in range(op_no): - if flow[j, i] == 1: + if flow[j, i] > 0: self.flow.append(self.iter_functions[j]) self.aver = None @@ -239,11 +221,21 @@ def create_feat_objects(params, feats): self.gen = gen self.prev_dir = dir self.alpha_dir = alpha_dir - self.support_obj = Support(self.params, self.dims, dir) + + # create or get initial support + if dir is None or not os.path.isfile(ut.join(dir, 'support.npy')): + init_support = [int(self.params['initial_support_area'][i] * self.dims[i]) for i in range(len(self.dims))] + center = devlib.full(init_support, 1) + self.support = dvut.pad_around(center, self.dims, 0) + else: + self.support = devlib.load(ut.join(dir, 'support.npy')) + if self.is_pc: self.pc_obj = ft.Pcdi(self.params, self.data, dir) # create the object even if the feature inactive, it will be empty - create_feat_objects(self.params, feats) + # If successful, it will return True, otherwise False + if not create_feat_objects(self.params, feats): + return -1 # for the fast GA the data needs to be saved, as it would be changed by each lr generation # for non-fast GA the Rec object is created in each generation with the initial data @@ -268,21 +260,27 @@ def create_feat_objects(params, feats): # the line below are for testing to set the initial guess to support # self.ds_image = devlib.full(self.dims, 1.0) + 1j * devlib.full(self.dims, 1.0) - self.ds_image *= self.support_obj.get_support() + self.ds_image *= self.support return 0 - def breed(self): - breed_mode = self.params['ga_breed_modes'][self.gen] - if breed_mode != 'none': - try: - self.ds_image = dvut.breed(breed_mode, self.alpha_dir, self.ds_image) - self.support_obj.params = dvut.shrink_wrap(self.ds_image, self.params['ga_sw_thresholds'][self.gen], - self.params['ga_sw_gauss_sigmas'][self.gen]) - except: - print('exception during breeding') - return -1 - return 0 + def clean_breeder(self): + self.breeder = None + devlib.clean_default_mem() + + + def breed(self): + def breed_retry(): + breed_mode = self.params['ga_breed_modes'][self.gen] + threshold = self.params['ga_sw_thresholds'][self.gen] + sigma = self.params['ga_sw_gauss_sigmas'][self.gen] + return self.breeder.breed(self, breed_mode, threshold, sigma) + try: + return breed_retry() + except Exception as e: + # retry + self.breeder = Breeder(self.alpha_dir) + return breed_retry() def iterate(self): @@ -291,15 +289,16 @@ def iterate(self): try: for f in self.flow: f() - except: + except Exception as error: + print(error) return -1 - print('iterate took ', (time.time() - start_t), ' sec') - if devlib.hasnan(self.ds_image): print('reconstruction resulted in NaN') return -1 + print('iterate took ', (time.time() - start_t), ' sec') + if self.aver is not None: ratio = self.get_ratio(devlib.from_numpy(self.aver), devlib.absolute(self.ds_image)) self.ds_image *= ratio / self.aver_iter @@ -309,66 +308,65 @@ def iterate(self): return 0 - def save_res(self, save_dir, only_image=False): + # center image's center_of_mass and sync support + self.ds_image, self.support = dvut.center_sync(self.ds_image, self.support) + from array import array if not os.path.exists(save_dir): os.makedirs(save_dir) - devlib.save(save_dir + '/image', self.ds_image) + devlib.save(ut.join(save_dir, 'image'), self.ds_image) + ut.save_tif(devlib.to_numpy(self.ds_image), ut.join(save_dir, 'image.tif')) if only_image: return 0 - devlib.save(save_dir + '/support', self.support_obj.get_support()) + devlib.save(ut.join(save_dir, 'support'), self.support) if self.is_pc: - devlib.save(save_dir + '/coherence', self.pc_obj.kernel) + devlib.save(ut.join(save_dir, 'coherence'), self.pc_obj.kernel) errs = array('f', self.errs) - with open(save_dir + "/errors.txt", "w+") as err_f: + with open(ut.join(save_dir, 'errors.txt'), 'w+') as err_f: err_f.write('\n'.join(map(str, errs))) - devlib.save(save_dir + '/errors', errs) + devlib.save(ut.join(save_dir, 'errors'), errs) metric = dvut.all_metrics(self.ds_image, self.errs) - with open(save_dir + "/metrics.txt", "w+") as f: + with open(ut.join(save_dir, 'metrics.txt'), 'w+') as f: f.write(str(metric)) - # for key, value in metric.items(): - # f.write(key + ' : ' + str(value) + '\n') return 0 - def get_metric(self, metric_type): return dvut.all_metrics(self.ds_image, self.errs) - # return dvut.get_metric(self.ds_image, self.errs, metric_type) def next(self): self.iter = self.iter + 1 - def lowpass_filter_trigger(self): + def lowpass_filter_operation(self): args = (self.data, self.iter, self.ds_image) - (self.iter_data, self.support_obj.support) = self.lowpass_filter_obj.apply_trigger(*args) + (self.iter_data, self.support) = self.lowpass_filter_obj.apply_trigger(*args) def reset_resolution(self): self.iter_data = self.data - def shrink_wrap_trigger(self): + def shrink_wrap_operation(self): args = (self.ds_image,) - self.support_obj.support = self.shrink_wrap_obj.apply_trigger(*args) + self.support = self.shrink_wrap_obj.apply_trigger(*args) - def phm_trigger(self): + def phc_operation(self): args = (self.ds_image,) - self.support_obj.support *= self.phm_obj.apply_trigger(*args) + self.support *= self.phc_obj.apply_trigger(*args) def to_reciprocal_space(self): self.rs_amplitudes = devlib.ifft(self.ds_image) - def new_func_trigger(self): + def new_func_operation(self): self.params['new_param'] = 1 - print('in new_func_trigger, new_param', self.params['new_param']) + print(f'in new_func_trigger, new_param {self.params["new_param"]}') - def pc_trigger(self): + def pc_operation(self): self.pc_obj.update_partial_coherence(devlib.absolute(self.rs_amplitudes)) def pc_modulus(self): @@ -387,24 +385,25 @@ def modulus(self): self.errs.append(error) self.rs_amplitudes *= ratio - def set_prev_pc_trigger(self): + def set_prev_pc(self): self.pc_obj.set_previous(devlib.absolute(self.rs_amplitudes)) def to_direct_space(self): self.ds_image_raw = devlib.fft(self.rs_amplitudes) def er(self): - self.ds_image = self.ds_image_raw * self.support_obj.get_support() + self.er_iter = True + self.ds_image = self.ds_image_raw * self.support def hio(self): + self.er_iter = False combined_image = self.ds_image - self.ds_image_raw * self.params['hio_beta'] - support = self.support_obj.get_support() - self.ds_image = devlib.where((support > 0), self.ds_image_raw, combined_image) + self.ds_image = devlib.where((self.support > 0), self.ds_image_raw, combined_image) def new_alg(self): - self.ds_image = 2.0 * (self.ds_image_raw * self.support_obj.get_support()) - self.ds_image_raw + self.ds_image = 2.0 * (self.ds_image_raw * self.support) - self.ds_image_raw - def twin_trigger(self): + def twin_operation(self): # TODO this will work only for 3D array, but will the twin be used for 1D or 2D? # com = devlib.center_of_mass(devlib.absolute(self.ds_image)) # sft = [int(self.dims[i] / 2 - com[i]) for i in range(len(self.dims))] @@ -421,7 +420,7 @@ def twin_trigger(self): else: self.ds_image[:, : half_y, :] = 0 - def average_trigger(self): + def average_operation(self): if self.aver is None: self.aver = devlib.to_numpy(devlib.absolute(self.ds_image)) self.aver_iter = 1 @@ -429,35 +428,101 @@ def average_trigger(self): self.aver = self.aver + devlib.to_numpy(devlib.absolute(self.ds_image)) self.aver_iter += 1 - def progress_trigger(self): - print('------iter', self.iter, ' error ', self.errs[-1]) + def progress_operation(self): + print(f'------iter {self.iter} error {self.errs[-1]}') def get_ratio(self, divident, divisor): - divisor = devlib.where((divisor != 0.0), divisor, 1.0) - ratio = divident / divisor + ratio = devlib.where((divisor > 1e-9), divident / divisor, 0.0) return ratio + def make_binary(self): + self.support = devlib.absolute(self.support) + self.support = devlib.where(self.support >= 0.9 * devlib.amax(self.support), 1, 0).astype("?") + + +def resample(data, matrix, size=None, plot=False): + s = data.shape + corners = [ + [0, 0, 0], + [s[0], 0, 0], [0, s[1], 0], [0, 0, s[2]], + [0, s[1], s[2]], [s[0], 0, s[2]], [s[0], s[1], 0], + [s[0], s[1], s[2]] + ] + new_corners = [np.linalg.inv(matrix)@c for c in corners] + new_shape = np.ceil(np.max(new_corners, axis=0) - np.min(new_corners, axis=0)).astype(int) + pad = np.max((new_shape - s) // 2) + data = devlib.pad(data, np.max([pad, 0])) + mid_shape = np.array(data.shape) + + offset = (mid_shape / 2) - matrix @ (mid_shape / 2) + + # The phasing data as saved is a magnitude, which is not smooth everywhere (there are non-differentiable points + # where the complex amplitude crosses zero). However, the intensity (magnitude squared) is a smooth function. + # Thus, in order to allow a high-order spline, we take the square root of the interpolated square of the image. + data = devlib.affine_transform(devlib.square(data), devlib.array(matrix), order=1, offset=offset) + data = devlib.sqrt(devlib.clip(data, 0)) + if plot: + arr = devlib.to_numpy(data) + plt.imshow(arr[np.argmax(np.sum(arr, axis=(1, 2)))]) + plt.title(np.linalg.det(matrix)) + plt.show() + + if size is not None: + data = pad_to_cube(data, size) + return data + + +def pad_to_cube(data, size): + shp = np.array([size, size, size]) // 2 + # Pad the array to the largest dimensions + shp1 = np.array(data.shape) // 2 + pad = shp - shp1 + pad[pad < 0] = 0 + data = devlib.pad(data, [(pad[0], pad[0]), (pad[1], pad[1]), (pad[2], pad[2])]) + # Crop the array to the final dimensions + shp1 = np.array(data.shape) // 2 + start, end = shp1 - shp, shp1 + shp + data = data[start[0]:end[0], start[1]:end[1], start[2]:end[2]] + return data + class Peak: """ Holds parameters related to a peak. """ - def __init__(self, dir_ornt): - (self.dir, self.orientation) = dir_ornt - - def set_data(self, G_0): - import tifffile as tf + def __init__(self, peak_dir): + geometry = ut.read_config(f"{peak_dir}/geometry") + self.hkl = geometry["peak_hkl"] + self.norm = np.linalg.norm(self.hkl) + print(f"Loading peak {self.hkl}") - self.g_vec = devlib.array([0, * self.orientation]) * G_0 + # Geometry parameters + self.g_vec = devlib.array([0, * self.hkl]) * 2*pi/geometry["lattice"] self.gdotg = devlib.array(devlib.dot(self.g_vec, self.g_vec)) - - fn = self.dir + '/phasing_data/data.tif' - data_np = tf.imread(fn.replace(os.sep, '/')) - data = devlib.from_numpy(data_np) - - # in the formatted data the max is in the center, we want it in the corner, so do fft shift - self.data = devlib.fftshift(devlib.absolute(data)) + self.size = geometry["final_size"] + self.rs_lab_to_det = np.array(geometry["rmatrix"]) / geometry["rs_voxel_size"] + self.rs_det_to_lab = np.linalg.inv(self.rs_lab_to_det) + self.ds_lab_to_det = self.rs_det_to_lab.T + self.ds_det_to_lab = self.rs_lab_to_det.T + + datafile = (Path(peak_dir) / "phasing_data/data.tif").as_posix() + self.full_data = devlib.absolute(devlib.from_numpy(tf.imread(datafile))) + self.full_data = devlib.moveaxis(self.full_data, 0, 2) + self.full_data = devlib.moveaxis(self.full_data, 0, 1) + self.full_size = np.max(self.full_data.shape) + self.full_data = pad_to_cube(self.full_data, self.full_size) + self.data = pad_to_cube(self.full_data, self.size) + + # Resample the diffraction patterns into the lab reference frame. + self.res_data = resample(self.full_data, self.rs_det_to_lab, self.size) + tf.imwrite((Path(peak_dir) / "phasing_data/data_resampled.tif").as_posix(), devlib.to_numpy(self.res_data)) + # self.res_mask = resample(devlib.array(np.ones(self.full_data.shape)), self.det_to_lab, self.size) + + # Do the fftshift on all the arrays + self.full_data = devlib.fftshift(self.full_data) + self.data = devlib.fftshift(self.data) + self.res_data = devlib.fftshift(self.res_data) class CoupledRec(Rec): @@ -482,38 +547,54 @@ class CoupledRec(Rec): __author__ = "Nick Porter" __all__ = [] - def __init__(self, params, peak_dir_orient): + def __init__(self, params, peak_dirs): super().__init__(params, None) - self.iter_functions = self.iter_functions + [self.switch_peaks] - - if "switch_peak_trigger" not in params: - params["switch_peak_trigger"] = [0, 10] - if "mp_max_weight" not in params: - params["mp_max_weight"] = 0.9 - if "mp_taper" not in params: - params["mp_taper"] = 0.75 + self.iter_functions = self.iter_functions + [self.switch_resampling_operation, self.switch_peaks_operation] - self.peak_objs = [Peak(dir_ornt) for dir_ornt in peak_dir_orient] + self.params["switch_peak_trigger"] = self.params.get("switch_peak_trigger", [0, 5]) + self.params["mp_max_weight"] = self.params.get("mp_max_weight", 0.9) + self.params["mp_taper"] = self.params.get("mp_taper", 0.75) + self.params["calc_strain"] = self.params.get("calc_strain", True) + self.peak_dirs = peak_dirs def init_dev(self, device_id): - self.dev = device_id + os.environ["CUDA_DEVICE_ORDER"] = "PCI_BUS_ID" if device_id != -1: + self.dev = device_id + if device_id != -1: + try: + devlib.set_device(device_id) + except Exception as e: + print(e) + print('may need to restart GUI') + return -1 + + # Allows for a control peak to be left out of the reconstruction for comparison + all_peaks = [Peak(d) for d in self.peak_dirs] + if "control_peak" in self.params: + control_peak = self.params["control_peak"] + self.peak_objs = [p for p in all_peaks if p.hkl != control_peak] try: - devlib.set_device(device_id) - except Exception as e: - print(e) - print('may need to restart GUI') - return -1 + self.ctrl_peak = next(p for p in all_peaks if p.hkl == control_peak) + self.ctrl_error = [] + except StopIteration: + self.ctrl_peak = None + self.ctrl_error = None + else: + self.peak_objs = all_peaks + self.ctrl_peak = None + self.ctrl_error = None - G_0 = 2*pi/self.params["lattice_size"] - for or_obj in self.peak_objs: - or_obj.set_data(G_0) + # Fast resampling means using the resampled data to reconstruct the object in a single common basis + # When this is set to False, the reconstruction will need to be resampled each time the peak is switched + self.fast_resample = True + self.lab_frame = True self.num_peaks = len(self.peak_objs) self.pk = 0 # index in list of current peak being reconstructed - self.data = self.peak_objs[self.pk].data + self.data = self.peak_objs[self.pk].res_data self.iter_data = self.data self.dims = self.data.shape @@ -540,58 +621,177 @@ def init(self, img_dir=None, alpha_dir=None, gen=None): def save_res(self, save_dir): from array import array - - self.shared_image = self.shared_image * self.support_obj.get_support()[:, :, :, None] + # J = self.calc_jacobian() + # s_xx = J[:,:,:,0,0][:,:,:,None] + # s_yy = J[:,:,:,1,1][:,:,:,None] + # s_zz = J[:,:,:,2,2][:,:,:,None] + # s_xy = ((J[:,:,:,0,1] + J[:,:,:,1,0]) / 2)[:,:,:,None] + # s_yz = ((J[:,:,:,1,2] + J[:,:,:,2,1]) / 2)[:,:,:,None] + # s_zx = ((J[:,:,:,2,0] + J[:,:,:,0,2]) / 2)[:,:,:,None] + # + sup = self.support[:, :, :, None] + self.shared_image = self.shared_image * sup if not os.path.exists(save_dir): os.makedirs(save_dir) + + if self.params["calc_strain"]: + J = self.calc_jacobian() + s_xx = J[:,:,:,0,0][:,:,:,None] + s_yy = J[:,:,:,1,1][:,:,:,None] + s_zz = J[:,:,:,2,2][:,:,:,None] + s_xy = ((J[:,:,:,0,1] + J[:,:,:,1,0]) / 2)[:,:,:,None] + s_yz = ((J[:,:,:,1,2] + J[:,:,:,2,1]) / 2)[:,:,:,None] + s_zx = ((J[:,:,:,2,0] + J[:,:,:,0,2]) / 2)[:,:,:,None] + final_rec = devlib.concatenate((self.shared_image, s_xx, s_yy, s_zz, s_xy, s_yz, s_zx, sup), axis=-1) + else: + s_00 = devlib.zeros(sup.shape) + final_rec = devlib.concatenate((self.shared_image, s_00, s_00, s_00, s_00, s_00, s_00, sup), axis=-1) + + devlib.save(save_dir + "/reconstruction", final_rec) devlib.save(save_dir + "/image", self.shared_image) devlib.save(save_dir + "/shared_density", self.shared_image[:, :, :, 0]) devlib.save(save_dir + "/shared_u1", self.shared_image[:, :, :, 1]) devlib.save(save_dir + "/shared_u2", self.shared_image[:, :, :, 2]) devlib.save(save_dir + "/shared_u3", self.shared_image[:, :, :, 3]) - devlib.save(save_dir + '/support', self.support_obj.get_support()) + devlib.save(save_dir + '/support', self.support) errs = array('f', self.errs) - with open(save_dir + "/errors.txt", "w+") as err_f: + with open(ut.join(save_dir, 'errors.txt'), 'w+') as err_f: err_f.write('\n'.join(map(str, errs))) - devlib.save(save_dir + '/errors', errs) + devlib.save(ut.join(save_dir, 'errors'), errs) metric = dvut.all_metrics(self.ds_image, self.errs) - with open(save_dir + "/metrics.txt", "w+") as f: + with open(ut.join(save_dir, 'metrics.txt'), 'w+') as f: f.write(str(metric)) + if self.ctrl_error is not None: + ctrl_error = np.array(self.ctrl_error) + np.save(save_dir + '/ctrl_error', ctrl_error) + np.savetxt(save_dir + '/ctrl_error.txt', ctrl_error) return 0 - def switch_peaks(self): + def get_phase_gradient(self): + """Calculate the phase gradient for a given peak""" + # Project onto the peak + if self.fast_resample: + self.iter_data = self.peak_objs[self.pk].res_data + else: + self.iter_data = self.peak_objs[self.pk].data + + self.to_working_image() + + # Iterate a few times to settle + for _ in range(25): + self.to_reciprocal_space() + self.modulus() + self.to_direct_space() + self.er() + + # Calculate the continuous phase gradient using Hofmann's method + # (https://doi.org/10.1103/PhysRevMaterials.4.013801) + gradient = [] + voxel_size = self.params["ds_voxel_size"] + dx0, dy0, dz0 = devlib.gradient(devlib.angle(self.ds_image), voxel_size) + dx1, dy1, dz1 = devlib.gradient(devlib.angle(self.ds_image * devlib.exp(1j*pi/2)), voxel_size) + dx2, dy2, dz2 = devlib.gradient(devlib.angle(self.ds_image * devlib.exp(-1j*pi/2)), voxel_size) + for stack in ([dx0, dx1, dx2], [dy0, dy1, dy2], [dz0, dz1, dz2]): + stack = devlib.array(stack) + index = devlib.argmin(stack**2, axis=0) + gradient.append(devlib.take_along_axis(stack, index[None, :, :, :], axis=0).squeeze(axis=0)) + + return gradient + + def calc_jacobian(self): + grad_phi = [] + for i, peak in enumerate(self.peak_objs): + self.prev_pk = self.pk + self.pk = i + grad_phi.append(self.get_phase_gradient()) + grad_phi = devlib.array(grad_phi) + G = devlib.array([peak.g_vec[1:] for peak in self.peak_objs]) + grad_phi = devlib.moveaxis(devlib.moveaxis(grad_phi, 0, -1), 0, -1) + cond = devlib.where(self.support, True, False) + final_shape = (grad_phi.shape[0], grad_phi.shape[1], grad_phi.shape[2], 3, 3) + ind = devlib.moveaxis(devlib.indices(cond.shape), 0, -1).reshape((-1, 3))[cond.ravel()].get() + jacobian = devlib.zeros(final_shape) + print("Calculating strain...") + for i, j, k in tqdm.tqdm(ind): + jacobian[i, j, k] = devlib.lstsq(G, grad_phi[i, j, k])[0] + return jacobian + + def switch_peaks_operation(self): self.to_shared_image() + self.get_control_error() self.pk = random.choice([x for x in range(self.num_peaks) if x not in (self.pk,)]) - self.iter_data = self.peak_objs[self.pk].data - self.to_working_image(self.peak_objs[self.pk]) + if self.fast_resample: + self.iter_data = self.peak_objs[self.pk].res_data + else: + self.iter_data = self.peak_objs[self.pk].data + self.to_working_image() def to_shared_image(self): + if not self.fast_resample: + if self.lab_frame: + self.lab_frame = False + else: + self.shift_to_center() + self.ds_image = resample(self.ds_image, self.peak_objs[self.pk].ds_det_to_lab, self.dims[0]) + self.support = resample(self.support, self.peak_objs[self.pk].ds_det_to_lab, + self.dims[0], plot=True) + self.make_binary() + beta = self.proj_weight[self.iter] - curr_peak = self.peak_objs[self.pk] - old_image = (devlib.dot(self.shared_image, curr_peak.g_vec) / curr_peak.gdotg)[:, :, :, None] * \ - curr_peak.g_vec + devlib.dot(self.shared_image, self.rho_hat)[:, :, :, None] * self.rho_hat - new_image = (devlib.angle(self.ds_image) / curr_peak.gdotg)[:, :, :, None] * curr_peak.g_vec + \ - devlib.absolute(self.ds_image)[:, :, :, None] * self.rho_hat + pk = self.peak_objs[self.pk] + old_image = (devlib.dot(self.shared_image, pk.g_vec) / pk.gdotg)[:, :, :, None] * pk.g_vec + \ + devlib.dot(self.shared_image, self.rho_hat)[:, :, :, None] * self.rho_hat + new_image = (devlib.angle(self.ds_image) / pk.gdotg)[:, :, :, None] * pk.g_vec + \ + devlib.absolute(self.ds_image)[:, :, :, None] * self.rho_hat * pk.norm self.shared_image = self.shared_image + beta * (new_image - old_image) - - - def to_working_image(self, peak_obj): - phi = devlib.dot(self.shared_image, peak_obj.g_vec) - rho = self.shared_image[:, :, :, 0] - self.ds_image = rho * devlib.exp(1j*phi) - - def progress_trigger(self): - ornt = self.peak_objs[self.pk].orientation - print(f'| iter {self.iter:>4} ' - f'| [{ornt[0]:>2}, {ornt[1]:>2}, {ornt[2]:>2}] ' - f'| err {self.errs[-1]:0.6f} ' - f'| max {self.shared_image[:, :, :, 0].max():0.5g}' - ) + if self.er_iter: + self.shared_image = self.shared_image * self.support[:, :, :, None] + self.shift_to_center() + + def to_working_image(self): + phi = devlib.dot(self.shared_image, self.peak_objs[self.pk].g_vec) + rho = self.shared_image[:, :, :, 0] / self.peak_objs[self.pk].norm + self.ds_image = rho * devlib.exp(1j * phi) + if not self.fast_resample: + self.ds_image = resample(self.ds_image, self.peak_objs[self.pk].ds_lab_to_det, self.dims[0]) + self.support = resample(self.support, self.peak_objs[self.pk].ds_lab_to_det, + self.dims[0]) + self.make_binary() + + def get_control_error(self): + if self.ctrl_peak is None: + return + phi = devlib.dot(self.shared_image, self.peak_objs[self.pk].g_vec) + rho = self.shared_image[:, :, :, 0] / self.peak_objs[self.pk].norm + img = devlib.absolute(devlib.ifft(rho * devlib.exp(1j * phi))) + lin_hist = devlib.histogram2d(self.ctrl_peak.res_data, img, log=False) + nmi = devlib.calc_nmi(lin_hist).get() + ehd = devlib.calc_ehd(lin_hist).get() + log_hist = devlib.histogram2d(self.ctrl_peak.res_data, img, log=True) + lnmi = devlib.calc_nmi(log_hist).get() + lehd = devlib.calc_ehd(log_hist).get() + self.ctrl_error.append([nmi, lnmi, ehd, lehd]) + + def progress_operation(self): + ornt = self.peak_objs[self.pk].hkl + prg = f'| iter {self.iter:>4} ' \ + f'| [{ornt[0]:>2}, {ornt[1]:>2}, {ornt[2]:>2}] ' \ + f'| err {self.errs[-1]:0.6f} ' \ + f'| max {self.shared_image[:, :, :, 0].max():0.5f} ' + if self.ctrl_error is not None: + prg += f"| NMI {self.ctrl_error[-1][0]:0.6f} " + prg += f"| LNMI {self.ctrl_error[-1][1]:0.6f} " + prg += f"| EHD {self.ctrl_error[-1][2]:0.6f} " + prg += f"| LEHD {self.ctrl_error[-1][3]:0.6f} " + print(prg) + + def switch_resampling_operation(self): + self.fast_resample = False def get_density(self): return self.shared_image[:, :, :, 0] @@ -602,12 +802,15 @@ def get_distortion(self): def flip(self): self.shared_image = devlib.flip(self.shared_image, axis=(0, 1, 2)) self.shared_image[:, :, :, 1:] *= -1 - self.support_obj.flip() + self.support = devlib.flip(self.support, axis=(0, 1, 2)) - def shift_to_center(self, ind, cutoff=None): - shift_dist = -devlib.array(ind) + (self.dims[0]//2) - self.shared_image = devlib.shift(self.shared_image, shift_dist, axis=(0, 1, 2)) - self.support_obj.support = devlib.shift(self.support_obj.support, shift_dist) + def shift_to_center(self): + ind = devlib.center_of_mass(self.shared_image[:, :, :, 0]) + shift_dist = (self.dims[0]//2) - devlib.round(devlib.array(ind)) + shift_dist = devlib.to_numpy(shift_dist).tolist() + self.shared_image = devlib.roll(self.shared_image, shift_dist, axis=(0, 1, 2)) + self.ds_image = devlib.roll(self.ds_image, shift_dist, axis=(0, 1, 2)) + self.support = devlib.roll(self.support, shift_dist, axis=(0, 1, 2)) def reconstruction(datafile, **kwargs): @@ -651,11 +854,11 @@ def reconstruction(datafile, **kwargs): used to calculate the Gaussian filter initial_support_area : list If the values are fractional, the support area will be calculated by multiplying by the data array dimensions. The support will be set to 1s to this dimensions centered. - phm_trigger : list + phc_trigger : list defines when to update support array using the parameters below by applaying phase constrain. - phm_phase_min : float + phc_phase_min : float point with phase below this value will be removed from support area - phm_phase_max : float + phc_phase_max : float point with phase over this value will be removed from support area pc_interval : int defines iteration interval to update coherence. @@ -675,33 +878,39 @@ def reconstruction(datafile, **kwargs): defines when to apply averaging. Negative start means it is offset from the last iteration. progress_trigger : list of int defines when to print info on the console. The info includes current iteration and error. - + debug : boolean + if in debug mode the verifier will not stop the progress, only print message """ + debug = kwargs.get('debug', False) error_msg = ver.verify('config_rec', kwargs) if len(error_msg) > 0: print(error_msg) - return + if not debug: + return if not os.path.isfile(datafile): - print('no file found', datafile) + print(f'no file found {datafile}') return if 'reconstructions' in kwargs and kwargs['reconstructions'] > 1: - print('Use cohere_core-ui package to run multiple reconstructions and GA. Proceding with single reconstruction.') - if 'processing' in kwargs: - pkg = kwargs['processing'] - else: - pkg = 'auto' + print('Use cohere_core-ui package to run multiple reconstructions. Processing single reconstruction.') + pkg = kwargs.get('processing', 'auto') if pkg == 'auto': try: import cupy as cp pkg = 'cp' except: - pkg = 'np' + try: + import torch + pkg = 'torch' + except: + pkg = 'np' if pkg == 'cp': devlib = importlib.import_module('cohere_core.lib.cplib').cplib elif pkg == 'np': devlib = importlib.import_module('cohere_core.lib.nplib').nplib + if pkg == 'torch': + devlib = importlib.import_module('cohere_core.lib.cplib').torchlib else: print('supporting cp and np processing') return @@ -715,10 +924,7 @@ def reconstruction(datafile, **kwargs): if worker.init_dev(device[0]) < 0: return - if 'continue_dir' in kwargs: - continue_dir = kwargs['continue_dir'] - else: - continue_dir = None + continue_dir = kwargs.get('continue_dir') ret_code = worker.init(continue_dir) if ret_code < 0: return diff --git a/cohere_core/controller/reconstruction_GA.py b/cohere_core/controller/reconstruction_GA.py index de73d64..78d55ce 100644 --- a/cohere_core/controller/reconstruction_GA.py +++ b/cohere_core/controller/reconstruction_GA.py @@ -11,15 +11,15 @@ This module controls a reconstructions using genetic algorithm (GA). Refer to cohere_core-ui suite for use cases. The reconstruction can be started from GUI x or using command line scripts x. """ - +import argparse import numpy as np import os -import cohere_core.controller.reconstruction_multi as multi import cohere_core.utilities.utils as ut -from multiprocessing import Process, Queue -import shutil +import cohere_core.utilities.ga_utils as gaut import importlib import cohere_core.controller.phasing as calc +from mpi4py import MPI +import datetime __author__ = "Barbara Frosik" @@ -28,110 +28,8 @@ __all__ = ['reconstruction'] -class Tracing: - def __init__(self, reconstructions, pars, dir): - self.init_dirs = [] - self.report_tracing = [] - - if 'init_guess' not in pars: - pars['init_guess'] = 'random' - if pars['init_guess'] == 'continue': - continue_dir = pars['continue_dir'] - for sub in os.listdir(continue_dir): - image, support, coh = ut.read_results(continue_dir + '/' + sub + '/') - if image is not None: - self.init_dirs.append(continue_dir + '/' + sub) - self.report_tracing.append([continue_dir + '/' + sub]) - if len(self.init_dirs) < reconstructions: - for i in range(reconstructions - len(self.init_dirs)): - self.report_tracing.append(['random' + str(i)]) - self.init_dirs = self.init_dirs + (reconstructions - len(self.init_dirs)) * [None] - elif pars['init_guess'] == 'AI_guess': - import cohere_core.controller.AI_guess as ai - - self.report_tracing.append(['AI_guess']) - for i in range(reconstructions - 1): - self.report_tracing.append(['random' + str(i)]) - # ai_dir = ai.start_AI(pars, datafile, dir) - # if ai_dir is None: - # return - self.init_dirs = [dir + '/results_AI'] + (reconstructions - 1) * [None] - else: - for i in range(reconstructions): - self.init_dirs.append(None) - self.report_tracing.append(['random' + str(i)]) - - - def set_map(self, map): - self.map = map - - - def append_gen(self, gen_ranks): - for key in gen_ranks: - self.report_tracing[self.map[key]].append(gen_ranks[key]) - - - def pretty_format_results(self): - """ - Takes in a list of report traces and formats them into human-readable tables. - Performs data conversion in 1 pass and formatting in a second to determine - padding and spacing schemes. - - Parameters - ---------- - none - - Returns - ------- - report_output : str - a string containing the formatted report - """ - col_gap = 2 - - num_gens = len(self.report_tracing[0]) - 1 - fitnesses = list(self.report_tracing[0][1][1].keys()) - - report_table = [] - report_table.append(['start'] + [f'generation {i}' for i in range(num_gens)]) - report_table.append([''] * len(report_table[0])) - - data_col_width = 15 - start_col_width = 15 - for pop_data in self.report_tracing: - report_table.append([str(pop_data[0])] + [str(ind_data[0]) for ind_data in pop_data[1:]]) - start_col_width = max(len(pop_data[0]), start_col_width) - - for fit in fitnesses: - fit_row = [''] - for ind_data in pop_data[1:]: - data_out = f'{fit} : {ind_data[1][fit]}' - data_col_width = max(len(data_out), data_col_width) - fit_row.append(data_out) - report_table.append(fit_row) - report_table.append([''] * len(report_table[0])) - - report_str = '' - for row in report_table: - report_str += row[0].ljust(start_col_width + col_gap) - report_str += (' ' * col_gap).join([cell.ljust(data_col_width) for cell in row[1:]]) + '\n' - - return report_str - - - def save(self, save_dir): - try: - report_str = self.pretty_format_results() - except Exception as e: - print(f'WARNING: Report formatting failed due to {type(e)}: {e}! Falling back to raw formatting.') - report_str = '\n'.join([str(l) for l in self.report_tracing]) - - with open(save_dir + '/ranks.txt', 'w+') as rank_file: - rank_file.write(report_str) - rank_file.flush() - - def set_lib(pkg): - global dvclib + global devlib if pkg == 'cp': devlib = importlib.import_module('cohere_core.lib.cplib').cplib elif pkg == 'np': @@ -141,173 +39,7 @@ def set_lib(pkg): calc.set_lib(devlib) -def set_ga_defaults(pars): - if 'reconstructions' not in pars: - pars['reconstructions'] = 1 - - if 'ga_generations' not in pars: - pars['ga_generations'] = 1 - - if 'init_guess' not in pars: - pars['init_guess'] = 'random' - - # check if pc feature is on - if 'pc' in pars['algorithm_sequence'] and 'pc_interval' in pars: - if not 'ga_gen_pc_start' in pars: - pars['ga_gen_pc_start'] = 0 - pars['ga_gen_pc_start'] = min(pars['ga_gen_pc_start'], pars['ga_generations']-1) - - if 'ga_fast' not in pars: - pars['ga_fast'] = False - - if 'ga_metrics' not in pars: - metrics = ['chi'] * pars['ga_generations'] - else: - metrics = pars['ga_metrics'] - if len(metrics) == 1: - metrics = metrics * pars['ga_generations'] - elif len(metrics) < pars['ga_generations']: - metrics = metrics + ['chi'] * (pars['ga_generations'] - len(metrics)) - pars['ga_metrics'] = metrics - - ga_reconstructions = [] - if 'ga_cullings' in pars: - worst_remove_no = pars['ga_cullings'] - if len(worst_remove_no) < pars['ga_generations']: - worst_remove_no = worst_remove_no + [0] * (pars['ga_generations'] - len(worst_remove_no)) - else: - worst_remove_no = [0] * pars['ga_generations'] - pars['worst_remove_no'] = worst_remove_no - # calculate how many reconstructions should continue - reconstructions = pars['reconstructions'] - for culling in worst_remove_no: - reconstructions = reconstructions - culling - if reconstructions <= 0: - return 'culled down to 0 reconstructions, check configuration' - ga_reconstructions.append(reconstructions) - pars['ga_reconstructions'] = ga_reconstructions - - if 'shrink_wrap_threshold' in pars: - sw_threshold = pars['shrink_wrap_threshold'] - else: - sw_threshold = .1 - if 'ga_sw_thresholds' in pars: - ga_sw_thresholds = pars['ga_sw_thresholds'] - if len(ga_sw_thresholds) == 1: - ga_sw_thresholds = ga_sw_thresholds * pars['ga_generations'] - elif len(ga_sw_thresholds) < pars['ga_generations']: - ga_sw_thresholds = ga_sw_thresholds + [sw_threshold] * (pars['ga_generations'] - len(ga_sw_thresholds)) - else: - ga_sw_thresholds = [sw_threshold] * pars['ga_generations'] - pars['ga_sw_thresholds'] = ga_sw_thresholds - - if 'sw_gauss_sigma' in pars: - sw_gauss_sigma = pars['sw_gauss_sigma'] - else: - sw_gauss_sigma = .1 - if 'ga_sw_gauss_sigmas' in pars: - ga_sw_gauss_sigmas = pars['ga_sw_gauss_sigmas'] - if len(ga_sw_gauss_sigmas) == 1: - ga_sw_gauss_sigmas = ga_sw_gauss_sigmas * pars['ga_generations'] - elif len(pars['ga_sw_gauss_sigmas']) < pars['ga_generations']: - ga_sw_gauss_sigmas = ga_sw_gauss_sigmas + [sw_gauss_sigma] * (pars['ga_generations'] - len(ga_sw_gauss_sigmas)) - else: - ga_sw_gauss_sigmas = [sw_gauss_sigma] * pars['ga_generations'] - pars['ga_sw_gauss_sigmas'] = ga_sw_gauss_sigmas - - if 'ga_breed_modes' not in pars: - ga_breed_modes = ['sqrt_ab'] * pars['ga_generations'] - else: - ga_breed_modes = pars['ga_breed_modes'] - if len(ga_breed_modes) == 1: - ga_breed_modes = ga_breed_modes * pars['ga_generations'] - elif len(ga_breed_modes) < pars['ga_generations']: - ga_breed_modes = ga_breed_modes + ['sqrt_ab'] * (pars['ga_generations'] - len(ga_breed_modes)) - pars['ga_breed_modes'] = ga_breed_modes - - if 'ga_lpf_sigmas' in pars: - pars['low_resolution_generations'] = len(pars['ga_lpf_sigmas']) - else: - pars['low_resolution_generations'] = 0 - - if pars['low_resolution_generations'] > 0: - if 'low_resolution_alg' not in pars: - pars['low_resolution_alg'] = 'GAUSS' - - print() - - return pars - - -def order_dirs(tracing, dirs, evals, metric_type): - """ - Orders results in generation directory in subdirectories numbered from 0 and up, the best result stored in the '0' subdirectory. The ranking is done by numbers in evals list, which are the results of the generation's metric to the image array. - - Parameters - ---------- - tracing : object - Tracing object that keeps fields related to tracing - dirs : list - list of directories where the reconstruction results files are saved - evals : list - list of evaluations of the results saved in the directories matching dirs list. The evaluations are dict - : - metric_type : metric type to be applied for evaluation - - Returns - ------- - list : - a list of directories where results are saved ordered from best to worst - dict : - evaluations of the best results - """ - # ranks keeps indexes of reconstructions from best to worst - # for most of the metric types the minimum of the metric is best, but for - # 'summed_phase' and 'area' it is opposite, so reversing the order - metric_evals = [e[metric_type] for e in evals] - ranks = np.argsort(metric_evals).tolist() - if metric_type == 'summed_phase' or metric_type == 'area': - ranks.reverse() - best_metrics = evals[ranks[0]] - - # Add tracing for the generation results - gen_ranks = {} - for i in range(len(evals)): - gen_ranks[ranks[i]] = (i, evals[ranks[i]]) - tracing.append_gen(gen_ranks) - - # find how the directories based on ranking, map to the initial start - prev_map = tracing.map - map = {} - for i in range(len(ranks)): - prev_seq = int(os.path.basename(dirs[ranks[i]])) - inx = prev_map[prev_seq] - map[i] = inx - prev_map.clear() - tracing.set_map(map) - - # order directories by rank - rank_dirs = [] - # append "_" to each result directory name - for i in range(len(ranks)): - dest = dirs[ranks[i]] + '_' + str(i) - src = dirs[ranks[i]] - shutil.move(src, dest) - rank_dirs.append(dest) - - # remove the number preceding rank from each directory name, so the directories are numbered - # according to rank - current_dirs = [] - for dir in rank_dirs: - last_sub = os.path.basename(dir) - parent_dir = os.path.dirname(dir).replace(os.sep, '/') - dest = parent_dir + '/' + last_sub.split('_')[-1] - shutil.move(dir, dest) - current_dirs.append(dest) - return current_dirs, best_metrics - - -def order_processes(tracing, proc_metrics, metric_type): +def order_ranks(tracing, proc_metrics, metric_type): """ Orders results in generation directory in subdirectories numbered from 0 and up, the best result stored in the '0' subdirectory. The ranking is done by numbers in evals list, which are the results of the generation's metric to the image array. @@ -327,19 +59,20 @@ def order_processes(tracing, proc_metrics, metric_type): dict : evaluations of the best results """ - proc_eval = [(key, proc_metrics[key][metric_type]) for key in proc_metrics.keys()] - ranked_proc = sorted(proc_eval, key=lambda x: x[1]) + rank_eval = [(key, proc_metrics[key][metric_type]) for key in proc_metrics.keys()] + ranked_proc = sorted(rank_eval, key=lambda x: x[1]) - # ranks keeps indexes of reconstructions from best to worstpro + # ranks keeps indexes of reconstructions from best to worst # for most of the metric types the minimum of the metric is best, but for - # 'summed_phase' and 'area' it is oposite, so reversing the order + # 'summed_phase' and 'area' it is opposite, so reversing the order if metric_type == 'summed_phase' or metric_type == 'area': ranked_proc.reverse() gen_ranks = {} for i in range(len(ranked_proc)): - pid = ranked_proc[i][0] - gen_ranks[pid] = (i, proc_metrics[pid]) + rk = ranked_proc[i][0] + gen_ranks[rk] = (i, proc_metrics[rk]) tracing.append_gen(gen_ranks) + return ranked_proc, proc_metrics[ranked_proc[0][0]] @@ -350,11 +83,27 @@ def cull(lst, no_left): return lst[0:no_left] +def write_log(rank: int, msg: str) -> None: + """ + Use this to force writes for debugging. PBS sometimes doesn't flush + std* outputs. MPI faults clobber greedy flushing of default python + logs. + """ + with open(f'{rank}.log', 'a') as log_f: + log_f.write(f'{datetime.datetime.now()} | {msg}\n') + + def reconstruction(lib, conf_file, datafile, dir, devices): """ Controls reconstruction that employs genetic algorith (GA). - This script is typically started with cohere_core-ui helper functions. The 'init_guess' parameter in the configuration file defines whether it is a random guess, AI algorithm determined (one reconstruction, the rest random), or starting from some saved state. It will set the initial guess accordingly and start GA algorithm. It will run multiple reconstructions for each generation in a loop. After each generation the best reconstruction, alpha is identified, and used for breeding. For each generation the results will be saved in g_x subdirectory, where x is the generation number, in configured 'save_dir' parameter or in 'results_phasing' subdirectory if 'save_dir' is not defined. + This script is typically started with cohere_core-ui helper functions. The 'init_guess' parameter in + the configuration file defines whether it is a random guess, AI algorithm determined + (one reconstruction, the rest random), or starting from some saved state. It will set the initial guess + accordingly and start GA algorithm. It will run multiple reconstructions for each generation in a loop. + After each generation the best reconstruction, alpha is identified, and used for breeding. + For each generation the results will be saved in g_x subdirectory, where x is the generation number, + in configured 'save_dir' parameter or in 'results_phasing' subdirectory if 'save_dir' is not defined. Parameters ---------- @@ -376,206 +125,195 @@ def reconstruction(lib, conf_file, datafile, dir, devices): list of GPUs available for this reconstructions """ + def save_metric(metric, file_name): + with open(file_name.replace(os.sep, '/'), 'w+') as f: + f.truncate(0) + linesep = os.linesep + for key, value in metric.items(): + f.write(f'{key} = {str(value)}{linesep}') + + comm = MPI.COMM_WORLD + size = comm.Get_size() + rank = comm.Get_rank() + pars = ut.read_config(conf_file) - pars = set_ga_defaults(pars) + if 'save_dir' in pars: + save_dir = pars['save_dir'] + else: + # the config_rec might be an alternate configuration with a postfix that will be included in save_dir + filename = conf_file.split('/')[-1] + save_dir = ut.join(dir, filename.replace('config_rec', 'results_phasing')) + alpha_dir = ut.join(save_dir, 'alpha') - if pars['ga_generations'] < 2: - print("number of generations must be greater than 1") - return + if rank == 0: + # create alpha dir and placeholder for the alpha's metrics + if not os.path.isdir(save_dir): + os.mkdir(save_dir) + if not os.path.isdir(alpha_dir): + os.mkdir(alpha_dir) - if pars['reconstructions'] < 2: - print("GA not implemented for a single reconstruction") - return + comm.Barrier() + + pars = gaut.set_ga_defaults(pars) + + reconstructions = size - if pars['ga_fast']: - reconstructions = min(pars['reconstructions'], int(len(devices)/2)) - else: - reconstructions = pars['reconstructions'] - print('GA starting', reconstructions, 'reconstructions') if 'ga_cullings' in pars: cull_sum = sum(pars['ga_cullings']) - if reconstructions - cull_sum < 2: - print("At least two reconstructions should be left after culling. Number of starting reconstructions is", reconstructions, "but ga_cullings adds to", cull_sum) - return - if pars['init_guess'] == 'AI_guess': - # run AI part first - import cohere_core.controller.AI_guess as ai - ai_dir = ai.start_AI(pars, datafile, dir) - if ai_dir is None: - return + set_lib(lib) - if 'save_dir' in pars: - save_dir = pars['save_dir'] - else: - filename = conf_file.split('/')[-1] - save_dir = dir + '/' + filename.replace('config_rec', 'results_phasing') - - # create alpha dir and placeholder for the alpha's metrics - alpha_dir = save_dir + '/alpha' - if not os.path.isdir(save_dir): - os.mkdir(save_dir) - if not os.path.isdir(alpha_dir): - os.mkdir(alpha_dir) - - tracing = Tracing(reconstructions, pars, dir) - - if pars['ga_fast']: # the number of processes is the same as available GPUs (can be same GPU if can fit more recs) - set_lib(lib) - - workers = [] - processes = {} - tracing_map = {} - init_dirs = {} - - for i in range(reconstructions): - workers.append(calc.Rec(pars, datafile)) - worker_qin = Queue() - worker_qout = Queue() - process = Process(target=workers[i].fast_ga, args=(worker_qin, worker_qout)) - process.start() - processes[process.pid] = [worker_qin, worker_qout] - init_dirs[process.pid] = tracing.init_dirs[i] - tracing_map[process.pid] = i + worker = calc.Rec(pars, datafile) + if rank == 0: + active_ranks = list(np.linspace(0, size-1, num=size, dtype=int)) + tracing = gaut.Tracing(reconstructions, pars, dir) + tracing_map = {r : r for r in active_ranks} tracing.set_map(tracing_map) - def handle_cmd(): - bad_processes = [] - for pid in processes: - worker_qout = processes[pid][1] - ret = worker_qout.get() + active = True + success = True + for g in range(pars['ga_generations']): + was_active = active + if g == 0: + ret = worker.init_dev(devices[rank]) + if ret < 0: + print(f'rank {rank} failed initializing device {devices[rank]}') + active = False + if active: + if pars['init_guess'] == 'AI_guess' and rank == 0: + # run AI part first + import cohere_core.controller.AI_guess as ai + ai_dir = ai.start_AI(pars, datafile, dir) + if ai_dir is None: + ret = worker.init(None, alpha_dir, g) + else: + ret = worker.init(ai_dir, alpha_dir, g) + else: + ret = worker.init(None, alpha_dir, g) if ret < 0: - worker_qin = processes[pid][0] - worker_qin.put('done') - bad_processes.append(pid) - # remove bad processes from dict (in the future we may reuse them) - for pid in bad_processes: - processes.pop(pid) - - for g in range(pars['ga_generations']): - print ('starting generation', g) - if g == 0: - for pid in processes: - worker_qin = processes[pid][0] - worker_qin.put(('init_dev', devices.pop())) - handle_cmd() - for pid in processes: - worker_qin = processes[pid][0] - worker_qin.put(('init', init_dirs[pid], alpha_dir, g)) - handle_cmd() - else: - for pid in processes: - worker_qin = processes[pid][0] - worker_qin.put(('init', None, alpha_dir, g)) - handle_cmd() - if g > 0: - for pid in processes: - worker_qin = processes[pid][0] - worker_qin.put('breed') - handle_cmd() - - for pid in processes: - worker_qin = processes[pid][0] - worker_qin.put('iterate') - handle_cmd() - # get metric, i.e the goodness of reconstruction from each run - proc_metrics = {} + print('failed init, check config') + break + else: + if active: + ret = worker.init(None, alpha_dir, g) + if ret < 0: + print(f'rank {rank} reconstruction failed, check algorithm sequence and triggers in configuration') + active = False + + if active: + ret = worker.breed() + worker.clean_breeder() + worker.breeder = None + if ret < 0: + active = False + comm.Barrier() + + if active: + ret = worker.iterate() + if ret < 0: + print(f'reconstruction for rank {rank} failed during iterations') + active = False + + if active: metric_type = pars['ga_metrics'][g] - for pid in processes: - worker_qin = processes[pid][0] - worker_qin.put(('get_metric', metric_type)) - for pid in processes: - worker_qout = processes[pid][1] - metric = worker_qout.get() - proc_metrics[pid] = metric + metric = worker.get_metric(metric_type) + else: + metric = None + + # send metric and rank to rank 0 + if rank > 0 and was_active: + comm.send(metric, dest=0) + + # rank 0 receives metrics from ranks and creates metrics dict + if rank == 0: + metrics = {} + if active: + metrics[0] = metric + elif was_active: + active_ranks.remove(0) + to_remove = [] + for r in active_ranks: + if r != 0: + metric = comm.recv(source=r) + if metric is None: + to_remove.append(r) + else: + metrics[r] = metric + for r in to_remove: + active_ranks.remove(r) + + comm.Barrier() + + if rank == 0: # order processes by metric - proc_ranks, best_metrics = order_processes(tracing, proc_metrics, metric_type) + proc_ranks, best_metrics = order_ranks(tracing, metrics, metric_type) + proc_ranks = [p[0] for p in proc_ranks] # cull culled_proc_ranks = cull(proc_ranks, pars['ga_reconstructions'][g]) - # remove culled processes from list (in the future we may reuse them) - for i in range(len(culled_proc_ranks), len(proc_ranks)): - pid = proc_ranks[i][0] - worker_qin = processes[pid][0] - worker_qin.put('done') - processes.pop(pid) - + # send out the ordered active ranks + for r in active_ranks: + if r != 0: + comm.send(culled_proc_ranks, dest=r) + # remove culled processes from active list + to_remove = proc_ranks[len(culled_proc_ranks) : len(proc_ranks)] + for r in to_remove: + if r in active_ranks: + active_ranks.remove(r) + elif active: + culled_proc_ranks = comm.recv(source=0) + if rank not in culled_proc_ranks: + active = False + + # check if this process is the best + if active and rank == culled_proc_ranks[0]: # compare current alpha and previous. If previous is better, set it as alpha. - if (g == 0 - or - best_metrics[metric_type] >= alpha_metrics[metric_type] and - (metric_type == 'summed_phase' or metric_type == 'area') - or - best_metrics[metric_type] < alpha_metrics[metric_type] and - (metric_type == 'chi' or metric_type == 'sharpness')): - pid = culled_proc_ranks[0][0] - worker_qin = processes[pid][0] - worker_qin.put(('save_res', alpha_dir, True)) - worker_qout = processes[pid][1] - ret = worker_qout.get() - alpha_metrics = best_metrics - - # save results, we may modify it later to save only some - gen_save_dir = save_dir + '/g_' + str(g) - if g == pars['ga_generations'] -1: - for i in range(len(culled_proc_ranks)): - pid = culled_proc_ranks[i][0] - worker_qin = processes[pid][0] - worker_qin.put(('save_res', gen_save_dir + '/' + str(i))) - for pid in processes: - worker_qout = processes[pid][1] - ret = worker_qout.get() - if len(processes) == 0: - break - for pid in processes: - worker_qin = processes[pid][0] - worker_qin.put('done') - else: # not fast GA - q = Queue() - prev_dirs = tracing.init_dirs - tracing.set_map({i:i for i in range(len(prev_dirs))}) - rec = multi - for g in range(pars['ga_generations']): - # delete previous-previous generation - if g > 1: - shutil.rmtree(save_dir + '/g_' + str(g-2)) - print ('starting generation', g) - gen_save_dir = save_dir + '/g_' + str(g) - metric_type = pars['ga_metrics'][g] - reconstructions = len(prev_dirs) - p = Process(target=rec.multi_rec, args=(lib, gen_save_dir, devices, reconstructions, pars, datafile, prev_dirs, metric_type, g, alpha_dir, q)) - p.start() - p.join() - - results = q.get() - evals = [] - temp_dirs = [] - for r in results: - eval = r[0] - if eval is not None: - evals.append(eval) - temp_dirs.append(r[1]) - - # results are saved in a list of directories - save_dir - # it will be ranked, and moved to temporary ranked directories - current_dirs, best_metrics = order_dirs(tracing, temp_dirs, evals, metric_type) - reconstructions = pars['ga_reconstructions'][g] - current_dirs = cull(current_dirs, reconstructions) - prev_dirs = current_dirs + if g == 0: + worker.save_res(alpha_dir, True) + save_metric(metric, ut.join(alpha_dir, 'alpha_metric')) + else: + def is_best(this_metric, alpha_metric, type): + if type == 'chi' or type == 'sharpness': + return this_metric[type] < alpha_metric[type] + elif type == 'summed_phase' or type == 'area': + return this_metric[type] > alpha_metric[type] - # compare current alpha and previous. If previous is better, set it as alpha. - # no need toset alpha for last generation - if (g == 0 - or - best_metrics[metric_type] >= alpha_metrics[metric_type] and - (metric_type == 'summed_phase' or metric_type == 'area') - or - best_metrics[metric_type] < alpha_metrics[metric_type] and - (metric_type == 'chi' or metric_type == 'sharpness')): - shutil.copyfile(current_dirs[0] + '/image.npy', alpha_dir + '/image.npy') - alpha_metrics = best_metrics - # remove the previous gen - shutil.rmtree(save_dir + '/g_' + str(pars['ga_generations'] - 2)) - - tracing.save(save_dir) - - print('done gen') \ No newline at end of file + # read the previous alpha metric + alpha_metric = ut.read_config(ut.join(alpha_dir, 'alpha_metric')) + + if is_best(metric, alpha_metric, metric_type): + worker.save_res(alpha_dir, True) + save_metric(metric, ut.join(alpha_dir, 'alpha_metric')) + + # save results, we may modify it later to save only some + gen_save_dir = ut.join(save_dir, f'g_{str(g)}') + if g == pars['ga_generations'] -1: + if active: + worker.save_res(ut.join(gen_save_dir, str(culled_proc_ranks.index(rank)))) + + if not active: + worker = None + devlib.clean_default_mem() + + comm.Barrier() + + if rank == 0 and success: + tracing.save(save_dir) + + +def main(): + import ast + + parser = argparse.ArgumentParser() + parser.add_argument("lib", help="lib") + parser.add_argument("conf_file", help="conf_file") + parser.add_argument("datafile", help="datafile") + parser.add_argument('dir', help='dir') + parser.add_argument('dev', help='dev') + + args = parser.parse_args() + dev = ast.literal_eval(args.dev) + reconstruction(args.lib, args.conf_file, args.datafile, args.dir, dev) + + +if __name__ == "__main__": + exit(main()) \ No newline at end of file diff --git a/cohere_core/controller/reconstruction_coupled.py b/cohere_core/controller/reconstruction_coupled.py index 5da8994..36846d0 100644 --- a/cohere_core/controller/reconstruction_coupled.py +++ b/cohere_core/controller/reconstruction_coupled.py @@ -6,13 +6,12 @@ """ cohere_core.reconstruction_coupled -================================= +================================== This module controls a multipeak reconstruction process. Refer to cohere_core-ui suite for use cases. The reconstruction can be started from GUI or using command line scripts, see :ref:`use`. """ -import numpy as np import os import importlib import cohere_core.controller.phasing as calc @@ -37,23 +36,8 @@ def set_lib(pkg, ndim=None): def rec_process(lib, pars, peak_dirs, dev, continue_dir): set_lib(lib) - # It is assumed that the calling script uses peak_dirs containing - # peak orientation. It is parsed here, and passed to the CoupledRec constractor as list of - # touples (, ) - peak_dir_orient = [] - packed_orientations = [(str(o[0]) + str(o[1]) + str(o[2])) for o in pars['orientations']] - # find the directory that matches the packed orientation - for dir in peak_dirs: - found = False - i = 0 - while i < len(packed_orientations) and not found: - if dir.endswith(packed_orientations[i]): - found = True - peak_dir_orient.append((dir, pars['orientations'][i])) - i += 1 - print(peak_dir_orient) - - worker = calc.CoupledRec(pars, peak_dir_orient) + + worker = calc.CoupledRec(pars, peak_dirs) if worker.init_dev(dev[0]) < 0: return worker.init(continue_dir) @@ -62,7 +46,7 @@ def rec_process(lib, pars, peak_dirs, dev, continue_dir): if 'save_dir' in pars: save_dir = pars['save_dir'] else: - save_dir = os.path.dirname(peak_dirs[0]) + '/results_phasing' + save_dir = ut.join(os.path.dirname(peak_dirs[0]), 'results_phasing') worker.save_res(save_dir) diff --git a/cohere_core/controller/reconstruction_multi.py b/cohere_core/controller/reconstruction_multi.py index 328558f..86fe8aa 100644 --- a/cohere_core/controller/reconstruction_multi.py +++ b/cohere_core/controller/reconstruction_multi.py @@ -13,32 +13,18 @@ """ import os -import numpy as np +import sys +import argparse import importlib import cohere_core.utilities.utils as ut import cohere_core.controller.phasing as calc -from multiprocessing.pool import ThreadPool as Pool -from multiprocessing import Process -from functools import partial -import threading +from mpi4py import MPI __author__ = "Barbara Frosik" __copyright__ = "Copyright (c) 2016, UChicago Argonne, LLC." __docformat__ = 'restructuredtext en' -__all__ = ['multi_rec', - 'reconstruction'] - - -class Devices: - def __init__(self, devices): - self.devices = devices - self.index = 0 - - def assign_gpu(self): - thr = threading.current_thread() - thr.gpu = self.devices[self.index] - self.index = self.index + 1 +__all__ = ['reconstruction'] def set_lib(pkg): @@ -52,126 +38,15 @@ def set_lib(pkg): calc.set_lib(devlib) -def single_rec_process(metric_type, gen, alpha_dir, rec_attrs): - """ - This function runs a single reconstruction process. - - Parameters - ---------- - proc : str - string defining library, supported 'np' and 'cp' - - pars : Object - Params object containing parsed configuration - - data : numpy array - data array - - req_metric : str - defines metric that will be used if GA is utilized - - dirs : list - tuple of two elements: directory that contain results of previous run or None, and directory where the results of this processing will be saved - - Returns - ------- - metric : float - a calculated characteristic of the image array defined by the metric - """ - worker, prev_dir, save_dir = rec_attrs - thr = threading.current_thread() - if worker.init_dev(thr.gpu) < 0: - metric = None - else: - ret_code = worker.init(prev_dir, alpha_dir, gen) - if ret_code == 0: - if gen is not None and gen > 0: - worker.breed() - - ret_code = worker.iterate() - if ret_code == 0: - worker.save_res(save_dir) - metric = worker.get_metric(metric_type) - if ret_code < 0: # bad reconstruction - metric = None - - return [metric, save_dir] - - -def multi_rec(lib, save_dir, devices, no_recs, pars, datafile, prev_dirs, metric_type='chi', gen=None, alpha_dir=None, q=None): - """ - This function controls the multiple reconstructions. - - Parameters - ---------- - lib : str - library acronym to use for reconstruction. Supported: - np - to use numpy - cp - to use cupy - - save_dir : str - a directory where the subdirectories will be created to save all the results for multiple reconstructions - - devices : list - list of GPUs available for this reconstructions - - no_recs : int - number of reconstructions - - pars : dict - parameters for reconstruction - - datafie : str - name of file containing data for reconstruction - - previous_dirs : list - directories that hols results of previous reconstructions if it is continuation or None(s) - - metric_type : str - a metric defining algorithm by which to evaluate the quality of reconstructed array - - gen : int - which generation is the reconstruction for - - q : queue - if provided the results will be queued - - Returns - ------- - None - """ - results = [] - - def collect_result(result): - results.append(result) - - set_lib(lib) - - workers = [calc.Rec(pars, datafile) for _ in range(no_recs)] - dev_obj = Devices(devices) - iterable = [] - save_dirs = [] - - for i in range(len(workers)): - save_sub = save_dir + '/' + str(i) - save_dirs.append(save_sub) - iterable.append((workers[i], prev_dirs[i], save_sub)) - func = partial(single_rec_process, metric_type, gen, alpha_dir) - with Pool(processes=len(devices), initializer=dev_obj.assign_gpu, initargs=()) as pool: - pool.map_async(func, iterable, callback=collect_result) - pool.close() - pool.join() - pool.terminate() - - if q is not None: - q.put(results[0]) - - def reconstruction(lib, conf_file, datafile, dir, devices): """ Controls multiple reconstructions, the reconstructions run concurrently. - This script is typically started with cohere_core-ui helper functions. The 'init_guess' parameter in the configuration file defines whether guesses are random, or start from some saved states. It will set the initial guesses accordingly and start phasing process, running each reconstruction in separate thread. The results will be saved in configured 'save_dir' parameter or in 'results_phasing' subdirectory if 'save_dir' is not defined. + This script is typically started with cohere_core-ui helper functions. The 'init_guess' parameter in the + configuration file defines whether guesses are random, or start from some saved states. It will set the + initial guesses accordingly and start phasing process, running each reconstruction in separate thread. + The results will be saved in configured 'save_dir' parameter or in 'results_phasing' subdirectory if + 'save_dir' is not defined. Parameters ---------- @@ -193,36 +68,80 @@ def reconstruction(lib, conf_file, datafile, dir, devices): list of GPUs available for this reconstructions """ - pars = ut.read_config(conf_file) + comm = MPI.COMM_WORLD + rank = comm.Get_rank() - if 'reconstructions' in pars: - reconstructions = pars['reconstructions'] - else: - reconstructions = 1 + pars = ut.read_config(conf_file) - prev_dirs = [] - if 'init_guess' not in pars: - pars['init_guess'] = 'random' - if pars['init_guess'] == 'continue': - continue_dir = pars['continue_dir'] - for sub in os.listdir(continue_dir): - image, support, coh = ut.read_results(continue_dir + '/' + sub + '/') - if image is not None: - prev_dirs.append(continue_dir + '/' + sub) - if len(prev_dirs) < reconstructions: - prev_dirs = prev_dirs + (reconstructions - len(prev_dirs)) * [None] - elif pars['init_guess'] == 'AI_guess': - print('multiple reconstruction do not support AI_guess initial guess') - return - else: - for _ in range(reconstructions): - prev_dirs.append(None) if 'save_dir' in pars: save_dir = pars['save_dir'] else: + # the config_rec might be an alternate configuration with a postfix that will be included in save_dir filename = conf_file.split('/')[-1] - save_dir = dir + '/' + filename.replace('config_rec', 'results_phasing') + save_dir = ut.join(dir, filename.replace('config_rec', 'results_phasing')) + if rank == 0: + if not os.path.isdir(save_dir): + os.mkdir(save_dir) + + comm.Barrier() + + if 'init_guess' in pars and pars['init_guess'] == 'AI_guess': + print('multiple reconstruction do not support AI_guess initial guess') + return + + prev_dir = None + if 'init_guess' in pars and pars['init_guess'] == 'continue': + if not 'continue_dir' in pars: + print('continue directory must be defined if initial guess is continue') + return + + continue_dir = pars['continue_dir'] + if os.path.isdir(continue_dir): + prev_dirs = os.listdir(continue_dir) + if len(prev_dirs) > rank: + prev_dir = ut.join(continue_dir, prev_dirs[rank]) + if not os.path.isfile(ut.join(prev_dir, 'image.npy')): + prev_dir = None + + set_lib(lib) + + # assume for now the size can accommodate the no_recs + worker = calc.Rec(pars, datafile) + ret = worker.init_dev(devices[rank]) + if ret < 0: + print(f'rank {rank} failed initializing device {devices[rank]}') + return + + ret = worker.init(prev_dir) + if ret < 0: + print(f'rank {rank} failed, reconstruction failed, check algorithm sequence and triggers in configuration') + return + + ret = worker.iterate() + if ret < 0: + print(f'reconstruction for rank {rank} failed during iterations') + return + + save_sub = ut.join(save_dir, str(rank)) + if not os.path.isdir(save_sub): + os.mkdir(save_sub) + worker.save_res(save_sub) + + +def main(): + import ast + + parser = argparse.ArgumentParser() + parser.add_argument("lib", help="lib") + parser.add_argument("conf_file", help="conf_file") + parser.add_argument("datafile", help="datafile") + parser.add_argument('dir', help='dir') + parser.add_argument('dev', help='dev') + + args = parser.parse_args() + dev = ast.literal_eval(args.dev) + reconstruction(args.lib, args.conf_file, args.datafile, args.dir, dev) + - p = Process(target=multi_rec, args=(lib, save_dir, devices, reconstructions, pars, datafile, prev_dirs)) - p.start() - p.join() +if __name__ == "__main__": + exit(main()) diff --git a/cohere_core/controller/reconstruction_populous_GA.py b/cohere_core/controller/reconstruction_populous_GA.py new file mode 100644 index 0000000..451ac72 --- /dev/null +++ b/cohere_core/controller/reconstruction_populous_GA.py @@ -0,0 +1,384 @@ +# ######################################################################### +# Copyright (c) , UChicago Argonne, LLC. All rights reserved. # +# # +# See LICENSE file. # +# ######################################################################### + +""" +cohere_core.reconstruction_GA +============================= + +This module controls a reconstructions using genetic algorithm (GA). +Refer to cohere_core-ui suite for use cases. The reconstruction can be started from GUI x or using command line scripts x. +""" + +import numpy as np +import os +import cohere_core.utilities.utils as ut +import cohere_core.utilities.ga_utils as gaut +from multiprocessing import Queue +import shutil +import importlib +import cohere_core.controller.phasing as calc +from multiprocessing.pool import ThreadPool as Pool +from multiprocessing import Process +from functools import partial +import threading + + +__author__ = "Barbara Frosik" +__copyright__ = "Copyright (c) 2016, UChicago Argonne, LLC." +__docformat__ = 'restructuredtext en' +__all__ = ['reconstruction'] + + +class Devices: + def __init__(self, devices): + self.devices = devices + self.index = 0 + + def assign_gpu(self): + thr = threading.current_thread() + thr.gpu = self.devices[self.index] + self.index = self.index + 1 + + +def set_lib(pkg): + global dvclib + if pkg == 'cp': + devlib = importlib.import_module('cohere_core.lib.cplib').cplib + elif pkg == 'np': + devlib = importlib.import_module('cohere_core.lib.nplib').nplib + elif pkg == 'torch': + devlib = importlib.import_module('cohere_core.lib.torchlib').torchlib + calc.set_lib(devlib) + + +def single_rec_process(metric_type, gen, alpha_dir, rec_attrs): + """ + This function runs a single reconstruction process. + + Parameters + ---------- + proc : str + string defining library, supported 'np' and 'cp' + + pars : Object + Params object containing parsed configuration + + data : numpy array + data array + + req_metric : str + defines metric that will be used if GA is utilized + + dirs : list + tuple of two elements: directory that contain results of previous run or None, and directory where the + results of this processing will be saved + + Returns + ------- + metric : float + a calculated characteristic of the image array defined by the metric + """ + worker, prev_dir, save_dir = rec_attrs + thr = threading.current_thread() + if worker.init_dev(thr.gpu) < 0: + print(f'reconstruction failed, device not initialized to {thr.gpu}') + metric = None + else: + ret_code = worker.init(prev_dir, alpha_dir, gen) + if ret_code < 0: + print('reconstruction failed, check algorithm sequence and triggers in configuration') + metric = None + else: + if gen is not None and gen > 0: + worker.breed() + + ret_code = worker.iterate() + if ret_code == 0: + worker.save_res(save_dir) + metric = worker.get_metric(metric_type) + else: # bad reconstruction + print('reconstruction failed during iterations') + metric = None + + return [metric, save_dir] + + +def multi_rec(save_dir, devices, no_recs, pars, datafile, prev_dirs, metric_type='chi', gen=None, alpha_dir=None, q=None): + """ + This function controls the multiple reconstructions. + + Parameters + ---------- + lib : str + library acronym to use for reconstruction. Supported: + np - to use numpy + cp - to use cupy + + save_dir : str + a directory where the subdirectories will be created to save all the results for multiple reconstructions + + devices : list + list of GPUs available for this reconstructions + + no_recs : int + number of reconstructions + + pars : dict + parameters for reconstruction + + datafie : str + name of file containing data for reconstruction + + previous_dirs : list + directories that hols results of previous reconstructions if it is continuation or None(s) + + metric_type : str + a metric defining algorithm by which to evaluate the quality of reconstructed array + + gen : int + which generation is the reconstruction for + + q : queue + if provided the results will be queued + + Returns + ------- + None + """ + results = [] + + def collect_result(result): + results.append(result) + + workers = [calc.Rec(pars, datafile) for _ in range(no_recs)] + dev_obj = Devices(devices) + iterable = [] + save_dirs = [] + + for i in range(len(workers)): + save_sub = ut.join(save_dir, str(i)) + save_dirs.append(save_sub) + iterable.append((workers[i], prev_dirs[i], save_sub)) + func = partial(single_rec_process, metric_type, gen, alpha_dir) + with Pool(processes=len(devices), initializer=dev_obj.assign_gpu, initargs=()) as pool: + pool.map_async(func, iterable, callback=collect_result) + pool.close() + pool.join() + pool.terminate() + + if q is not None: + if len(results) == 0: + q.put('failed') + else: + q.put(results[0]) + + +def order_dirs(tracing, dirs, evals, metric_type): + """ + Orders results in generation directory in subdirectories numbered from 0 and up, the best result stored in the '0' subdirectory. The ranking is done by numbers in evals list, which are the results of the generation's metric to the image array. + + Parameters + ---------- + tracing : object + Tracing object that keeps fields related to tracing + dirs : list + list of directories where the reconstruction results files are saved + evals : list + list of evaluations of the results saved in the directories matching dirs list. The evaluations are dict + : + metric_type : metric type to be applied for evaluation + + Returns + ------- + list : + a list of directories where results are saved ordered from best to worst + dict : + evaluations of the best results + """ + # ranks keeps indexes of reconstructions from best to worst + # for most of the metric types the minimum of the metric is best, but for + # 'summed_phase' and 'area' it is opposite, so reversing the order + metric_evals = [e[metric_type] for e in evals] + ranks = np.argsort(metric_evals).tolist() + if metric_type == 'summed_phase' or metric_type == 'area': + ranks.reverse() + best_metrics = evals[ranks[0]] + + # Add tracing for the generation results + gen_ranks = {} + for i in range(len(evals)): + gen_ranks[ranks[i]] = (i, evals[ranks[i]]) + tracing.append_gen(gen_ranks) + + # find how the directories based on ranking, map to the initial start + prev_map = tracing.map + map = {} + for i in range(len(ranks)): + prev_seq = int(os.path.basename(dirs[ranks[i]])) + inx = prev_map[prev_seq] + map[i] = inx + prev_map.clear() + tracing.set_map(map) + + # order directories by rank + rank_dirs = [] + # append "_" to each result directory name + for i in range(len(ranks)): + dest = f'{dirs[ranks[i]]}_{str(i)}' + src = dirs[ranks[i]] + shutil.move(src, dest) + rank_dirs.append(dest) + + # remove the number preceding rank from each directory name, so the directories are numbered + # according to rank + current_dirs = [] + for dir in rank_dirs: + last_sub = os.path.basename(dir) + parent_dir = os.path.dirname(dir) + dest = ut.join(parent_dir, last_sub.split('_')[-1]) + shutil.move(dir, dest) + current_dirs.append(dest) + return current_dirs, best_metrics + + +def cull(lst, no_left): + if len(lst) <= no_left: + return lst + else: + return lst[0:no_left] + + +def reconstruction(lib, conf_file, datafile, dir, devices): + """ + Controls reconstruction that employs genetic algorith (GA). + + This script is typically started with cohere_core-ui helper functions. The 'init_guess' parameter in the + configuration file defines whether it is a random guess, AI algorithm determined (one reconstruction, + the rest random), or starting from some saved state. It will set the initial guess accordingly and start + GA algorithm. It will run multiple reconstructions for each generation in a loop. After each generation + the best reconstruction, alpha is identified, and used for breeding. For each generation the results will + be saved in g_x subdirectory, where x is the generation number, in configured 'save_dir' parameter or in + 'results_phasing' subdirectory if 'save_dir' is not defined. + + Parameters + ---------- + lib : str + library acronym to use for reconstruction. Supported: + np - to use numpy, + cp - to use cupy, + + conf_file : str + configuration file name + + datafile : str + data file name + + dir : str + a parent directory that holds the reconstructions. For example experiment directory or scan directory. + + devices : list + list of GPUs available for this reconstructions + + """ + pars = ut.read_config(conf_file) + pars = gaut.set_ga_defaults(pars) + + if pars['ga_generations'] < 2: + print("number of generations must be greater than 1") + return + + if pars['reconstructions'] < 2: + print("GA not implemented for a single reconstruction") + return + + if pars['ga_fast']: + reconstructions = min(pars['reconstructions'], len(devices)) + else: + reconstructions = pars['reconstructions'] + + if 'ga_cullings' in pars: + cull_sum = sum(pars['ga_cullings']) + if reconstructions - cull_sum < 2: + print("At least two reconstructions should be left after culling. Number of starting reconstructions is", reconstructions, "but ga_cullings adds to", cull_sum) + return + + if pars['init_guess'] == 'AI_guess': + # run AI part first + import cohere_core.controller.AI_guess as ai + ai_dir = ai.start_AI(pars, datafile, dir) + if ai_dir is None: + return + + set_lib(lib) + + if 'save_dir' in pars: + save_dir = pars['save_dir'] + else: + # the config_rec might be an alternate configuration with a postfix that will be included in save_dir + filename = conf_file.split('/')[-1] + save_dir = ut.join(dir, filename.replace('config_rec', 'results_phasing')) + + # create alpha dir and placeholder for the alpha's metrics + alpha_dir = ut.join(save_dir, 'alpha') + if not os.path.isdir(save_dir): + os.mkdir(save_dir) + if not os.path.isdir(alpha_dir): + os.mkdir(alpha_dir) + + tracing = gaut.Tracing(reconstructions, pars, dir) + + q = Queue() + prev_dirs = tracing.init_dirs + tracing.set_map({i:i for i in range(len(prev_dirs))}) + for g in range(pars['ga_generations']): + # delete previous-previous generation + if g > 1: + shutil.rmtree(ut.join(save_dir, f'g_{str(g-2)}')) + print ('starting generation', g) + gen_save_dir = ut.join(save_dir, f'g_{str(g)}') + metric_type = pars['ga_metrics'][g] + reconstructions = len(prev_dirs) + p = Process(target=multi_rec, args=(gen_save_dir, devices, reconstructions, pars, datafile, prev_dirs, metric_type, g, alpha_dir, q)) + p.start() + p.join() + + results = q.get() + if results == 'failed': + print('reconstruction failed') + return + + evals = [] + temp_dirs = [] + for r in results: + eval = r[0] + if eval is not None: + evals.append(eval) + temp_dirs.append(r[1]) + + # results are saved in a list of directories - save_dir + # it will be ranked, and moved to temporary ranked directories + current_dirs, best_metrics = order_dirs(tracing, temp_dirs, evals, metric_type) + reconstructions = pars['ga_reconstructions'][g] + current_dirs = cull(current_dirs, reconstructions) + prev_dirs = current_dirs + + # compare current alpha and previous. If previous is better, set it as alpha. + # no need to set alpha for last generation + if (g == 0 + or + best_metrics[metric_type] >= alpha_metrics[metric_type] and + (metric_type == 'summed_phase' or metric_type == 'area') + or + best_metrics[metric_type] < alpha_metrics[metric_type] and + (metric_type == 'chi' or metric_type == 'sharpness')): + shutil.copyfile(ut.join(current_dirs[0], 'image.npy'), ut.join(alpha_dir, 'image.npy')) + alpha_metrics = best_metrics + # remove the previous gen + shutil.rmtree(ut.join(save_dir, 'g_' + str(pars['ga_generations'] - 2))) + + tracing.save(save_dir) + + print('done gen') \ No newline at end of file diff --git a/cohere_core/controller/reconstruction_single.py b/cohere_core/controller/reconstruction_single.py index 3823261..642fcbb 100644 --- a/cohere_core/controller/reconstruction_single.py +++ b/cohere_core/controller/reconstruction_single.py @@ -47,21 +47,32 @@ def rec_process(lib, pars, datafile, dev, continue_dir, save_dir): else: device = dev[0] if worker.init_dev(device) < 0: + print (f'reconstruction failed, device not initialized to {device}') return ret_code = worker.init(continue_dir) if ret_code < 0: + print ('reconstruction failed, check algorithm sequence and triggers in configuration') return + ret_code = worker.iterate() - if ret_code == 0: - worker.save_res(save_dir) + if ret_code < 0: + print ('reconstruction failed during iterations') + return + + worker.save_res(save_dir) + return def reconstruction(lib, conf_file, datafile, dir, dev=None): """ Controls single reconstruction. - This script is typically started with cohere-ui helper functions. The 'init_guess' parameter in the configuration file defines whether it is a random guess, AI algorithm determined, or starting from some saved state. It will set the initial guess accordingly and start phasing process. The results will be saved in configured 'save_dir' parameter or in 'results_phasing' subdirectory if 'save_dir' is not defined. + This script is typically started with cohere-ui helper functions. The 'init_guess' parameter in the + configuration file defines whether it is a random guess, AI algorithm determined, or starting from + some saved state. It will set the initial guess accordingly and start phasing process. The results + will be saved in configured 'save_dir' parameter or in 'results_phasing' subdirectory if 'save_dir' + is not defined. Parameters ---------- @@ -85,8 +96,7 @@ def reconstruction(lib, conf_file, datafile, dir, dev=None): """ pars = ut.read_config(conf_file) - if 'init_guess' not in pars: - pars['init_guess'] = 'random' + pars['init_guess'] = pars.get('init_guess', 'random') if pars['init_guess'] == 'continue': continue_dir = pars['continue_dir'] elif pars['init_guess'] == 'AI_guess': @@ -102,7 +112,7 @@ def reconstruction(lib, conf_file, datafile, dir, dev=None): save_dir = pars['save_dir'] else: filename = conf_file.split('/')[-1] - save_dir = dir + '/' + filename.replace('config_rec', 'results_phasing') + save_dir = ut.join(dir, filename.replace('config_rec', 'results_phasing')) p = Process(target=rec_process, args=(lib, pars, datafile, dev, continue_dir, save_dir)) diff --git a/cohere_core/data/alien_tools.py b/cohere_core/data/alien_tools.py index 1912cd9..9caf73c 100644 --- a/cohere_core/data/alien_tools.py +++ b/cohere_core/data/alien_tools.py @@ -110,7 +110,6 @@ def analyze_clusters(arr, labels, nz): # print("processing labels") # loop over the labels (clusters). label_counts[0] is the unique labels for n in range(1, nlabels): - # print(" %i %i "%(label_counts[0][n],label_counts[1][n]), end='\r') # the nth label from the first array of the label_counts tuple n_lab = label_counts[0][n] # the indicies of the points belonging to label n @@ -122,8 +121,6 @@ def analyze_clusters(arr, labels, nz): cluster_avg[cluspts] = np.sum(arr[cluspts]) / cluspts[0].size # compute average asym of each cluster and store in array. cluster_avg_asym[cluspts] = np.sum(asymmetry[cluspts]) / cluspts[0].size - # print(" %i %i %f %f "%(label_counts[0][n],label_counts[1][n],np.sum(asymmetry[cluspts]),cluspts[0].size), end='\n') - # print("largest clus size", cluster_size.max()) # compute relative cluster sizes to largest (main) cluster. rel_cluster_size = cluster_size / cluster_size.max() @@ -184,7 +181,7 @@ def save_arr(arr, dir, fname): """ if dir is not None: - full_name = dir + '/' + fname + full_name = ut.join(dir, fname) else: full_name = fname # save in the current dir tif.imsave(full_name, arr.transpose().astype(np.float32)) @@ -251,39 +248,21 @@ def auto_alien1(data, config, data_dir=None): data array with removed aliens """ data_dir = data_dir.replace(os.sep, '/') - if 'AA1_size_threshold' in config: - size_threshold = config['AA1_size_threshold'] - else: - size_threshold = 0.01 - if 'AA1_asym_threshold' in config: - asym_threshold = config['AA1_asym_threshold'] - else: - asym_threshold = 1.75 - if 'AA1_min_pts' in config: - min_pts = config['AA1_min_pts'] - else: - min_pts = 5 - if 'AA1_eps' in config: - eps = config['AA1_eps'] - else: - eps = 1.1 - if 'AA1_amp_threshold' in config: - threshold = config['AA1_amp_threshold'] - else: - threshold = 6 + size_threshold = config.get('AA1_size_threshold', 0.01) + asym_threshold = config.get('AA1_asym_threshold', 1.75) + min_pts = config.get('AA1_min_pts', 5) + eps = config.get('AA1_eps', 1.1) + threshold = config.get('AA1_amp_threshold', 6) if 'AA1_save_arrs' in config: save_arrs = config['AA1_save_arrs'] if save_arrs: - save_dir = data_dir + '/alien_analysis' + save_dir = ut.join(data_dir, 'alien_analysis') if not os.path.exists(save_dir): os.makedirs(save_dir) else: save_arrs = False - if 'AA1_expandcleanedsigma' in config: - expandcleanedsig = config['AA1_expandcleanedsigma'] - else: - expandcleanedsig = 0.0 + expandcleanedsig = config.get('AA1_expandcleanedsigma', 0.0) cuboid = crop_center(data) cuboid = np.where(cuboid >= threshold, cuboid, 0) @@ -323,7 +302,7 @@ def auto_alien1(data, config, data_dir=None): if (expandcleanedsig > 0): cuboid = np.where(cuboid > 0, 1.0, 0.0) sig = [expandcleanedsig, expandcleanedsig, 1.0] - cuboid = ut.gauss_conv_fft(cuboid, sig) + cuboid = np.gaussian_filter(cuboid, sig) no_thresh_cuboid = crop_center(data) cuboid = np.where(cuboid > 0.1, no_thresh_cuboid, 0.0) return cuboid @@ -382,7 +361,7 @@ def filter_aliens(data, config_map): return data = np.where((mask == 1), data, 0.0) else: - print('alien file does not exist ', alien_file) + print(f'alien file does not exist {alien_file}') else: print('alien_file parameter not configured') return data diff --git a/cohere_core/data/auto_prep.py b/cohere_core/data/auto_prep.py deleted file mode 100644 index e69de29..0000000 diff --git a/cohere_core/data/standard_preprocess.py b/cohere_core/data/standard_preprocess.py index 597022a..1fb137c 100644 --- a/cohere_core/data/standard_preprocess.py +++ b/cohere_core/data/standard_preprocess.py @@ -23,7 +23,7 @@ ] -def prep(datafile, **kwargs): +def prep(beamline_full_datafile_name, auto, **kwargs): """ This function formats data for reconstruction and saves it in data.tif file. The preparation consists of the following steps: - removing the alien: aliens are areas that are effect of interference. The area is manually set in a configuration file after inspecting the data. It could be also a mask file of the same dimensions that data. Another option is AutoAlien1 algorithm that automatically removes the aliens. @@ -35,11 +35,11 @@ def prep(datafile, **kwargs): Parameters ---------- - datafile : str - name of tif file containing raw data + beamline_full_datafile_name : str + full name of tif file containing beamline preprocessed data kwargs : keyword arguments - save_dir : str - Directory where results of reconstruction are saved as npy files. If not present, the reconstruction outcome will be save in the same directory where datafile is. + data_dir : str + directory where prepared data will be saved, default /phasing_data alien_alg : str Name of method used to remove aliens. Possible options are: ‘block_aliens’, ‘alien_file’, and ‘AutoAlien1’. The ‘block_aliens’ algorithm will zero out defined blocks, ‘alien_file’ method will use given file as a mask, and ‘AutoAlien1’ will use auto mechanism to remove aliens. Each of these algorithms require different parameters aliens : list @@ -68,36 +68,51 @@ def prep(datafile, **kwargs): Optional, enter center shift list the array maximum is centered before binning, and moved according to center_shift, [0,0,0] has no effect binning : list Optional, a list that defines binning values in respective dimensions, [1,1,1] has no effect. - + do_auto_binning : boolean + Optional, mandatory if auto_data is True. is True the auto binning wil be done, and not otherwise. + debug : boolean + It's a command line argument passed as parameter. If True, ignores verifier error. """ + debug = kwargs.pop("debug", False) er_msg = ver.verify('config_data', kwargs) if len(er_msg) > 0: # the error message is printed in verifier - return + if not debug: + return - datafile = datafile.replace(os.sep, '/') + beamline_full_datafile_name = beamline_full_datafile_name.replace(os.sep, '/') # The data has been transposed when saved in tif format for the ImageJ to show the right orientation - data = ut.read_tif(datafile) + data = ut.read_tif(beamline_full_datafile_name) + if debug: + print(f"Loaded array (max={int(data.max())}) as {beamline_full_datafile_name}") + prep_data_dir, beamline_datafile_name = os.path.split(beamline_full_datafile_name) if 'data_dir' in kwargs: data_dir = kwargs['data_dir'].replace(os.sep, '/') else: - data_dir, filename = os.path.split(datafile) + # assuming the directory structure and naming follows cohere-ui mechanism + data_dir = prep_data_dir.replace(os.sep, '/').replace('preprocessed_data', 'phasing_data') if 'alien_alg' in kwargs: data = at.remove_aliens(data, kwargs, data_dir) - if 'intensity_threshold' in kwargs: + if auto: + # the formula for auto threshold was found empirically, may be + # modified in the future if more tests are done + auto_threshold_value = 0.141 * data[np.nonzero(data)].mean() - 3.062 + intensity_threshold = auto_threshold_value + print(f'auto intensity threshold: {intensity_threshold}') + elif 'intensity_threshold' in kwargs: intensity_threshold = kwargs['intensity_threshold'] else: print('define amplitude threshold. Exiting') return # zero out the noise - prep_data = np.where(data <= intensity_threshold, 0.0, data) + data = np.where(data <= intensity_threshold, 0.0, data) # square root data - prep_data = np.sqrt(prep_data) + data = np.sqrt(data) if 'adjust_dimensions' in kwargs: crops_pads = kwargs['adjust_dimensions'] @@ -114,31 +129,55 @@ def prep(datafile, **kwargs): pair = crops_pads[2 * i:2 * i + 2] pairs.append(pair) - prep_data = ut.adjust_dimensions(prep_data, pairs) - if prep_data is None: + data = ut.adjust_dimensions(data, pairs) + if data is None: print('check "adjust_dimensions" configuration') return if 'center_shift' in kwargs: center_shift = kwargs['center_shift'] - prep_data = ut.get_centered(prep_data, center_shift) + data, shift = ut.center_max(data, center_shift) else: - prep_data = ut.get_centered(prep_data, [0, 0, 0]) - - if 'binning' in kwargs: + data, shift = ut.center_max(data, [0, 0, 0]) + + try: + # assuming the mask file is in directory of preprocessed data + mask = ut.read_tif(beamline_full_datafile_name.replace(beamline_datafile_name, 'mask.tif')) + mask = np.roll(mask, shift, tuple(range(mask.ndim))) + ut.save_tif(mask, ut.join(data_dir, 'mask.tif')) + except FileNotFoundError: + pass + + if auto and kwargs['do_auto_binning']: + # prepare data to make the oversampling ratio ~3 + wos = 3.0 + orig_os = ut.get_oversample_ratio(data) + # match oversampling to wos + wanted_os = [wos, wos, wos] + change = [np.round(f / t).astype('int32') for f, t in zip(orig_os, wanted_os)] + bins = [int(max([1, c])) for c in change] + print(f'auto binning size: {bins}') + data = ut.binning(data, bins) + elif 'binning' in kwargs: binsizes = kwargs['binning'] try: bins = [] for binsize in binsizes: bins.append(binsize) - filler = len(prep_data.shape) - len(bins) + filler = len(data.shape) - len(bins) for _ in range(filler): bins.append(1) - prep_data = ut.binning(prep_data, bins) + data = ut.binning(data, bins) except: print('check "binning" configuration') # save data - data_file = data_dir + '/data.tif' - ut.save_tif(prep_data, data_file) - print('data ready for reconstruction, data dims:', prep_data.shape) + data_file = ut.join(data_dir, 'data.tif') + ut.save_tif(data, data_file) + print(f'data ready for reconstruction, data dims: {data.shape}') + + # if auto save new config + if auto: + kwargs['binning'] = bins + kwargs['intensity_threshold'] = intensity_threshold + return kwargs \ No newline at end of file diff --git a/cohere_core/lib/aflib.py b/cohere_core/lib/aflib.py index d76c962..e6cf845 100644 --- a/cohere_core/lib/aflib.py +++ b/cohere_core/lib/aflib.py @@ -333,3 +333,52 @@ def gaussian_filter(arr, sigma, **kwargs): correction = arr_sum / af.sum(convag) convag *= correction return convag + + def diff(arr, axis=None, prepend=0): + raise NotImplementedError + + def gradient(arr, dx=1): + raise NotImplementedError + + def argmin(arr, axis=None): + raise NotImplementedError + + def take_along_axis(a, indices, axis): + raise NotImplementedError + + def moveaxis(arr, source, dest): + raise NotImplementedError + + def lstsq(A, B): + raise NotImplementedError + + def zeros(shape): + raise NotImplementedError + + def indices(dims): + raise NotImplementedError + + def concatenate(tup, axis=0): + raise NotImplementedError + + def amin(arr): + raise NotImplementedError + + def affine_transform(arr, matrix, order=3, offset=0): + raise NotImplementedError + + def pad(arr, padding): + raise NotImplementedError + + def histogram2d(meas, rec, n_bins=100, log=False): + raise NotImplementedError + + def calc_nmi(hgram): + raise NotImplementedError + + def calc_ehd(hgram): + raise NotImplementedError + + def clean_default_mem(): + pass + diff --git a/cohere_core/lib/cohlib.py b/cohere_core/lib/cohlib.py index 0e2ae8d..f473648 100644 --- a/cohere_core/lib/cohlib.py +++ b/cohere_core/lib/cohlib.py @@ -21,24 +21,30 @@ def __subclasshook__(cls, subclass): callable(subclass.from_numpy) and hasattr(subclass, 'dtype') and callable(subclass.dtype) and + hasattr(subclass, 'astype') and + callable(subclass.astype) and hasattr(subclass, 'size') and callable(subclass.size) and hasattr(subclass, 'hasnan') and callable(subclass.hasnan) and hasattr(subclass, 'copy') and callable(subclass.copy) and + hasattr(subclass, 'roll') and + callable(subclass.shift) and + hasattr(subclass, 'roll') and + callable(subclass.shift) and hasattr(subclass, 'fftshift') and callable(subclass.fftshift) and hasattr(subclass, 'ifftshift') and callable(subclass.ifftshift) and - hasattr(subclass, 'shift') and - callable(subclass.shift) and hasattr(subclass, 'fft') and callable(subclass.fft) and hasattr(subclass, 'ifft') and callable(subclass.ifft) and hasattr(subclass, 'fftconvolve') and callable(subclass.fftconvolve) and + hasattr(subclass, 'correlate') and + callable(subclass.correlate) and hasattr(subclass, 'where') and callable(subclass.where) and hasattr(subclass, 'dims') and @@ -98,13 +104,33 @@ def __subclasshook__(cls, subclass): hasattr(subclass, 'conj') and callable(subclass.conj) and hasattr(subclass, 'array_equal') and - callable(subclass.array_equal) or + callable(subclass.array_equal) and hasattr(subclass, 'cos') and - callable(subclass.cos) or + callable(subclass.cos) and hasattr(subclass, 'linspace') and - callable(subclass.linspace) or + callable(subclass.linspace) and hasattr(subclass, 'clip') and - callable(subclass.clip) or + callable(subclass.clip) and + hasattr(subclass, 'gradient') and + callable(subclass.gradient) and + hasattr(subclass, 'argmin') and + callable(subclass.argmin) and + hasattr(subclass, 'take_along_axis') and + callable(subclass.take_along_axis) and + hasattr(subclass, 'moveaxis') and + callable(subclass.moveaxis) and + hasattr(subclass, 'lstsq') and + callable(subclass.lstsq) and + hasattr(subclass, 'zeros') and + callable(subclass.zeros) and + hasattr(subclass, 'indices') and + callable(subclass.indices) and + hasattr(subclass, 'diff') and + callable(subclass.diff) and + hasattr(subclass, 'concatenate') and + callable(subclass.concatenate) and + callable(subclass, 'clean_default_mem') and + callable(subclass.clean_default_mem) or NotImplemented) @abc.abstractmethod @@ -143,6 +169,10 @@ def from_numpy(arr): def dtype(arr): pass + @abc.abstractmethod + def astype(arr): + pass + @abc.abstractmethod def size(arr): pass @@ -155,6 +185,14 @@ def hasnan(arr): def copy(arr): pass + @abc.abstractmethod + def roll(arr, sft, axis): + pass + + @abc.abstractmethod + def shift(arr, sft): + pass + @abc.abstractmethod def fftshift(arr): pass @@ -164,19 +202,19 @@ def ifftshift(arr): pass @abc.abstractmethod - def shift(arr, sft): + def fft(arr, norm='forward'): pass @abc.abstractmethod - def fft(arr): + def ifft(arr, norm='forward'): pass @abc.abstractmethod - def ifft(arr): + def fftconvolve(arr1, kernel): pass @abc.abstractmethod - def fftconvolve(arr1, arr2): + def correlate(arr1, arr2, mode='same', method='fft'): pass @abc.abstractmethod @@ -308,6 +346,70 @@ def cos(arr): def linspace(start, stop, num): pass + @abc.abstractmethod + def diff(arr, axis=None, prepend=0): + pass + @abc.abstractmethod def clip(arr, min, max=None): pass + + @abc.abstractmethod + def gradient(arr, dx=1): + pass + + @abc.abstractmethod + def argmin(arr, axis=None): + pass + + @abc.abstractmethod + def take_along_axis(a, indices, axis): + pass + + @abc.abstractmethod + def moveaxis(arr, source, dest): + pass + + @abc.abstractmethod + def lstsq(A, B): + pass + + @abc.abstractmethod + def zeros(shape): + pass + + @abc.abstractmethod + def indices(dims): + pass + + @abc.abstractmethod + def concatenate(tup, axis=0): + pass + + @abc.abstractmethod + def amin(arr): + pass + + @abc.abstractmethod + def affine_transform(arr, matrix, order=3, offset=0): + pass + + @abc.abstractmethod + def pad(arr, padding): + pass + + @abc.abstractmethod + def histogram2d(meas, rec, n_bins=100, log=False): + pass + + @abc.abstractmethod + def calc_nmi(hgram): + pass + + @abc.abstractmethod + def calc_ehd(hgram): + pass + + @abc.abstractmethod + def clean_default_mem(): + pass diff --git a/cohere_core/lib/cplib.py b/cohere_core/lib/cplib.py index df5b5ca..2a3c093 100644 --- a/cohere_core/lib/cplib.py +++ b/cohere_core/lib/cplib.py @@ -2,7 +2,9 @@ import cupy as cp import numpy as np import cupyx.scipy.ndimage -import cupyx.scipy.signal +import cupyx.scipy.stats as stats +import cupyx.scipy.ndimage as sc + class cplib(cohlib): def array(obj): @@ -28,11 +30,14 @@ def save(filename, arr): cp.save(filename, arr) def load(filename): - return cp.load(filename) + return cp.load(filename, allow_pickle=True) def dtype(arr): return arr.dtype + def astype(arr, dtype): + return arr.astype(dtype=dtype) + def size(arr): return arr.size @@ -50,25 +55,36 @@ def random(shape, **kwargs): rs = cp.random.RandomState(seed=seed) return cp.random.random(shape, dtype=cp.float32) + 1j * cp.random.random(shape, dtype=cp.float32) + def roll(arr, sft, axis=None): + if axis is None: + axis = list(range(len(sft))) + if type(sft) != list: + sft = [sft] + sft = [int(s) for s in sft] + return cp.roll(arr, sft, axis=axis) + + def shift(arr, sft): + return sc.fourier_shift(arr, sft) + def fftshift(arr): return cp.fft.fftshift(arr) def ifftshift(arr): return cp.fft.ifftshift(arr) - def shift(arr, sft, axis=None): - sft = [int(s) for s in sft] - return cp.roll(arr, sft, axis=axis) + def fft(arr, norm='forward'): + return cp.fft.fftn(arr, norm=norm) - def fft(arr): - return cp.fft.fftn(arr, norm='forward') + def ifft(arr, norm='forward'): + return cp.fft.ifftn(arr, norm=norm) - def ifft(arr): - return cp.fft.ifftn(arr, norm='forward') + def fftconvolve(arr1, kernel): + from cupyx.scipy import signal + return signal.fftconvolve(arr1, kernel, mode='same') - def fftconvolve(arr1, arr2): - return cupyx.scipy.ndimage.convolve(arr1, arr2) - # return cupyx.scipy.signal.aoconvolve(arr1, arr2, mode='same') + def correlate(arr1, arr2, mode='same', method='fft'): + from cupyx.scipy.signal import correlate + return correlate(arr1, arr2, mode, method) def where(cond, x, y): return cp.where(cond, x, y) @@ -88,8 +104,8 @@ def square(arr): def sum(arr, axis=None): sm = cp.sum(arr, axis) - if axis is None: - return sm.tolist() + # if axis is None: + # return sm.tolist() return sm def real(arr): @@ -148,14 +164,15 @@ def gaussian(shape, sigma, **kwargs): inarr = cp.zeros((n_el)) inarr[int(n_el / 2)] = 1.0 inarr = cp.reshape(inarr, shape) - gaussian = cupyx.scipy.ndimage.gaussian_filter(inarr, sigma) + gaussian = sc.gaussian_filter(inarr, sigma) return gaussian / cp.sum(gaussian) def gaussian_filter(arr, sigma, **kwargs): - return cupyx.scipy.ndimage.gaussian_filter(arr, sigma) + return sc.gaussian_filter(arr, sigma) def center_of_mass(inarr): - return cupyx.scipy.ndimage.center_of_mass(cp.absolute(inarr)) + t = sc.center_of_mass(cp.absolute(inarr)) + return t def meshgrid(*xi): return cp.meshgrid(*xi) @@ -166,15 +183,6 @@ def exp(arr): def conj(arr): return cp.conj(arr) - def cos(arr): - return cp.cos(arr) - - def linspace(start, stop, steps): - return cp.linspace(start, stop, steps) - - def diff(arr, axis=None, prepend=0): - return cp.diff(arr, axis=axis, prepend=prepend) - def array_equal(arr1, arr2): return cp.array_equal(arr1, arr2) @@ -186,3 +194,62 @@ def linspace(start, stop, num): def clip(arr, min, max=None): return cp.clip(arr, min, max) + + def diff(arr, axis=None, prepend=0): + return cp.diff(arr, axis=axis, prepend=prepend) + + def gradient(arr, dx=1): + return cp.gradient(arr, dx) + + def argmin(arr, axis=None): + return cp.argmin(arr, axis) + + def take_along_axis(a, indices, axis): + return cp.take_along_axis(a, indices, axis) + + def moveaxis(arr, source, dest): + return cp.moveaxis(arr, source, dest) + + def lstsq(A, B): + return cp.linalg.lstsq(A, B, rcond=None) + + def zeros(shape): + return cp.zeros(shape) + + def indices(dims): + return cp.indices(dims) + + def concatenate(tup, axis=0): + return cp.concatenate(tup, axis) + + def amin(arr): + return cp.amin(arr) + + def affine_transform(arr, matrix, order=3, offset=0): + return cupyx.scipy.ndimage.affine_transform(arr, matrix, order=order, offset=offset, prefilter=True) + + def pad(arr, padding): + return cp.pad(arr, padding) + + def histogram2d(meas, rec, n_bins=100, log=False): + norm = cp.max(meas) / cp.max(rec) + if log: + bins = cp.logspace(cp.log10(1), cp.log10(cp.max(meas)), n_bins+1) + else: + bins = n_bins + return cp.histogram2d(cp.ravel(meas), cp.ravel(norm*rec), bins)[0] + + def calc_nmi(hgram): + h0 = stats.entropy(np.sum(hgram, axis=0)) + h1 = stats.entropy(np.sum(hgram, axis=1)) + h01 = stats.entropy(np.reshape(hgram, -1)) + return (h0 + h1) / h01 + + def calc_ehd(hgram): + n = hgram.shape[0] * 1j + x, y = cp.mgrid[0:1:n, 0:1:n] + return cp.sum(hgram * cp.abs(x - y)) / cp.sum(hgram) + + def clean_default_mem(): + cp._default_memory_pool.free_all_blocks() + cp._default_pinned_memory_pool.free_all_blocks() diff --git a/cohere_core/lib/nplib.py b/cohere_core/lib/nplib.py index 4c5c6f9..cfc5bb1 100644 --- a/cohere_core/lib/nplib.py +++ b/cohere_core/lib/nplib.py @@ -1,6 +1,7 @@ from cohere_core.lib.cohlib import cohlib import numpy as np -from scipy.ndimage import convolve, gaussian_filter, center_of_mass +from scipy.ndimage import convolve, gaussian_filter, center_of_mass, shift +from scipy.signal import correlate class nplib(cohlib): @@ -31,6 +32,9 @@ def save(filename, arr): def dtype(arr): return arr.dtype + def astype(arr, dtype): + return arr.astype(dtype=dtype) + def size(arr): return arr.size @@ -51,24 +55,34 @@ def random(shape, **kwargs): return r #return rng.random(*shape).astype(float) + def roll(arr, sft, axis=None): + if type(sft) != list: + sft = [sft] + if axis is None: + axis = list(range(len(sft))) + sft = [int(s) for s in sft] + return np.roll(arr, sft, axis=axis) + + def shift(arr, sft): + return shift(arr, sft) + def fftshift(arr): return np.fft.fftshift(arr) def ifftshift(arr): return np.fft.ifftshift(arr) - def shift(arr, sft): - sft = [int(s) for s in sft] - return np.roll(arr, sft) - def fft(arr): return np.fft.fftn(arr, norm='forward') def ifft(arr): return np.fft.ifftn(arr, norm='forward') - def fftconvolve(arr1, arr2): - return convolve(arr1, arr2) + def fftconvolve(arr1, kernel): + return convolve(arr1, kernel) + + def correlate(arr1, arr2, mode='same', method='auto'): + return correlate(arr1, arr2, mode, method) def where(cond, x, y): return np.where(cond, x, y) @@ -186,6 +200,54 @@ def linspace(start, stop, num): def clip(arr, min, max=None): return np.clip(arr, min, max) + def diff(arr, axis=None, prepend=0): + return np.diff(arr, axis=axis, prepend=prepend) + + def gradient(arr, dx=1): + return np.gradient(arr, dx) + + def argmin(arr, axis=None): + return np.argmin(arr, axis) + + def take_along_axis(a, indices, axis): + return np.take_along_axis(a, indices, axis) + + def moveaxis(arr, source, dest): + return np.moveaxis(arr, source, dest) + + def lstsq(A, B): + return np.linalg.lstsq(A, B, rcond=None) + + def zeros(shape): + return np.zeros(shape) + + def indices(dims): + return np.indices(dims) + + def concatenate(tup, axis=0): + return np.concatenate(tup, axis) + + def amin(arr): + return np.amin(arr) + + def affine_transform(arr, matrix, order=3, offset=0): + raise NotImplementedError + + def pad(arr, padding): + raise NotImplementedError + + def histogram2d(meas, rec, n_bins=100, log=False): + raise NotImplementedError + + def calc_nmi(hgram): + raise NotImplementedError + + def calc_ehd(hgram): + raise NotImplementedError + + def clean_default_mem(): + pass + # a11 = np.array([0.1, 0.2, 0.3, 1.0, 1.2, 1.3]) # a12 = np.array([10.1, 10.2, 10.3, 11.0]) # print(np.convolve(a11,a12, mode='same')) diff --git a/cohere_core/lib/torchlib.py b/cohere_core/lib/torchlib.py index 2b88209..40c337c 100644 --- a/cohere_core/lib/torchlib.py +++ b/cohere_core/lib/torchlib.py @@ -4,8 +4,11 @@ from torch.nn.functional import conv1d import sys +import math + class torchlib(cohlib): + device = "cpu" # Interface def array(obj): return torch.Tensor(obj, device=torchlib.device) @@ -65,6 +68,13 @@ def from_numpy(arr): def dtype(arr): return arr.dtype + def astype(arr, dtype): + # this is kind of nasty, it does not understand 'int', so need to convert for the cases + # that will be used + if dtype == 'int32': + dtype = torch.int32 + return arr.type(dtype=dtype) + def size(arr): return torch.numel(arr) @@ -75,7 +85,26 @@ def copy(arr): return arr.clone() def random(shape, **kwargs): - return torch.rand(shape, device=torchlib.device) + arr = torch.rand(shape, device=torchlib.device) + print(arr.dtype) + # return torch.rand(shape, device=torchlib.device) + return arr + + def roll(arr, sft): + sft = [int(s) for s in sft] + dims = tuple([i for i in range(len(sft))]) + try: + return torch.roll(arr, sft, dims) + except Exception as e: + print('not supported error: ' + repr(e)) + + def shift(arr, sft): + sft = [int(s) for s in sft] + dims = tuple([i for i in range(len(sft))]) + try: + return torch.roll(arr, sft, dims) + except Exception as e: + print('not supported error: ' + repr(e)) def fftshift(arr): # if roll is implemented before fftshift @@ -97,14 +126,6 @@ def ifftshift(arr): except Exception as e: print('not supported error: ' + repr(e)) - def shift(arr, sft): - sft = [int(s) for s in sft] - dims = tuple([i for i in range(len(sft))]) - try: - return torch.roll(arr, sft, dims) - except Exception as e: - print('not supported error: ' + repr(e)) - def fft(arr): try: return torch.fft.fftn(arr, norm='forward') @@ -117,19 +138,21 @@ def ifft(arr): except Exception as e: print('not supported error: ' + repr(e)) - def fftconvolve(arr1, arr2): - shape1 = arr1.size() - shape2 = arr2.size() - pad1 = tuple([d//2 for d in shape2 for _ in (0, 1)]) - pad2 = tuple([d//2 for d in shape1 for _ in (0, 1)]) - parr1 = torch.nn.functional.pad(arr1, pad1) - parr2 = torch.nn.functional.pad(arr2, pad2) - conv = torch.fft.ifftn(torch.fft.fftn(parr1) * torch.fft.fftn(parr2)) - - for i in range(len(shape2)): - splitted = torch.split(conv, [shape2[i] //2, shape1[i], shape2[i] //2], dim=i) - conv = splitted[1] - return conv + def fftconvolve(arr1, kernel): + print('not supported yet in torch, use different library') + raise + # # kernel shape can be smaller than arr1 shape in each dim + # sh1 = list(arr1.size()) + # sh2 = list(kernel.size()) + # if sh1 != sh2: + # # the pad is added from last dim to first + # sh1.reverse() + # sh2.reverse() + # pad = [((sh1[i]-sh2[i])//2, sh1[i] - sh2[i] - (sh1[i]-sh2[i])//2) for i in range(len(sh1))] + # pad = tuple(sum(pad, ())) + # kernel = torch.nn.functional.pad(kernel, pad) + # conv = torch.fft.ifftn(torch.fft.fftn(arr1) * torch.fft.fftn(kernel)) + # return conv def where(cond, x, y): return torch.where(cond, x, y) @@ -209,15 +232,19 @@ def gaussian1(dim, sig): return gauss / gauss.sum() if isinstance(sigma, int) or isinstance(sigma, float): - sigma = [sigma] * len(dims) + sigmas = [sigma] * len(dims) + else: # assuming a list of values + sigmas = sigma + + sigmas = [dims[i] / (2.0 * math.pi * sigmas[i]) for i in range(len(dims))] last_axis = len(dims) - 1 last_dim = dims[-1] - gaussian = gaussian1(dims[last_axis], sigma[last_axis]).repeat(*dims[:-1], 1) + gaussian = gaussian1(dims[last_axis], sigmas[last_axis]).repeat(*dims[:-1], 1) for i in range(int(len(dims)) - 1): tdims = dims[:-1] tdims[i] = last_dim - temp = gaussian1(dims[i], sigma[i]) + temp = gaussian1(dims[i], sigmas[i]) temp = temp.repeat(*tdims, 1) temp = torch.swapaxes(temp, i, last_axis) gaussian *= temp @@ -225,14 +252,16 @@ def gaussian1(dim, sig): def gaussian_filter(arr, sigma, **kwargs): # will not work on M1 until the fft fuctions are supported + normalizer = torch.sum(arr) dims = list(arr.shape) - dist = torchlib.gaussian(dims, 1/sigma) + dist = torchlib.gaussian(dims, 1 / sigma) arr_f = torch.fft.ifftshift(torch.fft.fftn(torch.fft.ifftshift(arr))) filter = arr_f * dist filter = torch.fft.ifftshift(torch.fft.ifftn(torch.fft.ifftshift(filter))) filter = torch.real(filter) filter = torch.where(filter >= 0, filter, 0.0) - return filter + corr = normalizer/torch.sum(filter) + return filter * corr def center_of_mass(arr): normalizer = torch.sum(arr) @@ -261,10 +290,50 @@ def linspace(start, stop, num): def clip(arr, min, max=None): return torch.clip(arr, min, max) -# a1 = torch.Tensor([0.1, 0.2, 0.3, 1.0, 1.2, 1.3]) -# a2 = torch.Tensor([10.1, 10.2, 10.3, 11.0]) -# conv = torchlib.fftconvolve(a1,a2) -# print('torch conv', conv) -# print(conv.real) -# print(torch.abs(conv)) -# print(torch.nn.functional.conv1d(a1,a2)) + def diff(arr, axis=None, prepend=0): + raise NotImplementedError + + def gradient(arr, dx=1): + raise NotImplementedError + + def argmin(arr, axis=None): + raise NotImplementedError + + def take_along_axis(a, indices, axis): + raise NotImplementedError + + def moveaxis(arr, source, dest): + raise NotImplementedError + + def lstsq(A, B): + raise NotImplementedError + + def zeros(shape): + raise NotImplementedError + + def indices(dims): + raise NotImplementedError + + def concatenate(tup, axis=0): + raise NotImplementedError + + def amin(arr): + raise NotImplementedError + + def affine_transform(arr, matrix, order=3, offset=0): + raise NotImplementedError + + def pad(arr, padding): + raise NotImplementedError + + def histogram2d(meas, rec, n_bins=100, log=False): + raise NotImplementedError + + def calc_nmi(hgram): + raise NotImplementedError + + def calc_ehd(hgram): + raise NotImplementedError + + def clean_default_mem(): + pass diff --git a/cohere_core/utilities/config_errors_dict.py b/cohere_core/utilities/config_errors_dict.py index 2a32c24..a5ae7e5 100644 --- a/cohere_core/utilities/config_errors_dict.py +++ b/cohere_core/utilities/config_errors_dict.py @@ -89,7 +89,7 @@ 'AI_trained_model parameter is mandatory when init_guess is "AI algorithm"'], 'Reconstruction':['reconstructions parameter should be int', 'reconstructions parameter parsing error'], - 'Device':['device not a list of ints', + 'Device':['device should be a list of ints, or string "all", or dictionary for cluster configuration', 'device parameter parsing error'], 'Algorithmsequence':['algorithm_sequence should be string', 'algorithm_sequence configuration error, only numerical, digital, and the following characters: *, +, (, ), space, are allowed ', @@ -125,11 +125,11 @@ 'Shrinkwrapgausssigma':['sw_gauss_sigma should be float', 'sw_gauss_sigma should be a list of floats if multiple shrink wraps'], 'Phasesupporttrigger':['Trigger should be a list of int', - 'Each sub-trigger should be a list of int if multiple phase modulus'], - 'Phmphasemin':['phm_phase_min should be float', - 'phm_phase_min should be a list of floats if multiple phase modulus'], - 'Phmphasemax':['phm_phase_max should be float', - 'phm_phase_max should be a list of floats if multiple phase modulus'], + 'Each sub-trigger should be a list of int if multiple phase constrain'], + 'Phcphasemin':['phc_phase_min should be float', + 'phc_phase_min should be a list of floats if multiple phase constrain'], + 'Phcphasemax':['phc_phase_max should be float', + 'phc_phase_max should be a list of floats if multiple phase constrain'], 'Pcinterval':['pc_interval should be int'], 'Pctype':['pc_type parameter should be string', 'pc_type parameter can be configured "LUCY"', diff --git a/cohere_core/utilities/config_verifier.py b/cohere_core/utilities/config_verifier.py index 35c895e..27a0e2d 100644 --- a/cohere_core/utilities/config_verifier.py +++ b/cohere_core/utilities/config_verifier.py @@ -13,9 +13,9 @@ import os from cohere_core.utilities.config_errors_dict import * -from cohere_core.controller.op_flow import algs +import cohere_core.utilities.utils as ut -__author__ = "Barbara Frosik, Dave Cyl" +__author__ = "Dave Cyl" __copyright__ = "Copyright (c) 2016, UChicago Argonne, LLC." __docformat__ = 'restructuredtext en' __all__ = ['verify'] @@ -39,11 +39,11 @@ def ver_list_int(param_name, param_value): True if all elements are int, False otherwise """ if not issubclass(type(param_value), list): - print (param_name + ' is not a list') + print (f'{param_name} is not a list') return False for e in param_value: if type(e) != int: - print (param_name + ' should be list of integer values') + print (f'{param_name} should be list of integer values') return False return True @@ -66,11 +66,11 @@ def ver_list_float(param_name, param_value): True if all elements are float, False otherwise """ if not issubclass(type(param_value), list): - print (param_name + ' is not a list') + print (f'{param_name} is not a list') return False for e in param_value: if type(e) != float: - print (param_name + ' should be list of float values') + print (f'{param_name} should be list of float values') return False return True @@ -234,14 +234,16 @@ def parse_entry(ent): return '', sum([e[0] for e in seq]) - def verify_trigger(trigger, no_iter): + def verify_trigger(trigger, no_iter, trigger_name): + if type(trigger) != list: + return(f'{trigger_name} trigger type should be list') if len(trigger) == 0: - return ('empty trigger ' + str(trigger)) + return (f'empty {trigger_name} trigger {str(trigger)}') elif trigger[0] >= no_iter: - return ('trigger start ' + str(trigger[0]) + ' exceeds number of iterations ' + str(no_iter)) + return (f'{trigger_name} trigger start {str(trigger[0])} exceeds number of iterations {str(no_iter)}') if len(trigger) == 3: if trigger[2] >= no_iter: - return ('trigger end ' + str(trigger[2]) + ' exceeds number of iterations ' + str(no_iter)) + return (f'{trigger_name} trigger end {str(trigger[2])} exceeds number of iterations {str(no_iter)}') return '' @@ -261,7 +263,7 @@ def verify_trigger(trigger, no_iter): error_message = get_config_error_message(fname, config_map_file, config_parameter, config_error) print(error_message) return (error_message) - if not os.path.isfile(data_dir + '/data.tif') and not os.path.isfile(data_dir + '/data.npy'): + if not os.path.isfile(ut.join(data_dir, 'data.tif')) and not os.path.isfile(ut.join(data_dir, 'data.npy')): config_error = 2 error_message = get_config_error_message(fname, config_map_file, config_parameter, config_error) print(error_message) @@ -326,12 +328,27 @@ def verify_trigger(trigger, no_iter): config_parameter = 'Device' if 'device' in config_map: + def ver_dev(d): + if d == 'all' or ver_list_int('device', device): + return '' + else: + config_error = 0 + error_message = get_config_error_message(fname, config_map_file, config_parameter, config_error) + print(error_message) + return (error_message) + device = config_map['device'] - if not ver_list_int('device', device): - config_error = 0 - error_message = get_config_error_message(fname, config_map_file, config_parameter, config_error) - print(error_message) - return (error_message) + if issubclass(type(device), dict): + for d in device.values(): + error_message = ver_dev(d) + if len(error_message) > 0: + print(error_message) + return (error_message) + else: + error_message = ver_dev(device) + if len(error_message) > 0: + print(error_message) + return (error_message) config_parameter = 'Algorithmsequence' if 'algorithm_sequence' in config_map: @@ -506,8 +523,9 @@ def verify_trigger(trigger, no_iter): config_parameter = 'Twintrigger' if 'twin_trigger' in config_map: twin_trigger = config_map['twin_trigger'] - m = verify_trigger(twin_trigger, iter_no) + m = verify_trigger(twin_trigger, iter_no, 'twin') if len(m) > 0: + print(m) return(m) if not ver_list_int('twin_trigger', twin_trigger): config_error = 0 @@ -527,8 +545,9 @@ def verify_trigger(trigger, no_iter): config_parameter = 'Shrinkwraptrigger' if 'shrink_wrap_trigger' in config_map: if '.SW' not in config_map['algorithm_sequence']: - m = verify_trigger(config_map['shrink_wrap_trigger'], iter_no) + m = verify_trigger(config_map['shrink_wrap_trigger'], iter_no, 'shrink wrap') if len(m) > 0: + print(m) return(m) if not ver_list_int('shrink_wrap_trigger', config_map['shrink_wrap_trigger']): config_error = 0 @@ -591,53 +610,53 @@ def verify_trigger(trigger, no_iter): return (error_message) config_parameter = 'Phasesupporttrigger' - if 'phm_trigger' in config_map: - if '.PHM' not in config_map['algorithm_sequence']: - m = verify_trigger(config_map['phm_trigger'], iter_no) - print(m) + if 'phc_trigger' in config_map: + if '.PHC' not in config_map['algorithm_sequence']: + m = verify_trigger(config_map['phc_trigger'], iter_no, 'phase constrain') if len(m) > 0: + print(m) return(m) - if not ver_list_int('phm_trigger', config_map['phm_trigger']): + if not ver_list_int('phc_trigger', config_map['phc_trigger']): config_error = 0 error_message = get_config_error_message(fname, config_map_file, config_parameter, config_error) print(error_message) return (error_message) - config_parameter = 'Phmphasemin' - if 'phm_phase_min' in config_map: - phm_phase_min = config_map['phm_phase_min'] - if type(phm_phase_min) != float: + config_parameter = 'Phcphasemin' + if 'phc_phase_min' in config_map: + phc_phase_min = config_map['phc_phase_min'] + if type(phc_phase_min) != float: config_error = 0 error_message = get_config_error_message(fname, config_map_file, config_parameter, config_error) print(error_message) return (error_message) - config_parameter = 'Phmphasemax' - if 'phm_phase_max' in config_map: - phm_phase_max = config_map['phm_phase_max'] - if type(phm_phase_max) != float: + config_parameter = 'Phcphasemax' + if 'phc_phase_max' in config_map: + phc_phase_max = config_map['phc_phase_max'] + if type(phc_phase_max) != float: config_error = 0 error_message = get_config_error_message(fname, config_map_file, config_parameter, config_error) print(error_message) return (error_message) else: - for t in config_map['phm_trigger']: - if not ver_list_int('phm_trigger', t): + for t in config_map['phc_trigger']: + if not ver_list_int('phc_trigger', t): config_error = 1 error_message = get_config_error_message(fname, config_map_file, config_parameter, config_error) print(error_message) return (error_message) - config_parameter = 'Phmphasemin' - if 'phm_phase_min' in config_map: - if not ver_list_float('phm_phase_min', config_map['phm_phase_min']): + config_parameter = 'Phcphasemin' + if 'phc_phase_min' in config_map: + if not ver_list_float('phc_phase_min', config_map['phc_phase_min']): config_error = 1 error_message = get_config_error_message(fname, config_map_file, config_parameter, config_error) print(error_message) return (error_message) - config_parameter = 'Phmphasemax' - if 'phm_phase_max' in config_map: - if not ver_list_float('phm_phase_max', config_map['phm_phase_max']): + config_parameter = 'Phcphasemax' + if 'phc_phase_max' in config_map: + if not ver_list_float('phc_phase_max', config_map['phc_phase_max']): config_error = 1 error_message = get_config_error_message(fname, config_map_file, config_parameter, config_error) print(error_message) @@ -703,8 +722,9 @@ def verify_trigger(trigger, no_iter): config_parameter = 'Resolutiontrigger' if 'lowpass_filter_trigger' in config_map: if '.LPF' not in config_map['algorithm_sequence']: - m = verify_trigger(config_map['lowpass_filter_trigger'], iter_no) + m = verify_trigger(config_map['lowpass_filter_trigger'], iter_no, 'lowpass filter') if len(m) > 0: + print(m) return(m) if not ver_list_int('lowpass_filter_trigger', config_map['lowpass_filter_trigger']): config_error = 0 @@ -767,8 +787,9 @@ def verify_trigger(trigger, no_iter): config_parameter = 'Averagetrigger' if 'average_trigger' in config_map: - m = verify_trigger(config_map['average_trigger'], iter_no) + m = verify_trigger(config_map['average_trigger'], iter_no, 'average') if len(m) > 0: + print(m) return(m) if not ver_list_int('average_trigger', config_map['average_trigger']): config_error = 0 @@ -778,7 +799,7 @@ def verify_trigger(trigger, no_iter): config_parameter = 'Progresstrigger' if 'progress_trigger' in config_map: - m = verify_trigger(config_map['progress_trigger'], iter_no) + m = verify_trigger(config_map['progress_trigger'], iter_no, 'progress') if len(m) > 0: return(m) if not ver_list_int('progress_trigger', config_map['progress_trigger']): @@ -848,12 +869,7 @@ def ver_config_data(config_map): config_error = 0 error_message = get_config_error_message(fname, config_map_file, config_parameter, config_error) print(error_message) - return (error_message) - else: - config_error = 1 - error_message = get_config_error_message(fname, config_map_file, config_parameter, config_error) - print(error_message) - return (error_message) + return '' config_parameter = 'Alienalg' if 'alien_alg' in config_map: @@ -1329,5 +1345,7 @@ def verify(file_name, conf_map): return ver_config_disp(conf_map) elif file_name == 'config_instr': return ver_config_instr(conf_map) + elif file_name == 'config_mp': + return '' else: return ('verifier has no fumction to check config file named', file_name) diff --git a/cohere_core/utilities/dvc_utils.py b/cohere_core/utilities/dvc_utils.py index 5917b29..a2238c7 100644 --- a/cohere_core/utilities/dvc_utils.py +++ b/cohere_core/utilities/dvc_utils.py @@ -1,4 +1,33 @@ import math +import cmath +import cohere_core.utilities.utils as ut + +_docformat__ = 'restructuredtext en' +__all__ = ['set_lib', + 'arr_property', + 'pad_around', + 'center_sync', + 'gauss_conv_fft', + 'shrink_wrap', + 'shift_phase', + 'zero_phase', + 'conj_reflect', + 'cross_correlation', + 'check_get_conj_reflect', + 'get_metric', + 'all_metrics', + 'dftups', + 'dftregistration', + 'register_3d_reconstruction', + 'sub_pixel_shift', + 'align_arrays_subpixel', + 'align_arrays_pixel', + 'fast_shift', + 'correlation_err', + 'remove_ramp', + 'zero_phase_cc', + 'breed' + ] def set_lib(dlib): @@ -6,22 +35,25 @@ def set_lib(dlib): dvclib = dlib -def crop_center(arr, shape): - dims = dvclib.dims(arr) - principio = [] - finem = [] - for i in range(3): - principio.append(int((dims[i] - shape[i]) / 2)) - finem.append(principio[i] + shape[i]) - if len(shape) == 1: - cropped = arr[principio[0]: finem[0]] - elif len(shape) == 2: - cropped = arr[principio[0]: finem[0], principio[1]: finem[1]] - elif len(shape) == 3: - cropped = arr[principio[0]: finem[0], principio[1]: finem[1], principio[2]: finem[2]] - else: - raise NotImplementedError - return cropped +def arr_property(arr): + """ + Used only in development. Prints max value of the array and max coordinates. + + Parameters + ---------- + arr : ndarray + array to find max + """ + import numpy as np + arr1 = abs(arr) + print('norm', dvclib.sum(pow(abs(arr), 2))) + print('mean across all data', arr1.mean()) + print('std all data', arr1.std()) + print('sum',np.sum(arr)) + # print('mean non zeros only', arr1[np.nonzero(arr1)].mean()) + # print('std non zeros only', arr1[np.nonzero(arr1)].std()) + max_coordinates = list(dvclib.unravel_index(dvclib.argmax(arr1), arr.shape)) + print('max coords, value', max_coordinates, arr[max_coordinates[0], max_coordinates[1], max_coordinates[2]]) def pad_around(arr, shape, val=0): @@ -43,6 +75,34 @@ def pad_around(arr, shape, val=0): return padded +def center_sync(image, support): + """ + Shifts the image and support arrays so the center of mass is in the center of array. + Parameters + ---------- + image, support : ndarray, ndarray + image and support arrays to evaluate and shift + Returns + ------- + image, support : ndarray, ndarray + shifted arrays + """ + shape = image.shape + + # set center phase to zero, use as a reference + phi0 = math.atan2(image.flatten().imag[int(image.flatten().shape[0] / 2)], + image.flatten().real[int(image.flatten().shape[0] / 2)]) + image = image * cmath.exp(-1j * phi0) + + com = dvclib.center_of_mass(dvclib.astype(support, dtype='int32')) + # place center of mass in the center + sft = [shape[i] / 2.0 - com[i] + .5 for i in range(len(shape))] + image = dvclib.roll(image, sft) + support =dvclib.roll(support, sft) + + return image, support + + def gauss_conv_fft(arr, distribution): dims = dvclib.dims(arr) arr_sum = dvclib.sum(arr) @@ -66,7 +126,8 @@ def shrink_wrap(arr, threshold, sigma): subject array threshold : float - support is formed by points above this valuue + support is formed by points above this value + sigma: float sigmas : list sigmas in all dimensions @@ -207,14 +268,18 @@ def get_metric(image, errs, metric_type): calculated metric for a given type """ if metric_type == 'chi': - return errs[-1] + eval = errs[-1] elif metric_type == 'sharpness': - return dvclib.sum(dvclib.square(dvclib.square(dvclib.absolute(image)))) + eval = dvclib.sum(dvclib.square(dvclib.square(dvclib.absolute(image)))) elif metric_type == 'summed_phase': ph = shift_phase(image, 0) - return dvclib.sum(abs(ph)) + eval = dvclib.sum(abs(ph)) elif metric_type == 'area': - return dvclib.sum(shrink_wrap(image, .2, .5)) + eval = dvclib.sum(shrink_wrap(image, .2, .5)) + + if type(eval) != float: + eval = eval.item() + return eval def all_metrics(image, errs): @@ -232,10 +297,23 @@ def all_metrics(image, errs): calculated metrics for all types """ metrics = {} - metrics['chi'] = errs[-1] - metrics['sharpness'] = dvclib.sum(dvclib.square(dvclib.square(dvclib.absolute(image)))) - metrics['summed_phase'] = dvclib.sum(abs(shift_phase(image, 0))) - metrics['area'] = dvclib.sum(shrink_wrap(image, .2, .5)) + eval = errs[-1] + if type(eval) != float: + eval = eval.item() + metrics['chi'] = eval + eval = dvclib.sum(dvclib.square(dvclib.square(dvclib.absolute(image)))) + if type(eval) != float: + eval = eval.item() + metrics['sharpness'] = eval + eval = dvclib.sum(abs(shift_phase(image, 0))) + if type(eval) != float: + eval = eval.item() + metrics['summed_phase'] = eval + eval = dvclib.sum(shrink_wrap(image, .2, .5)) + if type(eval) != float: + eval = eval.item() + metrics['area'] = eval + return metrics @@ -317,7 +395,7 @@ def dftregistration(ref_arr, arr, usfac=2): c_c = pad_around(dvclib.fftshift(ref_arr) * dvclib.conj(dvclib.fftshift(arr)), large_shape, val=0j) # Compute crosscorrelation and locate the peak - c_c = dvclib.ifft(dvclib.ifftshift(c_c)) #TODO will not work for aflib + c_c = dvclib.ifft(dvclib.ifftshift(c_c)) max_coord = dvclib.unravel_index(dvclib.argmax(c_c), dvclib.dims(c_c)) if max_coord[0] > shape[0]: @@ -371,9 +449,9 @@ def register_3d_reconstruction(ref_arr, arr): shift_2, shift_1, shift_0 : float, float pixel shifts between images """ - r_shift_2, c_shift_2 = dftregistration(dvclib.fft(dvclib.sum(ref_arr, 2)), dvclib.fft(dvclib.sum(arr, 2)), 100) - r_shift_1, c_shift_1 = dftregistration(dvclib.fft(dvclib.sum(ref_arr, 1)), dvclib.fft(dvclib.sum(arr, 1)), 100) - r_shift_0, c_shift_0 = dftregistration(dvclib.fft(dvclib.sum(ref_arr, 0)), dvclib.fft(dvclib.sum(arr, 0)), 100) + r_shift_2, c_shift_2 = dftregistration(dvclib.fft(dvclib.sum(ref_arr, 2)), dvclib.fft(dvclib.sum(arr, 2)), 4) + r_shift_1, c_shift_1 = dftregistration(dvclib.fft(dvclib.sum(ref_arr, 1)), dvclib.fft(dvclib.sum(arr, 1)), 4) + r_shift_0, c_shift_0 = dftregistration(dvclib.fft(dvclib.sum(ref_arr, 0)), dvclib.fft(dvclib.sum(arr, 0)), 4) shift_2 = sum([r_shift_2, r_shift_1]) * 0.5 shift_1 = sum([c_shift_2, r_shift_0]) * 0.5 @@ -411,7 +489,7 @@ def sub_pixel_shift(arr, row_shift, col_shift, z_shift): return dvclib.ifft(Greg) -def align_arrays(ref_arr, arr): +def align_arrays_subpixel(ref_arr, arr): """ Shifts array to align with reference array. @@ -432,6 +510,128 @@ def align_arrays(ref_arr, arr): return sub_pixel_shift(arr, shift_2, shift_1, shift_0) +def remove_ramp(arr, ups=1): + """ + Smooths image array (removes ramp) by applaying math formula. + Parameters + ---------- + arr : ndarray + array to remove ramp + ups : int + upsample factor + Returns + ------- + ramp_removed : ndarray + smoothed array + """ + # the remove ramp is called (now) by visualization when arrays are on cpu + # thus using functions from utilities will be ok for now + import cohere_core.utilities.utils as ut + shape = arr.shape + new_shape = [dim * ups for dim in shape] + padded = ut.pad_center(arr, new_shape) + padded_f = dvclib.fftshift(dvclib.fft(dvclib.ifftshift(padded))) + com = dvclib.center_of_mass(dvclib.square(dvclib.absolute(padded_f))) + sft = [new_shape[i] / 2.0 - com[i] + .5 for i in range(len(shape))] + sub_pixel_shifted = sub_pixel_shift(padded_f, *sft) + ramp_removed_padded = dvclib.fftshift(dvclib.ifft(dvclib.fftshift(sub_pixel_shifted))) + ramp_removed = ut.crop_center(ramp_removed_padded, arr.shape) + + return ramp_removed + + +# supposedly this is faster than np.roll or scipy interpolation shift. +# https://stackoverflow.com/questions/30399534/shift-elements-in-a-numpy-array +def fast_shift(arr, shifty, fill_val=0): + """ + Shifts array by the shifty parameter. + Parameters + ---------- + arr : ndarray + array to shift + shifty : ndarray + an array of integers/shifts in each dimension + fill_val : float + + Returns + ------- + ndarray, shifted array + """ + dims = arr.shape + result = dvclib.full(dims, 1.0) + result *= fill_val + result_slices = [] + arr_slices = [] + for n in range(len(dims)): + if shifty[n] > 0: + result_slices.append(slice(shifty[n], dims[n])) + arr_slices.append(slice(0, -shifty[n])) + elif shifty[n] < 0: + result_slices.append(slice(0, shifty[n])) + arr_slices.append(slice(-shifty[n], dims[n])) + else: + result_slices.append(slice(0, dims[n])) + arr_slices.append(slice(0, dims[n])) + result_slices = tuple(result_slices) + arr_slices = tuple(arr_slices) + result[result_slices] = arr[arr_slices] + return result + + +def align_arrays_pixel(ref, arr): + """ + Aligns two arrays of the same dimensions with the pixel resolution. Used in reciprocal space. + Parameters + ---------- + ref : ndarray + array to align with + arr : ndarray + array to align + + Returns + ------- + ndarray : aligned array + """ + CC = dvclib.correlate(ref, arr, mode='same', method='fft') + err = correlation_err(ref, arr, CC) + CC_shifted = dvclib.ifftshift(CC) + shape = dvclib.array(CC_shifted.shape) + amp = dvclib.absolute(CC_shifted) + shift = dvclib.unravel_index(amp.argmax(), shape) + if dvclib.sum(shift) == 0: + return [arr, err] + intshift = dvclib.array(shift) + pixelshift = dvclib.where(intshift >= shape / 2, intshift - shape, intshift) + shifted_arr = fast_shift(arr, pixelshift) + return [shifted_arr, err] + + +def correlation_err(refarr, arr, CC=None): + """ + author: Paul Frosik + The method is based on "Invariant error metrics for image reconstruction" by J. R. Fienup. + Returns correlation error between two arrays. + + Parameters + ---------- + refarr : ndarray + reference array + arr : ndarray + array to measure likeness with reference array + + Returns + ------- + float, correlation error + """ + if CC is None: + CC = dvclib.correlate(refarr, arr, mode='same', method='fft') + CCmax = CC.max() + rfzero = dvclib.sum(dvclib.absolute(refarr) ** 2) + rgzero = dvclib.sum(dvclib.absolute(arr) ** 2) + e = abs(1 - CCmax * dvclib.conj(CCmax) / (rfzero * rgzero)) ** 0.5 + return e + + def zero_phase_cc(arr1, arr2): """ Sets avarage phase in array to reference array. @@ -472,7 +672,7 @@ def breed(breed_mode, alpha_dir, image): """ # load alpha from alpha dir - alpha = dvclib.load(alpha_dir + '/image.npy') + alpha = dvclib.load(ut.join(alpha_dir, 'image.npy')) if dvclib.array_equal(image, alpha): # it is alpha, no breeding return zero_phase(image) @@ -482,7 +682,7 @@ def breed(breed_mode, alpha_dir, image): beta = image beta = zero_phase(beta) alpha = check_get_conj_reflect(beta, alpha) - alpha = align_arrays(beta, alpha) + alpha = align_arrays_subpixel(beta, alpha) alpha = zero_phase(alpha) ph_alpha = dvclib.angle(alpha) beta = zero_phase_cc(beta, alpha) diff --git a/cohere_core/utilities/ga_utils.py b/cohere_core/utilities/ga_utils.py new file mode 100644 index 0000000..f1295c1 --- /dev/null +++ b/cohere_core/utilities/ga_utils.py @@ -0,0 +1,183 @@ +import os +import cohere_core.utilities.utils as ut + +class Tracing: + def __init__(self, reconstructions, pars, dir): + self.init_dirs = [] + self.report_tracing = [] + + pars['init_guess'] = pars.get('init_guess', 'random') + if pars['init_guess'] == 'continue': + continue_dir = pars['continue_dir'] + for sub in os.listdir(continue_dir): + image, support, coh = ut.read_results(ut.join(continue_dir, sub)) + if image is not None: + self.init_dirs.append(ut.join(continue_dir, sub)) + self.report_tracing.append([ut.join(continue_dir, sub)]) + if len(self.init_dirs) < reconstructions: + for i in range(reconstructions - len(self.init_dirs)): + self.report_tracing.append([f'random{str(i)}']) + self.init_dirs = self.init_dirs + (reconstructions - len(self.init_dirs)) * [None] + elif pars['init_guess'] == 'AI_guess': + self.report_tracing.append(['AI_guess']) + for i in range(reconstructions - 1): + self.report_tracing.append([f'random{str(i)}']) + self.init_dirs = [ut.join(dir, 'results_AI')] + (reconstructions - 1) * [None] + else: + for i in range(reconstructions): + self.init_dirs.append(None) + self.report_tracing.append([f'random{str(i)}']) + + + def set_map(self, map): + self.map = map + + + def append_gen(self, gen_ranks): + for key in gen_ranks: + self.report_tracing[self.map[key]].append(gen_ranks[key]) + + + def pretty_format_results(self): + """ + Takes in a list of report traces and formats them into human-readable tables. + Performs data conversion in 1 pass and formatting in a second to determine + padding and spacing schemes. + + Parameters + ---------- + none + + Returns + ------- + report_output : str + a string containing the formatted report + """ + col_gap = 2 + + num_gens = len(self.report_tracing[0]) - 1 + fitnesses = list(self.report_tracing[0][1][1].keys()) + + report_table = [] + report_table.append(['start'] + [f'generation {i}' for i in range(num_gens)]) + report_table.append([''] * len(report_table[0])) + + data_col_width = 15 + start_col_width = 15 + for pop_data in self.report_tracing: + report_table.append([str(pop_data[0])] + [str(ind_data[0]) for ind_data in pop_data[1:]]) + start_col_width = max(len(pop_data[0]), start_col_width) + + for fit in fitnesses: + fit_row = [''] + for ind_data in pop_data[1:]: + data_out = f'{fit} : {ind_data[1][fit]}' + data_col_width = max(len(data_out), data_col_width) + fit_row.append(data_out) + report_table.append(fit_row) + report_table.append([''] * len(report_table[0])) + + report_str = '' + linesep = os.linesep + for row in report_table: + report_str += row[0].ljust(start_col_width + col_gap) + report_str += (' ' * col_gap).join([cell.ljust(data_col_width) for cell in row[1:]]) + linesep + + return report_str + + + def save(self, save_dir): + try: + report_str = self.pretty_format_results() + except Exception as e: + print(f'WARNING: Report formatting failed due to {type(e)}: {e}! Falling back to raw formatting.') + report_str = '\n'.join([str(l) for l in self.report_tracing]) + + with open(ut.join(save_dir, 'ranks.txt'), 'w+') as rank_file: + rank_file.write(report_str) + rank_file.flush() + + + +def set_ga_defaults(pars): + pars['reconstructions'] = pars.get('reconstructions', 1) + pars['ga_generations'] = pars.get('ga_generations', 1) + pars['init_guess'] = pars.get('init_guess', 'random') + + # check if pc feature is on + if 'pc' in pars['algorithm_sequence'] and 'pc_interval' in pars: + if not 'ga_gen_pc_start' in pars: + pars['ga_gen_pc_start'] = 0 + pars['ga_gen_pc_start'] = min(pars['ga_gen_pc_start'], pars['ga_generations']-1) + + pars['ga_fast'] = pars.get('ga_fast', False) + + if 'ga_metrics' not in pars: + metrics = ['chi'] * pars['ga_generations'] + else: + metrics = pars['ga_metrics'] + if len(metrics) == 1: + metrics = metrics * pars['ga_generations'] + elif len(metrics) < pars['ga_generations']: + metrics = metrics + ['chi'] * (pars['ga_generations'] - len(metrics)) + pars['ga_metrics'] = metrics + + ga_reconstructions = [] + if 'ga_cullings' in pars: + worst_remove_no = pars['ga_cullings'] + if len(worst_remove_no) < pars['ga_generations']: + worst_remove_no = worst_remove_no + [0] * (pars['ga_generations'] - len(worst_remove_no)) + else: + worst_remove_no = [0] * pars['ga_generations'] + pars['worst_remove_no'] = worst_remove_no + # calculate how many reconstructions should continue + reconstructions = pars['reconstructions'] + for culling in worst_remove_no: + reconstructions = reconstructions - culling + if reconstructions <= 0: + return 'culled down to 0 reconstructions, check configuration' + ga_reconstructions.append(reconstructions) + pars['ga_reconstructions'] = ga_reconstructions + + sw_threshold = .1 + if 'ga_sw_thresholds' in pars: + ga_sw_thresholds = pars['ga_sw_thresholds'] + if len(ga_sw_thresholds) == 1: + ga_sw_thresholds = ga_sw_thresholds * pars['ga_generations'] + elif len(ga_sw_thresholds) < pars['ga_generations']: + ga_sw_thresholds = ga_sw_thresholds + [sw_threshold] * (pars['ga_generations'] - len(ga_sw_thresholds)) + else: + ga_sw_thresholds = [sw_threshold] * pars['ga_generations'] + pars['ga_sw_thresholds'] = ga_sw_thresholds + + sw_gauss_sigma = 1.0 + if 'ga_sw_gauss_sigmas' in pars: + ga_sw_gauss_sigmas = pars['ga_sw_gauss_sigmas'] + if len(ga_sw_gauss_sigmas) == 1: + ga_sw_gauss_sigmas = ga_sw_gauss_sigmas * pars['ga_generations'] + elif len(pars['ga_sw_gauss_sigmas']) < pars['ga_generations']: + ga_sw_gauss_sigmas = ga_sw_gauss_sigmas + [sw_gauss_sigma] * (pars['ga_generations'] - len(ga_sw_gauss_sigmas)) + else: + ga_sw_gauss_sigmas = [sw_gauss_sigma] * pars['ga_generations'] + pars['ga_sw_gauss_sigmas'] = ga_sw_gauss_sigmas + + if 'ga_breed_modes' not in pars: + ga_breed_modes = ['sqrt_ab'] * pars['ga_generations'] + else: + ga_breed_modes = pars['ga_breed_modes'] + if len(ga_breed_modes) == 1: + ga_breed_modes = ga_breed_modes * pars['ga_generations'] + elif len(ga_breed_modes) < pars['ga_generations']: + ga_breed_modes = ga_breed_modes + ['sqrt_ab'] * (pars['ga_generations'] - len(ga_breed_modes)) + pars['ga_breed_modes'] = ga_breed_modes + + if 'ga_lpf_sigmas' in pars: + pars['low_resolution_generations'] = len(pars['ga_lpf_sigmas']) + else: + pars['low_resolution_generations'] = 0 + + if pars['low_resolution_generations'] > 0: + pars['low_resolution_alg'] = pars.get('low_resolution_alg', 'GAUSS') + + return pars + diff --git a/cohere_core/utilities/host_utils.py b/cohere_core/utilities/host_utils.py new file mode 100644 index 0000000..e78ec9a --- /dev/null +++ b/cohere_core/utilities/host_utils.py @@ -0,0 +1,36 @@ +import argparse +import socket +import cohere_core.utilities as ut + + +def main(): + import ast + + parser = argparse.ArgumentParser() + parser.add_argument("devices", help="configured devices") + parser.add_argument("run_size", help="memory requirement for one job") + args = parser.parse_args() + + # devices parameter is dict with host names as keys + devs = ast.literal_eval(args.devices) + run_size = ast.literal_eval(args.run_size) + + # host name in configuration can be full host name or a first part of full name + # it needs to be matched with the hostname returned by socket utilities. + configured_hosts = devs.keys() + host_name = socket.gethostname() + if host_name in configured_hosts: + use_host_name = host_name + else: + use_host_name = host_name.split('.')[0] + if not use_host_name in configured_hosts: + return + + # The result available is a dictionary with key, value pars of GPU ID, number pf available jobs + available = ut.get_avail_gpu_runs(devs[use_host_name], run_size) + # the printed hostname and available devices will be received through pipe object by calling code + print([use_host_name, available]) + + +if __name__ == "__main__": + exit(main()) \ No newline at end of file diff --git a/cohere_core/utilities/utils.py b/cohere_core/utilities/utils.py index e8437a0..d495e18 100644 --- a/cohere_core/utilities/utils.py +++ b/cohere_core/utilities/utils.py @@ -16,30 +16,49 @@ import os import logging import stat +import scipy.ndimage as ndi +import GPUtil +from functools import reduce +import subprocess +import ast __author__ = "Barbara Frosik" __copyright__ = "Copyright (c), UChicago Argonne, LLC." __docformat__ = 'restructuredtext en' -__all__ = ['get_logger', +__all__ = ['join', + 'get_logger', 'read_tif', 'save_tif', 'read_config', - 'get_good_dim', - 'binning', - 'get_centered', - 'get_zero_padded_centered', - 'adjust_dimensions', - 'gaussian', - 'gauss_conv_fft', + 'write_config', 'read_results', + 'save_results', 'save_metrics', 'write_plot_errors', - 'save_results', - 'arr_property', + 'get_good_dim', + 'threshold_by_edge', + 'select_central_object', + 'get_central_object_extent', + 'get_oversample_ratio', + 'Resize', + 'binning', + 'crop_center', + 'pad_center', + 'center_max', + 'adjust_dimensions', + 'normalize', + 'get_avail_gpu_runs', + 'get_one_dev', + 'get_gpu_use', + 'estimate_no_proc' ] +def join(*args): + return os.path.join(*args).replace(os.sep, '/') + + def get_logger(name, ldir=''): """ Creates looger instance that will write to default.log file in a given directory. @@ -56,7 +75,7 @@ def get_logger(name, ldir=''): """ logger = logging.getLogger(name) logger.setLevel(logging.DEBUG) - log_file = ldir + '/default.log' + log_file = join(ldir, 'default.log') fh = logging.FileHandler(log_file) fh.setLevel(logging.DEBUG) formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s') @@ -79,8 +98,7 @@ def read_tif(filename): ndarray an array containing the data parsed from the file """ - ar = tf.imread(filename.replace(os.sep, '/')).transpose() - return ar + return tf.imread(filename.replace(os.sep, '/')).transpose() def save_tif(arr, filename): @@ -94,7 +112,7 @@ def save_tif(arr, filename): filename : str tif format file name """ - tf.imsave(filename.replace(os.sep, '/'), arr.transpose().astype(np.float32)) + tf.imsave(filename.replace(os.sep, '/'), np.abs(arr).transpose().astype(np.float32)) def read_config(config): @@ -111,8 +129,6 @@ def read_config(config): dict dictionary containing parsed configuration, None if the given file does not exist """ - import ast - config = config.replace(os.sep, '/') if not os.path.isfile(config): print(config, 'is not a file') @@ -124,7 +140,7 @@ def read_config(config): while line: # Ignore comment lines and move along line = line.strip() - if line.startswith('//') or line.startswith('#'): + if line.startswith('/') or line.startswith('#'): line = input.readline() continue elif "=" in line: @@ -139,6 +155,160 @@ def read_config(config): return param_dict +def write_config(param_dict, config): + """ + Writes configuration to a file. + Parameters + ---------- + param_dict : dict + dictionary containing configuration parameters + config : str + configuration name theparameters will be written into + """ + with open(config.replace(os.sep, '/'), 'w+') as cf: + cf.truncate(0) + linesep = os.linesep + for key, value in param_dict.items(): + if type(value) == str: + value = f'"{value}"' + cf.write(f'{key} = {str(value)}{linesep}') + + +def read_results(read_dir): + """ + Reads results and returns array representation. + + Parameters + ---------- + read_dir : str + directory to read the results from + Returns + ------- + ndarray, ndarray, ndarray (or None) + image, support, and coherence arrays + """ + try: + imagefile = join(read_dir, 'image.npy') + image = np.load(imagefile) + except: + image = None + + try: + supportfile = join(read_dir, 'support.npy') + support = np.load(supportfile) + except: + support = None + + try: + cohfile = join(read_dir, 'coherence.npy') + coh = np.load(cohfile) + except: + coh = None + + return image, support, coh + + +def save_results(image, support, coh, errs, save_dir, metric=None): + """ + Saves results of reconstruction. Saves the following files: image.np, support.npy, errors.npy, + optionally coherence.npy, plot_errors.py, graph.npy, flow.npy, iter_array.npy + + + Parameters + ---------- + image : ndarray + reconstructed image array + + support : ndarray + support array related to the image + + coh : ndarray + coherence array when pcdi feature is active, None otherwise + + errs : ndarray + errors "chi" by iterations + + save_dir : str + directory to write the files + + metrics : dict + dictionary with metric type keys, and metric values + """ + if not os.path.exists(save_dir): + os.makedirs(save_dir) + image_file = join(save_dir, 'image') + np.save(image_file, image) + support_file = join(save_dir, 'support') + np.save(support_file, support) + + errs_file = join(save_dir, 'errors') + np.save(errs_file, errs) + if not coh is None: + coh_file = join(save_dir, 'coherence') + np.save(coh_file, coh) + + write_plot_errors(save_dir) + + save_metrics(errs, save_dir, metric) + + +def save_metrics(errs, dir, metrics=None): + """ + Saves arrays metrics and errors by iterations in text file. + + Parameters + ---------- + errs : list + list of "chi" error by iteration + + dir : str + directory to write the file containing array metrics + + metrics : dict + dictionary with metric type keys, and metric values + """ + metric_file = join(dir, 'summary') + linesep = os.linesep + with open(metric_file, 'w+') as mf: + if metrics is not None: + mf.write(f'metric result{linesep}') + for key in metrics: + value = metrics[key] + mf.write(f'{key} = {str(value)}{linesep}') + mf.write(f'{linesep}errors by iteration{linesep}') + for er in errs: + mf.write(f'str({er}) ') + mf.close() + + +def write_plot_errors(save_dir): + """ + Creates python executable that draw plot of error by iteration. It assumes that the given directory + contains "errors.npy" file + + Parameters + ---------- + save_dir : str + directory containing errors.npy file + """ + plot_file = join(save_dir, 'plot_errors.py') + f = open(plot_file, 'w+') + f.write("#! /usr/bin/env python\n") + f.write("import matplotlib.pyplot as plt\n") + f.write("import numpy as np\n") + f.write("import sys\n") + f.write("import os\n") + f.write("current_dir = sys.path[0]\n") + f.write("errs = np.load(current_dir + '/errors.npy').tolist()\n") + f.write("errs.pop(0)\n") + f.write("plt.plot(errs)\n") + f.write("plt.ylabel('errors')\n") + f.write("plt.show()") + f.close() + st = os.stat(plot_file) + os.chmod(plot_file, st.st_mode | stat.S_IEXEC) + + def get_good_dim(dim): """ This function calculates the dimension supported by opencl library (i.e. is multiplier of 2, 3, or 5) and is closest to the given starting dimension. @@ -174,6 +344,103 @@ def is_correct(x): return new_dim +def threshold_by_edge(fp: np.ndarray) -> np.ndarray: + """ + Author: Yudong Yao + + :param fp: + :return: + """ + # threshold by left edge value + mask = np.ones_like(fp, dtype=bool) + mask[tuple([slice(1, None)] * fp.ndim)] = 0 + zero = 1e-6 + cut = np.max(fp[mask]) + binary = np.zeros_like(fp) + binary[(np.abs(fp) > zero) & (fp > cut)] = 1 + return binary + + +def select_central_object(fp: np.ndarray) -> np.ndarray: + """ + :param fp: + :return: + """ + # import scipy.ndimage as ndimage + zero = 1e-6 + binary = np.abs(fp) + binary[binary > zero] = 1 + binary[binary <= zero] = 0 + + # cluster by connectivity + struct = ndi.morphology.generate_binary_structure(fp.ndim, 1).astype("uint8") + label, nlabel = ndi.label(binary, structure=struct) + + # select largest cluster + select = np.argmax(np.bincount(np.ravel(label))[1:]) + 1 + + binary[label != select] = 0 + + fp[binary == 0] = 0 + return fp + + +def get_central_object_extent(fp: np.ndarray) -> list: + """ + :param fp: + :return: + """ + fp_cut = threshold_by_edge(np.abs(fp)) + need = select_central_object(fp_cut) + + # get extend of cluster + extent = [np.max(s) + 1 - np.min(s) for s in np.nonzero(need)] + return extent + + +def get_oversample_ratio(fp: np.ndarray) -> np.ndarray: + """ + :param fp: diffraction pattern + :return: oversample ratio in each dimension + """ + # autocorrelation + acp = np.fft.fftshift(np.fft.ifftn(np.abs(fp)**2.)) + aacp = np.abs(acp) + + # get extent + blob = get_central_object_extent(aacp) + + # correct for underestimation due to thresholding + correction = [0.025, 0.025, 0.0729][:fp.ndim] + + extent = [ + min(m, s + int(round(f * aacp.shape[i], 1))) + for i, (s, f, m) in enumerate(zip(blob, correction, aacp.shape)) + ] + + # oversample ratio + oversample = [ + 2. * s / (e + (1 - s % 2)) for s, e in zip(aacp.shape, extent) + ] + return np.round(oversample, 3) + + +def Resize(IN, dim): + """ + :param IN: + :param dim: + :return: + """ + ft = np.fft.fftshift(np.fft.fftn(IN)) / np.prod(IN.shape) + + pad_value = np.array(dim) // 2 - np.array(ft.shape) // 2 + pad = [[pad_value[0], pad_value[0]], [pad_value[1], pad_value[1]], + [pad_value[2], pad_value[2]]] + ft_resize = adjust_dimensions(ft, pad) + output = np.fft.ifftn(np.fft.ifftshift(ft_resize)) * np.prod(dim) + return output + + def binning(array, binsizes): """ This function does the binning of the array. The array is binned in each dimension by the corresponding binsizes elements. @@ -229,17 +496,43 @@ def get_centered(arr, center_shift): ndarray centered array """ - max_coordinates = list(np.unravel_index(np.argmax(arr), arr.shape)) - max_coordinates = np.add(max_coordinates, center_shift) - shape = arr.shape - centered = arr - for i in range(len(max_coordinates)): - centered = np.roll(centered, int(shape[i] / 2) - max_coordinates[i], i) + shift = (np.array(arr.shape)/2) - np.unravel_index(np.argmax(arr), arr.shape) + np.array(center_shift) + return np.roll(arr, shift.astype(int), tuple(range(arr.ndim))), shift.astype(int) - return centered - -def get_zero_padded_centered(arr, new_shape): +def crop_center(arr, new_shape): + """ + This function crops the array to the new size, leaving the array in the center. + The new_size must be smaller or equal to the original size in each dimension. + Parameters + ---------- + arr : ndarray + the array to crop + new_shape : tuple + new size + Returns + ------- + cropped : ndarray + the cropped array + """ + shape = arr.shape + principio = [] + finem = [] + for i in range(3): + principio.append(int((shape[i] - new_shape[i]) / 2)) + finem.append(principio[i] + new_shape[i]) + if len(shape) == 1: + cropped = arr[principio[0]: finem[0]] + elif len(shape) == 2: + cropped = arr[principio[0]: finem[0], principio[1]: finem[1]] + elif len(shape) == 3: + cropped = arr[principio[0]: finem[0], principio[1]: finem[1], principio[2]: finem[2]] + else: + raise NotImplementedError + return cropped + + +def pad_center(arr, new_shape): """ This function pads the array with zeros to the new shape with the array in the center. Parameters @@ -254,20 +547,42 @@ def get_zero_padded_centered(arr, new_shape): the zero padded centered array """ shape = arr.shape - pad = [] - c_vals = [] - for i in range(len(new_shape)): - pad.append((0, new_shape[i] - shape[i])) - c_vals.append((0.0, 0.0)) - arr = np.lib.pad(arr, (pad), 'constant', constant_values=c_vals) + centered = np.zeros(new_shape, arr.dtype) + if len(shape) == 1: + centered[: shape[0]] = arr + elif len(shape) == 2: + centered[: shape[0], : shape[1]] = arr + elif len(shape) == 3: + centered[: shape[0], : shape[1], : shape[2]] = arr - centered = arr for i in range(len(new_shape)): - centered = np.roll(centered, int((new_shape[i] - shape[i] + 1) / 2), i) + centered = np.roll(centered, (new_shape[i] - shape[i] + 1) // 2, i) return centered +def center_max(arr, center_shift): + """ + This function finds maximum value in the array, and puts it in a center of a new array. If center_shift is not zeros, + the array will be shifted accordingly. The shifted elements are rolled into the other end of array. + + Parameters + ---------- + arr : ndarray + the original array to be centered + + center_shift : list + a list defining shift of the center + + Returns + ------- + ndarray + centered array + """ + shift = (np.array(arr.shape)/2) - np.unravel_index(np.argmax(arr), arr.shape) + np.array(center_shift) + return np.roll(arr, shift.astype(int), tuple(range(arr.ndim))), shift.astype(int) + + def adjust_dimensions(arr, pads): """ This function adds to or subtracts from each dimension of the array elements defined by pad. If the pad is positive, the array is padded in this dimension. If the pad is negative, the array is cropped. @@ -388,142 +703,208 @@ def gauss_conv_fft(arr, sigmas): return convag -def read_results(read_dir): - """ - Reads results and returns array representation. +def normalize(vec): + return vec / np.linalg.norm(vec) - Parameters - ---------- - read_dir : str - directory to read the results from - Returns - ------- - ndarray, ndarray, ndarray (or None) - image, support, and coherence arrays - """ - read_dir = read_dir.replace(os.sep, '/') - try: - imagefile = read_dir + '/image.npy' - image = np.load(imagefile) - except: - image = None - try: - supportfile = read_dir + '/support.npy' - support = np.load(supportfile) - except: - support = None +def get_one_dev(ids): + """ + Returns GPU ID that is included in the configuration, is on a local node, and has the most available memory. - try: - cohfile = read_dir + '/coherence.npy' - coh = np.load(cohfile) - except: - coh = None + :param ids: list or string or dict + list of gpu ids, or string 'all' indicating all GPUs included, or dict by hostname + :return: int + selected GPU ID + """ + import socket - return image, support, coh + # if cluster configuration, look only at devices on local machine + if issubclass(type(ids), dict): # a dict with cluster configuration + ids = ids[socket.gethostname()] # configured devices on local host + os.environ["CUDA_DEVICE_ORDER"] = "PCI_BUS_ID" + gpus = GPUtil.getGPUs() + dev = -1 + max_mem = 0 + # select one with the highest availbale memory + for gpu in gpus: + if ids == 'all' or gpu.id in ids: + free_mem = gpu.memoryFree + if free_mem > max_mem: + dev = gpu.id + max_mem = free_mem + return dev + + +def get_avail_gpu_runs(devices, run_mem): + """ + Finds how many jobs of run_mem size can run on configured GPUs on local host. -def save_metrics(errs, dir, metrics=None): + :param devices: list or string + list of GPU IDs or 'all' if configured to use all available GPUs + :param run_mem: int + size of GPU memory (in MB) needed for one job + :return: dict + pairs of GPU IDs, number of available jobs """ - Saves arrays metrics and errors by iterations in text file. + import GPUtil - Parameters - ---------- - errs : list - list of "chi" error by iteration + os.environ["CUDA_DEVICE_ORDER"] = "PCI_BUS_ID" + gpus = GPUtil.getGPUs() + available = {} - dir : str - directory to write the file containing array metrics + for gpu in gpus: + if devices == 'all' or gpu.id in devices: + available[gpu.id] = gpu.memoryFree // run_mem - metrics : dict - dictionary with metric type keys, and metric values - """ - dir = dir.replace(os.sep, '/') - metric_file = dir + '/summary' - with open(metric_file, 'w+') as f: - if metrics is not None: - f.write('metric result\n') - for key in metrics: - value = metrics[key] - f.write(key + ' = ' + str(value) + '\n') - f.write('\nerrors by iteration\n') - for er in errs: - f.write(str(er) + ' ') - f.close() + return available -def write_plot_errors(save_dir): +def get_avail_hosts_gpu_runs(devices, run_mem): """ - Creates python executable that draw plot of error by iteration. It assumes that the given directory - contains "errors.npy" file + This function is called in a cluster configuration case, i.e. devices parameter is configured as dictionary of hostnames and GPU IDs (either list of the IDs or 'all' for all GPUs per host). + It starts mpi subprocess that targets each of the configured host. The subprocess returns tuples with hostname and available GPUs. The tuples are converted into dictionary and returned. - Parameters - ---------- - save_dir : str - directory containing errors.npy file + :param devices: + :param run_mem: + :return: """ - save_dir = save_dir.replace(os.sep, '/') - plot_file = save_dir + '/plot_errors.py' - f = open(plot_file, 'w+') - f.write("#! /usr/bin/env python\n") - f.write("import matplotlib.pyplot as plt\n") - f.write("import numpy as np\n") - f.write("import sys\n") - f.write("import os\n") - f.write("current_dir = sys.path[0]\n") - f.write("errs = np.load(current_dir + '/errors.npy').tolist()\n") - f.write("errs.pop(0)\n") - f.write("plt.plot(errs)\n") - f.write("plt.ylabel('errors')\n") - f.write("plt.show()") - f.close() - st = os.stat(plot_file) - os.chmod(plot_file, st.st_mode | stat.S_IEXEC) + hosts = ','.join(devices.keys()) + script = join(os.path.realpath(os.path.dirname(__file__)), 'host_utils.py') + command = ['mpiexec', '-n', str(len(devices)), '--host', hosts, 'python', script, str(devices), str(run_mem)] + result = subprocess.run(command, stdout=subprocess.PIPE, text=True).stdout + mem_map = {} + for entry in result.splitlines(): + host_devs = ast.literal_eval(entry) + mem_map[host_devs[0]] = host_devs[1] + return mem_map -def save_results(image, support, coh, errs, save_dir, metric=None): +def get_balanced_load(avail_runs, runs): """ - Saves results of reconstruction. Saves the following files: image.np, support.npy, errors.npy, - optionally coherence.npy, plot_errors.py, graph.npy, flow.npy, iter_array.npy - - - Parameters - ---------- - image : ndarray - reconstructed image array - - support : ndarray - support array related to the image - - coh : ndarray - coherence array when pcdi feature is active, None otherwise - - errs : ndarray - errors "chi" by iterations + This function distributes the runs proportionally to the GPUs availability. + If number of available runs is less or equal to the requested runs, the input parameter avail_runs becomes load. + The function also returns number of available runs. - save_dir : str - directory to write the files - - metrics : dict - dictionary with metric type keys, and metric values + :param avail_runs: dict + keys are GPU IDs, and values are available runs + for cluster configuration the keys are prepended with the hostnames + :param runs: int + number of requested jobs + :return: dict, int + a dictionary with the same structure as avail_runs input parameter, but with values indicating runs modified to achieve balanced distribution. """ - save_dir = save_dir.replace(os.sep, '/') - if not os.path.exists(save_dir): - os.makedirs(save_dir) - image_file = save_dir + '/image' - np.save(image_file, image) - support_file = save_dir + '/support' - np.save(support_file, support) + if len(avail_runs) == 0: + return {} - errs_file = save_dir + '/errors' - np.save(errs_file, errs) - if not coh is None: - coh_file = save_dir + '/coherence' - np.save(coh_file, coh) + # if total number of available runs is less or equal runs, return the avail_runs, + # and total number of available jobs + total_available = reduce((lambda x, y: x + y), avail_runs.values()) + if total_available <= runs: + return avail_runs, total_available - write_plot_errors(save_dir) - - save_metrics(errs, save_dir, metric) + # initialize variables for calculations + need_runs = runs + available = total_available + load = {} + + # add one run from each available + for k, v in avail_runs.items(): + if v > 0: + load[k] = 1 + avail_runs[k] = v - 1 + need_runs -= 1 + if need_runs == 0: + return load, runs + available -= 1 + + # use proportionally from available + distributed = 0 + ratio = need_runs / available + for k, v in avail_runs.items(): + if v > 0: + share = int(v * ratio) + load[k] = load[k] + share + avail_runs[k] = v - share + distributed += share + need_runs -= distributed + available -= distributed + + # need to add the few remaining + for k, v in avail_runs.items(): + if v > 0: + load[k] = load[k] + 1 + need_runs -= 1 + if need_runs == 0: + break + + return load, runs + + +def get_gpu_use(devices, no_jobs, job_size): + """ + Determines available GPUs that match configured devices, and selects the optimal distribution of jobs on available devices. If devices is configured as dict (i.e. cluster configuration) then a file "hosts" is created in the running directory. This file contains hosts names and number of jobs to run on that host. + Parameters + ---------- + devices : list or dict or 'all' + Configured parameter. list of GPU ids to use for jobs or 'all' if all GPUs should be used. If cluster configuration, then + it is dict with keys being host names. + no_jobs : int + wanted number of jobs + job_size : float + a GPU memory requirement to run one job + Returns + ------- + picked_devs : list or list of lists(if cluster conf) + list of GPU ids that were selected for the jobs + available jobs : int + number of jobs allocated on all GPUs + cluster_conf : boolean + True is cluster configuration + """ + + def unpack_load(load): + picked_devs = [] + for ds in [[k] * int(v) for k, v in load.items()]: + picked_devs.extend(ds) + return picked_devs + + if type(devices) != dict: # a configuration for local host + hostfile_name = None + avail_jobs = get_avail_gpu_runs(devices, job_size) + balanced_load, avail_jobs_no = get_balanced_load(avail_jobs, no_jobs) + picked_devs = unpack_load(balanced_load) + else: # cluster configuration + hosts_avail_jobs = get_avail_hosts_gpu_runs(devices, job_size) + avail_jobs = {} + # collapse the host dict into one dict by adding hostname in front of key (gpu id) + for k, v in hosts_avail_jobs.items(): + host_runs = {(f'{k}_{str(kv)}'): vv for kv, vv in v.items()} + avail_jobs.update(host_runs) + balanced_load, avail_jobs_no = get_balanced_load(avail_jobs, no_jobs) + + # un-collapse the balanced load by hosts + host_balanced_load = {} + for k, v in balanced_load.items(): + idx = k.rfind('_') + host = k[:idx] + if host not in host_balanced_load: + host_balanced_load[host] = {} + host_balanced_load[host].update({int(k[idx + 1:]): v}) + + # create hosts file and return corresponding picked devices + hosts_picked_devs = [(k, unpack_load(v)) for k, v in host_balanced_load.items()] + + picked_devs = [] + hostfile_name = f'hostfile_{os.getpid()}' + host_file = open(hostfile_name, mode='w+') + linesep = os.linesep + for h, ds in hosts_picked_devs: + host_file.write(f'{h}:{str(len(ds))}{linesep}') + picked_devs.append(ds) + host_file.close() + + return picked_devs, min(avail_jobs_no, no_jobs), hostfile_name def arr_property(arr): @@ -536,9 +917,13 @@ def arr_property(arr): array to find max """ arr1 = abs(arr) - print('norm', np.sum(pow(abs(arr), 2))) - max_coordinates = list(np.unravel_index(np.argmax(arr1), arr.shape)) - print('max coords, value', max_coordinates, arr[max_coordinates[0], max_coordinates[1], max_coordinates[2]]) + #print('norm', np.sum(pow(abs(arr), 2))) + print('mean across all data', arr1.mean()) + print('std all data', arr1.std()) + print('mean non zeros only', arr1[np.nonzero(arr1)].mean()) + print('std non zeros only', arr1[np.nonzero(arr1)].std()) + # max_coordinates = list(np.unravel_index(np.argmax(arr1), arr.shape)) + # print('max coords, value', max_coordinates, arr[max_coordinates[0], max_coordinates[1], max_coordinates[2]]) # https://stackoverflow.com/questions/51503672/decorator-for-timeit-timeit-method/51503837#51503837 @@ -546,6 +931,34 @@ def arr_property(arr): from time import time +def estimate_no_proc(arr_size, factor): + """ + Estimates number of processes the prep can be run on. Determined by number of available cpus and size + of array. + Parameters + ---------- + arr_size : int + size of array + factor : int + an estimate of how much memory is required to process comparing to array size + Returns + ------- + int + number of processes + """ + from multiprocessing import cpu_count + import psutil + + ncpu = cpu_count() + freemem = psutil.virtual_memory().available + nmem = freemem / (factor * arr_size) + # decide what limits, ncpu or nmem + if nmem > ncpu: + return ncpu + else: + return int(nmem) + + def measure(func): @wraps(func) def _time_it(*args, **kwargs): @@ -556,4 +969,4 @@ def _time_it(*args, **kwargs): end_ = int(round(time() * 1000)) - start print("Total execution time: {end_ if end_ > 0 else 0} ms") - return _time_it \ No newline at end of file + return _time_it diff --git a/docs/source/about.rst b/docs/source/about.rst index 1580950..d1db27f 100755 --- a/docs/source/about.rst +++ b/docs/source/about.rst @@ -8,8 +8,8 @@ The reconstruction has very good performance, in particular when utilizing GPU. A powerful features that can deliver good reconstruction result offered by the cohere package is genetic algorithm (GA) or AI guess. -The tools are supplemented by another package cohere-ui which contains users scripts ready to use and easy to modify. Refer to :ref:`use` for instructions on how to use the supplemental software. Combined together both of the packages: cohere and cohere-tools offer a full solution from beamline preprocessing data followed by standard preprocessing data, reconstruction, and beamline visualization. The project was implemented for the Advanced Photon Source beamline 34-IDC and thus the beamline data preparation and data visualization is customized for that specific beamline. The measurements and parameters that were used during the experiment are parsed specifically for the beamline setup. We offer the code, as the community would benefit from it when customizing for own case. +The tools are supplemented by another package cohere-ui which contains users scripts ready to use and easy to modify. Refer to :ref:`use` for instructions on how to use the supplemental software. Combined together both of the packages: cohere and cohere-ui offer a full solution from beamline preprocessing data followed by standard preprocessing data, reconstruction, and beamline visualization. The project was implemented for the Advanced Photon Source beamline 34-IDC and thus the beamline data preparation and data visualization is customized for that specific beamline. The measurements and parameters that were used during the experiment are parsed specifically for that beamline setup. We offer the code, as the community would benefit from it when customizing for own cases. -The cohere is written in Python as this choice offers simplicity of development, and thus the tool can be expanded by scintific community, +The cohere is written in Python as this choice offers simplicity of development, and thus the tool can be expanded by community. A legacy version (1.2.4) written in Python/c++ is still available in the legacy branch, and available from anaconda but there is no new development in this branch. \ No newline at end of file diff --git a/docs/source/conf.py b/docs/source/conf.py index 3ec0c55..8752fc6 100755 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -23,7 +23,7 @@ author = 'Barbara Frosik, Ross Harder' # The full version, including alpha/beta/rc tags -# release = '2.0' +release = '4.1-beta1' # -- General configuration --------------------------------------------------- @@ -62,6 +62,15 @@ # so a file named "default.css" will overwrite the builtin "default.css". #html_static_path = ['_static'] -autodoc_mock_imports = ['tifffile', 'scipy'] +autodoc_mock_imports = ['numpy', + 'scikit-learn', + 'tifffile', + 'tensorflow', + 'mpi4py', + 'GPUtil', + 'psutil', + 'tqdm', + 'scipy', + 'matplotlib'] master_doc = 'index' diff --git a/docs/source/config_data.rst b/docs/source/config_data.rst index 4b88724..4dd11b0 100755 --- a/docs/source/config_data.rst +++ b/docs/source/config_data.rst @@ -6,21 +6,21 @@ config_data Parameters ========== - data_dir: -| optional, used for specific cases. Driven by the scripts, the experiment directory contains four directories: conf, preprocessed_data, phasing_data, results_phasing, and results_viz. The formatted data is saved in the /data directory. If the data_dir parameter is configured then the data is saved in this directory. Warning: The script running reconstruction expects to read data from /data. +| Optional, defines directory that contains data.tif file with experiment data. Defaults to /data directory. | example: :: - data_dir = "/path/to/data_dir/data_dir" + data_dir = "/path/to/data_dir" - alien_alg: -| optional, name of method used to remove aliens. Possible options are: 'block_aliens', 'alien_file', and 'AutoAlien1'. The 'block_aliens' algorithm will zero out defined blocks, 'alien_file' method will use given file as a mask, and 'AutoAlien1' will use auto mechanism to remove aliens. Each of these algorithms require different parameters, explained below. +| Optional, name of method used to remove aliens. Possible options are: 'block_aliens', 'alien_file', and 'AutoAlien1'. The 'block_aliens' algorithm will zero out defined blocks, 'alien_file' method will use given file as a mask, and 'AutoAlien1' will use auto mechanism to remove aliens. Each of these algorithms require different parameters, explained below. | example: :: alien_alg = "AutoAlien1" - aliens: -| Needed when the 'block_aliens' method is configured. Used when the data contains regions with intensity produced by interference. The regions needs to be zeroed out. The aliens can be defined as regions each defined by coordinates of starting point, and ending point (i.e. [[xb0,yb0,zb0,xe0,ye0,ze0],[xb1,yb1,zb1,xe1,ye1,ze1],...[xbn,ybn,zbn,xen,yen,zen]] ). +| Needed when the 'block_aliens' method is configured. Used when the data contains regions with intensity produced by interference. The regions are zeroed out. The aliens can be defined as regions, each region defined by coordinates of starting point, and ending point (i.e. [[xb0,yb0,zb0,xe0,ye0,ze0],[xb1,yb1,zb1,xe1,ye1,ze1],...[xbn,ybn,zbn,xen,yen,zen]] ). | example: :: @@ -34,14 +34,14 @@ Parameters alien_file = "/path/to/mask_file/AlienImg.npy" - AA1_size_threshold: -| used in the 'AutoAliens1' method. If not given it will default to 0.01. The AutoAlien1 algorithm will calculate relative sizes of all clusters with respect to the biggest cluster. The clusters with relative size smaller than the given threshold will be possibly deemed aliens. It also depends on asymmetry. +| Used in the 'AutoAliens1' method. If not given it will default to 0.01. The AutoAlien1 algorithm will calculate relative sizes of all clusters with respect to the biggest cluster. The clusters with relative size smaller than the given threshold will be possibly deemed aliens. It also depends on asymmetry. | example: :: AA1_size_threshold = 0.01 - AA1_asym_threshold: -| used in the 'AutoAliens1' method. If not given it will default to 1.75. The AutoAlien1 algorithm will calculate average asymmetry of all clusters. The clusters with average asymmetry greater than the given threshold will be possibly deemed aliens. It also depends on relative size. +| Used in the 'AutoAliens1' method. If not given it will default to 1.75. The AutoAlien1 algorithm will calculate average asymmetry of all clusters. The clusters with average asymmetry greater than the given threshold will be possibly deemed aliens. It also depends on relative size. | example: :: @@ -55,56 +55,58 @@ Parameters AA1_min_pts = 5 - AA1_eps: -| used in the 'AutoAliens1' method. If not given it will default to 1.1. Used in the clustering algorithm. +| Used in the 'AutoAliens1' method. If not given it will default to 1.1. Used in the clustering algorithm. | example: :: AA1_eps = 1.1 - AA1_amp_threshold: -| mandatory in the 'AutoAliens1' method. Used to zero data points below that threshold. +| Mandatory in the 'AutoAliens1' method. Used to zero data points below that threshold. | example: :: AA1_amp_threshold = 6 - AA1_save_arrs -| used in the 'AutoAliens1' method, optional. If given and set to True multiple results of alien analysis will be saved in files. +| Used in the 'AutoAliens1' method, optional. If given and set to True multiple results of alien analysis will be saved in files. | example: :: AA1_save_arrs = True - AA1_expandcleanedsigma: -| used in the 'AutoAliens1' method, optional. If given the algorithm will apply last step of cleaning the data using the configured sigma. +| Used in the 'AutoAliens1' method, optional. If given the algorithm will apply last step of cleaning the data using the configured sigma. | example: :: AA1_expandcleanedsigma = 5.0 - intensity_threshold: -| mandatory, min data threshold. Intensity values below this are set to 0. The threshold is applied after removing aliens. +| Mandatory, data threshold. Intensity values below this value are set to 0. The threshold is applied after removing aliens. +| If auto_data is configured in main config file, this value will be overridden by calculated value. | example: :: intensity_threshold = 25.0 - adjust_dimensions: -| optional, a list of number to adjust the size at each side of 3D data. If number is positive, the array will be padded. If negative, cropped. The parameters correspond to [x left, x right, y left, y right, z left, z right] The final dimensions will be adjusted up to the good number for the FFT which also is compatible with opencl supported dimensions powers of 2 or a*2^n, where a is 3, 5, or 9 +| Optional, a list of numbers defining how to adjust the size at each side of 3D data. If number is positive, the array will be padded. If negative, cropped. The parameters correspond to [x left, x right, y left, y right, z left, z right]. The final dimensions will be adjusted up to the good number for the FFT such as: powers of 2 or a*2^n, where a is 3, 5, or 9. | example: :: adjust_dimensions = [13, 0, -65, -65, -65, -65] - center_shift: -| optional, enter center shift list the array maximum is centered before binning, and moved according to center_shift, [0,0,0] has no effect +| Optional, defines offset of max element from the array center. | example: :: center_shift = [0,0,0] - binning: -| optional, a list that defines binning values in respective dimensions, [1,1,1] has no effect +| Optional, a list that defines binning values in respective dimensions, [1,1,1] has no effect. +| If auto_data is configured in main config file, this list will be overridden by calculated values. | example: :: diff --git a/docs/source/config_disp.rst b/docs/source/config_disp.rst index 32e83ee..29be0bb 100755 --- a/docs/source/config_disp.rst +++ b/docs/source/config_disp.rst @@ -6,22 +6,36 @@ config_disp Parameters ========== - results_dir: -| optional, defaults to . A directory that has a tree, or leaf with reconstruction results. The results will be used as input to the visualization processing. There could be several results in the given directory tree, and all will be processed concurrently. +| Optional, defaults to . A directory that has a tree, or leaf with reconstruction results. The results will be used as input to the visualization processing. There could be several results in the given directory tree, and all will be processed concurrently. | example: :: results_dir = "/path/to/phasing_results_dir_tree" - rampups: -| optional, upsize when running ramp removal, default is 1. Expect long processing time for greater numbers. +| Optional, upsize when running ramp removal, default is 1. Expect long processing time for greater numbers. | example: :: rampups = 2 - crop: -| optional, defines a fraction of array at the center that will be subject of visualization. Defaults to [1.0, 1.0, 1.0], i.e. the size of the processed array. +| Optional, defines a fraction of array at the center that will be subject of visualization. Defaults to [1.0, 1.0, 1.0], i.e. the size of the processed array. | example: :: crop = [.5, .5, .5] + +- unwrap: +| Switch, if True the image.vts file will contain unwrapped phase in addition to phase and amplitude. +| example: +:: + + unwrap = False + +- make_twin: +| A switch to whether visualize twin image. +| example: +:: + + make_twin = False diff --git a/docs/source/config_instr.rst b/docs/source/config_instr.rst index 9bb9b32..cbaa605 100644 --- a/docs/source/config_instr.rst +++ b/docs/source/config_instr.rst @@ -2,19 +2,20 @@ config_instr ============ | The "config_instr" file defines parameters related to experiment and instrument used during the experiment, such as detector and diffractometer. -| The diffractometer and specfile are mandatory parameters for the 34-IDC beamline. The other parameters are obtained from the specfile. They can be overriden when included in this configuration file. +| The diffractometer and specfile are mandatory parameters for the 34-IDC beamline. The other parameters are obtained from the specfile. They can be overridden when included in this configuration file. Parameters ========== - diffractometer: -| mandatory, name of diffractometer used in experiment. +| Mandatory, name of diffractometer used in experiment. +| Must be defined for the beamline. | example: :: diffractometer = "34idc" - specfile: -| optional (but mandatory for aps_34idc), specfile recorded when the experiment was conducted. +| Optional (but mandatory for aps_34idc), specfile recorded when the experiment was conducted. | example: :: @@ -22,7 +23,7 @@ Parameters Parsed parameters ================= -| In a typical scenario at APS 34-idc beamline a spec file is generated during experiment and the parameters listed below are parsed from this file, so the user do not configure them. User may configure the following parameters if spec file is not configured in the main configuration or user wants to override the parsed parameters. +| In a typical scenario at APS 34-idc beamline a spec file is generated during experiment and the parameters listed below are parsed from this file, so the user do not configure them. User may override the parameters. - energy | example: :: diff --git a/docs/source/config_main.rst b/docs/source/config_main.rst index f8789da..323f9df 100755 --- a/docs/source/config_main.rst +++ b/docs/source/config_main.rst @@ -1,61 +1,70 @@ ===== cofig ===== -| The "config" is a main configuration defining the experiment and is applicable to all parts of processing: beamline preprocessing, standard preprocessing, phasing, and beamline visualization. +| The "config" file is a main configuration defining the experiment and is applicable to all parts of processing: beamline preprocessing, standard preprocessing, phasing, and beamline visualization. Parameters ========== - working_dir: -| mandatory, root directory of the experiment files +| Mandatory, root directory of the experiment files. | example: :: working_dir = "/path/to/the/experiment/" - experiment_id: -| mandatory, a string identifying the experiment +| Mandatory, a string identifying the experiment. | example: :: experiment_id = "ab" - scan: -| optional (but typically needed), string type encapsulating a single number, a range, or combination of numbers and ranges separated by comma, defining scans that will be read and combined to create data file. This parameter should not be set when running multi peak case. -| example: +| Optional (but typically needed), string type encapsulating a single number, a range, or combination of numbers and ranges separated by comma, defining scans that will be read and combined to create data file. This parameter should not be set when running multi-peak case. +| examples: :: + scan = "54" + scan = "2-7" scan = "2-7, 10, 15, 20-22" - beamline: -| optional, beamline where the experiment was conducted. If not configured, the beamline preprocessing and beamline visualization scripts are not available. +| Optional, beamline where the experiment was conducted. If not configured, the beamline preprocessing and beamline visualization scripts are not available. | example: :: beamline = "aps_34idc" - multipeak: -| optional, boolean parameter indicating whether it is multi peak case. Defaults to False. +| Optional, boolean parameter indicating whether it is multi-peak case. Defaults to False. | example: :: multipeak = True - separate_scans: -| in typical scenario the data from all scans in experiment are combined. If specified as separate scans, each scan will be processed separately and will have sub-experiment name containing scan index ex. "scan_9", where 9 is scan index +| Optional, deafults to False. In typical scenario the data from all scans in experiment are combined. If specified as separate scans, each scan will be processed separately and will have sub-experiment name containing scan index ex. "scan_9", where 9 is scan index. | example: :: separate_scans = False - separate_scan_ranges: -| in typical scenario the data from all scans in experiment are combined. If specified as separate scan ranges, each scan or scan range in the experiment will be processed separately and will have sub-experiment name containing scan index ex. "scan_9", where 9 is scan index, or "scan_10-15", where 10-15 is the scan range. The scans and scan ranges are defined in main configuration "config" file as scan parameter, and are part of experiment name. +| Optional, defaults to False. In typical scenario the data from all scans in experiment are combined. If specified as separate scan ranges, each scan or scan range in the experiment will be processed separately and will have sub-experiment name containing scan index, or scan index range, ex. "scan_9", where 9 is scan index, or "scan_10-15", where 10-15 is the scan range. The scans and scan ranges are defined in main configuration "config" file as scan parameter, and are part of experiment name. | example: :: separate_scan_ranges = True +- auto-data: +| Optional, boolean parameter indicating automated data preprocessing, which includes exclusion of outlier scans in multi-scan data and auto calculation of intensity threshold. Defaults to False. +| example: +:: + + auto_data = True + - converter_ver: -| optional, if not configured, it will be auto-updated by scripts to match the latest converter version. +| Mandatory after ver 3.0. If not configured, it will be auto-updated by scripts to match the latest converter version. | example: :: diff --git a/docs/source/config_mp.rst b/docs/source/config_mp.rst index a1aa72c..257d065 100644 --- a/docs/source/config_mp.rst +++ b/docs/source/config_mp.rst @@ -1,81 +1,78 @@ -========= +======== config_mp -========= -| Multi-peak reconstruction is marked by a parameter multipeak set to True in the config file. -| For the multi-peak reconstruction the conf directory must have config_mp file with the parameters listed below.. - +======== - scan: -| mandatory, string type encapsulating scans or ranges of scans containing data for each peak. The scans/scan ranges should be arranged in ascending order. +| Mandatory, string type encapsulating scans or ranges of scans containing data for each peak. The scans/scan ranges should be arranged in ascending order. | example: :: scan = "898-913,919-934,940-955,961-976" - orientations: -| mandatory, a list of lists, each inner list defining the orientation of a peak. +| Mandatory, a list of lists, each inner list defining the orientation of a peak. | example: :: orientations = [[-1, -1, 1], [0, 0, 2], [1, 1, 1], [2, 0, 0]] - hkl_in: -| mandatory, list with Miller indices representing the in-plane lattice vector. +| Mandatory, list with Miller indices representing the in-plane lattice vector. | example: :: hkl_in = [3.031127677370978, 12.31353345906843, 8.75104158816168] - hkl_out: -| mandatory, list with Miller indices representing the out-of-plane lattice vector. +| Mandatory, list with Miller indices representing the out-of-plane lattice vector. | example: :: hkl_out = [9.77805918769193, -8.402719849515048, 8.436553021703112] - twin_plane: -| mandatory, list with Miller indices of the twin plane. If there is not a twin plane, this should be the same as sample_axis +| Mandatory, list with Miller indices of the twin plane. If there is not a twin plane, this should be the same as sample_axis. | example: :: twin_plane = [1, -1, 1] - sample_axis: -| mandatory, axis of the sample. The data will be rotated so that the twin plane vector in the lattice frame (hkl) corresponds to this vector in the array frame (xyz). +| Mandatory, axis of the sample. The data will be rotated so that the twin plane vector in the lattice frame (hkl) corresponds to this vector in the array frame (xyz). | example: :: sample_axis = [0, 1, 0] - final_size: -| mandatory, a size in each dimension of the array holding reconstructed object. +| Mandatory, a size in each dimension of the array holding reconstructed object. | example: :: final_size = 180 - mp_max_weight: -| optional, a number between 0 and 1.0 specifying the initial coupling weight assigned to all peaks. if not provided, defaults to 1.0. +| Optional, a number between 0 and 1.0 specifying the initial coupling weight assigned to all peaks. if not provided, defaults to 1.0. | example: :: mp_max_weight = 1.0 - mp_taper: -| mandatory, fraction of the total iterations at which to begin tapering the coupling weight. Coupling weight will decrease for all peaks starting at mp_taper times the total number of iterations. +| Mandatory, fraction of the total iterations at which to begin tapering the coupling weight. Coupling weight will decrease for all peaks starting at mp_taper times the total number of iterations. | example: :: mp_taper = 0.6 - lattice_size: -| mandatory, lattice parameter of the reconstructing crystal. This is used to define the reciprocal lattice vectors, which are required for projecting to each peak. +| Mandatory, lattice parameter of the reconstructing crystal. This is used to define the reciprocal lattice vectors, which are required for projecting to each peak. | example: :: lattice_size = 0.4078 - switch_peak_trigger: -| mandatory, a trigger defining at which iteration to switch the peak +| Mandatory, a trigger defining at which iteration to switch the peak. | example: :: diff --git a/docs/source/config_prep.rst b/docs/source/config_prep.rst index 1188cc7..85f3d7a 100755 --- a/docs/source/config_prep.rst +++ b/docs/source/config_prep.rst @@ -7,42 +7,42 @@ Parameters ========== - data_dir: -| mandatory for the 34-IDC prep, the directory containing raw data. This directory will have several subdirectories each containing tif format data files. Each subdirectory represents separate scan and is numbered with the scan index +| Mandatory for the 34-IDC prep, the directory containing raw data. This directory will have several subdirectories each containing tif format data files. Each subdirectory represents separate scan and is numbered with the scan index. | example: :: data_dir = "/path/to/data/files/ADStaff19-1a" - darkfield_filename: -| optional, depending on the detector. Dark field file taken to clean bad pixels or correct background in CCD. Each detector needs to implement this correctly. +| Optional, depending on the detector. Dark field file taken to clean bad pixels or correct background in CCD. Each detector needs to implement this correctly. | example: :: darkfield_filename = "/path/to/darkfield_file/dark.tif" - whitefield_filename: -| optional, depending on the detector. White field file taken to correct over saturated pixels. Each detector needs to implement this correctly. +| Optional, depending on the detector. White field file taken to correct over saturated pixels. Each detector needs to implement this correctly. | example: :: whitefield_filename = "/path/to/whitefield_file/whitefield.tif" - roi: -| typically read from a spec file, but if not present there, must be configured. Defines the region of a detector that the data has been read from. Needed to determine corrections on the detector. The values in the list (refer to example below) are as follows: x coordinate starting point, x length, y coordinate starting point, y length. +| Typically read from a spec file. Defines the region of a detector that the data has been read from. The values in the list (refer to example below) are as follows: x coordinate starting point, x length, y coordinate starting point, y length. | example: :: roi = [0,256,0,256] - min_files: -| optional, defines a minimum number of raw data files in the scan directory. If number of files is less than minimum, the directory is not added to the data. +| Optional, defines a minimum number of raw data files in the scan directory. If number of files is less than minimum, the directory is not added to the data. | example: :: min_files = 80 - exclude_scans: -| a list containing scan indexes that will be excluded from preparation +| A list containing scan indexes that will be excluded from preparation. | example: :: @@ -55,10 +55,9 @@ Parameters Imult = 1000000 -- detector: -| optional, can be omitted if roi param is specified +- outliers_scans: +| This list is determined when running auto-data preprocessing. | example: :: - detector = "34idcTIM2" - + outliers_scans = [78,80] diff --git a/docs/source/config_rec.rst b/docs/source/config_rec.rst index e362f52..54a138c 100755 --- a/docs/source/config_rec.rst +++ b/docs/source/config_rec.rst @@ -3,87 +3,90 @@ ========== config_rec ========== -| The "config_rec" defines parameters used during reconstruction process. -| Parameters are grouped by function (feature) they are utilized for. There is a group of parameters that are not associated with any feature, and thus, we call them 'general'. -| The features are controlled by triggers. If trigger is defined then the feature is active and related parameters are used when processing the feature. Refer to :ref:`formula` page, for explanation of triggers, sub-triggers, and algorithm sequence. +| The "config_rec" file defines parameters used during reconstruction process. +| +| Parameters are grouped by trigger and corresponding parameters. There is a group of parameters that are not associated with any trigger, and thus, we call them 'general'. +| If trigger is defined then the feature is active and related parameters are used when processing the feature. Refer to :ref:`formula` page, for explanation of triggers, sub-triggers, and algorithm sequence. Parameters ========== General +++++++ - data_dir: -| optional, used for specific cases, default is /phasing_data. Driven by the scripts, the experiment directory contains four directories: conf, prep, data, and results. The standard preprocessed file, named data.tif should be in data_dir when starting reconstruction. If the data_dir parameter is configured then the data is read from this directory, otherwise from /data. This parameter is only accepted when running script from a command line. When using GUI, the data_dir is always default. +| Optional, defines a directory with file named "data.tif" that contains experiment data file ready for reconstruction. Default is /phasing_data. | example: :: data_dir = "/path/to/data_dir/data_dir" - save_dir: -| optional, used for specific cases, defualt is /results. Driven by the scripts, the experiment directory contains four directories: conf, preprocessed_data, phasing_data, results_phasing, and results_viz. The reconstruction results are saved in /results_phasing directory. If the save_dir parameter is configured then the reconstruction result is saved in this directory. +| Optional, used for specific cases, defualt is /results. Driven by the scripts, the experiment directory contains four directories: conf, preprocessed_data, phasing_data, results_phasing, and results_viz. The reconstruction results are saved in /results_phasing directory. If the save_dir parameter is configured then the reconstruction result is saved in this directory. | example: :: save_dir = "/path/to/save_dir/save_dir" - init_guess: -| optional, defines how the initial guess is set. Possible options are: 'random', 'continue', and 'AI_guess'. The choice "random" will generate random guess, "continue" will start from previously saved results, and "AI_guess" will run AI reconstruction that will be an initial guess saved at /results_AI. Each of these algorithms require different parameters, explained below. The default is 'random'. +| Optional, defines how the initial guess is set. Possible options are: 'random', 'continue', and 'AI_guess'. The choice "random" will generate random guess, "continue" will start from previously saved results, and "AI_guess" will run AI reconstruction that will be an initial guess saved at /results_AI. Each of these algorithms require different parameters, explained below. The default is 'random'. | example: :: init_guess = "random" - continue_dir: -| must be defined if the init_guess parameter is set to 'continue'. Directory from which initial guess, initial support are read for reconstruction continuation. If the directory contains multiple subdirectories with initial guesses, a thread will start for each subdirectory. +| Must be defined if the init_guess parameter is set to 'continue'. Directory from which initial guess, initial support are read for reconstruction continuation. If the directory contains multiple subdirectories with initial guesses, a thread will start for each subdirectory. | example: :: continue_dir = "/path/to/previous_results_dir/some_phasing_results_dir" - AI_trained_model: -| must be defined, if init_guess is "AI_guess". defines the file in hdf5 format that holds trained model. +| Must be defined, if init_guess is "AI_guess". Defines the file in hdf5 format that holds trained model. | example: :: AI_trained_model = "/path/to/trained/model/trained_model.hdf5" - reconstructions: -| optional, default is 1. Number of reconstructions to start with. Typically used when running genetic algorithm. +| Optional, default is 1. Number of reconstructions to start with. Typically used when running genetic algorithm. | example: :: reconstructions = 5 - processing: -| optional, the library used when running reconstruction. When the auto option is selected the program will use the best performing library that is available, in the following order: cupy, numpy. The cp option will utilize cupy, and np will utilize numpy, and af will leave selection to arrayfire. Default is auto. +| Optional, the library used when running reconstruction. When the auto option is selected the program will use the best performing library that is available, in the following order: cupy, torch, numpy. The cp option will utilize cupy, torch will utilize torch, and np will utilize numpy. Default is auto. | example: :: processing = "auto" - device: -| optional, IDs of the target devices (GPU) for each reconstruction thread. If not defined, the OS will select the GPU, but the processing will not be concurrent. Ignored when running cpu library. -| example: +| Optional, GPU IDs of the target devices for reconstruction(s) or 'all' if all available GPUs should be used. If not defined, the reconstruction process will run on CPU. For cluster configuration it is defined as dict with hosts names as keys and values as described. +| examples: :: device = [0,1,2,7] + device = 'all' + device = {'host1':'all', 'host2':[0,1,2,3,4]} - algorithm_sequence: -| mandatory, defines sequence of algorithms applied in each iteration during modulus projection and during modulus. The "*" charcter means repeat, and the "+" means add to the sequence. The sequence may contain single brackets defining a group that will be repeated by the preceding multiplier. The alphabetic entries: ER, ERpc, HIO, HIOpc define algorithms used in this iteration. The entries will invoke functions as follows: ER definition will invoke 'er' and 'modulus' functions, the ERpc will invoke 'er' and 'pc_modulus', HIO will invoke 'hio' and 'modulus', and HIOpc will invoke 'hio' and 'pc_modulus'. The pc_modulus is implementation of modulus with partial coherence correction. In second example the sequence contains subtriggers, explained in :ref:`formula` page. +| Mandatory, defines sequence of algorithms applied in each iteration during modulus projection and during modulus. The "*" character means repeat, and the "+" means add to the sequence. The sequence may contain single brackets defining a group that will be repeated by the preceding multiplier. The alphabetic entries: ER, ERpc, HIO, HIOpc define algorithms used in this iteration. The entries will invoke functions as follows: ER definition will invoke 'er' and 'modulus' functions, the ERpc will invoke 'er' and 'pc_modulus', HIO will invoke 'hio' and 'modulus', and HIOpc will invoke 'hio' and 'pc_modulus'. The pc_modulus is implementation of modulus with partial coherence correction. In second example the sequence contains sub-triggers, explained in :ref:`formula` page. | examples: :: algorithm_sequence = "2* (20*ER + 180*HIO) + 2* (20*ERpc + 180*HIOpc) + 20*ERpc" - algorithm_sequence = "20*ER.LPF0.PHM0 + 180*HIO.LPF1 + 2* (20*ER.SW0 + 180*HIO.SW1) + 20*ER.SW2" + algorithm_sequence = "20*ER.LPF0.PHC0 + 180*HIO.LPF1 + 2* (20*ER.SW0 + 180*HIO.SW1) + 20*ER.SW2" - hio_beta: -| optional, default is .9. A parameter used in hio algorithm. +| Optional, default is .9. A parameter used in hio algorithm. | example: :: hio_beta = .9 - initial_support_area: -| optional, defaults to [.5,.5,.5]. The list define dimensions of initial support area. If the values are fractional, the support area will be calculated by multiplying by the data array dimensions. The support array is centered. +| Optional, defaults to [.5,.5,.5]. The list define dimensions of initial support area. The support area is calculated by multiplying the fractions in this parameter by the corresponding data array dimensions. The support array is centered. | example: :: @@ -93,14 +96,14 @@ Twin ++++ - twin_trigger: -| optional, defines at which iteration to eliminate "twin", i.e. the image under reconstruction is trimmed by preserving a quadrant across x and y dimensions and zeroing the rest of the array. +| Defines at which iteration to eliminate "twin", i.e. the image under reconstruction is trimmed by preserving a quadrant across x and y dimensions and zeroing the rest of the array. | example: :: twin_trigger = [2] - twin_halves = [0, 0] -| optional, and only applied when twin_trigger is configured. Defines which quadrant of the array is preserved in x and y dimensions, defaults to (0, 0). +| Optional, defines which quadrant of the array is preserved in x and y dimensions, defaults to (0, 0). | Possible choices: [0, 0], [0, 1], [1, 0], [1,1] | example: :: @@ -109,11 +112,11 @@ Twin Shrink wrap +++++++++++ -| Support area is an array that defines region in which the image is meaningful. This area is recalculated at the trigger iteration shrinking along when the image develops. The calculations employ an algorithm defined here as shrink_wrap_type. +| Support area is an array that defines region in which the image is meaningful. This area is recalculated at the shrink wrap trigger iteration, shrinking along when the image develops. The calculations employ an algorithm defined here as shrink_wrap_type. - shrink_wrap_trigger: -| defines when to update support array using the parameters below. -| alternatively can be defined as list of sub-triggers. If sub-triggers are used, the parameters must be lists as well. +| Defines when to update support array using the parameters below. +| Alternatively can be defined as list of sub-triggers. If sub-triggers are used, the parameters must be lists as well. | examples: :: @@ -121,8 +124,7 @@ Shrink wrap shrink_wrap_trigger = [[10, 1],[0,5,100],[0,2]] # sub-triggers - shrink_wrap_type: -| optional, defaults to "GAUSS" which applies gaussian filter. Currently only "GAUSS" is supported. -| Currently only the GAUSS type is supported. +| Mandatory, defines type of shrink wrap. Currently only the "GAUSS" type is supported that applies gaussian filter to find support area. | examples: :: @@ -130,7 +132,7 @@ Shrink wrap shrink_wrap_type = [GAUSS, GAUSS, GAUSS] # sub-triggers - shrink_wrap_threshold: -| optional, defaults to 0.1. A threshold value used in the gaussian filter algorithm. +| Mandatory, defines a threshold value used in the gaussian filter algorithm. | examples: :: @@ -138,7 +140,7 @@ Shrink wrap shrink_wrap_threshold = [0.1, 0.11, .12] # sub-triggers - shrink_wrap_gauss_sigma: -| optional, defaults to 1.0. A sigma value used in the gaussian filter algorithm. +| Mandatory, defines a sigma value used in the gaussian filter algorithm. | examples: :: @@ -147,67 +149,67 @@ Shrink wrap Phase constrain +++++++++++++++ -| At the beginning iterations the support area is modified in respect to the phase. Support area will include only points with calculated phase intside of the defined bounds. +| At the beginning iterations the support area is modified in respect to the phase. Support area will include only points with calculated phase inside of the defined bounds. | Alternatively can be defined as list of sub-triggers. If sub-triggers are used, the parameters must be lists as well. -- phase_support_trigger: -| defines when to update support array using the parameters below by applying phase constrain. +- phc_trigger: +| Defines when to update support array using the parameters below by applying phase constrain. | examples: :: - phase_support_trigger = [0, 1, 310] - phase_support_trigger = [[0, 1, 310], [0,2]] # sub-triggers + phc_trigger = [0, 1, 310] + phc_trigger = [[0, 1, 310], [0,2]] # sub-triggers -- phm_phase_min: -| optional, defaults too -1.57. Defines lower bound phase. +- phc_phase_min: +| Mandatory, defines lower bound phase. | example: :: - phm_phase_min = -1.57 - phm_phase_min = [-1.5, -1.57] # sub-triggers + phc_phase_min = -1.57 + phc_phase_min = [-1.5, -1.57] # sub-triggers -- phm_phase_max: -| optional, defaults too 1.57. Defines upper bound phase. +- phc_phase_max: +| Mandatory, defines upper bound phase. | examples: :: - phm_phase_max = 1.57 - phm_phase_max = [1.5, 1.57] # sub-triggers + phc_phase_max = 1.57 + phc_phase_max = [1.5, 1.57] # sub-triggers Partial coherence +++++++++++++++++ -| Partial coherence triggers recalculation of coherence array for the amplitudes in reciprocal space. After coherence array is determined, it is used for convolution in subsequent iteration. The coherence array is updated as defined by the pc_interval. Partial coherence feature is active if the interval is defined and the algorithm sequence contains algorithm with partial coherence. +| Partial coherence operation initiates recalculation of coherence of data with respect to the amplitudes in reciprocal space. After coherence array is determined, it is used in convolution operation during modulus in subsequent iteration. The coherence array is updated periodically, as defined by the pc_interval. Partial coherence operation is active if the interval is defined and the algorithm sequence indicates "pc". - pc_interval: -| defines iteration interval between coherence update. +| Defines iteration interval between coherence update. | example: :: pc_interval = 50 - pc_type: -| mandatory, partial coherence algorithm. Currently "LUCY" is supported. +| Partial coherence algorithm. Currently "LUCY" is supported. | example: :: pc_type = "LUCY" - pc_LUCY_iterations: -| optional, defaults to 20. a number of iteration inside LUCY algorithm. +| Optional, defaults to 20. Defines number of iteration inside LUCY algorithm. | example: :: pc_LUCY_iterations = 20 - pc_normalize: -| optional, defaults to True. Internal. +| Optional, defaults to True. Internal. | example: :: pc_normalize = True - pc_LUCY_kernel: -| mandatory, coherence array area. +| Mandatory, coherence array area. | example: :: @@ -215,10 +217,10 @@ Partial coherence Lowpass Filter ++++++++++++++ -| When this feature is activated a lowpass Gaussian filter is applied on data with sigma being line-spaced over the range parameter. Simultaneously, the Gaussian type of shrink wrap is applied with the reverse sigma. The low resolution trigger is typically configured to be active at the first part of iterations and resume around the mid-run. +| When active, a lowpass Gaussian filter is applied on data, with iteration dependent sigma calculated by line-spacing the lowpass_filter_range parameter over trigger span iterations. Simultaneously, the Gaussian type of shrink wrap is applied with the reverse sigma. The low resolution trigger is typically configured to be active at the first part of iterations. - lowpass_filter_trigger: -| defines when to apply low resolution using the parameters below. Typically the last trigger is configured at half of total iterations. -| Alternatively can be defined as list of sub-triggers. If sub-triggers are used, the parameters must be lists as well. +| Defines when to apply lowpass filter operation using the parameters below. Typically the last trigger is configured at half of total iterations. +| Alternatively, it can be defined as list of sub-triggers. If sub-triggers are used, the parameters must be lists as well. | examples: :: @@ -226,7 +228,7 @@ Lowpass Filter lowpass_filter_trigger = [[0, 1], [0, 2, 100]] # sub-triggers - lowpass_filter_range: -| The range is linespaced over iterations to form a list of sigmas. The sigmas are used to apply gaussian lowpass filter over data. The inverse of sigma is used in shrink wrap.ue to last. If only one number is given, the last det will default to 1. +| The range is line-spaced over trigger iterations to form a list of iteration dependent sigmas. If only one number is given, the last sigma will default to 1. | examples: :: @@ -234,7 +236,7 @@ Lowpass Filter lowpass_filter_range = [[.7, .8], [.8, 1.0]] # sub-triggers - lowpass_filter_sw_threshold: -| during lowpass iterations a GAUSS type shrink wrap is applied with this threshold ans sigma calculated as reverse of low pass filter. +| During lowpass iterations a GAUSS type shrink wrap is applied with this threshold ans sigma calculated as reverse of low pass filter. | examples: :: @@ -245,7 +247,7 @@ averaging +++++++++ | When this feature is activated the amplitudes of the last several iterations are averaged. - average_trigger: -| defines when to apply averaging. Negative start means it is offset from the last iteration +| Defines when to apply averaging. Negative start means it is offset from the last iteration. | example: :: @@ -254,7 +256,7 @@ averaging progress ++++++++ - progress_trigger: -| defines when to print info on the console. The info includes current iteration and error +| Defines when to print info on the console. The info includes current iteration and error. | example: :: @@ -263,16 +265,17 @@ progress GA ++ - ga_generations: -| optional, number of generations. When defined, and the number is greater than 1, the genetic algorithm (GA) is activated +| Defines number of generations. When defined, and the number is greater than 1, the genetic algorithm (GA) is activated | example: :: ga_generations = 3 - ga_metrics: -| optional, a list of metrics that should be used to rank the reconstruction results for subsequent generations. If not defined, or shorter than number of generations, the metric defaults to "chi". -| supported: -| - 'chi': the last error calculated as norm(rs_amplitudes - data)/norm(data) +| Optional, a list of metrics that should be used to rank the reconstruction results for subsequent generations. If not defined, or shorter than number of generations, the metric defaults to "chi". +| If the list contains only one element, it will be used by all generations. +| Supported metrics: +| - 'chi': The last error calculated as norm(rs_amplitudes - data)/norm(data). | The smallest 'chi' value is the best. | - 'sharpness': sum(power(abs(image), 4)) | The smallest 'sharpness' value is the best. @@ -282,15 +285,17 @@ GA | - 'area': sum(support) | where support is calculated with shrink wrap using hardcoded threshold=.2 and sigma=.5 | The greatest 'area' value is the best. -| example: +| examples: :: ga_metrics = ["chi", "sharpness", "area"] + ga_metrics = ["chi"] - ga_breed_modes: -| optional, a list of breeding modes applied to breed consecutive generation. If not defined, or shorter that number of generations, the mode defaults to "sqrt_ab". +| Optional, a list of breeding modes applied to breed consecutive generation. If not defined, or shorter that number of generations, the mode defaults to "sqrt_ab". +| If the list contains only one element, it will be used by all generations. | Breeding starts with choosing alpha image. The rest of the images are crossed with alpha. Before the crossing, the image, called beta is aligned with alpha, and phases in both of the arrays are normalized to derive ph_alpha = angle(alpha), and ph_beta = angle(beta) -| supported: +| Supported modes: | - 'sqrt_ab': sqrt(abs(alpha) * abs(beta)) * exp(0.5j * (ph_beta + ph_alpha)) | - 'pixel_switch': where((cond > 0.5), beta, alpha); cond = random(shape(beta)) | - 'b_pa': abs(beta) * exp(1j * (ph_alpha)) @@ -305,41 +310,44 @@ GA | - 'max_ab_pa': max(abs(alpha), abs(beta)) * exp(1j * ph_alpha) | - 'avg_ab': 0.5 * (alpha + beta) | - 'avg_ab_pa: 0.5 * (abs(alpha) + abs(beta)) * exp(1j * (ph_alpha)) -| example: +| examples: :: ga_breed_modes = ["sqrt_ab", "pixel_switch", "none"] + ga_breed_modes = ["sqrt_ab"] - ga_cullings: -| optional, defines how many worst samples to remove in a breeding phase for each generation. If not defined for the generation, the culling defaults to 0. +| Optional, defines how many worst samples to remove in a breeding phase for each generation. If not defined for the generation, the culling defaults to 0. | example: :: ga_cullings = [2,1] - ga_shrink_wrap_thresholds: -| optional, a list of threshold values for each generation. The support is recalculated with this threshold after breeding phase. Defaults to configured value of support_threshold. +| Optional, a list of threshold values for each generation. The support is recalculated with this threshold after breeding phase. Defaults to configured value of support_threshold. +| If the list contains only one element, it will be used by all generations. | example: :: ga_shrink_wrap_thresholds = [.15, .1] - ga_shrink_wrap_gauss_sigmas: -| optional, a list of sigma values for each generation. The support is recalculated with this sigma after breeding phase. Defaults to configured value of support_sigma. +| Optional, a list of sigma values for each generation. The support is recalculated with this sigma after breeding phase. Defaults to configured value of support_sigma. +| If the list contains only one element, it will be used by all generations. | example: :: ga_shrink_wrap_gauss_sigmas = [1.1, 1.0] - ga_lowpass_filter_sigmas: -| optional, a list of sigmas that will be used in subsequent generations to calculate Gauss distribution in the space defined by the size of the data and applied it to the data. In the example given below this feature will be used in first two generations. +| Optional, a list of sigmas that will be used in subsequent generations to calculate Gaussian low-pass filter applied it to the data. In the example given below this feature will be used in first two generations. | example: :: ga_lowpass_filter_sigmas = [2.0, 1.5] - ga_gen_pc_start: -| optional, a number indicating at which generation the pcdi feature will start to be active. If not defined, and the pcdi feature is active, it will start at the first generation. +| Optional, a number indicating at which generation the partial coherence will start to be active. If not defined, and the pc feature is active, it will start at the first generation. | example: :: diff --git a/docs/source/define_alg_seq.rst b/docs/source/define_alg_seq.rst index 1398d2b..8699dca 100644 --- a/docs/source/define_alg_seq.rst +++ b/docs/source/define_alg_seq.rst @@ -1,33 +1,35 @@ .. _formula: ================= -execution formula +Execution formula ================= -| The "config_rec" defines parameters used during reconstruction process. -| Parameters are grouped by features. -| The features are controlled by triggers. If trigger is defined then the feature is active and related parameters are used when processing the feature. -| Triggers and algorithm_sequence parameter define which features are used in each iteration. Trigger can be defined as a single trigger or a list of sub-triggers. -| Algorithm_sequence define the algorithms used in iterations. The algorithm_sequence and triggers create a formula. +| The "config_rec" file defines parameters used during reconstruction process. +| The execution flow is defined by algorithm_sequence parameter and triggers. +| Algorithm_sequence defines the algorithms used in iterations. Algorithm is defined for each iteration. +| Triggers define iterations at which the corresponding triggered operations are active. +| The triggered operations use applicable parameters that are grouped by the corresponding trigger. +| If the trigger is not defined then the corresponding operation is not active during the execution. +| Trigger can be defined as a single trigger or a list of sub-triggers. | The trigger and algorithm sequence are explained in the sections below. Trigger ======= -| Trigger concludes at which iteration to apply the associated feature action. Trigger is defined as a list and can be configured for a single iteration, or multiple iterations. +| Trigger concludes at which iteration to apply the associated action. Trigger is defined as a list and can be configured for a single iteration, or multiple iterations. | examples: | [3] trigger at iteration 3 -| [20, 5] trigger starts at iteration 20, repeats every 5 iteration for the rest of run -| [20, 5, 40] trigger starts at iteration 20, repeats every 5 iteration until iteration 40 +| [20, 5] trigger starts at iteration 20, repeats every 5 iteration for the rest of run. +| [20, 5, 40] trigger starts at iteration 20, repeats every 5 iteration until iteration 40. Sub-Trigger =========== -| For the following features: shrink wrap, lowpass filter, and phm (phase modulus), the trigger can be defined as a list of sub-triggers. In this case the algorithm_sequence should include the sub-triggers in the sequence, for example: 30*ER.SW0, and the trigger should be defined as a list of sub-triggers. -| examples: +| For the following features: shrink wrap, lowpass filter, and phc (phase constrain), the trigger can be defined as a list of sub-triggers. In this case the algorithm_sequence should include the sub-triggers in the sequence, for example: 30*ER.SW0, and the trigger should be defined as a list of sub-triggers. +| example: | [[5, 1, 100], [0, 3], [5]] three sub-triggers are defined and the rules of the trigger apply for sub-trigger with the difference that the span of iterations in sub-trigger is the same as the algorithm it is attached to. In this example 30*ER.SW0, the start iteration is when the ER begins, and last iteration is when the ER ends. - -- algorithm_sequence: -| defines sequence of algorithms applied in each iteration during modulus projection and during modulus. The "*" character means repeat, and the "+" means add to the sequence. The sequence may contain single brackets defining a group that will be repeated by the preceding multiplier. The alphabetic entries: ER, ERpc, HIO, HIOpc define algorithms used in this iteration. The entries will invoke functions as follows: ER definition will invoke 'er' and 'modulus' functions, the ERpc will invoke 'er' and 'pc_modulus', HIO will invoke 'hio' and 'modulus', and HIOpc will invoke 'hio' and 'pc_modulus'. The pc_modulus is implementation of modulus with partial coherence correction. If defining ERpc or HIOpc the pcdi feature must be activated by configuring pc_interval. -| Algorithm sequence can include sub-triggers attached to the algorithm, for example 30*ER.SW0.PHM1. +Algorithm Sequence +================== +| This parameter defines sequence of algorithms applied in each iteration during modulus projection and during modulus. The "*" character means repeat, and the "+" means add to the sequence. The sequence may contain single brackets defining a group that will be repeated by the preceding multiplier. The alphabetic entries: ER, ERpc, HIO, HIOpc define algorithms used in this iteration. The entries will invoke functions as follows: ER definition will invoke 'er' and 'modulus' functions, the ERpc will invoke 'er' and 'pc_modulus', HIO will invoke 'hio' and 'modulus', and HIOpc will invoke 'hio' and 'pc_modulus'. The pc_modulus is implementation of modulus with partial coherence correction. If defining ERpc or HIOpc the pcdi (partial coherence)) must be activated by configuring pc_interval and pcdi parameters. +| Algorithm sequence can include sub-triggers attached to the algorithm, for example 30*ER.SW0.PHC1. Formula examples ================ @@ -38,24 +40,26 @@ Formula examples algorithm_sequence = "2* (20*ER + 180*HIO) + 2* (20*ERpc + 180*HIOpc) + 20*ERpc" config_rec: -| shrink_wrap_trigger = [1, 3] -| lowpass_trigger = [0, 2, 130] +| shrink_wrap_trigger = [1, 3] +| lowpass_trigger = [0, 2, 130] -| In this example the program will run twice 20 iterations with ER algorithm and modulus and 180 iterations with HIO algorithm and modulus -| followed by twice 20 iterations with ER algorithm and partial coherence modulus and 180 iterations with HIO algorithm and partial coherence modulus -| followed by 20 iterations with ER algorithm and partial coherence modulus -| In addition, based on the config_rec, every three iterations starting from first, for the rest of the run, the shrink wrap trigger will be applied -| and every two iterations, starting at the beginning, until iteration 130, the lowpass trigger will be applied. +| In this example the program will run: +| twice 20 iterations with ER algorithm and modulus and 180 iterations with HIO algorithm and modulus, +| followed by twice 20 iterations with ER algorithm and partial coherence modulus and 180 iterations with HIO algorithm and partial coherence modulus, +| followed by 20 iterations with ER algorithm and partial coherence modulus. +| In addition, based on the triggers configuration, every three iterations starting from the first, for the rest of the run, the shrink wrap operation will be applied, +| and every two iterations, starting at the beginning, until iteration 130, the lowpass operation will be applied. | example2: :: - algorithm_sequence = "20*ER.SW0.PHM0 + 180*HIO.SW1.PHM1 + 20*ER.SW2" + algorithm_sequence = "20*ER.SW0.PHC0 + 180*HIO.SW1.PHC1 + 20*ER.SW2" config_rec: -| shrink_wrap_trigger = [[1, 2], [0, 3], [0, 4]] -| phm_trigger = [[0, 2, 100], [0, 3, 100]] +| shrink_wrap_trigger = [[1, 2], [0, 3], [0, 4]] +| phc_trigger = [[0, 2, 100], [0, 3, 100]] -| In this example the program will run 20 iterations with ER algorithm and modulus. During this 20 iterations the first sub-trigger in shrink_wrap_trigger, defined as [1,2] will be applied, which will start at iteration one and proceed every second iteration. During these iterations first phase modulus sub-trigger from phm_trigger, defined as [0, 2 ,100] will be applied starting at the beginning iteration, repeating every other iteration, and will continue until the 20 ER iterations, even though it was configured for more iterations. -| The ER iterations will be followed by 180 iterations with HIO algorithm and modulus and sub-triggers. During this 180 iterations the second sub-trigger in shrink_wrap_trigger, defined as [0, 3] will be applied, which will start at first HIO iteration and proceed every third iteration. During these iterations second phase modulus sub-trigger from phm_trigger, defined as [0, 3 ,100] will be applied that will start at the beginning iteration, repeat every third iteration, and will continue until the 100 HIO iterations. +| In this example the program will run: +| 20 iterations with ER algorithm and modulus. During this 20 iterations the first sub-trigger in shrink_wrap_trigger, defined as [1,2] will be applied, which will start at iteration one and proceed every second iteration. During these iterations first phase constrain sub-trigger from phc_trigger, defined as [0, 2 ,100] will be applied starting at the beginning iteration, repeating every other iteration, and will continue until the 20 ER iterations, even though it was configured for more iterations. +| The ER iterations will be followed by 180 iterations with HIO algorithm and modulus and sub-triggers. During this 180 iterations the second sub-trigger in shrink_wrap_trigger, defined as [0, 3] will be applied, which will start at first HIO iteration and proceed every third iteration. During these iterations the second phase constrain sub-trigger from phc_trigger, defined as [0, 3 ,100] will be applied that will start at the beginning iteration, repeat every third iteration, and will continue until the 100 HIO iterations. | This sequence will be follow by 20 iterations with ER algorithm and modulus and shrink wrap sub-trigger, defined as [0,4]. diff --git a/docs/source/experiment_space.rst b/docs/source/experiment_space.rst index 3835627..4c34ead 100755 --- a/docs/source/experiment_space.rst +++ b/docs/source/experiment_space.rst @@ -2,7 +2,7 @@ Experiment Space ================ | In order to group the files associated with an experiment in a structured manner user scripts will create a space dedicated to the experiment. The space, experiment directory, will be a sub-directory in working directory. The name of experiment directory is descriptive as it contains an ID and scan ranges. -| ex: with experiment ID of "ABC", and scan "56-78", the experiment directory is ABC_56-78 +| ex: with experiment ID of "ABC", and scan "56-78", the experiment directory is ABC_56-78. | Below we show experiment directory structure for various use cases. Single reconstruction @@ -26,7 +26,7 @@ Single reconstruction - The experiment should be set up with "conf" subdirectory containing configuration files. There are scripts that create the initial experiment space. - The script "beamline_preprocess.py" creates "prep_data.tif" file in "preprocessed_data" subdirectory. This is a file ready to be formatted. - The script "standard_preprocess.py" reads the "preprocessed_data.tif" file, formats it, and saves the result in the "phasing_data" subdirectory as "data.tif" file. This file is ready for reconstruction. -- "run_reconstruction.py" script reads "data.tif" file and runs phasing. The results are stored in the "results_phasing" subdirectory in "image.npy" file, along with "support.npy, and "coherence.npy" if partial coherence feature is configured. +- The "run_reconstruction.py" script reads "data.tif" file and runs phasing. The results are stored in the "results_phasing" subdirectory in "image.npy" file, along with "support.npy, and "coherence.npy" if partial coherence feature is configured. - The "beamline_visualization.py" script loads the arrays from "results_phasing" directory, processes the image and saves it in the results_viz directory. | After running all the scripts the experiment will have the following files: @@ -45,6 +45,10 @@ Single reconstruction | \|--results_phasing | \|--image.npy | \|--support.npy +| \|--image.tif +| \|--errors.npy +| \|--errors.txt +| \|--metrics.txt | \|--results_viz | \|--image.vts | \|--support.viz @@ -59,12 +63,24 @@ Multiple reconstruction | \|--0 | \|--image.npy | \|--support.npy +| \|--image.tif +| \|--errors.npy +| \|--errors.txt +| \|--metrics.txt | \|--1 | \|--image.npy | \|--support.npy +| \|--image.tif +| \|--errors.npy +| \|--errors.txt +| \|--metrics.txt | \|--2 | \|--image.npy | \|--support.npy +| \|--image.tif +| \|--errors.npy +| \|--errors.txt +| \|--metrics.txt | \|--results_viz | \|--0 | \|--image.vts @@ -79,42 +95,34 @@ Multiple reconstruction Genetic Algorithm +++++++++++++++++ | Results of reconstruction when using GA are reflected in relevant directory structure. The "results" directory will have subdirectories reflecting the generation, and each generation subdirectory will have subdirectories reflecting the runs. The generation directory is a concatenation of "g", underscore, and the generation number. -| Below is an example of "results" directory structure when running two generations and three reconstructions: +| Below is an example of "results" directory structure when running two generations and three reconstructions. Only the last generation is saved.: | | \| | \|--results_phasing -| \|--g_0 -| \|--0 -| \--image.npy -| \|--support.npy -| \|--1 -| \|--image.npy -| \|--support.npy -| \|--2 -| \|--image.npy -| \|--support.npy -| \|--g_1 +| \|--g_2 | \|--0 | \|--image.npy | \|--support.npy +| \|--image.tif +| \|--errors.npy +| \|--errors.txt +| \|--metrics.txt | \|--1 | \|--image.npy | \|--support.npy +| \|--image.tif +| \|--errors.npy +| \|--errors.txt +| \|--metrics.txt | \|--2 | \|--image.npy | \|--support.npy +| \|--image.tif +| \|--errors.npy +| \|--errors.txt +| \|--metrics.txt | \|--results_viz -| \|--g_0 -| \|--0 -| \|--image.vts -| \|--support.vts -| \|--1 -| \|--image.vts -| \|--support.vts -| \|--2 -| \|--image.vts -| \|--support.vts -| \|--g_1 +| \|--g_2 | \|--0 | \|--image.vts | \|--support.vts @@ -145,6 +153,10 @@ Separate scans | \|--results_phasing | \|--image.npy | \|--support.npy +| \|--image.tif +| \|--errors.npy +| \|--errors.txt +| \|--metrics.txt | \|--results_viz | \|--image.vts | \|--support.vts @@ -156,6 +168,10 @@ Separate scans | \|--results_phasing | \|--image.npy | \|--support.npy +| \|--image.tif +| \|--errors.npy +| \|--errors.txt +| \|--metrics.txt | \|--results_viz | \|--image.vts | \|--support.vts @@ -182,18 +198,30 @@ Alternate configuration | \|--results_phasing | \|--image.npy | \|--support.npy +| \|--image.tif +| \|--errors.npy +| \|--errors.txt +| \|--metrics.txt | \|--results_viz | \|--image.vts | \|--support.viz | \|--results_phasing_aa | \|--image.npy | \|--support.npy +| \|--image.tif +| \|--errors.npy +| \|--errors.txt +| \|--metrics.txt | \|--results_viz_aa | \|--image.vts | \|--support.viz | \|--results_phasing_bb | \|--image.npy | \|--support.npy +| \|--image.tif +| \|--errors.npy +| \|--errors.txt +| \|--metrics.txt | \|--results_viz_bb | \|--image.vts | \|--support.viz diff --git a/docs/source/for_developers.rst b/docs/source/for_developers.rst index 663896e..eb679fd 100755 --- a/docs/source/for_developers.rst +++ b/docs/source/for_developers.rst @@ -5,14 +5,14 @@ develop Installation for development ============================ -The best practice is to create conda environment allocated for the development. +The best practice is to create conda environment dedicated for the development. :: - conda create -n python=3.x -c conda-forge + conda create -y -n -c conda-forge python=3.10 mayavi pyqt scikit-image xrayutilities conda activate -| Clone the latest cohere repository from GitHub. This will include the cohere-ui directory with all of the cohere-ui content, such running scripts and example. +| Clone the latest cohere repository from GitHub. This will include the cohere-ui directory with all of the cohere-ui content, such users scripts and example. :: @@ -27,7 +27,7 @@ The best practice is to create conda environment allocated for the development. git checkout -b pip install -e . -| Go to cohere-ui directory and checkout the Dev branch and create your own branch, then run setup.py and install_pkgs.sh/install_pkgs.bat scripts. The setup.py script modifies paths from relative to absolute in the provided example configuration. The install_pkgs script installs python packages xrayutilities, mayavi, and pyqt that are required to run the cohere-ui. During the installation user must interact with dialog to agree to the steps when installing the packages. +| Go to cohere-ui directory and checkout the Dev branch and create your own branch, then run setup.py. The setup.py script modifies paths from relative to absolute in the provided example configuration. :: @@ -35,37 +35,33 @@ The best practice is to create conda environment allocated for the development. git checkout Dev git checkout -b python setup.py - sh install_pkgs.sh # for Linux and OS_X - install_pkgs.bat # for Windows | If planning to use GPUs, install the packages/libraries that you wish to use. :: - conda install cupy -c conda-forge # if using cupy library + conda install -y cupy=12.2.0 -c conda-forge # if using cupy library pip install torch # if using torch Adding new trigger ================== The design allows to add a new feature in a standardized way. Typical feature is defined by a trigger and supporting parameters. The following modifications/additions need to be done to add a new feature: - - In cohere_core/controller/phasing.py, Rec constructor, insert a new function name to the self.iter_functions list in the correct order. + - In cohere_core/controller/phasing.py, Rec constructor, insert a new function name ending with '_operation' to the self.iter_functions list in the correct order. - Implement the new trigger function in cohere_core/controller/phasing.py, Rec class. - - In cohere_core/controller/op_flow.py add code to add the new trigger row into flow array. - - In cohere_core/controller/params.py add code to parse trigger and parameters applicable to the new feature. - - In utilities/config_verifier.py add code to verify added parameters + - In cohere_core/controller/phasing.py add code to set any new defaults when creating Rec object. + - In utilities/config_verifier.py add code to verify added parameters. -Adding new feature for sub-trigger -================================== +Adding new sub-trigger +====================== If the new feature will be used in a context of sub-triggers, in addition to the above steps, the following modifications/additions need to be done: - - In cohere_core/controller/op_flow.py add entry in the sub_triggers, where key is the trigger function name, and value is the feature mnemonics - - In cohere_core/controller/phasing.py, Rec.init function, create_feat_objects sub-function, add the new feature object, created the same as shrink_wrap_obj, and other features. + - In cohere_core/controller/op_flow.py add entry in the sub_triggers dictionary, where key is the arbitrary assigned mnemonics, and value is the trigger name. + - In cohere_core/controller/phasing.py, Rec.init function, create_feat_objects sub-function, add the new feature object, created the same way as shrink_wrap_obj, and other features. - In cohere_core/controller/phasing.py, Rec class add the trigger function. The code inside should call the trigger on the feature object with *args. - in cohere_core/controller/features.py add new feature class. + | The constructor factory function create should have a new lines to construct the new object. | The feature class should be subclass of Feature and - | have defined self.key = and - | have defined self.trig_name = - | should have implemented create_obj function that creates sub-object(s) + | should have implemented create_obj function that creates sub-object(s) and | should have defined the sub-object(s) class(es). The embedded class contains the apply_trigger function that has the trigger code. Some features can be configured to different types and therefore multiple classes can be defined. | | The easiest way to implement the feature is to copy one already implemented and modify. @@ -73,7 +69,7 @@ If the new feature will be used in a context of sub-triggers, in addition to the Adding new algorithm ==================== The algorithm sequence defines functions executed during modulus projection and during modulus. Adding new algorithm requires the following steps: - - In cohere_core/controller/op_flow.py add entry in the algs dictionary, where key is the mnemonic used in algorithm_sequence, and value is the tuple defining functions, ex: 'ER': ('er', 'modulus') + - In cohere_core/controller/op_flow.py add entry in the algs dictionary, where key is the mnemonic used in algorithm_sequence, and value is the tuple defining functions, ex: 'ER': ('to_reciprocal_space', 'modulus', 'to_direct_space', 'er') - In cohere_core/controller/phasing.py, Rec constructor, insert a new function name to the self.iter_functions list in the correct order. - In cohere_core/controller/phasing.py, Rec class add the new algorithm function(s). diff --git a/docs/source/how_to_use.rst b/docs/source/how_to_use.rst index 8efd3f9..5cf32a9 100755 --- a/docs/source/how_to_use.rst +++ b/docs/source/how_to_use.rst @@ -7,7 +7,7 @@ Using cohere with cohere-ui package Installing Scripts ################## | There are three ways to get the cohere-ui package onto your local computer, or remote computer. -| It is assumed the user activated the environment where cohere was installed. +| User should activate the environment where cohere was installed. :: @@ -28,7 +28,7 @@ Installing Scripts wget https://github.com/AdvancedPhotonSource/cohere-ui/archive/main.zip unzip main.zip -| After the package is installed run setup.py and install_pkgs.sh/install_pkgs.bat scripts. The setup.py script modifies paths from relative to absolute in the provided example configuration. The install_pkgs script installs python packages xrayutilities, mayavi, and pyqt that are required to run the cohere-ui. During the installation user must interact with dialog to agree to the steps when installing the packages.: +| After the package is installed run setup.py and install_pkgs.sh/install_pkgs.bat scripts. The setup.py script modifies paths from relative to absolute in the provided example configuration. The install_pkgs script installs python packages xrayutilities, mayavi, and pyqt that are required to run the cohere-ui. :: diff --git a/docs/source/index.rst b/docs/source/index.rst index 0f775bf..8bdd3a5 100755 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -8,7 +8,7 @@ Welcome to cohere documentation! about installation how_to_use - define_formula + define_alg_seq api_standard_preprocess api_phasing api_utils diff --git a/docs/source/installation.rst b/docs/source/installation.rst index b02c7ff..94c7469 100755 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -10,22 +10,25 @@ Conda Installation This will install latest official release. First you must have `Conda ` installed. -It is highly recommended to install the cohere_core package in conda environment. Supported python versions are 3.6 - 3.11 for Linux and Mac and 3.6 - 3.10 for Windows. +It is highly recommended to install the cohere_core package in conda environment. Supported python versions are 3.6 - 3.10. | To create and activate the environment run the commands below:: - conda create --name python=3.x -c conda-forge + conda create -y -n -c conda-forge python=3.10 conda activate -then run the installation command:: +then run the cohere installation command:: pip install cohere-core -The cohere_core package can be run as installed running on numpy library. Other available libraries are cupy and torch. +The cohere_core package can be run utilizing numpy library. Other available libraries are cupy and torch. User has to install the preferred library. -if using cupy library:: +If using cupy library:: + + conda -y install cupy=12.2.0 -c conda-forge +If torch is preferred, install with command:: - conda install cupy -c conda-forge -Cohere-core package does not install python packages used by user's scripts in cohere-ui package. If planning to use the scripts Refer to :ref:`use` page, section "Installing Scripts". + conda -y install -c conda-forge torch +The cohere-core package does not install python packages used by user's scripts in cohere-ui package. If planning to use the scripts Refer to :ref:`use` page, section "Installing Scripts". Note: The cupy installation on Windows may result in incompatible libraries, which makes the environment unusable. Run the repack.bat script from cohere-ui package and try running again. @@ -35,7 +38,7 @@ This will install the latest development. It might be not 100 percent tested. Th Create environment, activate it and clone cohere repository. It contains the cohere-ui submodule. Run the following commands:: - conda create --name python=3.x -c conda-forge + conda create --name -c conda-forge python=3.10 mayavi pyqt scikit-image xrayutilities conda activate git clone https://github.com/advancedPhotonSource/cohere --recurse-submodules cd cohere @@ -44,12 +47,13 @@ Create environment, activate it and clone cohere repository. It contains the coh cd cohere-ui git checkout Dev python setup.py - sh cohere-ui/install_pkgs.sh # for Linux and OS_X - cohere-ui/install_pkgs.bat # for Windows -For Windows make sure that numpy version is 1.23.5. The commands above will create conda environment and activate it, clone the packages, get the Dev branch, install, initialize. The install_pkgs script is interactive and the user must confirm the pacakages installations. -If using cupy library:: +For Windows make sure that numpy version is 1.23.5. The commands above will create conda environment and activate it, clone the packages, get the Dev branch, install, initialize. +| If using cupy library:: + + conda -y install cupy=12.2.0 -c conda-forge +| If using torch library:: - conda install cupy -c conda-forge + pip install torch After installation you may start using scripts from this directory, for example:: python cohere-scripts/cdi_window.py diff --git a/meta.yaml b/meta.yaml deleted file mode 100644 index dce1353..0000000 --- a/meta.yaml +++ /dev/null @@ -1,27 +0,0 @@ -package: - name: cohere_core - version: "4.0.0" - -source: - path: . - -build: - number: 0 - -requirements: - host: - - python=3.11 - - numpy - - run: - - python=3.11 - - numpy - - tensorflow - - scikit-learn - - tifffile - -about: - home: https://github.com/advancedPhotonSource/cohere - license: BSD - license_file: LICENSE - summary: AI, parallelized genetic algorithms and phase retrieval methods for Bragg CDI techniques. diff --git a/requirements.txt b/requirements.txt index 24ce15a..9c558e3 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1 +1 @@ -numpy +. diff --git a/setup.py b/setup.py index 441fdcc..52bcdb1 100644 --- a/setup.py +++ b/setup.py @@ -5,12 +5,16 @@ author = 'Barbara Frosik, Ross Harder', author_email = 'bfrosik@anl.gov', url='https://github.com/advancedPhotonSource/cohere', - version='4.0.0', + version='4.1-beta1', packages=setuptools.find_packages(), install_requires=['numpy', 'scikit-learn', 'tifffile', 'tensorflow', + 'mpi4py-mpich', + 'gputil', + 'psutil', + 'tqdm' ], classifiers=[ 'Intended Audience :: Science/Research', @@ -19,7 +23,6 @@ 'Programming Language :: Python :: 3.8', 'Programming Language :: Python :: 3.9', 'Programming Language :: Python :: 3.10', - 'Programming Language :: Python :: 3.11', ], )