From 5d14890e828d3a58877c642e9457b55eef42ad0f Mon Sep 17 00:00:00 2001 From: marco cammarata Date: Wed, 26 Apr 2017 09:04:01 +0200 Subject: [PATCH] storage and xray modules live now on their own --- storage.py | 298 ------------------------ xray/README.md | 24 -- xray/__init__.py | 12 - xray/azav.py | 465 ------------------------------------- xray/cell.py | 98 -------- xray/dataReduction.py | 261 --------------------- xray/example_main_salen.py | 356 ---------------------------- xray/example_main_tiox.py | 68 ------ xray/filters.py | 85 ------- xray/id9.py | 212 ----------------- xray/mask.py | 210 ----------------- xray/peaks.py | 32 --- xray/utils.py | 346 --------------------------- xray/xray.py | 49 ---- 14 files changed, 2516 deletions(-) delete mode 100644 storage.py delete mode 100644 xray/README.md delete mode 100644 xray/__init__.py delete mode 100644 xray/azav.py delete mode 100644 xray/cell.py delete mode 100644 xray/dataReduction.py delete mode 100644 xray/example_main_salen.py delete mode 100644 xray/example_main_tiox.py delete mode 100644 xray/filters.py delete mode 100644 xray/id9.py delete mode 100644 xray/mask.py delete mode 100644 xray/peaks.py delete mode 100644 xray/utils.py delete mode 100644 xray/xray.py diff --git a/storage.py b/storage.py deleted file mode 100644 index c933c69..0000000 --- a/storage.py +++ /dev/null @@ -1,298 +0,0 @@ -""" npy/npz/hdf5 file based storage; - this modules adds the possibility to dump and load objects in files and - a more convenient was of accessing the data via the .attributedict thanks - to the DataStorage class """ -import numpy as np -import os -import h5py -import collections -import logging -log = logging.getLogger(__name__) - -_array_cache = {} - -def unwrapArray(a,recursive=True,readH5pyDataset=True): - """ This function takes an object (like a dictionary) and recursively - unwraps it solving issues like: - * the fact that many objects are packaged as 0d array - This funciton has also some specific hack for handling h5py limits: - * handle the None python object - * numpy unicode ... - """ - # is h5py dataset convert to array - if isinstance(a,h5py.Dataset) and readH5pyDataset: a = a[...] - if isinstance(a,h5py.Dataset) and a.shape == (): a = a[...] - if isinstance(a,h5py.Group) and "IS_LIST_OF_ARRAYS" in a.attrs: - items = list(a.keys()) - items.sort() - a = np.asarray( [a[item][...] for item in items] ) - if isinstance(a,np.ndarray) and a.ndim == 0 : a = a.item() - if isinstance(a,np.ndarray) and a.dtype.char == "S": a = a.astype(str) - if recursive: - if "items" in dir(a): # dict, h5py groups, npz file - a = dict(a); # convert to dict, otherwise can't asssign values - for key,value in a.items(): a[key] = unwrapArray(value) - elif isinstance(a,(list,tuple)): - a = [unwrapArray(element) for element in a] - else: - pass - if isinstance(a,dict): a = DataStorage(a) - # restore None that cannot be saved in h5py - if isinstance(a,str) and a == "NONE_PYTHON_OBJECT": a = None - # h5py can't save numpy unicode - if isinstance(a,np.ndarray) and a.dtype.char == "S": a = a.astype(str) - return a - -def dictToH5Group(d,group,link_copy=True): - """ helper function that transform (recursive) a dictionary into an - hdf group by creating subgroups - link_copy = True, tries to save space in the hdf file by creating an internal link. - the current implementation uses memory though ... - """ - global _array_cache - for key,value in d.items(): - TOTRY = True - if isinstance(value,(list,tuple)): value = np.asarray(value) - if isinstance(value,dict): - group.create_group(key) - dictToH5Group(value,group[key],link_copy=link_copy) - elif value is None: - group[key] = "NONE_PYTHON_OBJECT" - elif isinstance(value,np.ndarray): - # take care of unicode (h5py can't handle numpy unicode arrays) - if value.dtype.char == "U": value = np.asarray([vv.encode('ascii') for vv in value]) - # check if it is list of array - elif isinstance(value,np.ndarray) and value.ndim == 1 and isinstance(value[0],np.ndarray): - group.create_group(key) - group[key].attrs["IS_LIST_OF_ARRAYS"] = True - for index,array in enumerate(value): dictToH5Group( { "index%010d"%index : array},group[key],link_copy=link_copy ); - TOTRY = False; # don't even try to save as generic call group[key]=value - if link_copy: - found_address = None - for address,(file_handle,array) in _array_cache.items(): - if np.array_equal(array,value) and group.file == file_handle: - log.info("Found array in cache, asked for %s/%s, found as %s"%(group.name,key,address)) - found_address = address - break - if found_address is not None: - value = group.file[found_address] - try: - if TOTRY: - group[key] = value - if link_copy: - log.info("Addind array %s to cache"%(group[key].name)) - _array_cache[ group[key].name ] = (group.file,value) - except Exception as e: - log.warning("Can't save %s, error was %s"%(key,e)) - # try saving everything else that is not dict or array - else: - try: - group[key] = value - except Exception as e: - log.error("Can't save %s, error was %s"%(key,e)) - - -def dictToH5(h5,d,link_copy=False): - """ Save a dictionary into an hdf5 file - TODO: add capability of saving list of array - h5py is not capable of handling dictionaries natively""" - h5 = h5py.File(h5,mode="w") -# group = h5.create_group("/") - dictToH5Group(d,h5["/"],link_copy=link_copy) - h5.close() - -def h5ToDict(h5,readH5pyDataset=True): - """ Read a hdf5 file into a dictionary """ - with h5py.File(h5,"r") as h: - ret = unwrapArray(h,recursive=True,readH5pyDataset=readH5pyDataset) - return ret - -def npzToDict(npzFile): - with np.load(npzFile) as npz: d = dict(npz) - d = unwrapArray(d,recursive=True) - return d - -def npyToDict(npyFile): - d = unwrapArray( np.load(npyFile).item() ,recursive=True) - return d - -def dictToNpz(npzFile,d): np.savez(npzFile,**d) -def dictToNpy(npyFile,d): np.save(npyFile,d) - -def objToDict(o,recursive=True): - """ convert a DictWrap to a dictionary (useful for saving); it should work for other objects too - TODO: this function does not catch a list of DataStorage instances like - objToDict( ( DataStorage(), DataStorage() ) ) - is not converted !! - """ - if "items" not in dir(o): return o - d = dict() - for k,v in o.items(): - try: - d[k] = objToDict( v ) - except Exception as e: - log.info("In objToDict, could not convert key %s to dict, error was"%\ - (k,e)) - d[k] = v - return d - - -def read(fname): - extension = os.path.splitext(fname)[1] - log.info("Reading storage file %s"%fname) - if extension == ".npz": - return DataStorage(npzToDict(fname)) - elif extension == ".npy": - return DataStorage(npyToDict(fname)) - elif extension == ".h5": - return DataStorage(h5ToDict(fname)) - else: - raise ValueError("Extension must be h5, npy or npz, it was %s"%extension) - -def save(fname,d,link_copy=True): - """ link_copy is used by hdf5 saving only, it allows to creat link of identical arrays (saving space) """ - # make sure the object is dict (recursively) this allows reading it - # without the DataStorage module - d = objToDict(d,recursive=True) - d['filename'] = fname - extension = os.path.splitext(fname)[1] - log.info("Saving storage file %s"%fname) - try: - if extension == ".npz": - return dictToNpz(fname,d) - elif extension == ".h5": - return dictToH5(fname,d,link_copy=link_copy) - elif extension == ".npy": - return dictToNpy(fname,d) - else: - raise ValueError("Extension must be h5, npy or npz, it was %s"%extension) - except Exception as e: - log.exception("Could not save %s"%fname) - - -class DataStorage(dict): - """ Storage for dict like object. - recursive : bool - recursively convert dict-like objects to DataStorage - It can save data to file (format npy,npz or h5) - - To initialize it: - - data = DataStorage( a=(1,2,3),b="add",filename='store.npz' ) - - # recursively by default - # data.a will be a DataStorage instance - data = DataStorage( a=dict( b = 1)) ); - - # data.a will be a dictionary - data = DataStorage( a=dict( b = 1),recursive=False ) - - # reads from file if it exists - data = DataStorage( 'mysaveddata.npz' ) ; - - DOES NOT READ FROM FILE (even if it exists)!! - data = DataStorage( filename = 'mysaveddata.npz' ); - - create empty storage (with default filename) - data = DataStorage() - """ - def __init__(self,*args,filename='data_storage.npz',recursive=True,**kwargs): -# self.filename = kwargs.pop('filename',"data_storage.npz") - self.filename = filename - self._recursive = recursive - # interpret kwargs as dict if there are - if len(kwargs) != 0: - fileOrDict = dict(kwargs) - elif len(kwargs)==0 and len(args)>0: - fileOrDict = args[0] - else: - fileOrDict = dict() - - - d = dict(); # data dictionary - if isinstance(fileOrDict,dict): - d = fileOrDict - elif isinstance(fileOrDict,str): - if os.path.isfile(fileOrDict): - d = read(fileOrDict) - else: - self.filename=fileOrDict - d = dict() - else: - raise ValueError("Invalid DataStorage definition") - - if recursive: - for k in d.keys(): - if not isinstance(d[k],DataStorage) and isinstance(d[k],dict): - d[k] = DataStorage(d[k]) - - # allow accessing with .data, .delays, etc. - for k,v in d.items(): setattr(self,k,v) - - # allow accessing as proper dict - self.update( **dict(d) ) - - def __setitem__(self, key, value): - #print("__setitem__") - setattr(self,key,value) - super().__setitem__(key, value) - - def __setattr__(self, key, value): - """ allows to add fields with data.test=4 """ - # check if attr exists is essential (or it fails when defining an instance) - if hasattr(self,"_recursive") and self._recursive and \ - isinstance(value,(dict,collections.OrderedDict)): value = DataStorage(value) - super().__setitem__(key, value) - super().__setattr__(key,value) - - def __delitem__(self, key): - delattr(self,key) - super().__delitem__(key) - - def __str__(self): - keys = list(self.keys()) - keys.sort() - return "DataStorage obj containing: %s" % ",".join(keys) - - def __repr__(self): - keys = list(self.keys()) - keys.sort() - if len(keys) == 0: return "Empty DataStorage" - nchars = max(map(len,keys)) - fmt = "%%%ds %%s" % (nchars) - s = ["DataStorage obj containing (sorted): ",] - for k in keys: - if k[0] == "_": continue - obj = self[k] - if ( (isinstance(obj,np.ndarray) and obj.ndim == 1) or isinstance(obj,(list,tuple))) and all( [isinstance(v,np.ndarray) for v in obj]): - value_str = "list of arrays, shapes " + ",".join([str(v.shape) for v in obj[:5]]) + " ..." - elif isinstance(obj,np.ndarray): - value_str = "array, size %s, type %s"% ("x".join(map(str,obj.shape)),obj.dtype) - elif isinstance(obj,DataStorage): - value_str = str(obj)[:50] - elif isinstance(obj,(str,DataStorage)): - value_str = obj[:50] - elif self[k] is None: - value_str = "None" - else: - value_str = str(self[k]) - if len(str(obj))>50: value_str += " ..." - s.append( fmt % (k,value_str) ) - return "\n".join(s) - - def keys(self): - keys = list(super().keys()) - keys = [k for k in keys if k != 'filename' ] - keys = [k for k in keys if k[0] != '_' ] - return keys - - def save(self,fname=None,link_copy=False): - """ link_copy: only works in hfd5 format - save space by creating link when identical arrays are found, - it slows down the saving (3 or 4 folds) but saves A LOT of space - when saving different dataset together (since it does not duplicate - internal pyfai matrices - """ - if fname is None: fname = self.filename - assert fname is not None - save(fname,self,link_copy=link_copy) diff --git a/xray/README.md b/xray/README.md deleted file mode 100644 index eed31d3..0000000 --- a/xray/README.md +++ /dev/null @@ -1,24 +0,0 @@ -# xray - -Different utilities to work with time resolved solution/powder scattering data -They should be extended to absorption soon - -The idea is an hierarchical structure of modules - -> user → beamline → workhorses - -## Workhorses -azav.py: uses pyFAI to convert images into 1d curves -dataReduction.py: works on the normalization, referencing and averaging. - -Both modules uses some common utilities in utils.py including file based storage for persistency and convenience (numpy or hdf5 based); -Indeed the data_storage class wraps around dictionaries allowing accessing with .attribute syntax. -In utils plotfuncitons are presents - -## beamline -this step is needed to read the diagnostic information (time delay, monitors, etc). An example for current id9 data collection macros is given (id9.py). Essentially it means creating a dictionary where each key will be added to the main data dictionary (that has keys like q, data, files, etc.) - -## user -at the user level, a rather simple macro like the one provided as an example (example_main_tiox.py) should be sufficient for most needs. - - diff --git a/xray/__init__.py b/xray/__init__.py deleted file mode 100644 index e57bba7..0000000 --- a/xray/__init__.py +++ /dev/null @@ -1,12 +0,0 @@ -from .. import storage -from . import utils -from . import mask -from . import peaks -from . import cell -from . import filters -from . import id9 -try: - from . import azav -except ImportError: - print("Can't import azav ...") -from . import dataReduction diff --git a/xray/azav.py b/xray/azav.py deleted file mode 100644 index 4aa81c2..0000000 --- a/xray/azav.py +++ /dev/null @@ -1,465 +0,0 @@ -from __future__ import print_function,division - -import logging -log = logging.getLogger(__name__) - - -import numpy as np -np.seterr(all='ignore') -import os -import collections -import glob -import pathlib -from . import storage -from . import utils -from . import filters -import re -import fabio -import pyFAI - -try: - import matplotlib.pyplot as plt -except ImportError: - log.warn("Can't import matplotlib !") - -def _read(fname): - """ read data from file using fabio """ - f = fabio.open(fname) - data = f.data - del f; # close file - return data - -def read(fnames): - """ read data from file(s) using fabio """ - if isinstance(fnames,str): - data = _read(fnames) - else: - # read one image to know img size - temp = _read(fnames[0]) - shape = [len(fnames),]+list(temp.shape) - data = np.empty(shape) - data[0] = temp - for i in range(1,len(fnames)): data[i] = _read(fnames[i]) - return data - -def ai_as_dict(ai): - """ ai is a pyFAI azimuthal intagrator""" - methods = dir(ai) - methods = [m for m in methods if m.find("get_") == 0] - names = [m[4:] for m in methods] - values = [getattr(ai,m)() for m in methods] - ret = dict( zip(names,values) ) - ret["detector"] = ai.detector.get_name() - return ret - -def ai_as_str(ai): - """ ai is a pyFAI azimuthal intagrator""" - s=[ "# Detector : %s" % ai.detector.name, - "# Pixel [um] : %.2fx%.2f" % (ai.pixel1*1e6,ai.pixel2*1e6), - "# Distance [mm] : %.3f" % (ai.dist*1e3), - "# Center [mm] : %.3f,%.3f" % (ai.poni1*1e3,ai.poni2*1e3), - "# Center [px] : %.3f,%.3f" % (ai.poni1/ai.pixel1,ai.poni2/ai.pixel2), - "# Wavelength [A] : %.5f" % (ai.wavelength*1e10), - "# rot[1,2,3] [rad]: %.3f,%.3f,%.3f" % (ai.rot1,ai.rot2,ai.rot3) ] - return "\n".join(s) - -def do1d(ai, imgs, mask = None, npt_radial = 600, method = 'csr',safe=True,dark=10., polCorr = 1,dezinger=None): - """ ai is a pyFAI azimuthal intagrator - it can be defined with pyFAI.load(ponifile) - dezinger: None or float (used as percentile of ai.separate) - mask: True are points to be masked out """ - # force float to be sure of type casting for img - if isinstance(dark,int): dark = float(dark); - if imgs.ndim == 2: imgs = (imgs,) - out_i = np.empty( ( len(imgs), npt_radial) ) - out_s = np.empty( ( len(imgs), npt_radial) ) - for _i,img in enumerate(imgs): - if dezinger is not None and dezinger > 0: - _,img=ai.separate(img,npt_rad=npt_radial,npt_azim=512,unit='q_A^-1', - method=method,mask=mask,percentile=dezinger) - q,i, sig = ai.integrate1d(img-dark, npt_radial, mask= mask, safe = safe,\ - unit="q_A^-1", method = method, error_model = "poisson", - polarization_factor = polCorr) - out_i[_i] = i - out_s[_i] = sig - return q,np.squeeze(out_i),np.squeeze(out_s) - -def do2d(ai, imgs, mask = None, npt_radial = 600, npt_azim=360,method = 'csr',safe=True,dark=10., polCorr = 1): - """ ai is a pyFAI azimuthal intagrator - it can be defined with pyFAI.load(ponifile) - mask: True are points to be masked out """ - # force float to be sure of type casting for img - if isinstance(dark,int): dark = float(dark); - if imgs.ndim == 2: imgs = (imgs,) - out = np.empty( ( len(imgs), npt_azim,npt_radial) ) - for _i,img in enumerate(imgs): - i2d,q,azTheta = ai.integrate2d(img-dark, npt_radial, npt_azim=npt_azim, - mask= mask, safe = safe,unit="q_A^-1", method = method, - polarization_factor = polCorr ) - out[_i] = i2d - return q,azTheta,np.squeeze(out) - -def getAI(poni=None,folder=None,**kwargs): - """ get AzimuthalIntegrator instance: - → if poni is a string, it is used as filename to read. - in this case if folder is given it is used (together with all its - subfolder) as search path (along with ./ and home folder) - → kwargs if present can be used to define (or override) parameters from files, - dist,xcen,ycen,poni1,poni2,rot1,rot2,rot3,pixel1,pixel2,splineFile, - detector,wavelength - """ - if isinstance(poni,pyFAI.azimuthalIntegrator.AzimuthalIntegrator): - ai = poni - elif isinstance(poni,str): - # look is file exists in cwd - if os.path.isfile(poni): - fname = poni - # if file does not exist look for one with that name around - else: - # build search paths - folders = [] - if folder is not None: - temp = os.path.abspath(folder) - path = pathlib.Path(temp) - folders = [ str(path), ] - for p in path.parents: folders.append(str(p)) - folders.append( "./" ) - folders.append( os.path.expanduser("~/") ) - # look for file - for path in folders: - fname = path + "/" + poni - if os.path.isfile(fname): - log.info("Found poni file %s",fname) - break - else: - log.debug("Could not poni file %s",fname) - ai = pyFAI.load(fname) - else: - ai = pyFAI.azimuthalIntegrator.AzimuthalIntegrator() - for par,value in kwargs.items(): setattr(ai,par,value) - # provide xcen and ycen for convenience (note: xcen changes poni2 - # and ycen changes poni1) - if 'xcen' in kwargs: ai.poni2 = kwargs['xcen'] * ai.pixel2 - if 'ycen' in kwargs: ai.poni1 = kwargs['ycen'] * ai.pixel1 - ai.reset(); # needed in case of overridden parameters - return ai - -g_mask_str = re.compile("(\w)\s*(<|>)\s*(\d+)") - -def _interpretMask(mask,shape=None): - """ - if mask is an existing filename, returns it - if mask is a string like [x|y] [<|>] int; - for example y>500 will dis-regard out for y>500 - """ - maskout = None - if isinstance(mask,str) and os.path.isfile(mask): - maskout = read(mask).astype(np.bool) - elif isinstance(mask,str) and not os.path.isfile(mask): - err_msg = ValueError("The string '%s' could not be interpreted as simple\ - mask; it should be something like x>10"%mask) - assert shape is not None - # interpret string - maskout = np.zeros(shape,dtype=bool) - match = g_mask_str.match(mask) - if match is None: raise err_msg - (axis,sign,lim) = match.groups() - if axis not in ("x","y"): raise err_msg - if sign not in (">","<"): raise err_msg - lim = int(lim) - idx = slice(lim,None) if sign == ">" else slice(None,lim) - if axis == 'y': - maskout[idx,:] = True - else: - maskout[:,idx] = True - elif isinstance(mask,np.ndarray): - maskout = mask.astype(np.bool) - elif mask is None: - assert shape is not None - maskout = np.zeros(shape,dtype=bool) - else: - maskout = None - raise ValueError("Could not interpret %s as mask input"%mask) - - if shape is not None and maskout.shape != shape: - raise ValueError("The mask shape %s does not match the shape given as\ - argument %s"%(maskout.shape,shape)) - return maskout - - -def interpretMask(masks,shape=None): - """ - if masks is a list of masks, eachone can be: - * an existing filename - * a string like [x|y] [<|>] int; - """ - if not isinstance( masks, (list,tuple,np.ndarray) ): - masks = (masks,) - masks = [_interpretMask(mask,shape) for mask in masks] - # put them all together - mask = masks[0] - for m in masks[1:]: - mask = np.logical_or(mask,m) - return mask - -def removeBackground(data,qlims=(0,10),max_iter=30,background_regions=[],force=False, - storageFile=None,save=True,**removeBkg): - """ similar function to the zray.utils one, this works on dataset created by - doFolder """ - idx = utils.findSlice(data.orig.q,qlims) - # see if there are some to do ... - if force: - idx_start = 0 - else: - idx_start = len(data.data) - if idx_start < len(data.orig.data): - _q,_data = utils.removeBackground(data.orig.q[idx],data.orig.data[:,idx], - max_iter=max_iter,background_regions=background_regions,**removeBkg) - data.q = _q - data.data = np.concatenate( (data.data,_data ) ) - data.err = np.concatenate( (data.err ,data.err[idx_start,idx] ) ) - if save: data.save(storageFile); # if None uses .filename - return data - - -def doFolder(folder,files='*.edf*',nQ = 1500,force=False,mask=None,dark=10, - norm='auto',save_pyfai=False,saveChi=True,poni='pyfai.poni', - storageFile='auto',save=True,logDict=None,dezinger=None,skip_first=0,last=None): - """ calc 1D curves from files in folder, returning a dictionary of stuff - nQ : number of Q-points (equispaced) - force : if True, redo from beginning even if previous data are found - if False, do only new files - mask : can be a list of [filenames|array of booleans|mask string] - pixels that are True are dis-regarded - saveChi: self-explanatory - dezinger: None or 0 to disable; good value is ~50. Needs good center and mask - logDict: dictionary(-like) structure. has to have 'file' key - save_pyfai: store all pyfai's internal arrays (~110 MB) - poni : could be: - → an AzimuthalIntegrator instance - → a filename that will be look for in - 1 'folder' first - 2 in ../folder - 3 in ../../folder - .... - n-1 in pwd - n in homefolder - → a dictionary (use to bootstrap an AzimuthalIntegrator using - AzimuthalIntegrator(**poni) - """ - if storageFile == 'auto': storageFile = folder + "/" + "pyfai_1d.h5" - - if os.path.isfile(storageFile) and not force: - saved = storage.DataStorage(storageFile) - log.info("Found %d images in storage file"%saved.data.shape[0]) - else: - saved = None - - files = utils.getFiles(folder,files) - if logDict is not None: - files = [f for f in files if utils.getBasename(f) in logDict['file'] ] - files = files[skip_first:last] - - - if saved is not None: - files = [f for f in files if f not in saved["files"]] - log.info("Will do azimuthal integration for %d files"%(len(files))) - - files = np.asarray(files) - basenames = np.asarray( [ utils.getBasename(file) for file in files] ) - - if len(files) > 0: - # which poni file to use: - ai = getAI(poni,folder) - - - shape = read(files[0]).shape - mask = interpretMask(mask,shape) - - data = np.empty( (len(files),nQ) ) - err = np.empty( (len(files),nQ) ) - for ifname,fname in enumerate(files): - img = read(fname) - q,i,e = do1d(ai,img,mask=mask,npt_radial=nQ,dark=dark,dezinger=dezinger) - data[ifname] = i - err[ifname] = e - if saveChi: - chi_fname = utils.removeExt(fname) + ".chi" - utils.saveTxt(chi_fname,q,np.vstack((i,e)),info=ai_as_str(ai),overwrite=True) - - if saved is not None: - files = np.concatenate( (saved["files"] ,basenames ) ) - data = np.concatenate( (saved["data"] ,data ) ) - err = np.concatenate( (saved["err"] ,err ) ) - theta_rad = utils.qToTheta(q,wavelength=ai.wavelength) - theta_deg = utils.qToTheta(q,wavelength=ai.wavelength,asDeg=True) - orig = dict(data=data.copy(),err=err.copy(),q=q.copy()) - ret = dict(q=q,folder=folder,files=files,data=data,err=err, - orig = orig, theta_rad = theta_rad, theta_deg=theta_deg, - pyfai=ai_as_dict(ai),pyfai_info=ai_as_str(ai),mask=mask) - if not save_pyfai: - ret['pyfai']['chia'] = None - ret['pyfai']['dssa'] = None - ret['pyfai']['q'] = None - ret['pyfai']['ttha'] = None - - - ret = storage.DataStorage(ret) - - # add info from logDict if provided - if logDict is not None: - ret['log']=logDict - # sometime saving is not necessary (if one has to do it after subtracting background - if storageFile is not None and save: ret.save(storageFile) - else: - ret = saved - return ret - - -def removeBackground(data,qlims=(0,10),max_iter=30,background_regions=[], - storageFile=None,save=True,**removeBkg): - """ similar function to the zray.utils one, this works on dataset created by - doFolder """ - idx = utils.findSlice(data.q,qlims) - data.q,data.data = utils.removeBackground(data.q[idx],data.data[:,idx], - max_iter=max_iter,background_regions=background_regions,**removeBkg) - data.err = data.err[idx] - if save: data.save(storageFile); # if None uses .filename - return data - -def _calc_R(x,y, xc, yc): - """ calculate the distance of each 2D points from the center (xc, yc) """ - return np.sqrt((x-xc)**2 + (y-yc)**2) - -def _chi2(c, x, y): - """ calculate the algebraic distance between the data points and the mean - circle centered at c=(xc, yc) """ - Ri = _calc_R(x, y, *c) - return Ri - Ri.mean() - -def leastsq_circle(x,y): - from scipy import optimize - # coordinates of the barycenter - center_estimate = np.nanmean(x), np.nanmean(y) - center, ier = optimize.leastsq(_chi2, center_estimate, args=(x,y)) - xc, yc = center - Ri = _calc_R(x, y, *center) - R = Ri.mean() - residu = np.sum((Ri - R)**2) - return xc, yc, R - -def find_center(img,psize=100e-6,dist=0.1,wavelength=0.8e-10,center=None,reference=None,**kwargs): - """ center is the initial centr (can be None) - reference is a reference position to be plot in 2D plots """ - plt.ion() - kw = dict( pixel1 = psize, pixel2 = psize, dist = dist,wavelength=wavelength ) - kw.update(kwargs) - ai = pyFAI.azimuthalIntegrator.AzimuthalIntegrator(**kw) - fig_img,ax_img = plt.subplots(1,1) - fig_pyfai,ax_pyfai = plt.subplots(1,1) - fig_pyfai = plt.figure(2) - temp= ax_img.imshow(img) - plt.sca(ax_img); # set figure to use for mouse interaction - temp.set_clim( *np.percentile(img,(2,95) ) ) - ans = "" - print("Enter 'end' when done") - while ans != "end": - if center is None: - print("Click on beam center:") - plt.sca(ax_img); # set figure to use for mouse interaction - center = plt.ginput()[0] - print("Selected center:",center) - ai.set_poni1(center[1]*psize) - ai.set_poni2(center[0]*psize) - q,az,i = do2d(ai,img) - mesh = ax_pyfai.pcolormesh(q,az,i) - mesh.set_clim( *np.percentile(i,(2,95) ) ) - ax_pyfai.set_title(str(center)) - if reference is not None: ax_pyfai.axvline(reference) - plt.pause(0.01) - plt.draw() - plt.draw_all() - ans=input("Enter to continue with clinking or enter xc,yc values ") - if ans == '': - center = None - else: - try: - center = list(map(float,ans.split(","))) - except Exception as e: - center = None - if center == []: center = None - print("Final values: (in pixels) %.3f %.3f"%(center[0],center[1])) - return ai - -def average(fileOrFolder,delays=slice(None),scale=1,norm=None,returnAll=False,plot=False, - showTrend=False): - data = storage.DataStorage(fileOrFolder) - if isinstance(delays,slice): - idx = np.arange(data.delays.shape[0])[delays] - elif isinstance(delays,(int,float)): - idx = data.delays == float(delays) - else: - idx = data.delays < 0 - if idx.sum() == 0: - print("No data with the current filter") - return None - i = data.data[idx] - q = data.q - if isinstance(norm,(tuple,list)): - idx = ( q>norm[0] ) & (qnorm[0] ) & (q=idx_ref_after-1: _ref += 1 - log.debug("For image %d : %d-%d"%(_i,idx_ref_before,idx_ref_after)) - # take care of the reference for the references ... - if len(idx_ref) > 2: - iref[idx_ref[0]] = i[idx_ref[1]] - iref[idx_ref[-1]] = i[idx_ref[-2]] - for _i in range(1,len(idx_ref)-1): - idx_ref_before = idx_ref[_i-1] - idx_ref_after = idx_ref[_i+1] - ref_before = i[idx_ref_before] - ref_after = i[idx_ref_after] - weight_before = float(idx_ref[_i]-idx_ref_before)/(idx_ref_after-idx_ref_before) - weight_after = 1-weight_before - # normal reference for an on chi, the weighted average - iref[idx_ref[_i]] = weight_before*ref_before + weight_after*ref_after - log.debug("For reference image %d : %d-%d"%(idx_ref[_i],idx_ref_before,idx_ref_after)) - else: - #print(idx_ref) - #print(iref[idx_ref]) - iref[idx_ref]=i[idx_ref[0]] - #print(iref[idx_ref]) - if useRatio: - i /= iref - else: - i -= iref - return i - -def averageScanPoints(scan,data,errAbs=None,isRef=None,lpower=None,useRatio=False,\ - funcForAveraging=np.nanmean): - """ given scanpoints in 'scan' and corresponding data in 'data' - average all data corresponding the exactly the same scanpoint. - If the values in scan are coming from a readback, rounding might be - necessary. - No normalization is done inside this function - if isRef is provided must be a boolean array of the same shape as 'scan' - is there is at least one scanpoint marked as True, the data are - subtracted/divided by the interpolated reference - if lpower is provided the data is divided by it (how it is done depends - if one uses the ratio or not - funcForAveraging: is usually np.nanmean or np.nanmedian. it can be any - function that support axis=0 as keyword argument -""" - data = data.astype(np.float) - avData = np.nanmedian( data , axis = 0 ) - - if isRef is None: isRef = np.zeros( data.shape[0], dtype=bool ) - assert data.shape[0] == isRef.shape[0] - - # subtract reference only is there is at least one - if isRef.sum()>0: - # create a copy (subtractReferences works in place) - diff = subtractReferences(data.copy(),np.argwhere(isRef), useRatio=useRatio) - avNeg = funcForAveraging(data[isRef],axis=0) - else: - diff = data - avNeg = np.zeros_like(avData) - - # normalize signal for laser intensity if provided - if lpower is not None: - lpower = utils.reshapeToBroadcast(lpower,data) - if useRatio is False: - diff /= lpower - else: - diff = (data-1)/lpower+1 - - scan_pos = np.unique(scan) - shape_out = [len(scan_pos),] + list(diff.shape[1:]) - ret = np.empty(shape_out) - err = np.empty(shape_out) - data_abs = np.empty(shape_out) - diffsInScanPoint = [] - chi2_0 = [] - for i,t in enumerate(scan_pos): - shot_idx = (scan == t) - - # select data for the scan point - diff_for_scan = diff[shot_idx] - #if errAbs is not None: - # noise = np.nanmean(errAbs[shot_idx],axis = 0) - #else: - noise = np.nanstd(diff_for_scan, axis = 0) - - # if it is the reference take only every second ... - if np.all( shot_idx == isRef ): - diff_for_scan = diff_for_scan[::2] - - diffsInScanPoint.append( diff_for_scan ) - - # calculate average - ret[i] = funcForAveraging(diff_for_scan,axis=0) - data_abs[i] = funcForAveraging(data[shot_idx],axis=0) - - # calculate chi2 of different repetitions - chi2 = np.power( (diff_for_scan - ret[i])/noise,2) - # sum over all axis but first - for _ in range(diff_for_scan.ndim-1): - chi2 = np.nansum( chi2, axis=-1 ) - - # store chi2_0 - chi2_0.append( chi2/ret[i].size ) - - # store error of mean - err[i] = noise/np.sqrt(shot_idx.sum()) - ret = dict(scan=scan_pos,data=ret,dataUnmasked=ret.copy(),err=err, - errUnmasked=err.copy(),chi2_0=chi2_0,diffsInScanPoint=diffsInScanPoint, - dataAbsAvNeg = avNeg, dataAsAbs=ret+avNeg,errAbs=errAbs, - dataAbsAvAll=avData,dataAbsAvScanPoint=data_abs,dataAbs=data.copy()) - ret = storage.DataStorage(ret) - return ret - - -def calcTimeResolvedSignal(scan,data,err=None,reference="min",monitor=None,q=None, - saveTxt=True,folder="./",**kw): - """ - reference: can be 'min', 'max', a float|integer or an array of booleans - q : is needed if monitor is a tuple|list - monitor : normalization vector (if it is interpreted a list it is - interpreted as q-range normalization) - saveTxt : will save txt outputfiles (diff_av_*) in folder - other keywords are passed to averageScanPoints - """ - if reference == "min": - isRef = (scan == scan.min()) - elif reference == "max": - isRef = (scan == scan.max()) - elif isinstance(reference,(float,int)): - isRef = (scan == reference) - else: - isRef = reference - # normalize if needed - if monitor is not None: - if isinstance(monitor,(tuple,list)): - assert q is not None - assert data.ndim == 2 - idx = (q>= monitor[0]) & (q<= monitor[1]) - monitor = np.nanmedian(data[:,idx],axis=1) - monitor = utils.reshapeToBroadcast(monitor,data) - data = data/monitor - if err is not None: err = err/monitor - ret = averageScanPoints(scan,data,errAbs=err,isRef=isRef,**kw) - if q is not None: ret["q"] = q - return ret - -def saveTxt(folder,data,delayToStr=True,basename='auto',info="",**kw): - """ data must be a DataStorage instance """ - # folder ends usually with sample/run so use the last two subfolders - # the abspath is needed in case we analyze the "./" - folder = os.path.abspath(folder); - if basename == 'auto': - basename = "_".join(folder.rstrip("/").split("/")[-2:]) + "_" - q = data.q if "q" in data else np.arange(data.data.shape[-1]) - # save one file with all average diffs - fname = "%s/%sdiff_av_matrix.txt" % (folder,basename) - utils.saveTxt(fname,q,data.data,headerv=data.scan,**kw) - # save error bars in the matrix form - fname = "%s/%sdiff_av_matrix_err.txt" % (folder,basename) - utils.saveTxt(fname,q,data.err,headerv=data.scan,**kw) - - for iscan,scan in enumerate(data.scan): - scan = utils.timeToStr(scan) if delayToStr else "%+10.5e" % scan - - # try retreiving info on chi2 - try: - chi2_0 = data.chi2_0[iscan] - info_delay = [ "# rep_num : chi2_0 , discarded by chi2masking ?", ] - for irep,value in enumerate(chi2_0): - info_delay.append( "# %d : %.3f" % (irep,value)) - if 'chi2' in data.masks: info_delay[-1] += " %s"%str(data.masks['chi2'][iscan][irep]) - info_delay = "\n".join(info_delay) - if info != '': info_delay = "%s\n%s" % (info,info_delay) - except AttributeError: - info_delay = info - - # save one file per timedelay with average diff (and err) - fname = "%s/%sdiff_av_%s.txt" % (folder,basename,scan) - if 'mask' in data: - tosave = np.vstack( (data.data[iscan],data.err[iscan], - data.dataUnmasked[iscan],data.errUnmasked[iscan] ) ) - columns = 'q diffmask errmask diffnomask errnomask'.split() - else: - tosave = np.vstack( (data.data[iscan],data.err[iscan] ) ) - columns = 'q diff err'.split() - utils.saveTxt(fname,q,tosave,info=info_delay,columns=columns) - - # save one file per timedelay with all diffs for given delay - fname = "%s/%sdiffs_%s.txt" % (folder,basename,scan) - utils.saveTxt(fname,q,data.diffsInScanPoint[iscan],info=info_delay,**kw) - - -def read_diff_av(folder,plot2D=False,save=None): - print("Never tested !!!") - basename = folder+"/"+"diff_av*" - files = glob.glob(basename) - files.sort() - if len(files) == 0: - print("No file found (basename %s)" % basename) - return None - temp = [os.path.basename(f[:-4]) for f in files] - delays = [f.split("_")[-1] for f in temp ] - diffav = collections.OrderedDict() - diffs = collections.OrderedDict() - for d,f in zip(delays,files): - data = np.loadtxt(f) - diffav[d]=data[:,1] - diffs[d] = np.loadtxt(folder+"/diffs_%s.dat"%d)[:,1:] - q =data[:,0] - t = np.asarray( [mc.strToTime(ti) for ti in delays] ) - if plot2D: - idx = t>0 - i = np.asarray( diffav.values() ) - plt.pcolor(np.log10(t[idx]),q,i[idx].T) - plt.xlabel(r"$\log_{10}(t)$") - plt.ylabel(r"q ($\AA^{-1}$)") - it=np.asarray(diffav.values()) - if save: - tosave = np.vstack( (q,it) ) - header = np.hstack( (len(it),t) ) - tosave = np.vstack( (header,tosave.T) ) - np.savetxt(folder + "/all_diffs_av_matrix.txt",tosave) - return q,it,diffs,t diff --git a/xray/example_main_salen.py b/xray/example_main_salen.py deleted file mode 100644 index afabf5e..0000000 --- a/xray/example_main_salen.py +++ /dev/null @@ -1,356 +0,0 @@ -from __future__ import print_function,division -import os -import collections - -# prepare logging -import logging -logfname = os.path.splitext(__file__)[0] + ".log" -logging.basicConfig(level=logging.DEBUG, - format='%(asctime)s %(name)-12s %(levelname)-8s %(message)s', - datefmt='%y-%m-%d %H:%M:%S', - filename=logfname, - filemode='w') -import numpy as np -import pylab as plt -import mcutils as mc -import mcutils.xray as xray -from mcutils.xray import id9 -#id9 = xray.id9 - -# use npz files (they can handle more stuff (list of arrays,unicode) than h5py) -g_default_ext = '.npz' - -id9.default_extension = g_default_ext -g_default_mask = '../masks/fesalen1_run8_careful.edf' -g_default_folder = "../fesalen/fesalen1/run%d" - -runsE = collections.OrderedDict() -runsT = collections.OrderedDict() -runsE[7] = 60; runsT[7] = 135 -runsE[8] = 60; runsT[8] = 135 -runsE[9] = 120; runsT[9] = 135 -runsE[10] = 240; runsT[10] = 135 -runsE[11] = 30; runsT[11] = 135 -runsE[12] = 45; runsT[12] = 135 -runsE[13] = 60; runsT[13] = 120 -runsE[14] = 120; runsT[14] = 120 -runsE[15] = 30; runsT[15] = 120 - - -runs135=(11,12,8,9,10) -runs120=(15,13,14) - -def prepareLog(): - """ It allows printing to terminal on top of logfile """ - # define a Handler which writes INFO messages or higher to the sys.stderr - console = logging.StreamHandler() - console.setLevel(logging.INFO) - # set a format which is simpler for console use - formatter = logging.Formatter('%(message)s') - # tell the handler to use this format - console.setFormatter(formatter) - # add the handler to the root logger (if needed) - if len(logging.getLogger('').handlers)==1: - logging.getLogger('').addHandler(console) - -prepareLog() - -def defineCalibrants(dmin=3): - c120 = xray.cell.Triclinic(a=10.297,b=10.716,c=13.955,alpha=72.22,beta=70.45,gamma=74.13) - p120,_ = c120.reflection_list() - c300 = xray.cell.Triclinic(a=10.485,b=11.015,c=14.280,alpha=74.61,beta=69.29,gamma=75.29) - p300,_ = c300.reflection_list() - return xray.storage.DataStorage( p120=p120, c120=c120,p300=p300,c300=c300) - -def readCalculations(folder="../fesalen/cif",force=False,plot=False): - npz = folder + "/fesalen_calc.npz" - if not os.path.isfile(npz) or force: - ret = dict() - for T in ("120LS","200HS","300HS"): - # calculated with lambda = 1 Ang - theta,i = np.loadtxt( folder+"/fesalen_%s_fwhm0.1.tsv"%T,unpack=True ) - theta = theta/180*np.pi - q = 4*np.pi/1*np.sin(theta/2) - key = "T%s"%T - ret[ key ] = dict( q = q, i = i ) - ret = xray.storage.DataStorage(ret) - ret.save(folder + "/fesalen_calc.npz" ) - ret = xray.storage.read(npz) - if plot: - T = list(ret.keys()) - T.sort() - for i,t in enumerate(T): - plt.plot( ret[t].q,ret[t].i, label=t)#,color=plt.cm.Blues((0.3+0.7*i/2)) ) - plt.grid() - plt.legend() - plt.xlabel(r"q ($\AA^{-1}$)") - return ret - - -def findCenter(cen=(484,479.5),dist=0.251): - files = xray.utils.getEdfFiles("../fesalen/fesalen1/run17/",nFiles=100) - img = xray.azav.read(files).mean(axis=0) - 10. - wavelength=12.398/18.*1e-10 - # 1,0,0 peak are the one close to azimuth 0 and are at 0.6596Ang-1 at 120K - # and 0.65098 at 300K (run 17 is @ RT) - xray.azav.find_center(img,dist=dist,psize=0.0001770834,center=cen, - reference=0.65098,wavelength=wavelength) - - - - -def removeBaseline(folder,qlims=(0.5,2.5),max_iter=30,force=False): - if isinstance(folder,int): folder = g_default_folder%folder - fname = folder + "/pyfai_1d_nobaseline"+g_default_ext - if not os.path.isfile(fname) or force: - logging.info("Calling azav with default parameters") - data = azav(folder,withBaseline=True) - # first peak (around 0.6) is weak and requires special care ... - data.q,data.data = xray.utils.removeBackground(data.q,data.data, - xlims=qlims,max_iter=max_iter,background_regions=[ (0,0.58) ]) - data.save(fname) - else: - data = xray.storage.read(fname) - return data - -def azav(folder,nQ=3000,force=False,saveChi=True,mask=g_default_mask, - withBaseline=True): - if isinstance(folder,int): folder = g_default_folder%folder - if withBaseline: - data = id9.doFolder_azav(folder,nQ=nQ,force=force,mask=mask, - saveChi=saveChi) - else: - data=removeBaseline(folder) - return data - -def datared(folder,monitor=(0.5,4),showPlot=True,errMask=5,chi2Mask=1.5, - qlims=(0.5,3),withBaseline=False,forceDatared=False, - select=slice(None),**kw): - if isinstance(folder,int): folder = g_default_folder%folder - data = azav(folder,withBaseline=withBaseline) - if withBaseline: - diffFile = folder + "/diffs" + g_default_ext - else: - diffFile = folder + "/diffs_nobaseline" + g_default_ext - - if not os.path.isfile(diffFile) or forceDatared: - data,diffs = id9.doFolder_dataRed(data,monitor=monitor, - errMask=errMask,chi2Mask=chi2Mask,qlims=qlims, - outStorageFile=diffFile) - else: - diffs = xray.storage.read(diffFile) - - if showPlot: - xray.utils.plotdiffs(diffs.q,diffs.data,t=diffs.scan,select=select, - absSignal=diffs.dataAbsAvNeg,absSignalScale=10) - plt.title(folder + " norm %s" % str(monitor)) - plt.figure() - xray.utils.plotdiffs(diffs.q,diffs.dataAsAbs,t=diffs.scan, - select=select) - plt.title(folder + " norm %s" % str(monitor)) - return data,diffs - - -def doall(folder,force=False,removeBaseline=True): - azav(folder,force=force) - return datared(folder) - -def _getAmplitudeFromSVD(i,comp=0): - u,s,v = np.linalg.svd(i) - return u[:,comp]*s[comp] - -def _getAmplitudeFromLinFit(i,ref=1): - amp = [mc.linFit(i[ref],ii) for ii in i] - return np.asarray(amp) - -def _getAmplitudeFromAbs(i,e): - amp = np.abs(i).mean(axis=1) - baseline = np.median(e) - return amp-baseline - -def _getAmplitudeFromPeak(i): - """ get max-min (to be less sensitive to baseline change """ - amp = np.max(i,axis=1)-np.min(i,axis=1) - return amp - -def anaAmplitue(run=12,plot=False,scaleWithE=False,withBaseline=False): - folder = g_default_folder % run - pars = dict( chi2Mask=1.5,errMask=3,monitor=(0.4,2.5),showPlot=False ) - data,diffs=datared(folder,**pars,withBaseline=withBaseline) - ranges = ( (0.66,0.72), (0.74,0.8), (0.82,0.88) ) #, (0.96,1.02 ) ) - nPlot = len(ranges) - if plot: fig,ax = plt.subplots(nPlot,1,sharex=True) - amp = collections.OrderedDict() - for num,r in enumerate(ranges): - idx = (diffs.q>r[0]) & (diffs.qr[0]) & (diffs.q470'): - return id9.doFolder_azav(folder,nQ=nQ,force=force,mask=mask,saveChi=saveChi) - -def datared(folder,monitor=(1,5),forceDatared=False,showPlot=True,**kw): - azavFile = folder + "/pyfai_1d" + g_default_ext - diffFile = folder + "/diffs" + g_default_ext - if not os.path.isfile(diffFile) or forceDatared: - data,diffs = id9.doFolder_dataRed(folder,monitor=monitor,**kw) - else: - data = xray.storage.read(azavFile) - diffs = xray.storage.read(diffFile) - if showPlot: - xray.utils.plotdiffs(diffs.q,diffs.data,t=diffs.scan, - absSignal=diffs.dataAbsAvAll,absSignalScale=30) - plt.plot(calc.lam.q,calc.lam.i/1000+0.1,label='calc') - plt.title(folder + " norm %s" % str(monitor)) - return data,diffs - -def doall(folder,force=False): - azav(folder,force=force) - return datared(folder,forceDatared=force) - -def anaAmplitue(run=6): - fname = "../tiox/tiox1/run%d/diffs.npz" % run - data = xray.storage.DataStorage(fname) - ranges = ( (1.75,1.85), (2.2,2.4), (3.25,3.4) ) - nPlot = len(ranges) - fig,ax = plt.subplots(nPlot,1,sharex=True) - for r,a in zip(ranges,ax): - idx = (data.q>r[0]) & (data.q threshold - log.debug("Removed %d zingers from %d curves"%(idx.sum(),len(curves))) - print("Removed %d zingers from %d curves"%(idx.sum(),len(curves))) - if idx.sum()>0: - curves[idx]=np.nan - #curves = np.ma.MaskedArray(data=curves,mask=idx) - return curves - -def filterOutlier(curves,errs=None,norm=None,threshold=10): - # normalize - if norm == 'auto': - norm = np.nanmean(curves,axis=1) - norm = utils.reshapeToBroadcast(n,curves) - elif norm is None: - norm = 1 - - curves = curves/norm - if errs is None: - errs = statsmodels.robust.mad(curves,axis=0) - else: - errs = errs/norm - - median = np.median(curves) - diff = np.abs(curves-median)/errs - chi2 = np.sum(diff**2)/len(curves) - idx = chi2 < threshold - return curves[idx] - -def chi2Filter(diffs,threshold=10): - """ Contrary to removeZingers, this removes entire curves """ - idx_mask = [] - for iscan in range(len(diffs.diffsInScanPoint)): - idx = diffs.chi2_0[iscan] > threshold - # expand along other axis (q ...) - #idx = utils.reshapeToBroadcast(idx,data.diffsInScanPoint[iscan]) - idx_mask.append(idx) - log.debug("Chi2 mask, scanpoint: %s, curves filtereout out %d/%d (%.2f%%)"%\ - (data.scan[iscan],idx.sum(),len(idx),idx.sum()/len(idx)*100) ) - print("Chi2 mask, scanpoint: %s, curves filtereout out %d/%d (%.2f%%)"%\ - (data.scan[iscan],idx.sum(),len(idx),idx.sum()/len(idx)*100) ) - - if "masks" not in data: data['masks'] = dict() - if "masks_pars" not in data: data['masks_pars'] = dict() - data['masks']['chi2'] = idx_mask - data['masks_pars']['chi2_threshold'] = threshold - return data - diff --git a/xray/id9.py b/xray/id9.py deleted file mode 100644 index 221b6ed..0000000 --- a/xray/id9.py +++ /dev/null @@ -1,212 +0,0 @@ -import logging -log = logging.getLogger(__name__) - -import os -import collections -import numpy as np -from . import azav -from . import dataReduction -from . import utils -from . import storage -from . import filters - -default_extension = ".h5" - -def _conv(x): - try: - x = float(x) - except: - x = np.nan - return x - -def _readDiagnostic(fname,retry=3): - ntry = 0 - while ntry1: log.warn("Found more than one *.log file that is not diagnostics.log: %s"%files) - return logfile - -def readLogFile(fnameOrFolder,subtractDark=False,skip_first=0,asDataStorage=True,last=None): - """ read id9 style logfile """ - if os.path.isdir(fnameOrFolder): - fname = findLogFile(fnameOrFolder) - else: - fname = fnameOrFolder - f = open(fname,"r") - lines = f.readlines() - f.close() - lines = [line.strip() for line in lines] - darks = {} - for line in lines: - if line.find("pd1 dark/sec")>=0: darks['pd1ic'] = _findDark(line) - if line.find("pd2 dark/sec")>=0: darks['pd2ic'] = _findDark(line) - if line.find("pd3 dark/sec")>=0: darks['pd3ic'] = _findDark(line) - if line.find("pd4 dark/sec")>=0: darks['pd4ic'] = _findDark(line) - for iline,line in enumerate(lines): - if line.lstrip()[0] != "#": break - data=np.genfromtxt(fname,skip_header=iline-1,names=True,comments="%",dtype=None,converters = {'delay': lambda s: _delayToNum(s)}) - if subtractDark: - for diode in ['pd1ic','pd2ic','pd3ic','pd4ic']: - if diode in darks: data[diode]=data[diode]-darks[diode]*data['timeic'] - data = data[skip_first:last] - if asDataStorage: - # rstrip("_") is used to clean up last _ that appera for some reason in file_ - data = storage.DataStorage( dict((name.rstrip("_"),data[name]) for name in data.dtype.names ) ) - data.file = data.file.astype(str) - return data - - - -def doFolder_azav(folder,nQ=1500,files='*.edf*',force=False,mask=None, - saveChi=True,poni='pyfai.poni',storageFile='auto',dark=9.9,dezinger=None, - qlims=(0,10),removeBack=False,removeBack_kw=dict(),skip_first=0, - last=None): - """ very small wrapper around azav.doFolder, essentially just reading - the id9 logfile or diagnostics.log """ - - try: - loginfo = readLogFile(folder,skip_first=skip_first,last=last) - except Exception as e: - log.warn("Could not read log file, trying to read diagnostics.log") - loginfo = readDiagnostic(folder) - if storageFile == 'auto' : storageFile = folder + "/" + "pyfai_1d" + default_extension - - data = azav.doFolder(folder,files=files,nQ=nQ,force=force,mask=mask, - saveChi=saveChi,poni=poni,storageFile=storageFile,logDict=loginfo,dark=dark,save=False,dezinger=dezinger) - #try: - # if removeBack is not None: - # _,data.data = azav.removeBackground(data,qlims=qlims,**removeBack_kw) - #except Exception as e: - # log.error("Could not remove background, error was %s"%(str(e))) - -# idx = utils.findSlice(data.q,qlims) -# n = np.nanmean(data.data[:,idx],axis=1) -# data.norm_range = qlims -# data.norm = n -# n = utils.reshapeToBroadcast(n,data.data) -# data.data_norm = data.data/n - - data.save(storageFile) - - - return data - - - -def doFolder_dataRed(azavStorage,monitor=None,funcForAveraging=np.nanmean, - qlims=None,outStorageFile='auto',reference='min'): - """ azavStorage if a DataStorage instance or the filename to read - monitor : normalization vector that can be given as - 1. numpy array - 2. a list (interpreted as q-range normalization) - 3. a string to look for as key in the log, e.g. - monitor="pd2ic" would reult in using - azavStorage.log.pd2ic - - """ - - if isinstance(azavStorage,storage.DataStorage): - data = azavStorage - folder = azavStorage.folder - elif os.path.isfile(azavStorage): - folder = os.path.dirname(azavStorage) - data = storage.DataStorage(azavStorage) - else: - # assume is just a folder name - folder = azavStorage - azavStorage = folder + "/pyfai_1d" + default_extension - data = storage.DataStorage(azavStorage) - - #assert data.q.shape[0] == data.data.shape[1] == data.err.shape[1] - if qlims is not None: - idx = (data.q>qlims[0]) & (data.qx2: x1,x2=x2,x1 - if y1>y2: y1,y2=y2,y1 - return (X>x1) & (Xy1) & (Y shape[1]-snapRange: snapped[0] = shape[1] - if snapped[1] < snapRange: snapped[1] = 0 - if snapped[1] > shape[0]-snapRange: snapped[1] = shape[0] - return tuple(snapped) - -def getPoints(N=1,shape=(100,100),snapRange=0): - if N<1: print('Right click cancels last point, middle click ends the polygon') - c = plt.ginput(N) - c = [ snap(point,shape,snapRange=snapRange) for point in c ] - if len(c) == 1: c = c[0] - return c - -def makeMaskGui(img,snapRange=60): - """ snapRange controls border snapping (in pixels, use <= 0 to disable """ - mask = MyMask(img) - ans='ok' - while (ans != 'done'): - plt.imshow(img) - plt.clim(np.percentile(img,(2,98))) - plt.imshow(mask.getMatplotlibMask()) - plt.pause(0.01) - ans = input("What's next p/P/c/C/r/R/done? (capitals = subtract)") - if ans == "c": - print("Adding circle, click on center then another point to define radius") - vertices = getPoints(N=2,shape=img.shape,snapRange=snapRange) - mask.addCircle(*vertices) - if ans == "C": - print("Subtracting circle, click on center then another point to define radius") - vertices = getPoints(N=2,shape=img.shape,snapRange=snapRange) - mask.subtractCircle(*vertices) - if ans == "r": - print("Adding rectangle, click on one corner and then on the opposite one") - vertices = getPoints(N=2,shape=img.shape,snapRange=snapRange) - mask.addRectangle(*vertices) - if ans == "R": - print("Subtracting rectangle, click on one corner and then on the opposite one") - vertices = getPoints(N=2,shape=img.shape,snapRange=snapRange) - mask.subtractRectangle(*vertices) - if ans == 'p': - print("Adding polygon") - vertices = getPoints(N=-1,shape=img.shape,snapRange=snapRange) - mask.addPolygon(*vertices) - if ans == 'P': - print("Subtracting polygon") - vertices = getPoints(N=-1,shape=img.shape,snapRange=snapRange) - mask.subtractPolygon(*vertices) - - plt.imshow(mask.getMatplotlibMask()) - plt.pause(0.01) - fname = input("Enter a valid filename (ext .edf or .npy) if you want to save the mask (empty otherwise)") - try: - if fname != '': - ext = os.path.splitext(fname)[1] - if ext == '.edf': - mask.save(fname) - elif ext == '.npy': - np.save(fname,mask.getMask()) - except Exception as e: - log.error("Error in saving mask") - log.error(e) - finally: - return mask - -def maskBorder(width,shape): - mask = np.zeros(shape,dtype=bool) - mask[ :width , : ] = True - mask[ -width: , : ] = True - mask[ : , :width ] = True - mask[ : , -width: ] = True - return mask - -def maskCenterLines(width,shape): - mask = np.zeros(shape,dtype=bool) - if isinstance(width,int): width = (width,width) - c0 = int(shape[0]/2) - c1 = int(shape[1]/2) - w0 = int(width[0]/2) - w1 = int(width[1]/2) - mask[ c0-w0:c0+w0 , : ] = True - mask[ : , c1-w1:c1+w1 ] = True - return mask - - -def test(shape=(1000,2000)): - mask = MyMask() - mask.addCircle(400,300,250) - mask.subtractCircle(400,300,150) - mask.addRectangle(350,250,1500,700) - plt.imshow( mask.getMask(shape) ) - return mask - - - -if __name__ == "__main__": - test() - plt.show() - ans=input("Enter to finish") diff --git a/xray/peaks.py b/xray/peaks.py deleted file mode 100644 index a1f42ad..0000000 --- a/xray/peaks.py +++ /dev/null @@ -1,32 +0,0 @@ -import lmfit -import numpy as np -import logging as log - -pv = lmfit.models.PseudoVoigtModel() - -def fitPeak(x,y,err=1,autorange=False): - if isinstance(err,np.ndarray): - if np.all(err==0): - err = 1 - log.warn("Asked to fit peak but all errors are zero, forcing them to 1") - elif np.isfinite(err).sum() != len(err): - idx = np.isfinite(err) - x = x[idx] - y = y[idx] - err = err[idx] - log.warn("Asked to fit peak but some errorbars are infinite or nans,\ - excluding those points") - if autorange: - # find fwhm - idx = np.ravel(np.argwhere( ynp.argmax(y)][0] - c = int( (p1+p2)/2 ) - dp = int( np.abs(p1-p2) ) - idx = slice(c-dp,c+dp) - x = x[idx] - y = y[idx] - pars = pv.guess(y,x=x) - ret = pv.fit(y,x=x,weights=1/err,params=pars) - return ret diff --git a/xray/utils.py b/xray/utils.py deleted file mode 100644 index d1f2efa..0000000 --- a/xray/utils.py +++ /dev/null @@ -1,346 +0,0 @@ -from __future__ import print_function,division - -import logging -log = logging.getLogger(__name__) # __name__ is "foo.bar" here - -import numpy as np -np.seterr(all='ignore') -import os -import glob -import pathlib -import re -import numbers -from . import storage as storage - -try: - import matplotlib.pyplot as plt -except ImportError: - log.warn("Can't import matplotlib !") - -_time_regex = re.compile( "(-?\d+\.?\d*(?:ps|ns|us|ms)?)") -_timeInStr_regex = re.compile("_(-?\d+\.?\d*(?:ps|ns|us|ms)?)") - -def getFiles(folder,basename="*.edf*",nFiles=None): - files = glob.glob(folder + "/" + basename) - files.sort() - if nFiles is not None: files = files[:nFiles] - return files - -def getEdfFiles(folder,**kw): - return getFiles(folder,basename="*.edf*",**kw) - -def getDelayFromString(string) : - match = _timeInStr_regex.search(string) - return match and match.group(1) or None - -_time_regex = re.compile("(-?\d+\.?\d*)((?:s|fs|ms|ns|ps|us)?)") - -def strToTime(delay) : - if isinstance(delay,bytes): delay = delay.decode('ascii') - _time2value = dict( fs = 1e-15, ps = 1e-12, ns = 1e-9, us = 1e-6, ms = 1e-3, s = 1) - - match = _time_regex.search(delay) - if match: - n,t = float(match.group(1)),match.group(2) - value = _time2value.get(t,1) - return n*value - else: - return None - -def timeToStr(delay,fmt="%+.0f"): - a_delay = abs(delay) - if a_delay >= 1: - ret = fmt % delay + "s" - elif 1e-3 <= a_delay < 1: - ret = fmt % (delay*1e3) + "ms" - elif 1e-6 <= a_delay < 1e-3: - ret = fmt % (delay*1e6) + "us" - elif 1e-9 <= a_delay < 1e-6: - ret = fmt % (delay*1e9) + "ns" - elif 1e-12 <= a_delay < 1e-9: - ret = fmt % (delay*1e12) + "ps" - elif 1e-15 <= a_delay < 1e-12: - ret = fmt % (delay*1e12) + "fs" - elif 1e-18 <= a_delay < 1e-15: - ret = fmt % (delay*1e12) + "as" - else: - ret = str(delay) +"s" - return ret - -def removeExt(fname): - """ special remove extension meant to work with compressed files.edf and .edf.gz files """ - if fname[-3:] == ".gz": fname = fname[:-3] - return os.path.splitext(fname)[0] - -def getBasename(fname): - return os.path.basename(removeExt(fname)); - -def findSlice(array,lims): - start = np.ravel(np.argwhere(array>lims[0]))[0] - stop = np.ravel(np.argwhere(array from mM to M - # 1e3 -> from L to cm3 - # so 1 mM = 1e-3M/L = 1e-6 M/cm3 - density_g_cm3 = density*1e-6*atom.mass - n = xsf.index_of_refraction(atom,density=density_g_cm3,wavelength=wavelength) - attenuation_length = xrayAttLenght( atom,density=density_g_cm3,energy=energy ) - # um-> m: 1e-6 - fraction_absorbed = 1.-np.exp(-length*1e-6/attenuation_length) - if verbose: - print("Fraction of x-ray photons absorbed by the sample:",fraction_absorbed) - ## detector ## - det_area = np.pi*det_radius**2 - det_solid_angle = det_area/(4*np.pi*det_dist**2) - if verbose: - print("Detector fraction of solid angle:",det_solid_angle) - det_att_len = xrayAttLenght(det_material,wavelength=atom.K_alpha) - det_abs = 1-np.exp(-det_thick*1e-6/det_att_len) - if verbose: - print("Detector efficiency (assuming E fluo = E_K_alpha):",det_abs) - eff = fraction_absorbed*det_solid_angle*det_abs - if verbose: - print("Overall intensity (as ratio to incoming one):",eff) - return eff