cleanup

README.md

@@ -0,0 +1,53 @@
+# PyFisher
+
+PyFisher is a python package for calculating Fisher matrices and forecasting parameter uncertainties for CMB experiments. 
+
+## Installation
+
+You will need to install a dependency called `orphics` that is not available on pip, conda, etc.
+https://github.com/msyriac/orphics
+
+The installation procedure for that dependency is the same as for this package. Repeat the following for each of `orphics` and `pyfisher` separately.
+
+1. Git clone the repository
+2. `cd` into the repository and run `pip install -e . --user`.
+
+The latter step just copies symbolic links to the relevant modules into a directory (managed by pip) that is in your python path.
+
+Once this is done, you should be able to do things like
+
+``
+import pyfisher.clFisher as clFish
+``
+
+from anywhere on your system.
+
+
+## Basic Usage
+
+1. Edit `input/params.ini`. It is documented. You can include or exclude BAO by changing `otherFishers`, and the details of lensing reconstruction by changing the `lensing` section. There are a whole bunch of sections for different experiment configurations too.
+2. See how `tests/testFisher.py` uses this ini file. Inside the root repo directory, run
+``
+python tests/testFisher.py <experiment section name> <lensing section name>
+``
+For example
+``
+python tests/testFisher.py S4-5m lensing
+``
+which will tell you what the mnu constraint and lensing S/N is. You can adapt this script to, for example, override one of the ini configurations and make plots of constraints and S/N varying that configuration.
+
+## Advanced Usage
+
+
+### Make derivatives
+
+Edit `input/makeDefaults.ini` and specify a fiducial cosmology, derivative step sizes and a prefix under which to save the derivatives to  `output/`. Run with
+
+`python bin/makeDerivs.py`
+
+### Re-make Planck or BAO Fisher matrices.
+
+Edit `input/fisherDebug.ini` and run
+```
+python bin/driver.py input/fisherDebug.ini
+```

bin/BAO.py

@@ -0,0 +1,113 @@
+'''
+Output BAO Fisher matrix into .csv file in output/
+'''
+import camb
+from camb import model, initialpower
+import numpy as np
+import sys
+import ConfigParser
+import itertools
+
+from makeDerivs import getHubbleCosmology
+
+def getBAOCamb(zrange,params,AccuracyBoost=False):
+    '''
+    Get BAO's fk=rs/D_V from CAMB. Return 1D array of fk for all redshifts z
+    '''
+    pars = camb.CAMBparams()
+    if AccuracyBoost:
+        pars.set_accuracy(AccuracyBoost=2.0, lSampleBoost=2.0, lAccuracyBoost=2.0)
+        #pars.set_accuracy(AccuracyBoost=3.0, lSampleBoost=3.0, lAccuracyBoost=3.0)
+    pars.set_cosmology(H0=params['H0'], ombh2=params['ombh2'], omch2=params['omch2'], tau=params['tau'],mnu=params['mnu'],nnu=params['nnu'],omk=params['omk'])
+    pars.set_dark_energy(w=params['w'])
+    pars.InitPower.set_params(As=params['As'],ns=params['ns'],r=params['r'])
+    #pars.set_for_lmax(lmax=int(params['lmax']), lens_potential_accuracy=1, max_eta_k=2*params['lmax'])
+    camb.set_z_outputs(zrange)
+    results = camb.get_results(pars)
+    fk = results.get_background_outputs() #results.get_BAO(zrange,pars)[:,0]
+    #print fk
+    fk = np.nan_to_num(fk[:,0])
+    return fk
+
+def main(argv):
+    verbose=True
+
+    # Read Config
+    iniFile = "input/makeDefaultsBAO.ini"
+    Config = ConfigParser.SafeConfigParser()
+    Config.optionxform = str
+    Config.read(iniFile)
+    out_pre = Config.get('general','output_prefix')
+    AccuracyBoost = Config.getboolean('general','AccuracyBoost')
+    paramList = []
+    fparams = {}
+    stepSizes = {}
+    for (key, val) in Config.items('camb'):
+        if ',' in val:
+            param, step = val.split(',')
+            paramList.append(key)
+            fparams[key] = float(param)
+            stepSizes[key] = float(step)
+        else:
+            fparams[key] = float(val)    
+    zrange = [float(x) for x in Config.get('DESI','redshift').split(',')]
+    sigmafk = [float(x)/1000. for x in Config.get('DESI','sigmafkx1000').split(',')]
+
+    # Uncomment if need to calculate Derivs
+
+    # Save fiducials
+    print "Calculating and saving fiducial cosmology..."
+    if not('H0' in fparams):
+        fparams['H0'] = getHubbleCosmology(theta=fparams['theta100'],params=fparams)
+    fidfk = getBAOCamb(zrange,fparams,AccuracyBoost=AccuracyBoost)
+    np.savetxt("output/"+out_pre+"_fidfk.csv",fidfk,delimiter=",")
+
+
+    # Calculate and save derivatives
+    for paramName in paramList:
+        h = stepSizes[paramName]        
+        print "Calculating forward difference for ", paramName
+        pparams = fparams.copy()
+        pparams[paramName] = fparams[paramName] + 0.5*h
+        if paramName=='theta100':
+            pparams['H0'] = getHubbleCosmology(theta=pparams['theta100'],params=pparams)
+        pfk = getBAOCamb(zrange,pparams,AccuracyBoost=AccuracyBoost)
+
+
+        print "Calculating backward difference for ", paramName
+        mparams = fparams.copy()
+        mparams[paramName] = fparams[paramName] - 0.5*h
+        if paramName=='theta100':
+            mparams['H0'] = getHubbleCosmology(theta=mparams['theta100'],params=mparams)
+        mfk = getBAOCamb(zrange,mparams,AccuracyBoost=AccuracyBoost)
+
+        dfk = (pfk-mfk)/h
+
+        np.savetxt("output/"+out_pre+"_dfk_"+paramName+".csv",dfk,delimiter=",")
+
+    # Calculate Fisher Matrix
+    paramCombs = itertools.combinations_with_replacement(paramList,2)
+    Fisher = np.zeros((len(paramList),len(paramList)))
+    dfks = {}
+    for paramName in paramList:
+        dfks[paramName] = np.loadtxt("output/"+out_pre+"_dfk_"+paramName+".csv",delimiter=",")
+
+    for param1,param2 in paramCombs:
+        if verbose: print "Parameter combination : ", param1,param2
+        i = paramList.index(param1)
+        j = paramList.index(param2)
+        Fz = 0.		
+        for k in range(0,len(zrange)):
+            dfk1 = dfks[param1][k]
+            dfk2 = dfks[param2][k]
+            Fz += dfk1*dfk2/sigmafk[k]**2.
+            #if verbose: print "dfk1,dfk2,sigmafk,Fz:",dfk1*dfk2/sigmafk[k]**2,Fz
+        Fisher[i,j] = Fz
+        Fisher[j,i] = Fz
+    if verbose:
+        print Fisher
+        print np.diagonal(Fisher)
+    np.savetxt("output/"+out_pre+"_Fisher.csv",Fisher,delimiter=",")
+
+if (__name__ == "__main__"):
+    main(sys.argv[1:])

bin/Cls_plot.py

@@ -0,0 +1,196 @@
+import numpy as np
+import sys,os
+import matplotlib.pyplot as plt
+from makeDerivs import getPowerCamb
+import ConfigParser
+import itertools
+
+colors = itertools.cycle(['b', 'g', 'r', 'c', 'm', 'y'])
+
+'''
+iniFile = os.environ['FISHER_DIR']+'/input/cambRed_test.ini'
+Config = ConfigParser.SafeConfigParser()
+Config.optionxform = str
+Config.read(iniFile)
+spec = Config.get('general','spec')
+out_pre = Config.get('general','output_prefix')
+AccuracyBoost = Config.getboolean('general','AccuracyBoost')
+fparams = {}
+for (key, val) in Config.items('camb'):
+    fparams[key] = float(val)
+Cls = getPowerCamb(fparams,spec,AccuracyBoost=AccuracyBoost)
+np.savetxt(os.environ['FISHER_DIR']+'/output/test_lensing_pycambCls.dat',Cls)
+
+log = True
+prefix = ''
+name = 'test1_lensing_scalCovCls'
+data = np.loadtxt('output/'+name+'.dat')
+ell=1
+
+n = float(len(data[0,:])-ell)
+a = int(round(np.sqrt(n),0))
+b = int(np.ceil(n/a))
+fig, ax = plt.subplots(a,b,figsize=(20,10))
+x = data[:,0]
+#x = np.arange(float(len(data)))
+i = ell
+
+for xpos in range(a):
+    for ypos in range(b):
+        ax[xpos,ypos].plot(x,data[:,i])
+        #ax[xpos,ypos].plot(x,data[:,i]*2*np.pi/(x*(x+1.0)))
+        
+        if i <= 3:
+            ax[xpos,ypos].plot(x,data[:,i]/2/np.pi*(x*(x+1.0)))
+        else:
+            #ax[xpos,ypos].plot(x,data[:,i]*(x*(x+1.0))**2)
+            ax[xpos,ypos].plot(x,data[:,i]*4/2/np.pi)
+        
+        if log:
+            ax[xpos,ypos].set_xscale('log')
+            #ax[xpos,ypos].set_yscale('log')
+            prefix = '(log)'
+        i+=1
+        if i>n:
+            break
+'''
+'''
+name1 = 'test_lensing_pycambCls'
+name2 = 'z1.0_scalCovCls'
+name22 = 'z1.0_lenspotentialCls'
+name3 = 'z1.0_lensedCls'
+
+data1 = np.loadtxt('output/'+name1+'.dat')
+data2 = np.loadtxt('output/'+name2+'.dat')
+data22 = np.loadtxt('output/'+name22+'.dat')
+data3 = np.loadtxt('output/'+name3+'.dat')
+
+fig = plt.figure()
+ax = fig.add_subplot(111)
+
+x1 = np.arange(float(len(data1)))
+y1 = data1[:,3]
+#y1 = data1[:,4]/x1**4
+#y1 = data1[:,5]*np.sqrt(7.4311e12)*2/np.sqrt(x1*(x1+1.))
+
+
+x2 = data2[:,0]
+y2 = data2[:,2]*2.*np.pi/(x2*(x2+1.))
+#y2 = data2[:,11]*2.*np.pi/4.*(x2*(x2+1.))
+#y2 = data2[:,3]*2.*np.pi/4.*np.sqrt(x2*(x2+1.))
+#y2 = data2[:,3]*2.*np.pi/2.
+
+x22 = data22[:,0]
+#y22 = data22[:,5]*2*np.pi/(x2*(x2+1))**2
+#y22 = data22[:,5]*2.*np.pi/4.
+y22 = data22[:,6]*2.*np.pi/4.
+y22 = data22[:,6]*2.*np.pi/2./np.sqrt(x22*(x22+1.))
+
+x3 = data3[:,0]
+#y3 = data3[:,4]/x3**4/7.4311e12/4
+#y3 = data3[:,4]/x3**4/7.4311e12/4*(x3*(x3+1))**2
+y3 = data3[:,4]
+
+
+for label in ['cc','cg','cs','gg','gs','ss']:
+    data = np.loadtxt('output/Das_'+label+'.csv',delimiter=',')
+    ax.plot(data[:,0],data[:,1],colors.next()+'--',label=label)
+
+names = ['July12_highAcc_Das_bias1_galaxy_fCls']#,'July12_highAcc_Das_bias2_galaxy_fCls']
+names = ['July12_highAcc_Das_galaxy_fCls','July12_highAcc_Das_bias1_galaxy_fCls']
+
+lss = itertools.cycle(['-','-.'])
+for name in names:
+    ls = lss.next()
+    data = np.loadtxt('output/'+name+'.csv',delimiter=',')
+    x = np.arange(float(len(data)))
+    #y = data[:,2]
+    
+    ax.plot(x,data[:,0],colors.next()+ls,label='$\kappa_{CMB} \kappa_{CMB}$')
+    ax.plot(x,data[:,1],colors.next()+ls,label='$\kappa_{CMB} \kappa_{gal}$')
+    ax.plot(x,data[:,2],colors.next()+ls,label='$\kappa_{CMB} \Sigma$')
+    ax.plot(x,data[:,3],colors.next()+ls,label='$\kappa_{gal} \kappa_{gal}$')
+    ax.plot(x,data[:,4],colors.next()+ls,label='$\kappa_{gal} \Sigma$')
+    ax.plot(x,data[:,5],colors.next()+ls,label='$\Sigma \Sigma$')
+    #ax.plot(x1,y1,label='pycamb')
+    #ax.plot(x2,y2,label='scalCov')
+    #ax.plot(x22,y22,label='lenspotentialCov')
+    #ax.plot(x3,y3,label='lensed')
+    #,x2,y2,x3,y3)
+
+ax.set_xscale('log')
+ax.set_yscale('log')
+
+#plt.legend(loc='lower left',ncol=2)
+plt.show()
+fig.suptitle(name+prefix,fontsize=20)
+fileName = 'output/'+name+prefix+'.png'
+#plt.savefig(fileName,format='png')
+#print "Saved figure ",fileName
+'''
+'''
+Cls = np.loadtxt('output/Nov3_highAcc_CDM_unlensed_axion_fCls.csv',delimiter=',')
+l = np.array(range(len(Cls)))
+imax = 3
+
+for i in range(imax):
+    if i == 5:
+        Cls[:,i] = Cls[:,i] * np.sqrt(l*(l+1.))
+    elif i == 4:
+        continue
+    else:
+        Cls[:,i] = Cls[:,i] * (l*(l+1.))
+print Cls
+
+lmin = 0
+fig = plt.figure()
+for i in range(imax):
+    ax = fig.add_subplot(2,3,i+1)
+    ax.plot(l[lmin:],Cls[lmin:,i])
+    ax.set_xscale('log')
+    if i in [4,5]:
+        ax.set_yscale('log')
+'''
+#Cls = np.loadtxt(os.environ['CMBPROJ_DIR']+'/data/TheorySpectra/ell28k_highacc_lensedCls.dat')
+Cls = np.loadtxt(os.environ['CMBPROJ_DIR']+'/data/TheorySpectra/Nov10_highAcc_CDM_lensedCls.dat')
+#Cls = np.loadtxt(os.environ['AXIONCAMB_DIR']+'Nov8_highAcc_CDM_scalCls.dat')
+print Cls
+lmin = 6000
+fig = plt.figure()
+j = 1
+for i in range(1,5):
+    ax = fig.add_subplot(2,4,j)
+    ax.plot(Cls[lmin:,0],Cls[lmin:,i])
+    #ax.set_xscale('log')
+    #ax.set_yscale('log')
+    j+=1
+    #vmax = max(Cls[lmin:,i])
+    #ax.set_ylim([-vmax/1.e3,vmax])
+#Cls = np.loadtxt('/home/nhnguyen/CAMB-May2016/Nov18_highAcc_CDM_lensedCls.dat')    
+Cls = np.loadtxt(os.environ['CMBPROJ_DIR']+'/data/TheorySpectra/ell28k_highacc_lensedCls.dat')
+#Cls = np.loadtxt(os.environ['CMBPROJ_DIR']+'/data/TheorySpectra/Nov9_highAcc_CDM_lensedCls.dat')
+for i in range(1,5):
+    ax = fig.add_subplot(2,4,j)
+    ax.plot(Cls[lmin:,0],Cls[lmin:,i])
+    #ax.set_xscale('log')
+    #ax.set_yscale('log')
+    j+=1
+print Cls
+'''
+i+=1
+Cls = np.loadtxt('output/Nov3_highAcc_CDM_unlensed_axion_fCls.csv',delimiter=',')
+x = np.array(range(len(Cls)))
+N = len(Cls)
+print x
+data = x.reshape([N,1])
+for j in range(4):
+    ax = fig.add_subplot(2,4,i)
+    y = Cls[:,j]*(x*(x+1.))/2./np.pi
+    ax.plot(x[lmin:],y[lmin:])
+    ax.set_xscale('log')
+    i+=1
+    data = np.hstack([data,y.reshape([N,1])])
+    #print y
+print data
+'''
+plt.show()