3-Reconstruct_template (no regularization) | mumott v 0.2.1

• This is tutorial notebook is inspired by the one uploaded on the Mumott Git repository (https://gitlab.com/liebi-group/software/mumott.git) e.g. mumott/tutorial/reconstruct_and_visualize.ipynb. It runs with the installation of mumott v 0.2.1 (should have been installed via pip install mumott v 0.2.1)

• More detailed information on Mumott usage can be found in this detailed notebook and it is recommended to take a closer look there for details

• The notebook runs on CPU resources, but requires that you allocate multiple nodes. An example to get those for instance in ra would be: salloc -p day -c 20 -mem-per-cpu 10000 --time 1-00:00:00. This gets you 20 cores for a day on ra.

Step 1 - Set environment parameters for Mumott

These 2 cells need to run in order. They are required to set parameters for the CPU multithreading for the numba package. The optimal number of threads is hardware- and system-dependent. For large systems, it may be necessary to restrict the number of threads to avoid using too much RAM. In most cases, 4-20 threads is a good choice.

import os
#--- User Input with the number of cores that you requested, 12-20 seems to work good for numba --------------------#
nthreads = 20

#--- End Input  ----------------------------------------------------------------------------------------------------#

#---- Do not modify ----#
os.environ["NUMBA_NUM_THREADS"] = f'{nthreads}'
os.environ["NUMBA_DEFAULT_NUM_THREADS"] = f'{nthreads}'
os.environ["OMP_NUM_THREADS"] = f'{nthreads}'

import numba

Step 1.1 - Please run this one after the cell above

#---- Do not modify ----#
numba.config.NUMBA_DEFAULT_NUM_THREADS = nthreads
numba.config.NUMBA_NUM_THREADS = nthreads
numba.set_num_threads(nthreads)

Step 1.2 - Package and libraries

#---- Do not modify ----#
%matplotlib ipympl
import h5py
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import display
from ipywidgets import interact
# Mumott related inports
from mumott.data_handling import DataContainer, TransformParameters
from mumott.output_handling import OutputHandler, LiveViewHandler, ProjectionViewer
from mumott.optimization import OptimizationParameters, Optimizer

Step 2 - Load data

• The correctly prepared data needs to be loaded first. For this, please use the script 2-Export_template.

• In Mumott this is performed using the path, name, and input type are defined as strings. The path can be either relative or absolute. The file name needs to include the file ending.

• The two allowed formats are 'mat' (Aligned CSAXS-type Matlab file with a projection structure containing necessary data and metadata), and 'h5' (which requires a specially formated hdf5 file as done in the 2-Export_template).

First load args from file

#---- Modify input for your sample! ----#
filepath_args = '/data/visitors/formax/20220566/2023031408/process/work_dir_Christian/big_brain/args_big_brain.pickle'
#--- End Input to load the exported args file ---------------------#

import pickle
## ---------- Don't change this if not needed ----------------- ##
with open(filepath_args,'rb') as fid:
    args = pickle.load(fid)

Fix path for data loading of Mumott

#---- Modify input for your sample! ----#
#path = '/data/visitors/formax/20220566/2023031408/process/reconstructions/big_brain/'
#path = '/data/visitors/formax/20220566/2023031408/process/reconstructions/9425L_4w_XHP_NT/'
path = os.path.join(args['work_directory'] , args['sample_name'])

input_type = 'h5'
name = f'dataset_q_0.016_0.020.{input_type}'#f'dataset_q_0.079_0.090.{input_type}' # f"dataset_q_{q_range[0]:.3f}_{q_range[1]:.3f}"

Step 2.1 - Initializing the DataContainer

Mumott is working with the object DataContainer. It will be created based on the input paths

#---- Do not modify ----#
try:
    data_container = DataContainer(data_path=path,
                                   data_filename=name,
                                   data_type=input_type)
except FileNotFoundError:
    print('No data file found!')

We can now take a quick look at the container object to see that the parameters make sense. Correct input values in case they are off, e.g. corrected for transmission is a candidate that might be off in case data has already been normalized during the Export

display(data_container)

DataContainer

Field	Size
Number of projections	419
Volume shape	[ 87 131 87]
Corrected for transmission	False
Angles in radians	True
Transformation applied	False
Reconstruction parameters generated	False

Up until now, no geometrical adjustments have been performed to the data. In the dataContainer, you should check and confirm the number of projections, the volume shape, and the fact that we have converted our angles into radians. Corrected for transmission is set to False because the code at the moment has no way of telling whether our data is transmission corrected or not, but we know and can adjust this if we like. Otherwise it can also remain there. We also see that Transformation applied is False, and this leads us to the next important step.

# Note, the projection viewer below can be used to also directly check the imported data
p = ProjectionViewer(data_container.stack.data, data_container.stack.dimensions, data_container.stack.detector_angles)

Step 2.2 - Transform the data container to the correct geometry

Data from different beamlines come in different coordinate systems and sorted in different ways, and they must be transformed into the coordinate system consistent with the one used by this method. There is some redundancy in this definition, but the three-dimensional coordinate system native to the code has the properties of: - principal_rotation_right_handed and secondary_rotation_right_handed - These simply denote whether the rotation and tilt, after considering data_sorting and data_index_origin, are right-handed. Setting to False will mean that the rotation or tilt value is multiplied by -1.

Generally, the interdependence of these parameters can be counterintuitive to work out and setting them correctly requires some trial and error, however, for a particular measurement setup and data processing pipeline, they only need to be worked out once. It can be very difficult to work them out correctly with a highly symmetrical sample, such as a cylinder. In this case, it is strongly recommended to work them out using a less symmetrical sample measured in the same system.

#---- Do not modify ----# ForMAX, loaded from args.
transform_parameters = TransformParameters(
    data_sorting = (0, 1, 2),#args['data_sorting'], #(1, 0, 2)
    data_index_origin = (0, 1),#args['data_index_origin'], #(0, 0),
    principal_rotation_right_handed = True, #args['principal_rotation_right_handed'], #True,
    secondary_rotation_right_handed = True, #args['secondary_rotation_right_handed'], #True,
    detector_angle_0 = (1,0), #args['detector_angle_0'], #(1, 0),
    detector_angle_right_handed = False, #args['detector_angle_right_handed'], #=False,
    offset_positive = (True, True)) #args['offset_positive']) #(True, True))

Perform transformation

#---- Do not modify ----#
data_container.transform(transform_parameters)

Step 2.3 - Visualize transformed data and check projections

The sample should appear as it was measured at the beamline. Furthermore, the sample should also rotate and tilt as expected. It is crucial to check these aspects in this step

#---- Do not modify ----#
p = ProjectionViewer(data_container.stack.data, data_container.stack.dimensions, data_container.stack.detector_angles)

Having specified the settings (in this case, the settings are simply the default ones), we can hand this object over to the data_container.transform() method.

Step 3 - Visualize transformed data and check projections

• IMPORTANT: The current version of Mumott has a bug reshaping the reconstruction volume. This needs to be manually adjusted in the next cell.

• For this, use the visualization from the projection viewer to restore the volume shape. In this case, the dimension in direction of the principal (tomography) axis needs to be restored

• X , Y , Z are dimensions |

• Z is direction of principal rotation axis (tomography rotation axis)

• X direction of the secondary rotation axis (tilt rotation axis)

• Y is beam direction

#### Bug in Mumott with reshaping the recon volume shape when transformations are applied  ##########
#---- User Input - Adjust this based on your sample ----#

print(data_container.stack.volume_shape)

# The volume shape needs to be restored here
data_container.stack.volume_shape = np.array([data_container.stack[0].diode.shape[0], *data_container.stack[0].diode.shape])

#---- End Input ----#

display(data_container)

[ 87 131  87]

DataContainer

Field	Size
Number of projections	419
Volume shape	[131 131 85]
Corrected for transmission	False
Angles in radians	True
Transformation applied	True
Reconstruction parameters generated	False

Step 3.1 - Create reconstruction parameters

• This is the input to the optimization algorithm and used for reconstructing the data. data_container.reconstruction_parameters

• In the case of having to change data (e.g. removing NaNs manually), this must be handled after the reconstruction parameters and by modifying them directly

• They are a flattened list, with the following ordering.. (to be added)

#---- Do not modify ----#
reconstruction_parameters = data_container.reconstruction_parameters
mask = reconstruction_parameters.projection_weights>0

print(f" A few checks of the data, having 0, NAN or Infs values in the data: \n"
      f" There are vaules of 0 in the masked reconstruction_parameters.data : {(reconstruction_parameters.data[mask]==0).any()} \n"
      f" There are Infs in the reconstruction_parameters.data : {np.isinf(reconstruction_parameters.data).any()} \n"
      f" There are Infs in the reconstruction_parameters.projection_weights : {np.isinf(reconstruction_parameters.projection_weights).any()} \n"
      f" There are Nans in the reconstruction_parameters.data : {np.isnan(reconstruction_parameters.data).any()} \n"
      f" There are Nans in the reconstruction_parameters.projection_weights : {np.isnan(reconstruction_parameters.projection_weights).any()}")

 A few checks of the data, having 0, NAN or Infs values in the data:
 There are vaules of 0 in the masked reconstruction_parameters.data : False
 There are Infs in the reconstruction_parameters.data : False
 There are Infs in the reconstruction_parameters.projection_weights : False
 There are Nans in the reconstruction_parameters.data : False
 There are Nans in the reconstruction_parameters.projection_weights : False

Step 3.2 Remove Nans in case it is needed

Example code below that should work is (need to be applied to both, data and weights):

mask = np.isnan(reconstruction_parameters.data)

reconstruction_parameters.data[mask] = np.nan_2_num(reconstruction_parameters.data[mask])

reconstruction_parameters.projection_weights[mask] = np.nan_2_num(reconstruction_parameters.data[mask])

for the mask array

mask = np.isnan(reconstruction_parameters.projection_weights)

reconstruction_parameters.data[mask] = np.nan_2_num(reconstruction_parameters.data[mask])

reconstruction_parameters.projection_weights[mask] = np.nan_2_num(reconstruction_parameters.data[mask])

Step 3.3 Check data via histogram of transmission and scattering

• Useful to also identify negative values in scattering that might be created during the background correction

• Since that might cause problems in the optimization, one may now directly access the reconstruction parameters and mask out those values! An example is given in the next cell

fig, [ax1, ax2] = plt.subplots(1,2, figsize = (8,3))
#print(reconstruction_parameters.projection_weights.shape)
#mask = reconstruction_parameters.data*reconstruction_parameters.projection_weights>0
mask = reconstruction_parameters.projection_weights>0
#vmin, vmax =
#print(mask.shape, (reconstruction_parameters.data*reconstruction_parameters.projection_weights).shape)
#print(np.max(reconstruction_parameters.data*reconstruction_parameters.projection_weights))
ax1.hist(np.log10((reconstruction_parameters.data[mask])),1000)
ax1.set_title('scattering')
ax2.hist(data_container._stack.diode, 1000)
ax2.set_title('diode')

Text(0.5, 1.0, 'diode')

Step 3.3.1 Remove outliers of the data

• Sometimes, scattering data can include some outliers. Althought it would be nice to understand why and where they come from, the most easy solution is to remove them from here.

• For this purpose, we look at the instensity distribution for the highest 5% of intensity

• Use threshold value or manually pick one from the plot in case you find a very high value towards 100.

threshold = 99.90

# ---- Do not modify this section ---- #
perc= np.percentile(reconstruction_parameters.data, np.linspace(95,100,100))

plt.figure()
plt.loglog(np.linspace(95,100,100), perc.T)

plt.axvline(threshold, ls = '--', color='r')
plt.show()

Step 3.3.2 Mask outliers

Use threshold or use manual value

thres = threshold# or e.g. 99.9

# ---- Do not modify this section ---- #
vmin, vmax = np.percentile(reconstruction_parameters.data, [0,thres])
# Set weights of these entries to 0
reconstruction_parameters.projection_weights[reconstruction_parameters.data>vmax] = 0

Double check if outliers were removed

# ---- Do not modify this section ---- #

mask = reconstruction_parameters.projection_weights > 0
perc= np.percentile(reconstruction_parameters.data[mask], np.linspace(95,100,100))

plt.figure()
plt.loglog(np.linspace(95,100,100), perc.T)

plt.axvline(threshold, ls = '--', color='r')

<matplotlib.lines.Line2D at 0x149bf0465d60>

Step 3.4 Setting Optimization Parameters

The OptimizationParameters object takes some inputs used to initialize some members of the ReconstructionParameters instance, as well as settings used by scipy.optimize.minimize; these are documented in the scipy.optimize documentation. The SIGTT-specific settings are:

1. integration_step_size This determines how large steps are taken to perform the integration necessary to compute projections; the optimal settings depend on various factors, such as how much variation there is within your sample and how many projections you have. The absolute maximum setting is 1, but generally a setting of 0.5 to 0.25 will yield better results and avoid artefacts. We will use 0.5.

2. maximum_order This is the maximum order of the reconstructed spherical polynomial of each voxel. This must be an even number, and should be no larger than 180 / detector_segment_width - 1. Since we have eight detector segments (segment width of 22.5 degrees), 6 is an appropriate setting. The total number of coefficients fitted for each order is (maximum_order / 2 + 1) * (maximum_order + 1), or 28 for a setting of 6.

3. initial_value The initial value of the isotropic component. The anisotropic components will be set to a fraction of this value. This will apply to all voxels, so it should not be set too large.

optimization_parameters = OptimizationParameters(
    reconstruction_parameters=reconstruction_parameters,
    integration_step_size=0.5,
    maximum_order=6,
    initial_value=1e-5,
    minimize_args=dict(method='L-BFGS-B'),
    minimize_options=dict(maxiter=25,  # for final reconstructions 40-50, initial ones at 25 is good
                          ftol=1e-5, # convergence criteria, input on useful values will help, 1e-5 typically good start, decreasing may help with convergence
                          gtol=1e-5, # tolerance for the largest absolute value in the gradient, convergence criteria for gradient
                          disp=1, # display updates
                          maxls=10, # maximum number of line search for stepsize
                          maxcor=3)) # how many terms in its approx of the Hess matrix (gradient matrix), momentum parameter of gradient

Increasing maximum order ...

Step 3.4.1 Setting save data - Creating the output Handler

• To create output, we need to create an OutputHandler object. We want to store our output in a folder called output. This can be relative to our execution directory of the current script, or given a basepath for the output, e.g. “…/reconstructions/sample_name/.” based on preferences

• It will overwrite existing files in case the argument is set to True. Enclosing folder gets overwritten?

output_handler = OutputHandler(
    reconstruction_parameters=reconstruction_parameters,
    base_output_path = path, # Change to absolute path of where output folder should be generated
    create_enclosing_folder = True, # overwrite enclosing folder, should be True to automatically create folders
    overwrite_existing_files = True, # Overwrite existing files, should be True to automatically create folders
    enclosing_folder_name = 'output', # filename for output folder
    output_file_name = f'{name}_output') # filename appendix, name is the dataset name as from the import

Output directory already exists!

Step 3.4.2 LiveViewHandler - Feedback on running reconstructions

We would like live feedback on our reconstruction, so we create a LiveViewHandler object (to get live feedback in a notebook, one can run it with %matplotlib widget). Apart from the reconstruction_parameters argument, the others are equal to the default arguments.

1. shown_orders shown_orders=[2, 4, 6] means that we will get a separate plot of the order-2, 4 and 6 contributions to the overall standard deviation of the spherical function.

2. plane plane=1 means that we will see a cut of the plane orthogonal to the Y-direction. cut_index='middle' means that we will see a central cut;

3. orientation orientation='transversal' means that we expect the orientation to be defined by fiber-like or great circle symmetry in the spherical polynomial, as opposed to 'longitudinal' or polar symmetry.

live_view_handler = LiveViewHandler(
    reconstruction_parameters=reconstruction_parameters,
    shown_orders=[0, 2],
    plane=0,
    cut_index= 'middle',
    orientation='transversal')

Step 3.4.3 Init the optimizer

Finally, we need to create the Optimizer object and pass on to it the OptimizationParameters instance (which also contain our ReconstructionParameters instance), as well as the LiveViewHandler instance.

#---- Do not modify ----#
optimizer = Optimizer(
    optimization_parameters=optimization_parameters,
    live_view_handler=live_view_handler)

Initializing optimizer...

Step 4 - Start optimization

We then need to explicitly run the optimization, and it is possible to e.g. modify the OptimizationParameters instance and resume the reconstruction again if we are not satisfied after the first run. We will get feedback on the optimization progress as well as live updates in the plot created by the LiveViewHandler instance.

#---- Do not modify ----#
optimizer.run_optimization()

Running optimization...
Updating plot...
Calculating residual...
End time forward: 9.56
Total projection time forward: 9.36
Residual norm: 1.0009e+08
End time adjoint: 13.36
Total projection time adjoint: 11.55
Gradient norm: 1.9266e+00
RUNNING THE L-BFGS-B CODE

           * * *

Machine precision = 2.220D-16
 N =     40843180     M =            3

At X0         0 variables are exactly at the bounds

At iterate    0    f=  1.00088D+08    |proj g|=  1.92661D+00
Iteration wall time: 33.81
Total wall time: 33.81
Updating plot...
Calculating residual...
End time forward: 10.51
Total projection time forward: 10.31
Residual norm: 1.0009e+08
End time adjoint: 12.88
Total projection time adjoint: 11.14
Gradient norm: 1.9266e+00
Iteration wall time: 25.25
Total wall time: 59.06
Updating plot...
Calculating residual...
End time forward: 10.11
Total projection time forward: 9.90
Residual norm: 1.0008e+08
End time adjoint: 13.34
Total projection time adjoint: 11.58
Gradient norm: 1.9265e+00
Iteration wall time: 25.32
Total wall time: 84.38
Updating plot...
Calculating residual...
End time forward: 10.92
Total projection time forward: 10.73
Residual norm: 1.0005e+08
End time adjoint: 13.35
Total projection time adjoint: 11.61
Gradient norm: 1.9262e+00
Iteration wall time: 26.16
Total wall time: 110.54
Updating plot...
Calculating residual...
End time forward: 10.25
Total projection time forward: 10.07
Residual norm: 9.9919e+07
End time adjoint: 13.00
Total projection time adjoint: 11.13
Gradient norm: 1.9248e+00
Iteration wall time: 25.07
Total wall time: 135.61
Updating plot...
Calculating residual...
End time forward: 10.15
Total projection time forward: 9.96
Residual norm: 9.9413e+07
End time adjoint: 14.24
Total projection time adjoint: 12.46
Gradient norm: 1.9194e+00
Iteration wall time: 26.18
Total wall time: 161.78
Updating plot...
Calculating residual...
End time forward: 10.77
Total projection time forward: 10.58
Residual norm: 9.8133e+07
End time adjoint: 14.15
Total projection time adjoint: 12.36
Gradient norm: 1.9056e+00

At iterate    1    f=  9.81331D+07    |proj g|=  1.90559D+00
Iteration wall time: 28.82
Total wall time: 190.60
Updating plot...
Calculating residual...
End time forward: 11.64
Total projection time forward: 11.44
Residual norm: 1.6468e+07
End time adjoint: 15.59
Total projection time adjoint: 13.64
Gradient norm: 2.6238e-01

At iterate    2    f=  1.64683D+07    |proj g|=  2.62379D-01
Iteration wall time: 31.74
Total wall time: 222.34
Updating plot...
Calculating residual...
End time forward: 11.04
Total projection time forward: 10.83
Residual norm: 1.2532e+07
End time adjoint: 14.40
Total projection time adjoint: 12.64
Gradient norm: 2.6171e-01

At iterate    3    f=  1.25316D+07    |proj g|=  2.17808D-01
Iteration wall time: 30.54
Total wall time: 252.88
Updating plot...
Calculating residual...
End time forward: 10.90
Total projection time forward: 10.71
Residual norm: 9.0713e+06
End time adjoint: 14.68
Total projection time adjoint: 12.79
Gradient norm: 2.2456e-01

At iterate    4    f=  9.07134D+06    |proj g|=  9.49750D-02
Iteration wall time: 30.72
Total wall time: 283.60
Updating plot...
Calculating residual...
End time forward: 11.31
Total projection time forward: 11.14
Residual norm: 7.8776e+06
End time adjoint: 13.82
Total projection time adjoint: 12.05
Gradient norm: 2.0078e-01

At iterate    5    f=  7.87761D+06    |proj g|=  9.75840D-02
Iteration wall time: 30.14
Total wall time: 313.74
Updating plot...
Calculating residual...
End time forward: 10.84
Total projection time forward: 10.67
Residual norm: 7.1433e+06
End time adjoint: 13.99
Total projection time adjoint: 12.20
Gradient norm: 1.7035e-01

At iterate    6    f=  7.14325D+06    |proj g|=  6.65929D-02
Iteration wall time: 29.87
Total wall time: 343.61
Updating plot...
Calculating residual...
End time forward: 12.03
Total projection time forward: 11.87
Residual norm: 6.7695e+06
End time adjoint: 13.69
Total projection time adjoint: 11.88
Gradient norm: 1.5184e-01

At iterate    7    f=  6.76947D+06    |proj g|=  8.42317D-02
Iteration wall time: 30.68
Total wall time: 374.29
Updating plot...
Calculating residual...
End time forward: 11.65
Total projection time forward: 11.48
Residual norm: 6.3173e+06
End time adjoint: 13.74
Total projection time adjoint: 11.93
Gradient norm: 1.3724e-01

At iterate    8    f=  6.31731D+06    |proj g|=  6.55218D-02
Iteration wall time: 30.45
Total wall time: 404.74
Updating plot...
Calculating residual...
End time forward: 11.30
Total projection time forward: 11.14
Residual norm: 5.9555e+06
End time adjoint: 13.58
Total projection time adjoint: 11.53
Gradient norm: 1.0716e-01

At iterate    9    f=  5.95553D+06    |proj g|=  7.35670D-02
Iteration wall time: 30.52
Total wall time: 435.25
Updating plot...
Calculating residual...
End time forward: 12.12
Total projection time forward: 11.90
Residual norm: 5.6530e+06
End time adjoint: 13.84
Total projection time adjoint: 11.91
Gradient norm: 1.0410e-01

At iterate   10    f=  5.65301D+06    |proj g|=  4.56420D-02
Iteration wall time: 31.56
Total wall time: 466.82
Updating plot...
Calculating residual...
End time forward: 11.86
Total projection time forward: 11.69
Residual norm: 5.4469e+06
End time adjoint: 14.27
Total projection time adjoint: 12.29
Gradient norm: 9.1587e-02

At iterate   11    f=  5.44693D+06    |proj g|=  3.90272D-02
Iteration wall time: 31.79
Total wall time: 498.61
Updating plot...
Calculating residual...
End time forward: 11.88
Total projection time forward: 11.69
Residual norm: 5.2405e+06
End time adjoint: 14.16
Total projection time adjoint: 12.25
Gradient norm: 7.8012e-02

At iterate   12    f=  5.24050D+06    |proj g|=  4.75804D-02
Iteration wall time: 31.43
Total wall time: 530.04
Updating plot...
Calculating residual...
End time forward: 12.00
Total projection time forward: 11.79
Residual norm: 5.1052e+06
End time adjoint: 14.03
Total projection time adjoint: 12.11
Gradient norm: 7.4624e-02

At iterate   13    f=  5.10520D+06    |proj g|=  7.46236D-02
Iteration wall time: 31.36
Total wall time: 561.40
Updating plot...
Calculating residual...
End time forward: 12.46
Total projection time forward: 12.24
Residual norm: 4.9611e+06
End time adjoint: 14.23
Total projection time adjoint: 12.23
Gradient norm: 5.0312e-02

At iterate   14    f=  4.96108D+06    |proj g|=  3.45490D-02
Iteration wall time: 32.28
Total wall time: 593.69
Updating plot...
Calculating residual...
End time forward: 11.77
Total projection time forward: 11.58
Residual norm: 4.8522e+06
End time adjoint: 13.72
Total projection time adjoint: 11.70
Gradient norm: 5.3350e-02

At iterate   15    f=  4.85215D+06    |proj g|=  2.58477D-02
Iteration wall time: 31.74
Total wall time: 625.43
Updating plot...
Calculating residual...
End time forward: 13.16
Total projection time forward: 13.01
Residual norm: 4.7455e+06
End time adjoint: 14.48
Total projection time adjoint: 12.47
Gradient norm: 5.5942e-02

At iterate   16    f=  4.74554D+06    |proj g|=  2.90913D-02
Iteration wall time: 33.20
Total wall time: 658.63
Updating plot...
Calculating residual...
End time forward: 11.66
Total projection time forward: 11.48
Residual norm: 4.5907e+06
End time adjoint: 14.16
Total projection time adjoint: 12.31
Gradient norm: 4.8686e-02

At iterate   17    f=  4.59067D+06    |proj g|=  4.74387D-02
Iteration wall time: 31.15
Total wall time: 689.78
Updating plot...
Calculating residual...
End time forward: 12.24
Total projection time forward: 12.06
Residual norm: 4.5604e+06
End time adjoint: 14.01
Total projection time adjoint: 11.88
Gradient norm: 5.8901e-02

At iterate   18    f=  4.56041D+06    |proj g|=  5.86664D-02
Iteration wall time: 32.18
Total wall time: 721.96
Updating plot...
Calculating residual...
End time forward: 11.88
Total projection time forward: 11.68
Residual norm: 4.4159e+06
End time adjoint: 16.23
Total projection time adjoint: 12.41
Gradient norm: 4.9168e-02

At iterate   19    f=  4.41595D+06    |proj g|=  2.56139D-02
Iteration wall time: 33.63
Total wall time: 755.59
Updating plot...
Calculating residual...
End time forward: 11.27
Total projection time forward: 11.06
Residual norm: 4.3536e+06
End time adjoint: 14.21
Total projection time adjoint: 12.07
Gradient norm: 4.3051e-02

At iterate   20    f=  4.35361D+06    |proj g|=  1.89310D-02
Iteration wall time: 31.27
Total wall time: 786.86
Updating plot...
Calculating residual...
End time forward: 11.69
Total projection time forward: 11.53
Residual norm: 4.2511e+06
End time adjoint: 13.93
Total projection time adjoint: 11.96
Gradient norm: 3.5202e-02

At iterate   21    f=  4.25110D+06    |proj g|=  2.86460D-02
Iteration wall time: 30.85
Total wall time: 817.71
Updating plot...
Calculating residual...
End time forward: 11.63
Total projection time forward: 11.44
Residual norm: 4.1729e+06
End time adjoint: 13.96
Total projection time adjoint: 12.04
Gradient norm: 5.4647e-02

At iterate   22    f=  4.17290D+06    |proj g|=  5.46471D-02
Iteration wall time: 30.86
Total wall time: 848.57
Updating plot...
Calculating residual...
End time forward: 12.00
Total projection time forward: 11.82
Residual norm: 4.0716e+06
End time adjoint: 14.09
Total projection time adjoint: 12.19
Gradient norm: 2.9866e-02

At iterate   23    f=  4.07161D+06    |proj g|=  2.98662D-02
Iteration wall time: 31.43
Total wall time: 880.00
Updating plot...
Calculating residual...
End time forward: 12.03
Total projection time forward: 11.86
Residual norm: 3.9895e+06
End time adjoint: 13.53
Total projection time adjoint: 11.65
Gradient norm: 2.2883e-02

At iterate   24    f=  3.98955D+06    |proj g|=  1.96340D-02
Iteration wall time: 31.39
Total wall time: 911.39
Updating plot...
Calculating residual...
End time forward: 11.59
Total projection time forward: 11.39
Residual norm: 3.9059e+06
End time adjoint: 14.63
Total projection time adjoint: 12.51
Gradient norm: 2.7173e-02

At iterate   25    f=  3.90594D+06    |proj g|=  2.55145D-02

           * * *

Tit   = total number of iterations
Tnf   = total number of function evaluations
Tnint = total number of segments explored during Cauchy searches
Skip  = number of BFGS updates skipped
Nact  = number of active bounds at final generalized Cauchy point
Projg = norm of the final projected gradient
F     = final function value

           * * *

   N    Tit     Tnf  Tnint  Skip  Nact     Projg        F
*****     25     31 425558     0 *****   2.551D-02   3.906D+06
  F =   3905937.5014295108

STOP: TOTAL NO. of ITERATIONS REACHED LIMIT

Step 4.1 - Check results

• This is a quick and dirty hacked way to access the data for the reconstruction along different planes. Visualization might be a bit ackward (partly duplicated plot), but it works fine for the moment: 0: yz plane ; 1: xz plane ; 2: xy plane

• You can scroll through the tomogram and see if results match your expectations

• For details on representation, you can type LiveViewHandler? in a new tab - the widget/scrolling option is the hacked part

# ----- modify here which which plane to look at ----- #
plane_chosen = 0
data_container.stack.volume_shape

#---- Do not modify ----#
# Hacked version to access parameters from the live_view_handler of getting information for slices along different planes through the reconstruction
live_view_handler = LiveViewHandler(
    reconstruction_parameters=reconstruction_parameters,
    shown_orders=[0, 2],
    plane=plane_chosen,
    cut_index= data_container.stack.volume_shape[plane_chosen]//2,
    orientation='longitudinal')
#display(live_view_handler.figure)

@interact(slice_cut=(2,data_container.stack.volume_shape[plane_chosen]-1))
def plot(slice_cut):
    live_view_handler._cut_index = slice_cut
    live_view_handler.update_plots()
    #display(live_view_handler.figure)

We now run the save_output method of our OutputHandler instance.

Step 5 - Save output to file

This stores the results of the reconstruction to a file. It takes the input parameters as given within the initialization of the output_handler

#---- Do not modify ----#
# Please add here a printout where it saves the data to
#print(f"Saving data to file {os.path.join(output_handler._base_output_path, output_handler._output_file_name)}")
output_handler.save_output()

/das/home/appel_c/.local/lib/python3.9/site-packages/mumott/output_handling/output_handler.py:222: RuntimeWarning: divide by zero encountered in divide
  return np.nan_to_num(self._std_of_amplitude() / self._mean_of_amplitude(),

Saving coefficients to h5 ...
Saving rms_of_amplitude to h5 ...
Saving mean_of_amplitude to h5 ...
Saving std_of_amplitude to h5 ...
Saving relative_std_of_amplitude to h5 ...
Saving eigenvectors to h5 ...
Saving eigenvalues to h5 ...
Saving order_amplitude to h5 ...
Saving rank_2_tensor to h5 ...
Saving synthetic_data to h5 ...
reconstruction_type
algorithm

print(output_handler._output_file_name)

dataset_q_0.030_0.037.h5_output