extend fission avg mapping to consider spectrum uncertainties

gmapi/mappings/cross_section_fission_average_map.py

@@ -1,4 +1,5 @@
 import numpy as np
+from scipy.sparse import csr_matrix
 from .basic_integral_maps import (basic_integral_propagate,
         basic_integral_of_product_propagate, get_basic_integral_of_product_sensmats)
 from .helperfuns import return_matrix
@@ -168,7 +169,7 @@ def __modern_compute(self, datatable, refvals, what):
 
         priormask = (datatable['REAC'].str.match('MT:1-R1:', na=False) &
                      datatable['NODE'].str.match('xsid_', na=False))
-        priormask = np.logical_or(priormask, datatable['NODE'].str.fullmatch('fis_modern', na=False))
+        priormask = np.logical_or(priormask, datatable['NODE'].str.fullmatch('fis_modern|fis_errors', na=False))
         priortable = datatable[priormask]
 
         expmask = self.is_responsible(datatable)
@@ -185,23 +186,30 @@ def __modern_compute(self, datatable, refvals, what):
         if not np.isclose(normfact, 1., rtol=1e-4):
             raise ValueError('fission spectrum not normalized')
 
+        priortable_fiserror = priortable[priortable['NODE'].str.fullmatch(
+                                 'fis_errors', na=False)]
+        ensfis_error = priortable_fiserror['ENERGY'].to_numpy()
+        S = self.__make_covmapping(ensfis, ensfis_error, valsfis)
+
         for curexp in expids:
             exptable_red = exptable[exptable['NODE'].str.fullmatch(curexp, na=False)]
             if len(exptable_red) != 1:
                 raise IndexError('None or more than one rows associated with a ' +
                         'fission average, which must not happen!')
             curreac = exptable_red['REAC'].values[0]
             preac = curreac.replace('MT:6-', 'MT:1-')
-            priortable_red = priortable[priortable['REAC'].str.fullmatch(preac, na=False)]
+            priortable_reac = priortable[priortable['REAC'].str.fullmatch(preac, na=False)]
             # abbreviate some variables
-            ens1 = priortable_red['ENERGY'].to_numpy()
-            vals1 = refvals[priortable_red.index]
-            idcs1red = priortable_red.index
+            ens1 = priortable_reac['ENERGY'].to_numpy()
+            vals1_a = refvals[priortable_reac.index]
+            vals1_b = S @ refvals[priortable_fiserror.index]
+            idcs1red_a = priortable_reac.index
+            idcs1red_b = priortable_fiserror.index
             idcs2red = exptable_red.index
             if what == 'propagate':
 
                 curval = basic_integral_of_product_propagate(
-                        [ens1, ensfis], [vals1, valsfis],
+                        [ens1, ensfis], [vals1_a, valsfis+vals1_b],
                         ['lin-lin', 'lin-lin'], zero_outside=True,
                         maxord=16, rtol=1e-6)
                 curval = float(curval) * normfact
@@ -211,17 +219,23 @@ def __modern_compute(self, datatable, refvals, what):
 
             elif what == 'jacobian':
                 sensvecs = get_basic_integral_of_product_sensmats(
-                        [ens1, ensfis], [vals1, valsfis],
+                        [ens1, ensfis], [vals1_a, valsfis+vals1_b],
                         ['lin-lin', 'lin-lin'], zero_outside=True,
                         maxord=16, rtol=1e-6)
-                # because we assume that the fission spectrum is constant
-                # we can ignore the sensitivity to it
-                sensvec = np.ravel(sensvecs[0])
-                sensvec *= normfact
 
-                idcs1 = concat([idcs1, idcs1red])
-                idcs2 = concat([idcs2, np.full(len(idcs1red), idcs2red, dtype=int)])
-                coeff = concat([coeff, sensvec])
+                sensvec_a = np.ravel(sensvecs[0])
+                sensvec_a *= normfact
+
+                tmp = np.ravel(sensvecs[1])
+                tmp *= normfact
+                sensvec_b = np.ravel(tmp.reshape(1,-1) @ S)
+
+                # forward propagation
+
+                idcs1 = concat([idcs1, idcs1red_a, idcs1red_b])
+                idcs2 = concat([idcs2, np.full(len(idcs1red_a)+len(idcs1red_b),
+                                               idcs2red, dtype=int)])
+                coeff = concat([coeff, sensvec_a, sensvec_b])
 
             else:
                 raise ValueError('what must be either "propagate" or "jacobian"')
@@ -236,3 +250,25 @@ def __modern_compute(self, datatable, refvals, what):
             retdic['propvals'] = propvals
 
         return retdic
+
+
+    def __make_covmapping(self, ens, cov_ens, scl=None):
+        if not np.all(np.sort(cov_ens) == cov_ens):
+            raise ValueError('cov_ens must be sorted')
+        if scl is None:
+            scl = np.full(len(ens), 1.)
+        en_idcs = np.searchsorted(cov_ens, ens)
+        assert en_idcs[0] == 0
+        en_idcs[0] = 1
+        en_idcs -= 1
+        i_list = []
+        j_list = []
+        v_list = []
+        n = len(ens)
+        m = len(cov_ens)
+        for i in range(n):
+            i_list.append(i)
+            j_list.append(en_idcs[i])
+            v_list.append(scl[i])
+        S = csr_matrix((v_list, (i_list, j_list)), shape=(n, m), dtype='d')
+        return S