add support for triangular matrices to MultivariateNormalDistribution…

… and simplify loading of FFs (#247)

* [MultivariateNormalDistribution] add support for triangular matrices

* simplify loading of form factor parameters

* add tests with triangular input matrices

flavio/classes.py

@@ -148,6 +148,10 @@ def add_constraint(self, parameters, constraint):
 
         Note that if there already exists a constraint, it will be removed."""
         for num, parameter in enumerate(parameters):
+            try: # check if parameter object already exists
+                p = Parameter[parameter]
+            except: # otherwise, create a new one
+                p = Parameter(parameter)
             # remove constraint if there is one
             if parameter in self._parameters:
                 self.remove_constraint(parameter)

flavio/physics/bdecays/formfactors/b_p/bcl_parameters.py

@@ -7,17 +7,6 @@
 def load_parameters(filename, constraints):
     f = pkgutil.get_data('flavio.physics', filename)
     ff_dict = yaml.safe_load(f)
-    for parameter_name in ff_dict['parameters']:
-        try: # check if parameter object already exists
-            p = Parameter[parameter_name]
-        except: # otherwise, create a new one
-            p = Parameter(parameter_name)
-        else: # if parameter exists, remove existing constraints
-            constraints.remove_constraint(parameter_name)
     covariance = np.outer(ff_dict['uncertainties'], ff_dict['uncertainties'])*ff_dict['correlation']
-    if not np.allclose(covariance, covariance.T):
-        # if the covariance is not symmetric, it is assumed that only the values above the diagonal are present.
-        # then: M -> M + M^T - diag(M)
-        covariance = covariance + covariance.T - np.diag(np.diag(covariance))
     constraints.add_constraint(ff_dict['parameters'],
             MultivariateNormalDistribution(central_value=ff_dict['central_values'], covariance=covariance) )

flavio/physics/bdecays/formfactors/b_p/bcl_parameters_lmvd.py

@@ -22,17 +22,5 @@ def load_parameters(filename, constraints):
             except:
                 raise ValueError('No f or b in observable name')
         observables_renamed.append(o[:index]+ ' BCL ' + o[index:])
-
-    for parameter_name in observables_renamed:
-        try: # check if parameter object already exists
-            p = Parameter[parameter_name]
-        except: # otherwise, create a new one
-            p = Parameter(parameter_name)
-        else: # if parameter exists, remove existing constraints
-            constraints.remove_constraint(parameter_name)
-    if not np.allclose(covariance, covariance.T):
-        # if the covariance is not symmetric, it is assumed that only the values above the diagonal are present.
-        # then: M -> M + M^T - diag(M)
-        covariance = covariance + covariance.T - np.diag(np.diag(covariance))
     constraints.add_constraint(observables_renamed,
             MultivariateNormalDistribution(central_value=central_values, covariance=covariance) )

flavio/physics/bdecays/formfactors/b_p/bsz_parameters.py

@@ -46,8 +46,6 @@ def load_parameters(filename, process, constraints):
             _tex_ff = tex_ff[parameter_name.split(' ')[-1].split('_')[-1]]
             p.tex = r'$' + _tex_a + r'^{' + _tex_ff + r'}$'
             p.description = r'BSZ form factor parametrization coefficient $' + _tex_a + r'$ of $' + _tex_ff + r'$'
-        else:  # if parameter exists, remove existing constraints
-            constraints.remove_constraint(parameter_name)
     [central, unc, corr] = get_ffpar(filename)
     constraints.add_constraint(parameter_names,
         MultivariateNormalDistribution(central_value=central, covariance=np.outer(unc, unc)*corr))

flavio/physics/bdecays/formfactors/b_v/bsz_parameters.py

@@ -48,8 +48,6 @@ def load_parameters(filename, process, constraints):
             _tex_ff = tex_ff[parameter_name.split(' ')[-1].split('_')[-1]]
             p.tex = r'$' + _tex_a + r'^{' + _tex_ff + r'}$'
             p.description = r'BSZ form factor parametrization coefficient $' + _tex_a + r'$ of $' + _tex_ff + r'$'
-        else:  # if parameter exists, remove existing constraints
-            constraints.remove_constraint(parameter_name)
     [central, unc, corr] = get_ffpar(filename)
     constraints.add_constraint(parameter_names,
         MultivariateNormalDistribution(central_value=central, covariance=np.outer(unc, unc)*corr))

flavio/physics/bdecays/formfactors/b_v/lattice_parameters.py

@@ -33,13 +33,6 @@ def load_parameters(file_res, file_cov, process, constraints):
           + np.array([[ cov_dict.get((m,k),0) for m in keys_sorted] for k in keys_sorted])
           - np.diag([ cov_dict[(k,k)] for k in keys_sorted]) )
     parameter_names = [implementation_name + ' ' + coeff_name for coeff_name in keys_sorted]
-    for parameter_name in parameter_names:
-        try: # check if parameter object already exists
-            p = Parameter[parameter_name]
-        except: # otherwise, create a new one
-            p = Parameter(parameter_name)
-        else: # if parameter exists, remove existing constraints
-            constraints.remove_constraint(parameter_name)
     constraints.add_constraint(parameter_names,
             MultivariateNormalDistribution(central_value=res, covariance=cov ))
 

flavio/physics/bdecays/formfactors/lambdab_12/lattice_parameters.py

@@ -62,8 +62,6 @@ def load_parameters(file_res, file_cov, process, constraints):
             _tex_ff = tex_ff[parameter_name.split(' ')[-1].split('_')[-1]]
             p.tex = r'$' + _tex_a + r'^{' + _tex_ff + r'}$'
             p.description = r'SSE form factor parametrization coefficient $' + _tex_a + r'$ of $' + _tex_ff + r'$'
-        else: # if parameter exists, remove existing constraints
-            constraints.remove_constraint(parameter_name)
     constraints.add_constraint(parameter_names,
             MultivariateNormalDistribution(central_value=res, covariance=cov ))
 

flavio/statistics/probability.py

@@ -1325,6 +1325,10 @@ class MultivariateNormalDistribution(ProbabilityDistribution):
     If the covariance matrix is not specified, standard_deviation and the
     correlation matrix have to be specified.
 
+    If the covariance or correlation matrix is not symmetric, it is assumed that
+    only the values above the diagonal are present and the missing values are
+    filled in by reflecting the upper triangle at the diagonal.
+
     Methods:
 
     - get_random(size=None): get `size` random numbers (default: a single one)
@@ -1347,8 +1351,15 @@ def __init__(self, central_value, covariance=None,
 
         - central_value: vector of means, shape (n)
         - covariance: covariance matrix, shape (n,n)
+        - standard_deviation: vector of standard deviations, shape (n)
+        - correlation: correlation matrix, shape (n,n)
         """
         if covariance is not None:
+            covariance = np.array(covariance)
+            if not np.allclose(covariance, covariance.T):
+                # if the covariance is not symmetric, it is assumed that only the values above the diagonal are present.
+                # then: M -> M + M^T - diag(M)
+                covariance = covariance + covariance.T - np.diag(np.diag(covariance))
             self.covariance = covariance
             self.standard_deviation = np.sqrt(np.diag(self.covariance))
             self.correlation = self.covariance/np.outer(self.standard_deviation,
@@ -1366,7 +1377,12 @@ def __init__(self, central_value, covariance=None,
                     n_dim = len(central_value)
                     self.correlation = np.eye(n_dim) + (np.ones((n_dim, n_dim))-np.eye(n_dim))*float(correlation)
                 else:
-                    self.correlation = np.array(correlation)
+                    correlation = np.array(correlation)
+                    if not np.allclose(correlation, correlation.T):
+                        # if the correlation is not symmetric, it is assumed that only the values above the diagonal are present.
+                        # then: M -> M + M^T - diag(M)
+                        correlation = correlation + correlation.T - np.diag(np.diag(correlation))
+                    self.correlation = correlation
             self.covariance = np.outer(self.standard_deviation,
                                        self.standard_deviation)*self.correlation
         super().__init__(central_value, support=np.array([

flavio/statistics/test_probability.py

@@ -20,6 +20,28 @@ def test_multiv_normal(self):
         self.assertAlmostEqual(num_lpdf, ana_lpdf, delta=1e-6)
         self.assertEqual(len(pdf.get_random(10)), 10)
 
+    def test_multiv_normal_triangular_cov(self):
+        c = np.array([1e-3, 2])
+        cov = np.array([[(0.2e-3)**2, 0.2e-3*0.5*0.3],[0.2e-3*0.5*0.3, 0.5**2]])
+        cov_triangular = cov.copy()
+        cov_triangular[1,0] = 0
+        pdf_from_triangular_cov = MultivariateNormalDistribution(c, cov_triangular)
+        pdf = MultivariateNormalDistribution(c, cov)
+        np.testing.assert_equal(pdf.covariance, pdf_from_triangular_cov.covariance)
+        np.testing.assert_equal(pdf.correlation, pdf_from_triangular_cov.correlation)
+
+    def test_multiv_normal_triangular_corr(self):
+        c = np.array([1e-3, 2])
+        cov = np.array([[(0.2e-3)**2, 0.2e-3*0.5*0.3],[0.2e-3*0.5*0.3, 0.5**2]])
+        std = np.sqrt(np.diag(cov))
+        corr = cov / np.outer(std, std)
+        corr_triangular = corr.copy()
+        corr_triangular[1,0] = 0
+        pdf_from_triangular_corr = MultivariateNormalDistribution(c, standard_deviation=std, correlation=corr_triangular)
+        pdf = MultivariateNormalDistribution(c, standard_deviation=std, correlation=corr)
+        np.testing.assert_equal(pdf.covariance, pdf_from_triangular_corr.covariance)
+        np.testing.assert_equal(pdf.correlation, pdf_from_triangular_corr.correlation)
+
     def test_normal(self):
         d = NormalDistribution(2, 0.3)
         self.assertEqual(d.cdf(2), 0.5)
@@ -491,7 +513,7 @@ def test_repr(self):
         self.assertEqual(repr(GeneralGammaUpperLimit(limit=1e-9, confidence_level=0.95, counts_total=15, counts_background=10, background_std=0.2)),
                          fsp + 'GeneralGammaUpperLimit(limit=1e-09, confidence_level=0.95, counts_total=15, counts_background=10, background_std=0.2)')
         self.assertEqual(repr(MultivariateNormalDistribution([1., 2], [[2, 0.1], [0.1, 2]])),
-                         fsp + 'MultivariateNormalDistribution([1.0, 2], [[2, 0.1], [0.1, 2]])')
+                         fsp + 'MultivariateNormalDistribution([1.0, 2], [[2.  0.1]\n [0.1 2. ]])')
         self.assertEqual(repr(NumericalDistribution([1., 2], [3, 4.])),
                          fsp + 'NumericalDistribution([1.0, 2], [3, 4.0])')
         self.assertEqual(repr(GaussianKDE([1, 2, 3], 0.1)),