Merge pull request #130 from andreicuceu/bb-fix

bug fix in indices of broad band parameters

tests/full_configs/main.ini

@@ -67,18 +67,18 @@ par_sigma_smooth = 2.
 per_sigma_smooth = 2.
 BB-lyalya_lyalya-0-broadband_sky-scale-sky = 0.00896
 BB-lyalya_lyalya-0-broadband_sky-sigma-sky = 32.7
-BB-lyalya_lyalya-0 add post r,mu (0,0) = 0.1
-BB-lyalya_lyalya-0 add post r,mu (0,2) = 0.1
-BB-lyalya_lyalya-0 add post r,mu (0,4) = 0.1
-BB-lyalya_lyalya-0 add post r,mu (0,6) = 0.1
-BB-lyalya_lyalya-0 add post r,mu (1,0) = 0.1
-BB-lyalya_lyalya-0 add post r,mu (1,2) = 0.1
-BB-lyalya_lyalya-0 add post r,mu (1,4) = 0.1
-BB-lyalya_lyalya-0 add post r,mu (1,6) = 0.1
-BB-lyalya_lyalya-0 add post r,mu (2,0) = 0.1
-BB-lyalya_lyalya-0 add post r,mu (2,2) = 0.1
-BB-lyalya_lyalya-0 add post r,mu (2,4) = 0.1
-BB-lyalya_lyalya-0 add post r,mu (2,6) = 0.1
+BB-lyalya_lyalya-0 add post r,mu (0,0) = 0.001
+BB-lyalya_lyalya-0 add post r,mu (0,2) = 0.001
+BB-lyalya_lyalya-0 add post r,mu (0,4) = 0.001
+BB-lyalya_lyalya-0 add post r,mu (0,6) = 0.001
+BB-lyalya_lyalya-0 add post r,mu (1,0) = 0.001
+BB-lyalya_lyalya-0 add post r,mu (1,2) = 0.001
+BB-lyalya_lyalya-0 add post r,mu (1,4) = 0.001
+BB-lyalya_lyalya-0 add post r,mu (1,6) = 0.001
+BB-lyalya_lyalya-0 add post r,mu (2,0) = 0.001
+BB-lyalya_lyalya-0 add post r,mu (2,2) = 0.001
+BB-lyalya_lyalya-0 add post r,mu (2,4) = 0.001
+BB-lyalya_lyalya-0 add post r,mu (2,6) = 0.001
 
 [control]
 sampler = True

tests/test_vega.py

@@ -11,11 +11,11 @@ def test_vega_new():
     vega = VegaInterface('full_configs/main.ini')
 
     loglik = vega.log_lik()
-    assert isclose(loglik, -58967.38295296009)
+    assert isclose(loglik, -8766.997113218942)
 
     vega.minimize()
 
-    assert isclose(vega.bestfit.fmin.fval, 100380.65140993778)
+    assert isclose(vega.bestfit.fmin.fval, 0.6392334022603812)
 
 
 def test_vega_old():

vega/correlation_func.py

@@ -531,10 +531,21 @@ def broadband(self, bb_term, params):
             for j in r2_powers:
                 bb_params.append(params['{} ({},{})'.format(
                         bb_term['name'], i, j)])
-        bb_params = np.array(bb_params).reshape(-1, r_max - r_min + 1)
 
-        corr = (bb_params[:, :, None, None] * r1**r1_powers[:, None, None] *
-                r2**r2_powers[None, :, None]).sum(axis=(0, 1, 2))
+        # the first dimension of bb_params is that of r1 power indices
+        # the second dimension of bb_params is that of r2 power indices
+        bb_params = np.array(bb_params).reshape(r_max - r_min + 1, -1)
+
+        # we are summing 3D array along 2 dimensions.
+        # first dimension is the data array (2500 for the standard rp rt grid of 50x50)
+        # second dimension is the first dimension of bb_params  = r1 power indices
+        # third dimension is the sectond dimension of bb_params = r2 power indices
+        # dimensions addressed in an array get the indices ':'
+        # dimensions not addressed in an array get the indices 'None'
+        # we sum over the second and third dimensions which are powers of r
+        corr = (bb_params[None,:, :] * r1[:,None,None]**r1_powers[None, :, None] *
+                r2[:,None,None]**r2_powers[None, None, :]).sum(axis=(1, 2))
+
         return corr
 
     def compute_qso_radiation(self, params):