Merge pull request #18 from jorgenem/master

Going public with version 0.2

README.md

@@ -1,8 +1,8 @@
-# oslo_method_python
+# OMpy
 
-NB! This repo is currently in testing. Most things do not work correctly. Use at your own risk!
+NB! OMpy is currently in beta. Use at your own risk!
 
-This is a python implementation of the Oslo method. The goal is that this will become a complete Oslo method ecosystem within Python, which can also be interfaced with other Python packages to do e.g. Bayesian inference.
+This is `ompy`, a python implementation of the Oslo method. The goal is that this will become a complete Oslo method ecosystem within Python, which can also be interfaced with other Python packages to do e.g. Bayesian inference.
 
 
 ## How to use this package
@@ -68,11 +68,82 @@ An implementation of the first generation method present in Guttormsen, Ramsøy
 Like the unfolding function, please be on the lookout for bugs...
 
 ### The MatrixAnalysis class
-The MatrixAnalysis class is a convenience wrapper for the `unfold()` and `first_generation_method()` functions, along with some utility functions to use on the spectra. You do not have to use it except if you want to do error propagation. In that case, the parameters to use for unfolding and first generation are stored in the MatrixAnalysis instance that you pass to the ErrorPropagation instance, to ensure that the ensemble of perturbed copies of the spectra are treated identically.
+The steps described above make up the core matrix manipulation routines of `OMpy`. We therefore provide a class which helps to streamline the process, called `MatrixAnalysis`.
+```
+ma = ompy.MatrixAnalysis()
+# Load raw matrix:
+ma.raw.load(fname_raw)
+# Unfold the raw matrix and place it
+# in the variable ma.unfolded:
+ma.unfold(*args)
+# Apply first generation method 
+# and get result in ma.firstgen:
+ma.first_generation_method(*args)
+```
+The class methods `ma.unfold()` and `ma.first_generation_method()` accept the same input arguments as their respective stand-alone functions. The advantage of the `MatrixAnalysis` class is that all parameter choices for the unfolding and first-generation method are stored within the class instance. For this reason, the `MatrixAnalysis` class is necessary when one wants to do uncertainty propagation, to ensure that the same settings are applied to each member of the ensemble. 
 
 ## The ErrorPropagation class
-Propagate statistical errors through the Oslo method by making an ensemble of perturbed copies of the input spectrum. You can choose between gaussian or poisson errors. 
+To perform the uncertainty propagation, instantiate the class `ErrorPropagation` with the instance of `MatrixAnalysis` as an input argument:
+```
+ep = ompy.ErrorPropagation(ma)
+```
+Then, an ensemble of `N_ensemble_members` is generated by the command
+```
+ep.generate_ensemble(N_ensemble_members)
+```
+The standard deviations matrices for `raw`, `unfolded` and `firstgen` are automatically calculated and put in the variables `ep.std_raw`, `ep.std_unfolded` and `ep.std_firstgen`, respectively.
+
+## The FitRhoT class
+To fit $\rho$ and $\mathcal{T}$ to the first-generation matrix (also requiring `ep.std_firstgen`), do
+```
+# Instantiate the fitting class:
+fitting = ompy.FitRhoT(ma.firstgen,
+                       ep.std_firstgen,
+                       bin_width_out,
+                       Ex_min, Ex_max,
+                       Eg_min
+                       )
+# Perform the fit:
+fitting.fit()
+# The fitted functions are available 
+# as Vector instances:
+fitting.rho, fitting.T
+```
+Here, `bin_width_out` gives the energy bin width of the resulting fit, and the other arguments determine the region of the first-generation matrix to fit to. To get the uncertainty on the fitted quantities, one can run the fitting for each member in the ensemble:
+```
+# Allocate arrays to store each ensemble
+# member fit:
+rho_ens = np.zeros((N_ensemble_fit,
+                    len(rho.vector)))
+T_ens = np.zeros((N_ensemble_fit,
+                  len(T.vector)))
+# As a trick, we copy the instance ma
+# and replace its matrix every iteration:
+import copy
+ma_curr = copy.deepcopy(ma)
+# Loop through all and perform fit:
+for i_ens in range(N_ensemble_fit):
+    ma_curr.firstgen.matrix = \
+        ep.firstgen_ensemble[i_ens, :, :]
+    fitter_curr = ompy.FitRhoT(
+                    ma_curr.firstgen,
+                    ep.std_firstgen,
+                    bin_width_out,
+                    Ex_min, Ex_max,
+                    Eg_min
+                    )
+    rho_curr= fitter_curr.rho
+    T_curr = fitter_curr.T
+    rho_ens[i_ens, :] = rho_curr.vector
+    T_ens[i_ens, :] = T_curr.vector
+```
+
+## Adding new features
+OMpy is written with modularity in mind. We want it to be as easy as possible for the user to add custom functionality and interface OMpy with other Python packages. For example, it may be of interest to try other unfolding algorithms than the one presently implemented. To achieve this, one just has to write a wrapper function that has the same input and output structure as the function `unfold()`, found in the file `ompy/unfold.py`, and replace the calls to `unfold()` by the custom function. We hope to make such modifications even easier in a future version.
+
+It is our hope and goal that `OMpy` will be used, and we are happy to provide support. Feedback and suggestions are also very welcome. We encourage users who implement new features to share them by opening a pull request in the Github repository.
+
+
 
-[Write something about the probablity theory and assumptions. Also about how the spectra are stored and how the variance matrix is calculated at the end. Should also have some examples throughout.]
 
 <sup>1</sup>It stands for alfa-natrium, and comes from the fact that the outgoing projectile hitting the excitation energy detector SiRi used to be exclusively alpha particles and that the gamma-detector array CACTUS consisted of NaI scintillators.

oslo_method_python/__init__.py → ompy/__init__.py

@@ -8,3 +8,4 @@
 from .first_generation_method import *
 from .error_propagation import *
 from .rhosig import *
+from .spinfunctions import *

oslo_method_python/error_propagation.py → ompy/error_propagation.py

@@ -187,17 +187,22 @@ def generate_ensemble(self,
         # axs_fg[0,0].set_title("fg")
         # END TESTING
 
-        # Allocate a cube array to store all firstgen ensemble members. We need them to make the firstgen variance matrix.
+        # Allocate a cube array to store all ensemble members. We need them to make the std matrix.
         raw_ensemble = np.zeros((N_ensemble_members,
                                 matrix_analysis.raw.matrix.shape[0],
                                 matrix_analysis.raw.matrix.shape[1])
                                 )
+        unfolded_ensemble = np.zeros((N_ensemble_members,
+                                      matrix_analysis.unfolded.matrix.shape[0],
+                                      matrix_analysis.unfolded.matrix.shape[1])
+                                     )
         firstgen_ensemble = np.zeros((N_ensemble_members,
                                       matrix_analysis.firstgen.matrix.shape[0],
                                       matrix_analysis.firstgen.matrix.shape[1])
                                      )
 
-        # Loop over and generate the random perturbations, then unfold and first-generation-method them:
+        # Loop over and generate the random perturbations, then unfold and
+        # first-generation-method them:
         for i in range(N_ensemble_members):
             if verbose:
                 print("=== Begin ensemble member ", i, " ===", flush=True)
@@ -206,9 +211,9 @@ def generate_ensemble(self,
 
             # Allocate matrix_analysis instance for current ensemble member:
             ma_curr = copy.deepcopy(matrix_analysis)
-            # Check if the current ensemble member exists on file, create it if not:
+            # Check if the current ensemble member exists on file, create it if
+            # not:
             if os.path.isfile(fname_raw_current) and not purge_files:
-                # data_raw_ensemblemember, tmp, tmp, tmp = read_mama_2D(fname_raw_current)
                 ma_curr.raw.load(fname_raw_current, suppress_warning=True)
                 if verbose:
                     print("Read raw matrix from file", flush=True)
@@ -274,6 +279,9 @@ def generate_ensemble(self,
                 # if verbose:
                     # print("unfolded_ensemblemember.shape =", unfolded_ensemblemember.shape, flush=True)
 
+            # Store unfolded ensemble member in memory:
+            unfolded_ensemble[i, :, :] = ma_curr.unfolded.matrix
+
             # === Extract first generation spectrum ===:
             # Ex_max = 12000 # keV - maximum excitation energy
             # dE_gamma = 1000 # keV - allow gamma energy to exceed excitation energy by this much, to account for experimental resolution
@@ -322,6 +330,14 @@ def generate_ensemble(self,
                          )
         fname_raw_std = os.path.join(folder, "raw_std.m")
         std_raw.save(fname_raw_std)
+        # unfolded:
+        unfolded_ensemble_std = np.std(unfolded_ensemble, axis=0)
+        std_unfolded = Matrix(unfolded_ensemble_std,
+                         matrix_analysis.unfolded.E0_array,
+                         matrix_analysis.unfolded.E1_array
+                         )
+        fname_unfolded_std = os.path.join(folder, "unfolded_std.m")
+        std_unfolded.save(fname_unfolded_std)
         # firstgen:
         firstgen_ensemble_std = np.std(firstgen_ensemble, axis=0)
         std_firstgen = Matrix(firstgen_ensemble_std,
@@ -336,7 +352,14 @@ def generate_ensemble(self,
         # END TESTING
 
         self.std_firstgen = std_firstgen
+        self.std_unfolded = std_unfolded
         self.std_raw = std_raw
         # TODO add self.std_unfolded
 
+        # Also store a list containing all ensemble members of firstgen:
+        self.firstgen_ensemble = firstgen_ensemble
+        # self.firstgen_ensemble = []
+        # for i_ens in range(N_ensemble_members):
+
+
         # return std_firstgen

oslo_method_python/first_generation_method.py → ompy/first_generation_method.py

@@ -276,6 +276,11 @@ def first_generation_method(matrix_in,
         Eg_mesh, n_f) , np.power(Ef_mesh, 2) * np.exp(2 * np.sqrt(a_f * Ef_mesh / 1000))), 0)
     W_old = div0(W_old, W_old.sum(axis=1).reshape(N_Exbins, 1))
 
+    # Update 20190404: Tested using a box instead of Fermi gas, seems to work just as well.
+    # TODO: Should probably use this instead after double-checking it.
+    # W_old = np.ones_like(Eg_mesh)
+    # W_old = div0(W_old, W_old.sum(axis=1).reshape(N_Exbins, 1))
+
     mask_W = make_mask(Ex_array, Ex_array, Ex_array[0], Ex_array[
                        0] + dE_gamma, Ex_array[-1], Ex_array[-1] + dE_gamma)