Skip to content

GitLab

Commit 3a7e1658 authored by myron's avatar myron
Browse files

restructured response fitting

parent d7574d9a
Hide whitespace changes
Inline Side-by-side
nonregressions/volumes/interpola_1/batch_extraction_esynth1.py
View file @ 3a7e1658
......@@ -10,36 +10,29 @@ BATCH_PARALLELISM = 4
import os
def main():
peaks_shifts = h5py.File("peaks_positions_for_analysers.h5","r")["peaks_positions"][()]
assert( len(peaks_shifts) == 72)
os.system("xz -dk mask.h5.xz")
Enominal = np.median(peaks_shifts)
peaks_shifts-= Enominal
datadir = "/data/id20/inhouse/data/run3_20/run3_es949"
os.system("mkdir results")
roi_target_path = "results/myrois.h5:/ROIS"
# roi_target_path = "rr.h5:/extracted/ROI_AS_SELECTED"
filter_path = "mask.h5:/FILTER_MASK/filter"
roi_scan_num = [245,246,247]
reference_scan_list = [245, 246, 247]
roi_scan_num = [245,246,247]
monitor_column = "izero/0.000001"
first_scan_num = 651
Ydim = 25
Zdim = 10
Edim = 7
monitor_column = "izero/0.000001"
signals_target_file = "results/signals.h5"
interpolated_signals_target_file = "results/interpolated_signals.h5"
extracted_reference_target_file = "results/reference.h5"
peaks_shifts = h5py.File("peaks_positions_for_analysers.h5","r")["peaks_positions"][()]
assert( len(peaks_shifts) == 72)
reference_scan_list = [245, 246, 247]
Enominal = np.median(peaks_shifts)
peaks_shifts-= Enominal
datadir = "/data/id20/inhouse/data/run3_20/run3_es949"
# If reference_clip is not None, then a smaller part of the reference scan is considered
# This may be usefule to obtain smaller volumes containing the interesting part
......@@ -51,18 +44,41 @@ def main():
## in the reference scan for each position there is a spot with a maximum. We set zero the background which is further than
## such radius from the maximum
isolate_spot_by = 6
#### For the fit of the response function based on reference scans
response_fit_options = dict( [
["niter_optical" , 100],
["beta_optical" , 0.1],
["niter_global" , 3 ]
])
steps_to_do = {
"do_step_make_roi": False,
"do_step_sample_extraction": False,
"do_step_interpolation": False,
"do_extract_reference_scan": False,
"do_step_make_reference": False,
"do_step_scalar_products": False,
"do_step_finalise_for_fit": False
"do_fit_reference_response": False,
"do_resynthetise_reference": False,
"do_step_finalise_for_fit": False
}
os.system("mkdir results")
roi_target_path = "results/myrois.h5:/ROIS"
signals_target_file = "results/signals.h5"
interpolated_signals_target_file = "results/interpolated_signals.h5"
extracted_reference_target_file = "results/reference.h5"
response_target_file = "results/response.h5"
tools_sequencer( peaks_shifts = peaks_shifts ,
datadir = datadir ,
......@@ -84,7 +100,11 @@ def main():
reference_scan_list = reference_scan_list,
reference_clip = reference_clip,
extracted_reference_target_file = extracted_reference_target_file ,
isolate_spot_by = isolate_spot_by
isolate_spot_by = isolate_spot_by,
response_target_file = response_target_file,
response_fit_options = response_fit_options
)
......@@ -298,6 +318,77 @@ def interpolate( peaks_shifts, interp_file_str, interp_file_target_str):
fScan[str(k)]["monitor_divider"] = 1.0
def synthetise_response(scan_address=None, target_address=None, response_fit_options = None
):
inputstring = """
superR_fit_responses :
foil_scan_address : "{scan_address}"
nref : 7 # the number of subdivision per pixel dimension used to
# represent the optical response function at higher resolution
niter_optical : {niter_optical} # the number of iterations used in the optimisation of the optical
# response
beta_optical : {beta_optical} # The L1 norm factor in the regularisation
# term for the optical functions
pixel_dim : 1 # The pixel response function is represented with a
# pixel_dim**2 array
niter_pixel : 10 # The number of iterations in the pixel response optimisation
# phase. A negative number stands for ISTA, positive for FISTA
beta_pixel : 0.0 # L1 factor for the pixel response regularisation
## The used trajectories are always written whith the calculated response
## They can be reloaded and used as initialization(and freezed with do_refine_trajectory : 0 )
## Uncomment the following line if you want to reload a set of trajectories
## without this options trajectories are initialised from the spots drifts
##
# reload_trajectories_file : "response.h5"
filter_rois : 0
######
## The method first find an estimation of the foil scan trajectory on each roi
## then, based on this, obtain a fit of the optical response function
## assuming a flat pixel response. Finally the pixel response is optimised
##
## There is a final phase where a global optimisation
## is done in niter_global steps.
##
## Each step is composed of optical response fit, followed by a pixel response fit.
## If do_refine_trajectory is different from zero, the trajectory is reoptimised at each step
##
niter_global : {niter_global}
## if do_refine_trajectory=1 the start and end point of the trajectory are free
## if =2 then the start and end point are forced to a trajectory which is obtained
## from a reference scan : the foil scan may be short, then one can use the scan of
## an object to get another one : key *trajectory_reference_scan_address*
##
do_refine_trajectory : 1
## optional: only if do_refine_trajectory = 2
trajectory_reference_scansequence_address : "demo_newrois.h5:/ROI_FOIL/images/scans/"
trajectory_threshold : 0.1
## if the pixel response function is forced to be symmetrical
simmetrizza : 1
## where the found responses are written
target_file : {target_address}
# target_file : "fitted_responses.h5"
"""
s=inputstring.format( scan_address=scan_address ,
target_address=target_address,
niter_optical = response_fit_options[ "niter_optical"],
beta_optical=response_fit_options["beta_optical"],
niter_global=response_fit_options["niter_global"]
)
process_input( s , exploit_slurm_mpi = 1, stop_omp = True)
def tools_sequencer( peaks_shifts = None,
datadir = None,
......@@ -326,7 +417,8 @@ def tools_sequencer( peaks_shifts = None,
isolate_spot_by = None,
reference_scan_list = None,
extracted_reference_target_file = None ,
response_target_file = None,
response_fit_options = None
) :
......@@ -368,9 +460,7 @@ def tools_sequencer( peaks_shifts = None,
extracted_reference_target_file = extracted_reference_target_file,
isolate_spot_by = isolate_spot_by,
reference_scan_list = reference_scan_list
)
)
if reference_clip is not None:
clip1, clip2= reference_clip
......@@ -388,8 +478,7 @@ def tools_sequencer( peaks_shifts = None,
mat = target_group[k][dsn][()]
del target_group[k][dsn]
target_group[k][dsn] = mat[clip1:clip2]
ftarget.close()
ftarget.close()
ftarget = h5py.File( extracted_reference_target_file ,"r+")
ftarget["references/scans/ScansSum"] = ftarget["references/scans/Scan%03d"% reference_scan_list[0] ]
......@@ -405,5 +494,17 @@ def tools_sequencer( peaks_shifts = None,
ftarget.close()
if(steps_to_do["do_fit_reference_response"]): # of course we need the REFERENCE SCAN
synthetise_response(
scan_address= extracted_reference_target_file +":references/scans/ScansSum" ,
target_address = response_target_file +":/FIT",
response_fit_options = response_fit_options
)
main()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment