batch_extraction_esynth1.py 17.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
#  newfactors e' gia pronto , dopo verifica fuzionalita di roi_sel in ximages
import numpy as np
import h5py
import glob
import json

BATCH_PARALLELISM = 1

import os
myron's avatar
myron committed
10
11
12
13
14
15
def main():
    peaks_shifts = np.array(
                     [  6471.983002,
                        6471.318152,
                        6470.612314,
                        6471.217828,
myron's avatar
myron committed
16
17
                     ]) - 6470.6

myron's avatar
myron committed
18
19
20
21
22
23
24
    datadir =  "/data/scisofttmp/mirone/Loic1/data/"
    filter_path =  datadir + "../mask.h5:/filter"
    filename_template = "HC-D_2Dmap_%04d"
    data_path_template =  datadir + filename_template + ".nxs:/root.spyc.config1d_RIXS_0024/scan_data/data_03"
    monitor_path_template = datadir + filename_template +"_monitor.hd5:/monitor"

    # energy_custom_grid = np.array([6746.0,6754.0,6755.5,6756.2,6757.5,6759.3,6762.5,6770.0,6790.5])
myron's avatar
myron committed
25
26
27
28
    
    # energy_custom_grid = np.array([6744.0,6754.0,6755.5,6756.2,6757.5,6759.3,6762.5,6770.0,6792])
    
    energy_custom_grid = np.array([6744.0,6754.0,    6754.8,   6755.5,6756.2,6757.5,     6758.63,         6759.3,6762.5,6770.0,6792])
myron's avatar
myron committed
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

    energy_exp_grid  = h5py.File( datadir +(filename_template%1)+".nxs" ,"r")["/root.spyc.config1d_RIXS_0024/scan_data/actuator_1_1"][()]
    
    scan_interval   =  [1,476]   # from 1 to 475 included 

    Ydim = 25
    Zdim = 19
    Edim = 9

    roi_target_path          = "results/myrois.h5:/ROIS"
    reference_target_file    = "results/response.h5"
    signals_target_file      = "results/extracted.h5" 
    scalarprods_target_file  = "results/scalar_prods.h5"

    steps_to_do = {
        "do_step_make_roi":                      False,
        "do_step_make_reference":                False,
        "do_step_sample_extraction":             False,
        "do_step_scalar_products":               False, 
        "do_step_interpolation_coefficients":    True,
myron's avatar
myron committed
49
        "do_step_finalise_for_fit":              True
myron's avatar
myron committed
50
51
52
    }


myron's avatar
myron committed
53
54
    use_custom_components = None
    # use_custom_components = "components.h5"
myron's avatar
myron committed
55
56
57

    scans_to_use_for_roi = list(range(1,20))

myron's avatar
myron committed
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
    
    tools_sequencer(  peaks_shifts          = peaks_shifts          ,
                      datadir               = datadir               ,
                      filter_path           = filter_path           ,
                      filename_template     = filename_template     ,
                      data_path_template    = data_path_template    ,
                      energy_custom_grid    = energy_custom_grid    ,
                      energy_exp_grid       = energy_exp_grid       ,
                      monitor_path_template = monitor_path_template ,
                      scan_interval         = scan_interval         ,
                      scans_to_use_for_roi  = scans_to_use_for_roi  ,
                      Ydim                  = Ydim                  ,
                      Zdim                  = Zdim                  ,
                      Edim                  = Edim                  ,
                      
                      roi_target_path          = roi_target_path         ,
                      reference_target_file    = reference_target_file   , 
                      signals_target_file      = signals_target_file     , 
                      scalarprods_target_file  = scalarprods_target_file ,

                      
                      steps_to_do = steps_to_do,

                      use_custom_components = use_custom_components

    ) 



87
88
89
90
91
92
93
94
95
def process_input(s, go=0, exploit_slurm_mpi = 0, stop_omp = False):
    open("input_tmp_%d.par"%go, "w").write(s)
    background_activator = ""
    if (go % BATCH_PARALLELISM ):
        background_activator = "&"

    prefix=""
    if stop_omp:
        prefix = prefix +"export OMP_NUM_THREADS=1 ;"
myron's avatar
myron committed
96

97
98
        
    if (  exploit_slurm_mpi==0  ):
myron's avatar
myron committed
99
        comando = (prefix +"mpirun -n 1 XRS_swissknife  input_tmp_%d.par  %s"%(go, background_activator))
100
    elif (  exploit_slurm_mpi>0  ):
myron's avatar
myron committed
101
        comando = (prefix + "mpirun XRS_swissknife  input_tmp_%d.par  %s"%(go, background_activator) )
102
    else:
myron's avatar
myron committed
103
104
105
106
107
        comando = (prefix + "mpirun -n %d XRS_swissknife  input_tmp_%d.par  %s"%(abs( exploit_slurm_mpi  ), go, background_activator) )

    res = os.system( comando )

    assert (res==0) , " something went wrong running command : " + comando
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133

def select_rois(  data_path_template=None, filter_path=None, roi_target_path=None, scans_to_use=None   ):
    
    inputstring = """

    create_rois_galaxies :
       expdata     :  {data_path_template}
       filter_path  :  {filter_path}
       roiaddress :       {roi_target_path}      # the target destination for rois
       scans  : {scans_to_use}

    """ .format(data_path_template = data_path_template,
                filter_path = filter_path,
                roi_target_path =  roi_target_path,
                scans_to_use = scans_to_use
    ) 
    process_input( inputstring , exploit_slurm_mpi = 0 )

def extract_sample_givenrois(
        roi_path = None,
        data_path_template = None,
        monitor_path_template = None,
        scan_interval = None,
        Ydim = None,
        Zdim = None,
        Edim = None,
134
        signals_target_file = None
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
         ):

    inputstring = """
        loadscan_2Dimages_galaxies :

             roiaddress : {roi_path}

             expdata  :  {data_path_template}

             monitor_address : {monitor_path_template}

             scan_interval    :  {scan_interval}

             Ydim : {Ydim}
             Zdim : {Zdim}
             Edim : {Edim}

152
             signalfile : {signals_target_file}
153
154
155
156
157
158
159
160

    """.format( roi_path = roi_path,
                data_path_template = data_path_template,
                monitor_path_template = monitor_path_template,
                scan_interval = scan_interval,
                Ydim = Ydim,
                Zdim = Zdim,
                Edim = Edim,
161
                signals_target_file = signals_target_file) 
162
163
164
165
166
    

    process_input( inputstring, exploit_slurm_mpi = 0)


myron's avatar
myron committed
167
168


myron's avatar
myron committed
169
def InterpInfo_Esynt_components(peaks_shifts , energy_exp_grid = None,    custom_ene_list =  None, custom_components = None ):
myron's avatar
myron committed
170

myron's avatar
myron committed
171
        components = h5py.File( custom_components ,"r")["components"] [()]
myron's avatar
myron committed
172
173
174
175
176
177
178
179
180
181
        info_dict = {}

        
        for i_interval in range(len(components)):
            info_dict[str(i_interval)]      =   {}            
            info_dict[str(i_interval)]["E"] = custom_ene_list[ i_interval   ]
            info_dict[str(i_interval)]["coefficients"]={}

            for i_n in range(len(energy_exp_grid)):
                info_dict[str(i_interval)]["coefficients"  ][  str(i_n)  ]={}
myron's avatar
myron committed
182
                for roi_num, de in enumerate(     peaks_shifts   ):
myron's avatar
myron committed
183
184
185
186
187
188
189
190
191
192
                    info_dict[str(i_interval)]["coefficients"  ][ str(i_n) ][ str(roi_num) ] = 0


        for ic in range(len(components)):
                
            for i_interval in range(len(custom_ene_list)-1):
                cE1 =  custom_ene_list[ i_interval   ]
                cE2 =  custom_ene_list[ i_interval+1 ]
                for i_ene, t_ene  in enumerate( energy_exp_grid)   :

myron's avatar
myron committed
193
                    for roi_num, de in enumerate(     peaks_shifts   ):
myron's avatar
myron committed
194
195
                        if  t_ene+de < cE1 or t_ene+de > cE2:
                            continue
myron's avatar
myron committed
196
                        alpha = (cE2-(t_ene+de) )/(cE2+cE1)
myron's avatar
myron committed
197
198
199
200
201
202
203
204
205

                        info_dict[str(ic)]["coefficients"  ][  str(i_ene)  ][  str(roi_num)  ]   += alpha * components[ic][ i_interval ]

                        
                        info_dict[str(ic)]["coefficients"  ][  str(i_ene)  ][  str(roi_num)  ]   += (1-alpha)*components[ic][ i_interval+1  ]

        return info_dict    

    
myron's avatar
myron committed
206
def  InterpInfo_Esynt( peaks_shifts , energy_exp_grid = None,    custom_ene_list =  None):
207

myron's avatar
myron committed
208
    info_dict = {"energy_exp_grid":list(energy_exp_grid), "de_list": list(peaks_shifts)}
209
210
211
212
213
214
215
216
217
218
219
220

    N_custom = len(custom_ene_list)
    N_data   = len( energy_exp_grid ) 
    
    for i_interval in range(len(custom_ene_list)):
        info_dict[str(i_interval)]      =   {}            
        info_dict[str(i_interval)]["E"] = custom_ene_list[ i_interval   ]
        info_dict[str(i_interval)]["coefficients"]={}

        for i_n in range(len(energy_exp_grid)):

            info_dict[str(i_interval)]["coefficients"  ][  str(i_n)  ]={}
myron's avatar
myron committed
221
            for roi_num, de in enumerate(     peaks_shifts   ):
222
223
224
225
226
227
228
229
230
231
                info_dict[str(i_interval)]["coefficients"  ][ str(i_n) ][ str(roi_num) ] = 0


    for i_interval in range( N_custom -1):
        
        cE1 =  custom_ene_list[ i_interval   ]
        cE2 =  custom_ene_list[ i_interval+1 ]
        
        for i_ene, t_ene  in enumerate( energy_exp_grid)   :
            
myron's avatar
myron committed
232
            for roi_num, de in enumerate(     peaks_shifts   ):
233
234
235
236
237
238
239
240
241
242
243
                if  t_ene+de < cE1 or t_ene+de > cE2:
                    continue
                alpha = (cE2-(t_ene+de) )/(cE2-cE1)

                info_dict[str(i_interval)]["coefficients"  ][  str(i_ene)  ][  str(roi_num)  ]   = alpha
                info_dict[str(i_interval+1)]["coefficients"][  str(i_ene)  ][  str(roi_num)  ] = 1-alpha
                
    return info_dict    



myron's avatar
myron committed
244
    def __init__(self, peaks_shifts, interp_file, source,  custom_ene_list = None):
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
        
        volum_list = list(interp_file[source].keys())
        scan_num_list = np.array([ int( t.split("_") [1]) for t in volum_list])
        
        ene_list      = np.array([    interp_file[source][vn]["scans"]["Scan%03d"%sn ]["motorDict"]["energy"][()] for vn,sn in zip(volum_list, scan_num_list   )   ])

        print ( " ecco la scannumlist " , scan_num_list)
        print (" ecco ene_list", ene_list)
        
        
        self.volum_list    =  volum_list
        self.scan_num_list =  scan_num_list
        self.ene_list      =  ene_list

        order = np.argsort(  self.ene_list    )
        
        self.ene_list  = self.ene_list [order]

        if custom_ene_list is None:
            self.custom_ene_list      = self.ene_list
        else:
            self.custom_ene_list   = custom_ene_list

        self.scan_num_list  = self.scan_num_list [order]
        self.volum_list  = [ self.volum_list [ii]  for ii in order  ] 
        
        self.interp_file=interp_file
        self.source= source
myron's avatar
myron committed
273
        self.peaks_shifts=peaks_shifts
myron's avatar
myron committed
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288

# interpolation_infos_file = "interpolation_infos.json"

# info_dict={}
# for i in range(NC):
#     dizio = {}
#     info_dict[str(i)] = {"coefficients":dizio}
#     c = model.components_[i]
#     np = len(c)
#     for j in range(np):
#         dizio[str(j)] = float(c[j])

# json.dump(info_dict,open( interpolation_infos_file,"w"), indent=4)


289
290
291
292
        
    def interpola_Esynt(self,  roi_sel=roi_sel ):
        print ( " ECCO I DATI ")
        print (  self.ene_list  ) 
myron's avatar
myron committed
293
        print (  self.peaks_shifts   )
294
295
296
297
298
299
300
301
302
303
304
305
306
307

        info_dict = {}

        for i_intervallo in range(len(self.custom_ene_list)):
            info_dict[str(i_intervallo)]      =   {}            
            info_dict[str(i_intervallo)]["E"] = self.custom_ene_list[ i_intervallo   ]
            info_dict[str(i_intervallo)]["coefficients"]={}
            for t_vn, t_sn, t_ene in list(zip(self.volum_list,  self.scan_num_list, self.ene_list    )):
                info_dict[str(i_intervallo)]["coefficients"  ][  t_vn  ]={}

        for i_intervallo in range(len(self.custom_ene_list)-1):
            cE1 =  self.custom_ene_list[ i_intervallo   ]
            cE2 =  self.custom_ene_list[ i_intervallo+1 ]
            for t_vn, t_sn, t_ene in list(zip(self.volum_list,  self.scan_num_list, self.ene_list    ))[0:]:
myron's avatar
myron committed
308
                for roi_num, de in enumerate(     self.peaks_shifts   ):
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
                    if roi_num not in roi_sel:
                        continue
                    if  t_ene+de < cE1 or t_ene+de > cE2:
                        continue
                    alpha = (cE2-(t_ene+de) )/(cE2-cE1)

                    info_dict[str(i_intervallo)]["coefficients"  ][  str(t_vn)  ][  str(roi_num)  ]   = alpha
                    info_dict[str(i_intervallo+1)]["coefficients"][  str(t_vn)  ][  str(roi_num)  ] = 1-alpha

        return info_dict    

    
def get_reference( roi_path = None, reference_target_file = None ):
    inputstring = """

    loadscan_2Dimages_galaxies_foilscan :
        roiaddress : {roi_path}           
        expdata    :  None  
        signalfile : {reference_target_file}  # Target file for the signals
    """ .format( roi_path = roi_path,  reference_target_file = reference_target_file ) 

330
    process_input( inputstring , exploit_slurm_mpi = 0) 
331
332
333
334
335
336
337
338





def get_scalars( iE = None,
                 signals_file      = None,
                 reference_file    = None,
339
                 target_file    = None
340
341
342
343
344
345
):

    inputstring = """
    superR_scal_deltaXimages_Esynt :
         sample_address : {signals_file}:/E{iE}
         delta_address : {reference_file}:/Scan0
346
         load_factors_from : 
347
         nbin : 1
348
349
         target_address : {target_file}:/{iE}/scal_prods
    """ . format( iE = iE,
350
351
                  signals_file      = signals_file      ,
                  reference_file    = reference_file    ,
352
                  target_file    = target_file,
353
354
355
356
357
358
359
360
361
362
363
    )
    process_input( inputstring, exploit_slurm_mpi = 0)    


    
def get_volume_Esynt( scalarprods_file = None, interpolation_file = None):
    
    os.system("mkdir DATASFORCC")

    inputstring = """    
     superR_getVolume_Esynt :
364
365
        scalprods_address : {scalarprods_file}:/
        target_address : {scalarprods_file}:/data_for_volumes
366
367
368
369
370
        dict_interp : {interpolation_file}
        debin : [1, 1]
        output_prefix : DATASFORCC/test0_
    """.format( scalarprods_file = scalarprods_file ,
                interpolation_file = interpolation_file
371
    )     
372
    process_input( inputstring,  exploit_slurm_mpi = 0)
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407

    
def myOrder(tok):
    if("volume" not  in tok):
        tokens = tok.split("_")
        print( tokens)
        return int(tokens[1])*10000+ int(tokens[2])
    else:
        return 0
    
def reshuffle(   volumefile  = "volumes.h5",   nick = None    ):

    h5file_root = h5py.File( volumefile ,"r+" )
    h5file = h5file_root[nick]
    scankeys = list( h5file.keys())
    scankeys.sort(key=myOrder)
    print( scankeys) 
    
    volumes = []
    for k in scankeys:
        if k[:1]!="_":
            continue
        print( k)
        if "volume" in h5file[k]:
            volumes.append( h5file[k]["volume"]  )
    # volume = np.concatenate(volumes,axis=0)
    volume = np.array(volumes)
    if "concatenated_volume" in h5file:
        del h5file["concatenated_volume"]
    h5file["concatenated_volume"]=volume
    h5file_root.close()
    
## THE FOLLOWING PART IS THE RELEVANT ONE


myron's avatar
myron committed
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
def tools_sequencer(  peaks_shifts = None,
                      datadir = None,
                      filter_path = None,
                      filename_template = None,
                      data_path_template = None,
                      energy_custom_grid = None,
                      energy_exp_grid  = None,
                      monitor_path_template = None,
                      scan_interval = None,
                      scans_to_use_for_roi = None,
                      Ydim = None,
                      Zdim = None,
                      Edim = None, 
                      roi_target_path          = None,
                      reference_target_file    = None,
                      signals_target_file      = None,
                      scalarprods_target_file  = None,
                      steps_to_do = None,
                      use_custom_components = None
) :
428
429


myron's avatar
myron committed
430
431
432
433
    
    if(steps_to_do["do_step_make_roi"]):   # ROI selection and reference scan
        select_rois(  data_path_template = data_path_template, filter_path = filter_path, roi_target_path = roi_target_path, scans_to_use=scans_to_use_for_roi   )
    roi_path = roi_target_path
434
435


myron's avatar
myron committed
436
437
438
    if(steps_to_do["do_step_make_reference"]):  
        get_reference( roi_path = roi_path ,    reference_target_file = reference_target_file   )
    reference_file = reference_target_file
439
440
441



myron's avatar
myron committed
442
443
444
445
446
    if(steps_to_do["do_step_sample_extraction"]): # SAMPLE extraction
        extract_sample_givenrois(
            roi_path              = roi_path               ,
            data_path_template    = data_path_template     ,
            monitor_path_template = monitor_path_template  ,
447

myron's avatar
myron committed
448
449
450
451
452
453
454
            scan_interval      = scan_interval      ,
            Ydim               = Ydim               ,
            Zdim               = Zdim               ,
            Edim               = Edim               ,
            signals_target_file = signals_target_file 
        )
    signals_file = signals_target_file    
455
456
457
458




myron's avatar
myron committed
459
460
461
462
463
464
465
466
    if(steps_to_do["do_step_scalar_products"]):
        os.system("rm %s"%scalarprods_target_file)
        for iE in range(Edim) :
            get_scalars(  iE = iE,
                          signals_file = signals_file,
                          reference_file = reference_file,
                          target_file =  scalarprods_target_file
            )
467

myron's avatar
myron committed
468
    scalarprods_file = scalarprods_target_file
469
470


myron's avatar
myron committed
471
472
473
474
475
476
477
478
479
480
481
482
483
484
    interpolation_infos_file = "interpolation_infos.json"
    if(steps_to_do["do_step_interpolation_coefficients"]):    # INTERPOLATION  ESYNTH
        if use_custom_components is None:        
            info_dict = InterpInfo_Esynt(  peaks_shifts ,
                                           energy_exp_grid = energy_exp_grid,
                                           custom_ene_list =  energy_custom_grid
            )
        else:
            info_dict = InterpInfo_Esynt_components(  peaks_shifts,
                                                      energy_exp_grid = energy_exp_grid,
                                                      custom_ene_list =  energy_custom_grid,
                                                      custom_components = use_custom_components
            )
        json.dump(info_dict,open( interpolation_infos_file,"w"), indent=4)
485
486
487
488



    
myron's avatar
myron committed
489
490
491
492
    # ### ESYNTH
    if(steps_to_do["do_step_finalise_for_fit"]):
        get_volume_Esynt( scalarprods_file = scalarprods_file,
                          interpolation_file = interpolation_infos_file)
493
494
495



myron's avatar
myron committed
496
main()