XRS_swissknife.py 155 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:]
myron's avatar
myron committed
375
376
    if( len(groupname) and groupname[0:1] !="/" ):
        groupname = "/"+groupname
377
378
379
    return filename, groupname 


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

385
386
#     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
387
388


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

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

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

400
401
402
403
404
#          #############################################################
#          # 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
405

406
407
408
409
410
411
#          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
412
413


414
415
#     #
#     """
Alessandro Mirone's avatar
Alessandro Mirone committed
416

417
418
#     roiaddress=None
#     roiaddress = mydata["roiaddress"]
Alessandro Mirone's avatar
Alessandro Mirone committed
419

420
421
422
423
424
#     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
425

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


432
#     reader.set_roiObj(roiob)
Alessandro Mirone's avatar
Alessandro Mirone committed
433

434
435
436
437
#     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
438

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

    

446
#     reader.getspectrum()
Alessandro Mirone's avatar
Alessandro Mirone committed
447

448
449
450
451
452
453
454
455
456
457
#     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
458

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

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

464
465
466
467




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

472
473
#     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
474

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

477
#     """
478

Alessandro Mirone's avatar
Alessandro Mirone committed
479

480
#     reader = xrs_imaging.oneD_imaging( "bidon"  , "bidon",  "bidon"  , "bidon")
481

Alessandro Mirone's avatar
Alessandro Mirone committed
482

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

486
#     reader.load_state_hdf5( filename, groupname)
487

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

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

497
498
499
500
#         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
501
    
502
503
#     first_name = scan_names[0]
#     roi_n=0
504
    
505
506
507
508
509
510
511
512
513
#     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
514
            
515
#     Volume  = np.zeros((     shape[0],  shape[1]   ,     len(scan_names)    ))
516

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

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

523
524
525
526
527
528
529
#         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
530
531


532

533
534
535
536
# 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
537
538


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


546
547
#     isolevel = mydata['isolevel']
#     opacity = mydata['opacity']
548
549

    
550
#     view_Volume_myavi_(Volume,  isolevel, opacity)
551
552
553
    


554
555
556
557
558
559
560
561
# 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()
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577


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" 
    #
    """
578
579
    allowed_keys =["bariA","bariB","target",   ]
    check_allowed_keys(mydata, allowed_keys)
Alessandro Mirone's avatar
Alessandro Mirone committed
580
    
581
582
583
584
585
586
587
588
589
    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
590
    print( bariA_filename, bariA_groupname)
myron's avatar
myron committed
591
    print( " OPENIN FILE ", bariA_filename  , " FOR RECENTERING ")
592
593
594
595
596
597
598
599
    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]
600
    offs = {}
601
602

    chiavi = filterRoiList(h5A.keys(), prefix="")
603
    for c in chiavi:
604
605
606
607
        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]
608
        offs[c] = [[bBy, bAy ],[bB,bA]]
609
610
611
612
        
    if h5B_f is not h5A_f:
        h5B_f.close()
    h5A_f.close()
613
614
615


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

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

Alessandro Mirone's avatar
Alessandro Mirone committed
632

633
# def     sum_scans2maps(mydata):
Alessandro Mirone's avatar
Alessandro Mirone committed
634

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

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

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


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

    
657
658
#     scans_infos = []
#     signals    = []
Alessandro Mirone's avatar
Alessandro Mirone committed
659

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

663
664
#     roi_names = list(rois.keys())
#     roi_list = [  rois[k] for k in roi_names  ]
Alessandro Mirone's avatar
Alessandro Mirone committed
665
    
666
667
668
669
#     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
670
        
671
672
#         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
673

674
#         othermotorsname = [name for name in  scan.allmotors()  if name != Motor_Variable]
675
        
676
677
678
679
680
681
682
683
684
685
#         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
686
687


688
689
#         # print signal
#         # print " Appendo " , signal
Alessandro Mirone's avatar
Alessandro Mirone committed
690
691

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

696
#         # print " DONE scan " , i
Alessandro Mirone's avatar
Alessandro Mirone committed
697
        
698
699
700
701
702
703
704
705
706
707
#     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
708
            
709
710
711
#             done[kscan]=1
#             kinfos = scans_infos[kscan]
#             kM, kOM,  kV = kinfos
Alessandro Mirone's avatar
Alessandro Mirone committed
712
            
713
714
715
716
717
718
#             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
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
747
748
#                     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
749
        
750
#     h5 =  h5.require_group(target_subdir)
Alessandro Mirone's avatar
Alessandro Mirone committed
751
            
752
753
754
#     for kscan in    list(synthes.keys() ):
#         if DONE[kscan] :
#             continue
755
        
756
#         kM, kV, kOM, kdata  = synthes[kscan] 
Alessandro Mirone's avatar
Alessandro Mirone committed
757

758
#         # print " kdata.shape" , kdata.shape
Alessandro Mirone's avatar
Alessandro Mirone committed
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
784
785
#         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
786
                        
787
#                         myscans.append(odata)
Alessandro Mirone's avatar
Alessandro Mirone committed
788

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



798
799
800
801
802
803
804
805
#         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
806
            
807
808
#         ts = "scan_"+str(int(s1)+kscan)
#         h5t = h5.require_group(ts)
Alessandro Mirone's avatar
Alessandro Mirone committed
809

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

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

815
816
#         ly = myM[0]
#         hy = myM[-1]
817

818
819
820
821
822
823
824
825
826
#         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]
827
828
829


              
830
831
832
#         for i,rn in enumerate(roi_names):
#             tok = myscans[:,:,i]
#             h5t["signal_"+rn] = tok
Alessandro Mirone's avatar
Alessandro Mirone committed
833
            
834
835
836
837
#             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)
838
            
839
840
841
842
843
#         s=s+"pylab.plt.show()\n"
#         h5t["python_plot"]=s
                                                                        
#     h5f.flush()
#     h5f.close()
Alessandro Mirone's avatar
Alessandro Mirone committed
844
845
846

        
        
847
848
849
850
851
852
853
854
855
# """
#      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
# """  
856

857
def loadscan_2Dimages(mydata):
858
859
860
861
862
863
864
    """
    **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 ) ::


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

869
         scan_interval    :  [372,423]   # from 372 to 422 included
870

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

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

880
         sumto1D  : 1  # OPTIONAL, default 0
881

882
883
         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

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

892

893
    
894
895
896
897
898
899
900
901
902
903
904
905
906
    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")

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

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

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

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

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

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

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

myron's avatar
myron committed
971

972
    todo_list = []
myron's avatar
myron committed
973
974
975
976
977
978
979
980
981
982
983
984
    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)]
        


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

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

999
1000
            h5f = h5py.File(recenterings_confirmed_filename,"a")
        else: