XRS_swissknife.py 154 KB
Newer Older
1
2
3
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
4
import scipy
myron's avatar
myron committed
5
import json
6
### __doc__="""
7

8
9
10

generality_doc = 1
"""
Alessandro MIRONE's avatar
Alessandro MIRONE committed
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
Documentation of XRS_swissknife
-------------------------------

The script is called in this way ::

  XRS_swissknife yourinput.yaml

The input file is in *yaml* format. In *yaml* format each line introduces an item
and the indentation expresses the hierarchy.
An example is ::

  Fancy_Reduction :
        parameterTom : 3.14
        parameterJerry : False
        choicesBob     : [1,2,3]

In this example we create an item called *Fancy_Reduction* which contains three subitems.

The *XRS_swissknife* expects that for each operation that you want to apply, you provide an item
named as the operation, and the associated subitems that provide that values for the parameters.

*XRS_swissknife* will do what you want provided you respect the proper indentation. A thing which helps
is using emacs and activating the *python mode*, because python uses the same indentation principle
to structure the code. 

Each processing item has the additional, optional, key **active**.
This key can be set to **0** or **1** to desactivate or not (default is **1**, active) the processing.
Here a desactivation example ::

  Fancy_Reduction :
        active : 0
        parameterTom : 3.14
        parameterJerry : False
        choicesBob     : [1,2,3]

The following documentation has been created automatically, for each functionality, based on the documentation string
written in the code for the fonctionality itself.  You can also write mathematical expressions :  :math:`\\int{ x dx} = \\frac { x^2 }{ 2 }`
and even include graphics.

"""
51
import collections
52

53
54
55
try:
    from mayavi import mlab
except:
Alessandro Mirone's avatar
Alessandro Mirone committed
56
    print( " WAS not able to load mayavi, some feature might be missing ")
57

myron's avatar
myron committed
58
import string
59

60
import numpy as np
61
import math
62

myron's avatar
myron committed
63
import yaml
64
65
from yaml import load, dump

66
import numbers
alex mirone's avatar
alex mirone committed
67
68
69
import re
import yaml
import yaml.resolver
Alessandro Mirone's avatar
Alessandro Mirone committed
70
71
import PyMca5.PyMcaIO.specfilewrapper as SpecIO
import fabio
72
from six import u
Alessandro Mirone's avatar
Alessandro Mirone committed
73

myron's avatar
myron committed
74
75
76

from silx.gui  import qt as Qt
## from  PyQt4 import Qt, QtCore
myron's avatar
myron committed
77
78
from XRStools.roiNmaSelectionGui   import roiNmaSelectionWidget
from XRStools import roiSelectionWidget
79
80
81
82

import h5py
import sys

alex mirone's avatar
alex mirone committed
83
84
85
86
yaml.resolver.Resolver
Resolver = yaml.resolver.Resolver
Resolver.add_implicit_resolver(
        u'tag:yaml.org,2002:float',
87
        re.compile(u(r"""^(?:[-+]?(?:[0-9][0-9_]*)(\.[0-9_]*)?(?:[eE][-+]?[0-9]+)?
alex mirone's avatar
alex mirone committed
88
89
90
                    |\.[0-9_]+(?:[eE][-+][0-9]+)?
                    |[-+]?[0-9][0-9_]*(?::[0-5]?[0-9])+\.[0-9_]*
                    |[-+]?\.(?:inf|Inf|INF)
91
                    |\.(?:nan|NaN|NAN))$"""), re.X),
alex mirone's avatar
alex mirone committed
92
93
        list(u'-+0123456789.'))

94
95
96
97
98
99
try:
    from mpi4py import MPI
    myrank = MPI.COMM_WORLD.Get_rank()
    nprocs = MPI.COMM_WORLD.Get_size()
    procnm = MPI.Get_processor_name()
    comm = MPI.COMM_WORLD
Alessandro Mirone's avatar
Alessandro Mirone committed
100
    print( "MPI LOADED , nprocs = ", nprocs)
101
102
103
104
except:
    class FakeComm:
        def Barrier(self):
            pass
105
        def allreduce(self,number, operation):
106
            assert ( isinstance(number, numbers.Number)    )
107
            return number
108
109
110
        def Get_size(self):
            return 1

111
112
    myrank=0
    nprocs = 1
113
    comm = FakeComm()
Alessandro Mirone's avatar
Alessandro Mirone committed
114
    print( "no MPI LOADED , nprocs = ", nprocs)
115

116
117
118
119
120
121
122
123


def stripROI(t):
    if "ROI" in t:
        return t[3:]
    else:
        return t

124
def filterRoiList(l, strip=False, prefix="ROI"):
125
126
    
    if strip:
127
        result =  [ str(int(stripROI(t))) for t in l if ( t not in [ "motorDict", "response"] and t[:len(prefix)]==prefix )]
128
    else:
129
        result =  [ t for t in l  if ( t not in [ "motorDict", "response"] and  t[:len(prefix)]==prefix ) ]
myron's avatar
myron committed
130
131
132
133
134


    result  =      sorted(  result     , key = lambda x:  int(''.join(filter(str.isdigit, str(x) ))) )
    ### result.sort()
    
135
    return result
136
137
138

def filterScanList(l, prefix="Scan"):
    result =  [ t for t in l  if ( t not in [ "motorDict", "response"] and t[:len(prefix)]==prefix ) ]
myron's avatar
myron committed
139
140
141
142

    result = sorted(  result     , key = lambda x:  int(''.join(filter(str.isdigit, str(x) ))) )
    # result.sort()
    
143
    return result
144
    
145

146
147
148
def checkNoParallel( routineName):
    if nprocs>1:
        msg = " ERROR : %s feature not yet parallel : relaunch it with 1 process only "
Alessandro Mirone's avatar
Alessandro Mirone committed
149
        print( msg)
150
        raise Exception( msg)
151

Alessandro Mirone's avatar
Alessandro Mirone committed
152
153
154
155
156
# try:
#     from yaml import CLoader as Loader, CDumper as Dumper
# except ImportError:
#     from yaml import Loader, Dumper            

157

158
159
160
161
162
163
164
165
166
167
168
169
170
171


nprocs = comm.Get_size()

if nprocs>1:
    # circumvent issus with mpi4py not stopping mpirun
    def excepthook(type, value, traceback):
        res = sys.__excepthook__(type, value, traceback)
        comm.Abort(1)
        return res
    sys.excepthook=excepthook



172
import os
myron's avatar
myron committed
173
174
175
176
177
178
179
from XRStools import xrs_rois
from XRStools import roifinder_and_gui
from XRStools import xrs_scans
from XRStools import xrs_read
from XRStools import rixs_read
from XRStools import theory
from XRStools import extraction
180

myron's avatar
myron committed
181
182
from XRStools import xrs_prediction
from XRStools import xrs_imaging
183
184


myron's avatar
myron committed
185
186
from XRStools import xrs_imaging
from XRStools import superr
187

Alessandro Mirone's avatar
Alessandro Mirone committed
188
189
190
191
192
193
194
195
196
#################################################################
##  THIS redefinition of yaml is used to keep the entry ordering
## when accessing yamlData keys
##
#yaml_anydict.py
import yaml
from yaml.representer import Representer
from yaml.constructor import Constructor, MappingNode, ConstructorError

myron's avatar
myron committed
197
198
from XRStools import fit_spectra
from XRStools import reponse_percussionelle
199

200
201


myron's avatar
myron committed
202
def check_allowed_keys(  mydata, allowed_keys  ) :
203
204
    for k in mydata.keys():
        if not k in allowed_keys:
myron's avatar
myron committed
205
            raise ValueError( (" key "+str(k) +" not in allowed keys :"+ str(allowed_keys))  )
206
    
Alessandro Mirone's avatar
Alessandro Mirone committed
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
def dump_anydict_as_map( anydict):
    yaml.add_representer( anydict, _represent_dictorder)
def _represent_dictorder( self, data):
    return self.represent_mapping('tag:yaml.org,2002:map', data.items() )

class Loader_map_as_anydict( object):
    'inherit + Loader'
    anydict = None      #override
    @classmethod        #and call this
    def load_map_as_anydict( klas):
        yaml.add_constructor( 'tag:yaml.org,2002:map', klas.construct_yaml_map)

    'copied from constructor.BaseConstructor, replacing {} with self.anydict()'
    def construct_mapping(self, node, deep=False):
        if not isinstance(node, MappingNode):
            raise ConstructorError(None, None,
                    "expected a mapping node, but found %s" % node.id,
                    node.start_mark)
        mapping = self.anydict()
        for key_node, value_node in node.value:
            key = self.construct_object(key_node, deep=deep)
            try:
                hash(key)
            except TypeError as exc:
                raise ConstructorError("while constructing a mapping", node.start_mark,
                        "found unacceptable key (%s)" % exc, key_node.start_mark)
            value = self.construct_object(value_node, deep=deep)
            mapping[key] = value
        return mapping

    def construct_yaml_map( self, node):
        data = self.anydict()
        yield data
        value = self.construct_mapping(node)
        data.update(value)


244
245
class myOrderedDict (collections.OrderedDict):
    def __setitem__(self,a,b):
246
        if type(a)==type("") and a in self:
myron's avatar
myron committed
247
            self[a+"_TAG_this_key_is_given_twice"]=b 
248
        else:
249
            ## print super(myOrderedDict, self)
250
251
252
253
            super(myOrderedDict, self).__setitem__(a,b )


def cleaned(key):
myron's avatar
myron committed
254
    while key[-28:]=="_TAG_this_key_is_given_twice":
255
256
257
258
259
260
        key=key[:-8]
    return key
  

            
dictOrder = myOrderedDict
Alessandro Mirone's avatar
Alessandro Mirone committed
261
262
263
264
265

class Loader( Loader_map_as_anydict, yaml.Loader):
    anydict = dictOrder
Loader.load_map_as_anydict()
dump_anydict_as_map( dictOrder)
266
267


Alessandro Mirone's avatar
Alessandro Mirone committed
268
269
270
271
272
273
##
## END of yaml redefinition
###############################

def check_libre( h5 , groupname   ) :

274
275
276
    print("check_libre DESACTIVATED. RESULTS CAN BE OVERWRITTEN")
    return

Alessandro Mirone's avatar
Alessandro Mirone committed
277
278
    if type(h5)==type(""):
        h5f = h5py.File(h5, "r"  )
279
        h5f.close()
Alessandro Mirone's avatar
Alessandro Mirone committed
280
        if groupname in h5f:
281
            msg=(("ERROR: %s key already present in the hdf5 file %s. Erase it before if you dont need it.\n" % (groupname, h5))*10 )
Alessandro Mirone's avatar
Alessandro Mirone committed
282
            print( msg)
283
            raise Exception( msg)
Alessandro Mirone's avatar
Alessandro Mirone committed
284
285
        else:
            if groupname in h5:
286
                msg = (("ERROR: %s key already present in the hdf5 file. Erase it before if you dont need it.\n"%groupname)*10 )
Alessandro Mirone's avatar
Alessandro Mirone committed
287
                print( msg)
288
                raise Exception( msg)
Alessandro Mirone's avatar
Alessandro Mirone committed
289
290
291



Alessandro Mirone's avatar
Alessandro Mirone committed
292
inputtext=""
293

myron's avatar
myron committed
294
295
296
297
298
299
300
301
def yamlFileToDic(fn):
    global  inputtext
    filename = fn
    inputfile = open(filename,"r")
    yamlData = load(inputfile, Loader=Loader)
    return yamlData
    

302
def main():
Alessandro Mirone's avatar
Alessandro Mirone committed
303
    global  inputtext
304
305
    filename = sys.argv[1]

Alessandro Mirone's avatar
Alessandro Mirone committed
306
307
308
    inputfile = open(filename,"r")
    yamlData = load(inputfile, Loader=Loader)
    inputtext = open(filename,"r").read()
309

Alessandro Mirone's avatar
Alessandro Mirone committed
310
    for key in list(yamlData.keys()):
311
      
Alessandro Mirone's avatar
Alessandro Mirone committed
312
        mydata = yamlData[key]
313
        if isinstance(mydata,dict) and    "active"  in mydata :
Alessandro Mirone's avatar
Alessandro Mirone committed
314
315
            if mydata["active"]==0:
                continue
316

317
318
319
        if key != "help":
            swissknife_operations[cleaned(key)](  mydata )
        else:
320
321
322
323

            if cleaned(key)  not in parallelised_operations:
                checkNoParallel( cleaned(key)  )
            
324
            swissknife_operations[cleaned(key)]( yamlData  )
325
326
327
328
329
330
    
    # workbench_file = yamlData["workbench_file"]


def help(yamlData):
    """
Alessandro MIRONE's avatar
Alessandro MIRONE committed
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
    **help**

    Displays doc on the operations. In the input file ::

       help :

    will trigger printing of all the available operation names ::

       help :
           create_rois
           load_scans

    will print the help on *create_rois* and the help about *load_scans*.
    By the way, it is the same that you can read here because the *sphinx* doc-generation tool
    reads the same docstrings contained in the code.
     
347
    """
Alessandro Mirone's avatar
Alessandro Mirone committed
348
    print( " HELP " *15)
349
    if yamlData ["help"] is None:
350
        print( """
351
352
353
354
              Printing all the function names
              To get help on a specific function:

              help : "functionName"
355
        """)
myron's avatar
myron committed
356
        for key,func in swissknife_operations.items():
Alessandro Mirone's avatar
Alessandro Mirone committed
357
            print( " FUNCTION : " , key)
358
359
360

    else:
        func = swissknife_operations[ yamlData ["help"]]
Alessandro Mirone's avatar
Alessandro Mirone committed
361
362
        print( "---------------------------------------")
        print( func.__doc__)
363

Alessandro Mirone's avatar
Alessandro Mirone committed
364

Alessandro MIRONE's avatar
Alessandro MIRONE committed
365

366
367
def split_hdf5_address(dataadress):
 
Alessandro MIRONE's avatar
Alessandro MIRONE committed
368
369
    pos = dataadress.rfind(":")
    if ( pos==-1):
370
        raise Exception( """
371
roiaddress   must be given in the form  roiaddress : "myfile.hdf5:/path/to/hdf5/group"
Alessandro MIRONE's avatar
Alessandro MIRONE committed
372
but : was not found
373
""")
Alessandro Mirone's avatar
Alessandro Mirone committed
374
    filename, groupname = dataadress[:pos], dataadress[pos+1:]
375
376
377
    return filename, groupname 


Alessandro Mirone's avatar
Alessandro Mirone committed
378
    
379
380
381
# def load_scans(mydata):
#     """
#     **load_scans**
Alessandro Mirone's avatar
Alessandro Mirone committed
382

383
384
#     This command harvest the selected signals.
#     the instructions on the scans to be taken must be in the form( as submembers ofload_scans ) ::
Alessandro Mirone's avatar
Alessandro Mirone committed
385
386


387
388
389
#      load_scans :
#          roiaddress :  "hdf5filename:nameofroigroup"  # the same given in create_rois
#          expdata    :  "absolutepathtoaspecfile"  # this points to a spec file
christoph's avatar
christoph committed
390

391
392
393
394
#          elastic_scans    : [623]
#          fine_scans       : [626,630,634,638,642]
#          n_loop           : 4
#          long_scan        : 624
Alessandro Mirone's avatar
Alessandro Mirone committed
395

396
#          signaladdress : "nameofsignalgroup"  # Target group for writing Relative to ROI (and in the same file)!!!!
christoph's avatar
christoph committed
397

398
399
400
401
402
#          #############################################################
#          # OPTIONALS
#          #
#          order : [0,1,2,3,4,5]  #  list of integers (0-5) which describes the order of modules in which the 
#                                 #    ROIs were defined (default is VD, VU, VB, HR, HL, HB; i.e. [0,1,2,3,4,5])
christoph's avatar
christoph committed
403

404
405
406
407
408
409
#          rvd : -41              # mean tth angle of HL module (default is 0.0)
#          rvu : 85               # mean tth angle of HR module (default is 0.0)
#          rvb : 121.8            # mean tth angle of HB module (default is 0.0)
#          rhl : 41.0             # mean tth angle of VD module (default is 0.0)
#          rhr : 41.0             # mean tth angle of VU module (default is 0.0)
#          rhb : 121.8            # mean tth angle of VB module (default is 0.0)
Alessandro Mirone's avatar
Alessandro Mirone committed
410
411


412
413
#     #
#     """
Alessandro Mirone's avatar
Alessandro Mirone committed
414

415
416
#     roiaddress=None
#     roiaddress = mydata["roiaddress"]
Alessandro Mirone's avatar
Alessandro Mirone committed
417

418
419
420
421
422
#     filename, groupname = split_hdf5_address (roiaddress)
#     file= h5py.File(filename,"r")
#     rois = {}
#     shape=xrs_rois.load_rois_fromh5(file[groupname],rois)
#     file.close()
Alessandro Mirone's avatar
Alessandro Mirone committed
423

424
425
#     roiob = xrs_rois.roi_object()
#     roiob.load_rois_fromMasksDict(rois ,  newshape = shape, kind="zoom")
Alessandro Mirone's avatar
Alessandro Mirone committed
426
    
427
#     reader = xrs_read.read_id20(mydata["expdata"] , monitorcolumn='kapraman')
Alessandro Mirone's avatar
Alessandro Mirone committed
428
429


430
#     reader.set_roiObj(roiob)
Alessandro Mirone's avatar
Alessandro Mirone committed
431

432
433
434
435
#     elastic_scans   =  mydata["elastic_scans"][:] 
#     fine_scans      =  mydata["fine_scans"][:] 
#     n_loop          =  mydata["n_loop"] 
#     long_scan       =  mydata["long_scan"]
Alessandro Mirone's avatar
Alessandro Mirone committed
436

437
438
439
440
#     reader.loadelasticdirect(elastic_scans)
#     reader.loadloopdirect(fine_scans,n_loop)
#     print( " LUNGO " )
#     reader.loadlongdirect(long_scan)
Alessandro Mirone's avatar
Alessandro Mirone committed
441
442
443

    

444
#     reader.getspectrum()
Alessandro Mirone's avatar
Alessandro Mirone committed
445

446
447
448
449
450
451
452
453
454
455
#     reader.geteloss()
#     reader.gettths(
#         rvd   = gvord(mydata,"rvd",0.0) ,
#         rvu   = gvord(mydata,"rvu",0.0) ,
#         rvb   = gvord(mydata,"rvb",0.0) ,
#         rhl   = gvord(mydata,"rhl",0.0) ,
#         rhr   = gvord(mydata,"rhr",0.0) ,
#         rhb   = gvord(mydata,"rhb",0.0) ,
#         order = gvord(mydata,"order",        [0,1,2,3,4,5])
#         )
Alessandro Mirone's avatar
Alessandro Mirone committed
456

457
458
#     groupname = groupname+"/"+ mydata["signaladdress"]
#     check_libre( filename , groupname   ) 
Alessandro Mirone's avatar
Alessandro Mirone committed
459

460
#     reader.save_state_hdf5( filename, groupname , comment = inputtext )
Alessandro Mirone's avatar
Alessandro Mirone committed
461

462
463
464
465




466
467
468
# def volume_from_2Dimages(mydata):
#     """
#     imagesaddress :  "test_imaging.hdf5:/ROI_A/images"  # where the data have been saved
469

470
471
#     scan_interval    :  [372,375]                    # OPTIONAL : can be shorter then the scans effectively present in the file
#     roi_n            : 0           # OPTIONAL. if not given, the first non empty found roi. Starts from 0
472

473
#     imagesaddress : "myfile.hdf5:/path/to/hdf5/data"  # OPTIONAL. the target destination for volume. if not given mayavi is launched on the fly.
474

475
#     """
476

Alessandro Mirone's avatar
Alessandro Mirone committed
477

478
#     reader = xrs_imaging.oneD_imaging( "bidon"  , "bidon",  "bidon"  , "bidon")
479

Alessandro Mirone's avatar
Alessandro Mirone committed
480

481
482
#     imagesaddress = mydata["imagesaddress"]
#     filename, groupname = split_hdf5_address(imagesaddress)
Alessandro Mirone's avatar
Alessandro Mirone committed
483

484
#     reader.load_state_hdf5( filename, groupname)
485

486
487
488
#     scan_names = list( reader.twoDimages.keys() )
#     scan_ids = map(int, [name[4:] for name in scan_names  ] )
#     order = np.argsort(scan_ids)
489

490
491
492
493
#     if not  ('scan_interval') in mydata :
#         scan_names = [  scan_names[id] for id in order ] 
#     else:
#         scan_interval = mydata['scan_interval']
494

495
496
497
498
#         print( order)
#         print( scan_names)
#         print( scan_interval)
#         scan_names = [  scan_names[id] for id in order if scan_ids[id]>=scan_interval[0] and scan_ids[id]<scan_interval[1]  ] 
Alessandro Mirone's avatar
Alessandro Mirone committed
499
    
500
501
#     first_name = scan_names[0]
#     roi_n=0
502
    
503
504
505
506
507
508
509
510
511
#     if not  ('roi_n' in mydata ):
#         while(1):
#             shape  = reader.twoDimages[first_name][roi_n].matrix.shape
#             if shape != (0,) : 
#                 break
#             roi_n+=1
#     else:
#         roi_n = mydata["roi_n"]
#         shape  = reader.twoDimages[first_name][roi_n].matrix.shape
512
            
513
#     Volume  = np.zeros((     shape[0],  shape[1]   ,     len(scan_names)    ))
514

515
516
#     for i,scanname  in enumerate(scan_names):
#         Volume[:,:,i] = reader.twoDimages[scanname][roi_n].matrix
517

518
519
#     if   ('volumeaddress') in mydata :
#         filename, groupname = split_hdf5_address( mydata['volumeaddress'] )
520

521
522
523
524
525
526
527
#         h5=h5py.File(filename,'a')
#         check_libre( h5 , groupname   ) 
#         h5[groupname] = Volume
#         h5.close()
#         h5=None
#     else:
#        view_Volume_myavi_(Volume) 
Alessandro Mirone's avatar
Alessandro Mirone committed
528
529


530

531
532
533
534
# def view_Volume_myavi(mydata):
#     """
#     volume_address : "myfile.hdf5:/path/to/hdf5/group"  # the target destination for volume. 
#     """
Alessandro Mirone's avatar
Alessandro Mirone committed
535
536


537
538
539
540
541
#     filename, groupname = split_hdf5_address( mydata['volume_address'] )
#     h5=h5py.File(filename,'r')
#     Volume = h5[groupname] [:]
#     h5.close()
#     h5=None
542
543


544
545
#     isolevel = mydata['isolevel']
#     opacity = mydata['opacity']
546
547

    
548
#     view_Volume_myavi_(Volume,  isolevel, opacity)
549
550
551
    


552
553
554
555
556
557
558
559
# def view_Volume_myavi_(V,  isolevel, opacity) :
#     print( " IN view ")
#     src = mlab.pipeline.scalar_field(V)
#     mlab.pipeline.iso_surface(src, contours=[V.min()+isolevel*V.ptp(), ], opacity=opacity)
#     mlab.show()
#     src = mlab.pipeline.scalar_field(V)
#     mlab.pipeline.volume(src,vmin=1000.0, vmax=2000.0)
#     mlab.show()
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575


def calculate_recenterings(mydata):
    """
    **calculate_recenterings**

    calculates offsets to go from baricenter A to baricenter B, for all the rois
    in 
 

      calculate_recenterings:
         bariA : "demo_rois.h5:/ROI_AS_SELECTED/images_broad/scans/scan342"
         bariB : "demo_rois.h5:/ROI_AS_SELECTED/energy_scan/scans/scan237"
         target: "recenterings.h5:/recenterings4rois" 
    #
    """
576
577
    allowed_keys =["bariA","bariB","target",   ]
    check_allowed_keys(mydata, allowed_keys)
Alessandro Mirone's avatar
Alessandro Mirone committed
578
    
579
580
581
582
583
584
585
586
587
    bariA = mydata["bariA"]
    bariA_filename, bariA_groupname = split_hdf5_address( bariA )
    
    bariB = mydata["bariB"]
    bariB_filename, bariB_groupname = split_hdf5_address( bariB )

    target = mydata["target"]
    target_filename, target_groupname = split_hdf5_address( target )

Alessandro Mirone's avatar
Alessandro Mirone committed
588
    print( bariA_filename, bariA_groupname)
myron's avatar
myron committed
589
    print( " OPENIN FILE ", bariA_filename  , " FOR RECENTERING ")
590
591
592
593
594
595
596
597
    h5A_f = h5py.File(bariA_filename,"r")
    h5A = h5A_f [bariA_groupname]
    if bariB_filename == bariA_filename :
        h5B_f =  h5A_f
    else:
        h5B_f=h5py.File(bariB_filename,"r")
        
    h5B = h5B_f[bariB_groupname]
598
    offs = {}
599
600

    chiavi = filterRoiList(h5A.keys(), prefix="")
601
    for c in chiavi:
602
603
604
605
        bA = h5A[c]["barix"][()] + h5A[c]["cornerpos"][:][1]
        bB = h5B[c]["barix"][()] + h5B[c]["cornerpos"][:][1]
        bAy = h5A[c]["bariy"][()] + h5A[c]["cornerpos"][:][0]
        bBy = h5B[c]["bariy"][()] + h5B[c]["cornerpos"][:][0]
606
        offs[c] = [[bBy, bAy ],[bB,bA]]
607
608
609
610
        
    if h5B_f is not h5A_f:
        h5B_f.close()
    h5A_f.close()
611
612
613


    if os.path.exists(target_filename):
614
        # check_libre( target_filename ,  target_groupname  )
615
        h5f = h5py.File(target_filename,"a")
616
617
618
619
        if target_groupname in h5f:
            del h5f[target_groupname]

        
620
621
    else:
        h5f = h5py.File(target_filename,"w")
622
623
624
    h5f.require_group(  target_groupname  )
    h5 = h5f[target_groupname]
    for c in chiavi :
625
       h5[c] = np.array(  offs[c]    )
626
627
628
629
    h5f.flush()
    h5f.close()
    h5f = None

Alessandro Mirone's avatar
Alessandro Mirone committed
630

631
# def     sum_scans2maps(mydata):
Alessandro Mirone's avatar
Alessandro Mirone committed
632

633
634
#     roiaddress=None
#     roiaddress = mydata["mask_file"]
Alessandro Mirone's avatar
Alessandro Mirone committed
635
  
636
#     filename, groupname = split_hdf5_address( roiaddress)
Alessandro Mirone's avatar
Alessandro Mirone committed
637

638
639
640
641
#     file= h5py.File(filename,"r")
#     rois = {}
#     shape=xrs_rois.load_rois_fromh5(file[groupname],rois)
#     file.close()
Alessandro Mirone's avatar
Alessandro Mirone committed
642
 
643
644
645
#     specfile_name = mydata["spec_file"]
#     Scan_Variable = mydata["Scan_Variable"]
#     Motor_Variable = mydata["Motor_Variable"]
Alessandro Mirone's avatar
Alessandro Mirone committed
646
647

    
648
#     specfile = SpecIO.Specfile( specfile_name )
Alessandro Mirone's avatar
Alessandro Mirone committed
649
650


651
652
#     dirname  = os.path.dirname(   specfile_name )   
#     basename = os.path.basename(   specfile_name )  
Alessandro Mirone's avatar
Alessandro Mirone committed
653
654

    
655
656
#     scans_infos = []
#     signals    = []
Alessandro Mirone's avatar
Alessandro Mirone committed
657

658
659
#     s1 = int(mydata["first_scan"])
#     s2 = int(mydata["last_scan"])
Alessandro Mirone's avatar
Alessandro Mirone committed
660

661
662
#     roi_names = list(rois.keys())
#     roi_list = [  rois[k] for k in roi_names  ]
Alessandro Mirone's avatar
Alessandro Mirone committed
663
    
664
665
666
667
#     for i in range(s1,s2+1):
#         # print " SCAN lettura " , i
#         scan = specfile.select(str(i))
#         scan_data = scan.data()
Alessandro Mirone's avatar
Alessandro Mirone committed
668
        
669
670
#         scan_themotor   =  scan.motorpos( Motor_Variable  )
#         scan_othermotors = [ scan.motorpos( name   )  for name in  scan.allmotors()  if name != Motor_Variable   ]
Alessandro Mirone's avatar
Alessandro Mirone committed
671

672
#         othermotorsname = [name for name in  scan.allmotors()  if name != Motor_Variable]
673
        
674
675
676
677
678
679
680
681
682
683
#         labels = scan.alllabels()
#         scan_variable = scan_data[ labels.index( Scan_Variable )   , :]
#         scan_ccdnos  = scan_data[ labels.index( "ccdno" )   , :].astype("i")
#         signal = []
#         for no in scan_ccdnos:
#             print( " opening image ", os.path.join(   dirname   ,  "edf", basename+"_"+str(no)+".edf"))
#             data = fabio.open(  os.path.join(   dirname   ,  "edf", basename+"_"+str(no)+".edf" )   ).data 
#             tok = [ (data[corner[0]:corner[0]+mask.shape[0],   corner[1]:corner[1]+mask.shape[1]]*mask).sum() for corner, mask in roi_list  ]
#             signal.append(tok)
#             # print " OK "
Alessandro Mirone's avatar
Alessandro Mirone committed
684
685


686
687
#         # print signal
#         # print " Appendo " , signal
Alessandro Mirone's avatar
Alessandro Mirone committed
688
689

    
690
691
692
#         signals.append(np.array(signal))
#         # print "OK " 
#         scans_infos.append( [scan_themotor, scan_othermotors,  scan_variable  ]  )
Alessandro Mirone's avatar
Alessandro Mirone committed
693

694
#         # print " DONE scan " , i
Alessandro Mirone's avatar
Alessandro Mirone committed
695
        
696
697
698
699
700
701
702
703
704
705
#     done = np.zeros( len(scans_infos) ,"i")
#     DONE={}
#     synthes = {}
#     for kscan in range(len(scans_infos) ):
#         # print " kscan " , kscan
#         DONE[kscan]=0
#         if done[kscan]:
#             continue
#         else:
#             res  = np.array(signals[kscan])
Alessandro Mirone's avatar
Alessandro Mirone committed
706
            
707
708
709
#             done[kscan]=1
#             kinfos = scans_infos[kscan]
#             kM, kOM,  kV = kinfos
Alessandro Mirone's avatar
Alessandro Mirone committed
710
            
711
712
713
714
715
716
#             for oscan in  range(len(scans_infos)) :
#                 # print " oscan " , oscan
#                 if done[oscan]:
#                     continue
#                 else:
#                     oinfos = scans_infos[oscan]
Alessandro Mirone's avatar
Alessandro Mirone committed
717
                    
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
#                     oM, oOM,  oV = oinfos

#                     if kM==oM:
#                         if True or (np.abs(np.array(kOM)-np.array(oOM)).sum()==0.0):
#                             print( " SONO UGUALI " )
#                             if len(oV)== len(kV) :
# #                                   if  np.abs(kV-oV).sum()==0.0:
#                                     print( " AGGIUNGO " )
#                                     res = res+np.array(signals[oscan])
#                                     done[oscan]=1


#             print( " AGGIUNGO " , kM, len(kV), len(kOM))
#             synthes[kscan] = [  kM, kV, kOM, res  ]
#             done[kscan]    = 1


#     # print " QUI " 

#     target = mydata["scans_file"]
#     target_filename, target_groupname = split_hdf5_address( target)
#     if os.path.exists(target_filename):
#         h5f = h5py.File(target_filename,"a")
#     else:
#         h5f = h5py.File(target_filename,"w")
#     h5 = h5f.require_group(target_groupname)
#     target_subdir = "collection_"+str(s1)+"_"+str(s2)
#     if target_subdir in  h5:
#         del h5[ target_subdir]
Alessandro Mirone's avatar
Alessandro Mirone committed
747
        
748
#     h5 =  h5.require_group(target_subdir)
Alessandro Mirone's avatar
Alessandro Mirone committed
749
            
750
751
752
#     for kscan in    list(synthes.keys() ):
#         if DONE[kscan] :
#             continue
753
        
754
#         kM, kV, kOM, kdata  = synthes[kscan] 
Alessandro Mirone's avatar
Alessandro Mirone committed
755

756
#         # print " kdata.shape" , kdata.shape
Alessandro Mirone's avatar
Alessandro Mirone committed
757
758
         
        
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
#         myscans = [kdata]
#         myM     = [kM    ]
#         myV     = [kV    ]
#         DONE[kscan] = 1

#         for oscan in    list(synthes.keys() ):
#             if DONE[oscan]:
#                 continue
#             else:
#                 oM, oV, oOM, odata  = synthes[oscan]
#                 # print " COMPARO "
#                 # diff =  np.array(kOM)-np.array(oOM)
#                 # for nn ,dd in zip(othermotorsname, diff    ):
#                 #     if dd >0.0:
#                 #         print nn, dd
#                 if True or (np.abs(np.array(kOM)-np.array(oOM)).sum()==0.0):
#                     print( " OM equal " ,  len(kV) , len(oV))
#                     if( len(kV) == len(oV)):
#                       # print np.abs(kV-oV)
#                       # if np.abs(kV-oV).sum()==0.0:

#                         assert( kM!=oM  )
#                         DONE[oscan]=1

#                         # print " odata.shape " ,  odata.shape
Alessandro Mirone's avatar
Alessandro Mirone committed
784
                        
785
#                         myscans.append(odata)
Alessandro Mirone's avatar
Alessandro Mirone committed
786

787
#                         # print "  myscans.shape " , np.array(myscans).shape
Alessandro Mirone's avatar
Alessandro Mirone committed
788
                        
789
790
791
792
#                         myM.append(oM)
#                         myV.append(oV)
#         DONE[kscan]=1 
#         myscans = np.array(myscans)
Alessandro Mirone's avatar
Alessandro Mirone committed
793
794
795



796
797
798
799
800
801
802
803
#         if len(myM)>1:
#             order = np.argsort(myM)
#             myM     = np.array(myM)[order]
#             myscans = np.array(myscans)[order]
#             myV     = np.array(myV)
#         else:
#             myM     = np.array(myM)
#             myV     = np.array(myV)
Alessandro Mirone's avatar
Alessandro Mirone committed
804
            
805
806
#         ts = "scan_"+str(int(s1)+kscan)
#         h5t = h5.require_group(ts)
Alessandro Mirone's avatar
Alessandro Mirone committed
807

808
809
#         h5t["Variable"] = myV
#         h5t["Motor"]    = myM
Alessandro Mirone's avatar
Alessandro Mirone committed
810

811
#         s="f, axarr = pylab.plt.subplots(%d, sharex=True)\n"%len(roi_names)
812

813
814
#         ly = myM[0]
#         hy = myM[-1]
815

816
817
818
819
820
821
822
823
824
#         if (  np.abs(myV[-1,:]-myV[0,:]).sum()==0.0 ):
#               xlabel =  "An. Energy"
#               lx = myV[0,0]
#               hx = myV[0,-1]
#         else:
#               xlabel =  "Energy Loss"
#               diff = myM[0]- myV[0,:]
#               lx =  diff[0]
#               hx =  diff[-1]
825
826
827


              
828
829
830
#         for i,rn in enumerate(roi_names):
#             tok = myscans[:,:,i]
#             h5t["signal_"+rn] = tok
Alessandro Mirone's avatar
Alessandro Mirone committed
831
            
832
833
834
835
#             s=s+"axarr[%d].imshow( %s ,extent=( %e ,%e ,%e ,%e) ) \n"%(i, "self.signal_"+rn, lx,hx,ly,hy)
#             s=s+"axarr[%d].set_title('%s')\n"%(i,rn)
#             if i==0:
#                 s=s+"axarr[%d].set_xlabel('%s')\n"%(i,xlabel)
836
            
837
838
839
840
841
#         s=s+"pylab.plt.show()\n"
#         h5t["python_plot"]=s
                                                                        
#     h5f.flush()
#     h5f.close()
Alessandro Mirone's avatar
Alessandro Mirone committed
842
843
844

        
        
845
846
847
848
849
850
851
852
853
# """
#      sum_scans2maps :
#  spec_file : /mntdirect/_data_visitor/hc2892/id20/run7_hc2892/rixs
#  first_scan : 127
#  last_scan : 136
#  Scan_Variable : Anal Energy
#  Motor_Variable : energy
#  scans_file : /scisoft/users/mirone/WORKS/matlabID20/rawdata/scansfile.h5:SCAN
# """  
854

855
def loadscan_2Dimages(mydata):
856
857
858
859
860
861
862
    """
    **loadscan_2Dimages**

    This command harvest the selected signals.
    the instructions on the scans to be taken must be in the form( as submembers ofload_scans ) ::


863
      loadscan_2Dimages :
864
865
866
         roiaddress :  "hdf5filename:nameofroigroup"  # the same given in create_rois
         expdata    :  "absolutepathtoaspecfile"  # this points to a spec file

867
         scan_interval    :  [372,423]   # from 372 to 422 included
868

Alessandro Mirone's avatar
Alessandro Mirone committed
869
         signaladdress : "nameofsignalgroup"  # Target group for writing Relative to nameofroigroup/ (and in the same file)!!!!
Alessandro Mirone's avatar
Alessandro Mirone committed
870
                                              # unless it is in the format filename:groupname
871
872

         energycolumn : 'sty'        # OPTIONAL, if not give defaults to sty
873
         monitorcolumn : 'kapraman'  # OPTIONAL , default is kapraman. If the key is not found in spec file, then no normalisation will be done
874
875
                                     # You can also write kapraman/1000 in this case the monito will be divided by 1000
                                     # (or the other number that you write)
876
         edfName    :  'edfprefix'   # OPTIONAL, if not given autonmatically determined
877

878
         sumto1D  : 1  # OPTIONAL, default 0
879

880
881
         isolateSpot : 0   # is different from zero selects on each image  ( when sumto1d=0 ) ROI  the spot region and sets to zero outside a radius =  isolateSpot

882
883
         # the following defaults to None
         recenterings : "recenterings.h5:/recenterings4rois"
884
885
886
 
    #
    """
myron's avatar
myron committed
887
    allowed_keys = ["roiaddress", 'monitorcolumn', 'recenterings', 'recenterings_confirmed', 'energycolumn', 'edfName', 'isolateSpot', "expdata","scan_interval","scan_list","signaladdress","save_also_roi", 'sumto1D',]
888
889
    check_allowed_keys(mydata, allowed_keys)

890

891
    
892
893
894
895
896
897
898
899
900
901
902
903
904
    roiaddress=None
    roiaddress = mydata["roiaddress"]
  
    filename, groupname = split_hdf5_address( roiaddress)

    file= h5py.File(filename,"r")
    rois = {}
    shape=xrs_rois.load_rois_fromh5(file[groupname],rois)
    file.close()

    roiob = xrs_rois.roi_object()
    roiob.load_rois_fromMasksDict(rois ,  newshape = shape, kind="zoom")

905
    monitor_divider = 1.0
906
    
907
    if  ('monitorcolumn') in mydata :
908
909
910
911
        monitorcolumn = mydata['monitorcolumn']
    else:
        monitorcolumn = 'kapraman'

912
913
    pos = monitorcolumn.find("/")
    if pos != -1:
914
        monitor_divider = float(monitorcolumn[pos+1:])
915
916
        monitorcolumn = monitorcolumn[:pos]

Alessandro Mirone's avatar
Alessandro Mirone committed
917
918
    print( "monitorcolumn : " , monitorcolumn)
    print( "monitor_divider : " , monitor_divider)
919
        
920
    is_by_refinement =False
921
    if  ('recenterings') in mydata :
922
923
        recenterings = mydata['recenterings']
        recenterings_filename, recenterings_groupname = split_hdf5_address( recenterings )
924
925
        h5f = h5py.File(recenterings_filename,"r")
        h5 = h5f[recenterings_groupname]
926
        recenterings= {}
927
        chiavi = filterRoiList(h5.keys(), prefix="")
928
        for c in chiavi:
Alessandro Mirone's avatar
Alessandro Mirone committed
929
            recenterings[int(c)]= h5[c][:]
930
931
            if recenterings[int(c)].shape == (2,2):
                if nprocs>1:
932
                    raise Exception("When using recentering with refinement parallelism cannote be used")
933
934
935
936
937
938
                is_by_refinement = True
        h5f.close()
        if is_by_refinement:
            recenterings_confirmed = mydata['recenterings_confirmed']
            recenterings_confirmed_filename, recenterings_confirmed_groupname = split_hdf5_address( recenterings_confirmed )
           
939
    else:
940
        recenterings = None
941

942
        
943
    if  ('energycolumn') in mydata :
944
945
946
947
       energycolumn  = mydata['energycolumn']
    else:
       energycolumn  = 'sty'

948
    if  ('edfName') in mydata :
949
950
951
952
        edfName = mydata['edfName']
    else:
       edfName  = None

953
    if  ('sumto1D') in mydata :
954
955
        sumto1D = mydata['sumto1D']
    else:
956
        sumto1D  = 0
957
        
958
    if  ('isolateSpot') in mydata :
959
960
961
        isolateSpot = mydata['isolateSpot']
    else:
        isolateSpot  = 0
962

963
    
964
    reader = xrs_imaging.oneD_imaging( mydata["expdata"]  ,monitorcolumn = monitorcolumn , monitor_divider=monitor_divider,
965
966
                                       energycolumn = energycolumn  , edfName = edfName, sumto1D = sumto1D,
                                       recenterings=recenterings)
967
968
    reader.set_roiObj(roiob)

myron's avatar
myron committed
969

970
    todo_list = []
myron's avatar
myron committed
971
972
973
974
975
976
977
978
979
980
981
982
    if "scan_interval" in mydata:
        scan_interval = mydata["scan_interval"]
        ninterval = len(scan_interval)//2
        for i in range(ninterval):
            todo_list = todo_list + list(range(    int(scan_interval[2*i]),  int(scan_interval[2*i+1])) ) #   *scan_interval[2*i :2*i+2])
    else:
        scan_list = mydata["scan_list"]
        for i in scan_list:
            todo_list = todo_list + [int(i)]
        


983
    mytodo = np.array_split(todo_list, nprocs) [myrank]
myron's avatar
myron committed
984
    print(  " Process ", myrank, " is going to read the following scans ", mytodo)
myron's avatar
myron committed
985
986
987
    maxvalue=0.0
    if(len(mytodo)):
        maxvalue = reader.loadscan_2Dimages( list(mytodo) ,scantype=energycolumn, isolateSpot = isolateSpot)
alex's avatar
alex committed
988
989
    if nprocs>1:
        maxvalue = comm.allreduce(maxvalue, op=MPI.MAX)
Alessandro Mirone's avatar
Alessandro Mirone committed
990

991
992
    if is_by_refinement :
        if nprocs>1:
993
            raise Exception("When using recentering with refinement parallelism cannote be used")
994
995
        if os.path.exists(recenterings_confirmed_filename):
            check_libre(  recenterings_confirmed_filename    , recenterings_confirmed_groupname       )
myron's avatar
myron committed
996

997
998
999
1000
            h5f = h5py.File(recenterings_confirmed_filename,"a")
        else:
            h5f = h5py.File(recenterings_confirmed_filename,"w")
        h5 = h5f.require_group(recenterings_confirmed_groupname)