Merge pull request #25 from simonsobs/fisher

added fisher sampler

Author: mabitbol

bbpower/compsep.py

@@ -388,6 +388,19 @@ def chi2(par):
         res=minimize(chi2, self.params.p0, method="Powell")
         return res.x
 
+    def fisher(self):
+        """
+        Evaluate Fisher matrix
+        """
+        import numdifftools as nd
+        def lnprobd(p):
+            l = self.lnprob(p)
+            if l == -np.inf:
+                l = -1E100
+            return l
+        fisher = - nd.Hessian(lnprobd)(self.params.p0)
+        return fisher
+
     def singlepoint(self):
         """
         Evaluate at a single point
@@ -417,6 +430,15 @@ def run(self):
                      chain=sampler.chain,         
                      names=self.params.p_free_names)
             print("Finished sampling")
+        elif self.config.get('sampler')=='fisher':
+            fisher = self.fisher()
+            cov = np.linalg.inv(fisher)
+            for i,(n,p) in enumerate(zip(self.params.p_free_names,
+                                         self.params.p0)):
+                print(n+" = %.3lE +- %.3lE" % (p, np.sqrt(cov[i, i])))
+            np.savez(self.get_output('param_chains'),
+                     fisher=fisher,
+                     names=self.params.p_free_names)
         elif self.config.get('sampler')=='maximum_likelihood':
             sampler = self.minimizer()
             chi2 = -2*self.lnprob(sampler)

examples/generate_SO_maps.py

@@ -2,7 +2,7 @@
 import matplotlib.pyplot as plt
 from noise_calc import Simons_Observatory_V3_SA_noise,Simons_Observatory_V3_SA_beams
 import sys
-import pymaster as nmt
+#import pymaster as nmt
 import healpy as hp
 import os
 import sacc
@@ -212,6 +212,11 @@ def read_camb(fname):
         bpw_model[b1,1,b2,1,:]+=bpw_sync_bb*seds[b1,1]*seds[b2,1]
         bpw_model[b1,0,b2,0,:]+=bpw_dust_ee*seds[b1,2]*seds[b2,2]
         bpw_model[b1,1,b2,1,:]+=bpw_dust_bb*seds[b1,2]*seds[b2,2]
+np.savez("c_ells_sky",
+         ls = ells_bpw,
+         cls_ee = bpw_model[:,0,:,0,:],
+         cls_bb = bpw_model[:,1,:,1,:])
+exit(1)         
 tracers=[]
 for b in range(nfreqs):
     T=sacc.Tracer("band%d"%(b+1),'CMBP',

examples/generate_SO_spectra.py

@@ -27,7 +27,7 @@ def get_output_params(do_phase=False,do_angle=False):
 # Choose here whether to include the effects of
 #  - A frequency-dependent polarization angle (do_phase=True)
 #  - A non-zero constant polarization angle (do_angle=True)
-prefix_out,phase_nu,angles=get_output_params(do_phase=True,do_angle=True)
+prefix_out,phase_nu,angles=get_output_params(do_phase=False,do_angle=False)
 
 
 #CMB spectrum
@@ -123,7 +123,9 @@ def rotate(self,cl,transpose=False):
 temp_dust=19.6
 nu0_dust=353.
 
-Alens=1.
+prefix_out+="_2y_Al1p0"
+Alens=1.0
+nyears=2.
 
 #Bandpowers
 dell=10
@@ -205,6 +207,7 @@ def rotate_cells_mat(mat1, mat2, cls):
 nell=np.zeros([nfreqs,lmax+1])
 _,nell[:,2:],_=Simons_Observatory_V3_SA_noise(sens,knee,ylf,fsky,lmax+1,1)
 nell*=cl2dl[None,:]
+nell*=5./nyears
 n_bpw=np.sum(nell[:,None,:]*windows[None,:,:],axis=2)
 bpw_freq_noi=np.zeros_like(bpw_freq_sig)
 for ib,n in enumerate(n_bpw):