Detector class bug fixes (#743)

* increase numba requirement

* Fixes for test notebook

* Update detectors class and add test notebook

* Make it possible to weight events

* adjust hypersurfaces for weighted events

* Update hist.py

* update multi det notebook

* Add update_param_values function for Detector class

* add test script for Detectors class

* Add Detectors test and (somewhat) fix DistributionMaker test

* only change bin uncertainties if they exist

* make genie sys stage more generic

* more flexibility when splitting interactions

* Handle __builtins__ both as dict and as module

* cythonize

* move cython import

* remove cython (seems like I do not know what I am doing)

* NuFit 5.2 and small bug fixes

* Detectors class bug fix

* Add prior for deltacp

* initialize params after changing selection

* move param init into function

---------

Co-authored-by: Jan Weldert <jweldert@cobalt01.icecube.wisc.edu>

Author: JanWeldert

pisa/analysis/analysis.py

@@ -375,17 +375,21 @@ def update_param_values_detector(
     """
     Modification of the update_param_values function to use with the Detectors class.
     """
+    assert hypo_maker.__class__.__name__ == "Detectors", "hypo_maker is not Detectors class"   
+    
+    if isinstance(params, Param): params = ParamSet(params)
+        
     for distribution_maker in hypo_maker:
-        update_param_values(distribution_maker, params)
-
-    if isinstance(params, Param): params = ParamSet(params) # just for the following
-
-    for p in params.names: # now update params with det_names inside
-        for i, det_name in enumerate(hypo_maker.det_names):
-            if det_name in p:
-                cp = deepcopy(params[p])
-                cp.name = cp.name.replace('_'+det_name, "")
-                update_param_values(hypo_maker._distribution_makers[i], cp)
+        ps = deepcopy(params)
+        for p in ps.names:
+            if distribution_maker.detector_name in p:
+                p_name = p.replace('_'+distribution_maker.detector_name, "")
+                if p_name in ps.names:
+                    ps.remove(p_name)
+                ps[p].name = p_name
+        update_param_values(distribution_maker, ps, 
+                            update_nominal_values, update_range, update_is_fixed)
+    hypo_maker.init_params()
 
 # TODO: move this to a central location prob. in utils
 class Counter(object):
@@ -1166,7 +1170,10 @@ def _fit_octants(self, data_dist, hypo_maker, metric, external_priors_penalty,
         # Copy the fitted parameter values from the best fit case into the hypo maker's
         # parameter values.
         # Also reinstate the original parameter range for the angle
-        update_param_values(hypo_maker, best_fit_info.params.free, update_range=True)
+        if hypo_maker.__class__.__name__ == "Detectors":
+            update_param_values_detector(hypo_maker, best_fit_info.params.free, update_range=True)
+        else:
+            update_param_values(hypo_maker, best_fit_info.params.free, update_range=True)
 
         return best_fit_info
 
@@ -1310,7 +1317,10 @@ def _fit_grid_scan(self, data_dist, hypo_maker, metric,
                     mod_param.is_fixed = True
                 # It is important not to use hypo_maker.update_params(mod_param) here
                 # because we don't want to overwrite the memory reference!
-                update_param_values(hypo_maker, mod_param, update_is_fixed=True)
+                if hypo_maker.__class__.__name__ == "Detectors":
+                    update_param_values_detector(hypo_maker, mod_param, update_is_fixed=True)
+                else:
+                    update_param_values(hypo_maker, mod_param, update_is_fixed=True)
             new_fit_info = self.fit_recursively(
                 data_dist, hypo_maker, metric, external_priors_penalty,
                 local_fit_kwargs["method"], local_fit_kwargs["method_kwargs"],
@@ -1338,7 +1348,10 @@ def _fit_grid_scan(self, data_dist, hypo_maker, metric,
         best_fit_result = all_fit_results[best_idx]
 
         if do_refined_fit:
-            update_param_values(hypo_maker, best_fit_result.params.free)
+            if hypo_maker.__class__.__name__ == "Detectors":
+                update_param_values_detector(hypo_maker, best_fit_result.params.free)
+            else:
+                update_param_values(hypo_maker, best_fit_result.params.free)
             # the params stored in the best fit may come from a grid point where
             # parameters were fixed, so we free them up again
             for param in originally_free:
@@ -1507,9 +1520,14 @@ def _fit_ranges(self, data_dist, hypo_maker, metric,
                 logging.info(f"parameter with modified range:\n{mod_param}")
             # use update_param_values instead of hypo_maker.update_params so that we
             # don't overwrite the internal memory reference
-            update_param_values(
-                hypo_maker, mod_param, update_range=True, update_nominal_values=True
-            )
+            if hypo_maker.__class__.__name__ == "Detectors":
+                update_param_values_detector(
+                    hypo_maker, mod_param, update_range=True, update_nominal_values=True
+                )
+            else:
+                update_param_values(
+                    hypo_maker, mod_param, update_range=True, update_nominal_values=True
+                )
             fit_result = self.fit_recursively(
                 data_dist, hypo_maker, metric, external_priors_penalty,
                 local_fit_kwargs["method"], local_fit_kwargs["method_kwargs"],
@@ -1536,10 +1554,16 @@ def _fit_ranges(self, data_dist, hypo_maker, metric,
         # set the values of all parameters in the hypo_maker to the best fit values
         # without overwriting the memory reference.
         # Also reset ranges and nominal values that we might have changed above!
-        update_param_values(
-            hypo_maker, best_fit_result.params.free,
-            update_range=True, update_nominal_values=True
-        )
+        if hypo_maker.__class__.__name__ == "Detectors":
+            update_param_values_detector(
+                hypo_maker, best_fit_result.params.free,
+                update_range=True, update_nominal_values=True
+            )
+        else:
+            update_param_values(
+                hypo_maker, best_fit_result.params.free,
+                update_range=True, update_nominal_values=True
+            )
         return best_fit_result
 
     def _fit_staged(self, data_dist, hypo_maker, metric,
@@ -1565,9 +1589,14 @@ def _fit_staged(self, data_dist, hypo_maker, metric,
         for i, fit_kwargs in enumerate(local_fit_kwargs):
             logging.info(f"Beginning fit {i+1} / {len(local_fit_kwargs)}")
             if best_fit_params is not None:
-                update_param_values(
-                    hypo_maker, best_fit_params.free, update_nominal_values=True
-                )
+                if hypo_maker.__class__.__name__ == "Detectors":
+                    update_param_values_detector(
+                        hypo_maker, best_fit_params.free, update_nominal_values=True
+                    )
+                else:
+                    update_param_values(
+                        hypo_maker, best_fit_params.free, update_nominal_values=True
+                    )
             best_fit_info = self.fit_recursively(
                 data_dist, hypo_maker, metric, external_priors_penalty,
                 fit_kwargs["method"], fit_kwargs["method_kwargs"],
@@ -1588,9 +1617,14 @@ def _fit_staged(self, data_dist, hypo_maker, metric,
         best_fit_info._rehash()
         # Make sure that the hypo_maker has its params also at the best fit point
         # with the original nominal parameter values.
-        update_param_values(
-            hypo_maker, best_fit_info.params.free, update_nominal_values=True
-        )
+        if hypo_maker.__class__.__name__ == "Detectors":
+            update_param_values_detector(
+                hypo_maker, best_fit_info.params.free, update_nominal_values=True
+            )
+        else:
+            update_param_values(
+                hypo_maker, best_fit_info.params.free, update_nominal_values=True
+            )
         return best_fit_info
 
     def _fit_scipy(self, data_dist, hypo_maker, metric,
@@ -1918,6 +1952,8 @@ def annealing_callback(x, f, context):
         rescaled_pvals = optimize_result.pop('x')
         rescaled_pvals = np.where(flip_x0, 1 - rescaled_pvals, rescaled_pvals)
         hypo_maker._set_rescaled_free_params(rescaled_pvals) # pylint: disable=protected-access
+        if hypo_maker.__class__.__name__ == "Detectors":
+            update_param_values_detector(hypo_maker, hypo_maker.params.free) #updates values for ALL detectors
 
         # Get the best-fit metric value
         metric_val = sign * optimize_result.pop('fun')
@@ -1944,6 +1980,8 @@ def annealing_callback(x, f, context):
             # Reset to starting value of the fit, rather than nominal values because
             # the nominal value might be out of range if this is inside an octant check.
             hypo_maker._set_rescaled_free_params(x0)
+            if hypo_maker.__class__.__name__ == "Detectors":
+                update_param_values_detector(hypo_maker, hypo_maker.params.free) #updates values for ALL detectors
 
         # TODO: other metrics
         fit_info = HypoFitResult(
@@ -2115,6 +2153,8 @@ def loss_func(x):
         # values from [0,1] back to physical range)
         rescaled_pvals = np.array(m.values)
         hypo_maker._set_rescaled_free_params(rescaled_pvals) # pylint: disable=protected-access
+        if hypo_maker.__class__.__name__ == "Detectors":
+            update_param_values_detector(hypo_maker, hypo_maker.params.free) #updates values for ALL detectors
 
         # Get the best-fit metric value
         metric_val = sign * m.fval
@@ -2153,6 +2193,8 @@ def loss_func(x):
             # Reset to starting value of the fit, rather than nominal values because
             # the nominal value might be out of range if this is inside an octant check.
             hypo_maker._set_rescaled_free_params(x0)
+            if hypo_maker.__class__.__name__ == "Detectors":
+                update_param_values_detector(hypo_maker, hypo_maker.params.free) #updates values for ALL detectors
 
         # TODO: other metrics
         fit_info = HypoFitResult(
@@ -2286,6 +2328,8 @@ def loss_func(x, grad):
         # values from [0,1] back to physical range)
         rescaled_pvals = xopt
         hypo_maker._set_rescaled_free_params(rescaled_pvals) # pylint: disable=protected-access
+        if hypo_maker.__class__.__name__ == "Detectors":
+            update_param_values_detector(hypo_maker, hypo_maker.params.free) #updates values for ALL detectors
 
         # Get the best-fit metric value
         metric_val = sign * opt.last_optimum_value()
@@ -2322,6 +2366,8 @@ def loss_func(x, grad):
 
         if self.blindness > 1:  # only at stricter blindness level
             hypo_maker._set_rescaled_free_params(x0)
+            if hypo_maker.__class__.__name__ == "Detectors":
+                update_param_values_detector(hypo_maker, hypo_maker.params.free) #updates values for ALL detectors
 
         # TODO: other metrics
         fit_info = HypoFitResult(
@@ -2403,6 +2449,8 @@ def ineq_func(x, grad):
                     if grad.size > 0:
                         raise RuntimeError("gradients not supported")
                     hypo_maker._set_rescaled_free_params(x)
+                    if hypo_maker.__class__.__name__ == "Detectors":
+                        update_param_values_detector(hypo_maker, hypo_maker.params.free) #updates values for ALL detectors
                     # In NLOPT, the inequality function must stay negative, while in
                     # scipy, the inequality function must stay positive. We keep with
                     # the scipy convention by flipping the sign.
@@ -2518,6 +2566,8 @@ def _minimizer_callable(self, scaled_param_vals, hypo_maker, data_dist,
         scaled_param_vals = np.where(flip_x0, 1 - scaled_param_vals, scaled_param_vals)
         # Set param values from the scaled versions the minimizer works with
         hypo_maker._set_rescaled_free_params(scaled_param_vals) # pylint: disable=protected-access
+        if hypo_maker.__class__.__name__ == "Detectors":
+            update_param_values_detector(hypo_maker, hypo_maker.params.free) #updates values for ALL detectors
 
         # Get the map set
         try:
@@ -2548,7 +2598,6 @@ def _minimizer_callable(self, scaled_param_vals, hypo_maker, data_dist,
         #
         try:
             if hypo_maker.__class__.__name__ == "Detectors":
-                update_param_values_detector(hypo_maker, hypo_maker.params) #updates values for ALL detectors
                 metric_val = 0
                 for i in range(len(hypo_maker.distribution_makers)):
                     data = data_dist[i].metric_total(expected_values=hypo_asimov_dist[i],

pisa/core/detectors.py

@@ -91,6 +91,8 @@ def __init__(self, pipelines, label=None, set_livetime_from_data=True, profile=F
                     n += 1
             if n < 2:
                 raise NameError('Shared param %s only a free param in less than 2 detectors.' % sp)
+        
+        self.init_params()
 
     def __repr__(self):
         return self.tabulate(tablefmt="presto")
@@ -147,30 +149,33 @@ def get_outputs(self, **kwargs):
         return outputs
 
     def update_params(self, params):
-        for distribution_maker in self:
-            distribution_maker.update_params(params)
-
-        #if None in self.det_names: return # No detector names
-
-        if isinstance(params,Param): params = ParamSet(params) # just for the following
-
-        for p in params.names: # now update params with det_names inside
-            for i, det_name in enumerate(self.det_names):
-                if det_name in p:
-                    cp = deepcopy(params[p])
-                    cp.name = cp.name.replace('_'+det_name, "")
-                    self._distribution_makers[i].update_params(cp)
+        if isinstance(params,Param): params = ParamSet(params)
+            
+        for distribution_maker in self.distribution_makers:
+            ps = deepcopy(params)
+            for p in ps.names:
+                if distribution_maker.detector_name in p:
+                    p_name = p.replace('_'+distribution_maker.detector_name, "")
+                    if p_name in ps.names:
+                        ps.remove(p_name)
+                    ps[p].name = p_name
+            distribution_maker.update_params(ps)
+        self.init_params()
 
     def select_params(self, selections, error_on_missing=True):
         for distribution_maker in self:
             distribution_maker.select_params(selections=selections, error_on_missing=error_on_missing)
+        self.init_params()
 
     @property
     def distribution_makers(self):
         return self._distribution_makers
 
     @property
     def params(self):
+        return self._params
+    
+    def init_params(self):
         """Returns a ParamSet including all params of all detectors. First the shared params
         (if there are some), then all the "single detector" params. If two detectors use a
         parameter with the same name (but not shared), the name of the detector is added to the
@@ -196,7 +201,7 @@ def params(self):
                     params.extend(changed_param)
                 else:
                     params.extend(param)
-        return params
+        self._params = params
 
     @property
     def shared_param_ind_list(self):
@@ -221,7 +226,7 @@ def param_selections(self):
         selections = None
         for distribution_maker in self:
             if selections != None and sorted(distribution_maker.param_selections) != selections:
-                raise ('Different param_selections for different detectors.')
+                raise AssertionError('Different param_selections for different detectors.')
             selections = sorted(distribution_maker.param_selections)
         return selections
 

pisa_examples/resources/settings/osc/nufitv52.cfg

@@ -19,10 +19,12 @@ theta23_ih.3s_range = [39.9, 51.5] * units.deg
 theta23_ih.range = [30.0, 60.0] * units.deg
 
 deltacp_nh = 232 units.deg
+deltacp_nh_gauss = 232 +/- 36 units.deg
 deltacp_nh.3s_range = [144,350] * units.deg
 deltacp_nh.range = [0,360] * units.deg
 
 deltacp_ih = 276 units.deg
+deltacp_ih_gauss = 276 +/- 29 units.deg
 deltacp_ih.3s_range = [194,344] * units.deg
 deltacp_ih.range = [0,360] * units.deg
 
@@ -39,4 +41,4 @@ deltam31_nh.3s_range = [0.002427, 0.002590] * units.eV**2
 deltam31_ih = -0.002412 units.eV**2
 deltam31_ih_gauss = -0.002412 +/- 0.000028 * units.eV**2
 deltam31_ih.3s_range = [-0.002496, -0.002332] * units.eV**2
-deltam31_ih.range = [-0.007, -0.001] * units.eV**2
+deltam31_ih.range = [-0.007, -0.001] * units.eV**2