FisherIM-cross-fsig8.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Author: Alkistis Pourtsidou, ICG Portsmouth \n",
    "\n",
    "Using part of http://camb.readthedocs.io/en/latest/CAMBdemo.html\n",
    "\n",
    "To run this notebook you need to have CAMB and the CAMB python package installed. \n",
    "In order to install the CAMB python package on your computer follow the instructions in http://camb.readthedocs.io/en/latest/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import sys, platform, os\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import camb\n",
    "from camb import model, initialpower"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import scipy\n",
    "from scipy.interpolate import interp1d\n",
    "from __future__ import division"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy import integrate\n",
    "from scipy import linalg\n",
    "\n",
    "pi=np.pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "font = {'size'   : 16, 'family':'STIXGeneral'}\n",
    "axislabelfontsize='x-large'\n",
    "plt.rc('font', **font)\n",
    "plt.rcParams['legend.fontsize']='small'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Fiducial cosmological parameters\n",
    "c=3e5\n",
    "hubble=0.678\n",
    "omegab=0.022*pow(hubble,-2)\n",
    "omegac=0.119*pow(hubble,-2)\n",
    "om0=omegac+omegab\n",
    "H00=100*hubble\n",
    "Ass=2.14e-9\n",
    "nss = 0.968\n",
    "\n",
    "gamma=0.545\n",
    "\n",
    "print om0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Set up the fiducial cosmology\n",
    "pars = camb.CAMBparams()\n",
    "#Set cosmology\n",
    "pars.set_cosmology(H0=H00, ombh2=omegab*pow(hubble,2), omch2=omegac*pow(hubble,2),omk=0,mnu=0)\n",
    "pars.set_dark_energy() #LCDM (default)\n",
    "pars.InitPower.set_params(ns=nss, r=0, As=Ass)\n",
    "pars.set_for_lmax(2500, lens_potential_accuracy=0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#calculate results for these parameters\n",
    "results = camb.get_results(pars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Get matter power spectrum at z=0: P(k,z=0)\n",
    "\n",
    "#Not non-linear corrections couples to smaller scales than you want\n",
    "pars.set_matter_power(redshifts=[0.], kmax=5.0)\n",
    "\n",
    "#Non-Linear spectra (Halofit)\n",
    "pars.NonLinear = model.NonLinear_both\n",
    "results.calc_power_spectra(pars)\n",
    "kh, z, pk = results.get_matter_power_spectrum(minkh=1e-4, maxkh=5.0, npoints = 200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Construct P(k,z=0) interpolating function, in units of Mpc (no h)\n",
    "Pkz0 = interp1d(kh*hubble, pk[0]/pow(hubble,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Redshift bins\n",
    "zlist = np.arange(0.7,1.45,0.1)\n",
    "ztest = zlist[0]\n",
    "Nzbins = len(zlist)\n",
    "\n",
    "#mean number density of galaxies\n",
    "factor = pow(hubble,3)*1e-4\n",
    "nbarlist = [17.5, 19, 18, 17, 15, 13, 12, 10]\n",
    "\n",
    "print zlist\n",
    "print \"ztest =\", ztest\n",
    "print \"Number of redshift bins =\", Nzbins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Define E(z) = H(z)/H0\n",
    "def Ez(zc):\n",
    "    return np.sqrt(1-om0+om0*pow(1+zc,3))\n",
    "\n",
    "#Define the comoving distance\n",
    "def drdz(zp):\n",
    "    return (c/H00)/Ez(zp)\n",
    "def rcom(zc):\n",
    "    return scipy.integrate.romberg(drdz,0,zc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Define the growth function in LCDM\n",
    "def get_growth(zz):\n",
    "    omz=om0*pow(1+zz,3)/(om0*pow(1+zz,3)+1-om0)\n",
    "    return pow(omz,gamma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Get the growth factor \n",
    "def Dg_dz(zz):\n",
    "    return get_growth(zz)/(1+zz)\n",
    "def Dgz(zc):\n",
    "    ans = scipy.integrate.romberg(Dg_dz, 0.0, zc)\n",
    "    return np.exp(-ans)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#fiducial bHI from Bull et al 2015\n",
    "def bHI(zc):\n",
    "    return 0.67+0.18*zc+0.05*pow(zc,2)\n",
    "\n",
    "#fiducial bgal\n",
    "def bg(zc):\n",
    "    return np.sqrt(1+zc)\n",
    "\n",
    "#fiducial OmHI Mario's fit\n",
    "def OmHI(zc):\n",
    "    return 0.00048+0.00039*zc-0.000065*pow(zc,2)\n",
    "\n",
    "#mean brightness temperature [mK] Mario's fit\n",
    "def Tb(zc):\n",
    "    return 0.0559+0.2324*zc-0.024*pow(zc,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Construct PHIg(k,μ,z) \n",
    "def PHIg(kk,mu,zc):    \n",
    "    return Tb(zc)*(bHI(zc)+get_growth(zc)*mu**2)*(bg(zc)+get_growth(zc)*mu**2)*pow(Dgz(zc),2)*Pkz0(kk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Construct P_HI(k,μ,z) [mK^2]\n",
    "def PHI(kk,mu,zc):\n",
    "    return pow(Tb(zc),2)*pow(bHI(zc),2)*pow(1+(get_growth(zc)/bHI(zc))*mu**2,2)*pow(Dgz(zc),2)*Pkz0(kk)\n",
    "\n",
    "#Construct Pgg(k,μ,z) \n",
    "def Pgg(kk,mu,zc):    \n",
    "    return pow(bg(zc),2)*pow(1+(get_growth(zc)/bg(zc))*mu**2,2)*pow(Dgz(zc),2)*Pkz0(kk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#SKA1 noise specs\n",
    "Ndishes=190\n",
    "Ddish=15*100 #cm\n",
    "Nbeams=1\n",
    "\n",
    "def thetab(zc):\n",
    "    return 21*(1+zc)/Ddish\n",
    "\n",
    "def omegapix(zc):\n",
    "    return 1.13*pow(thetab(zc),2)\n",
    "\n",
    "Area=7000.0 #deg^2\n",
    "omegatot = Area*pow(pi/180,2)\n",
    "ttotal = 5000*60*60 #4000 hours\n",
    "\n",
    "def fc(zc):\n",
    "    return 1420.4/(1+zc)\n",
    "\n",
    "def Tsky(zc):\n",
    "    return 60*pow(300/fc(zc),2.55)*1e3\n",
    "\n",
    "#receiver temperature\n",
    "Trcv = 25.0*1e3 #mK\n",
    "\n",
    "def tobs(zc):\n",
    "    return ttotal*(omegapix(zc)/omegatot)*Ndishes*Nbeams\n",
    "\n",
    "Dzbin = 0.1\n",
    "dfpix = 50*1e3 #Hz\n",
    "midfreq = 1420.4e6 #Hz\n",
    "\n",
    "def dzpix(zc):\n",
    "    return pow(1+zc,2)*dfpix/midfreq\n",
    "def sigpix(zc,Tsys):\n",
    "    return Tsys/np.sqrt(dfpix*tobs(zc)) \n",
    "def dVpixdz(zz):    \n",
    "    return c*pow(rcom(zz),2)/(H00*Ez(zz))\n",
    "def Vpix(zc):\n",
    "    return omegapix(zc)*scipy.integrate.romberg(dVpixdz,zc-dzpix(zc)/2,zc+dzpix(zc)/2)\n",
    "\n",
    "def Wsq(kk,mu,zc):\n",
    "    #add very small offset to avoid division by zero\n",
    "    return 1e-20+np.exp(-pow(kk,2)*(1-mu**2)*pow(rcom(zc),2)*pow(thetab(zc),2)/(8*np.log(2)))\n",
    "\n",
    "def Pnoise(kk,mu,zc,Tsys):\n",
    "    return pow(sigpix(zc,Tsys),2)*Vpix(zc)*pow(Wsq(kk,mu,zc),-1.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#survey (bin) volume [Mpc^3]\n",
    "def dVsurdz(zz):    \n",
    "    return omegatot*c*pow(rcom(zz),2)/(H00*Ez(zz))\n",
    "    \n",
    "def Vsur(zc):\n",
    "    return scipy.integrate.romberg(dVsurdz,zc-Dzbin/2,zc+Dzbin/2)\n",
    "\n",
    "def Pshot(zc):\n",
    "    return 1/nbar\n",
    "\n",
    "def Veff(kk,mu,zc):\n",
    "    return Vsur(zc)*(pow(PHIg(kk,mu,zc),2)/(pow(PHIg(kk,mu,zc),2)+(PHI(kk,mu,zc)+Pnoise(kk,mu,zc,Tsys))\n",
    "                                         *(Pgg(kk,mu,zc)+Pshot(zc))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def kmin(zc):\n",
    "    return 2*pi*pow(Vsur(zc),-1/3)\n",
    "\n",
    "def kmax(zc):\n",
    "    return 0.14*pow(1+zc,2/(2+nss)) #non-linear cutoff (Smith et al 2003)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Npar = 4\n",
    "params = [\"0:fsig8\",\"1:bsig8\",\"2:DA\",\"3:H\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Fisher matrix derivatives\n",
    "\n",
    "def dlnP_dlnfsig8(kk,mu,zc):\n",
    "    return ((mu**2*get_growth(zc)/(bg(zc)+mu**2*get_growth(zc)))+\n",
    "           (mu**2*get_growth(zc)/(bHI(zc)+mu**2*get_growth(zc))))\n",
    "\n",
    "def dlnP_dlnbHIsig8(kk,mu,zc):\n",
    "    return bHI(zc)/(bHI(zc)+mu**2*get_growth(zc))\n",
    "\n",
    "def dlnP_dlnDA(kk,mu,zc):\n",
    "    dk = (kmax(zc)-kmin(zc))/400\n",
    "    return (-2.0+2*mu**2*(1-mu**2)*get_growth(zc)/(bHI(zc)+mu**2*get_growth(zc))\n",
    "            +2*mu**2*(1-mu**2)*get_growth(zc)/(bg(zc)+mu**2*get_growth(zc))\n",
    "            -kk*(1-mu**2)*(1/Pkz0(kk))*(Pkz0(kk+dk)-Pkz0(kk-dk))/(2*dk))\n",
    "\n",
    "def dlnP_dlnH(kk,mu,zc):\n",
    "    dk = (kmax(zc)-kmin(zc))/400\n",
    "    return (1.0+2*mu**2*(1-mu**2)*get_growth(zc)/(bHI(zc)+mu**2*get_growth(zc))\n",
    "            +2*mu**2*(1-mu**2)*get_growth(zc)/(bg(zc)+mu**2*get_growth(zc))\n",
    "            +kk*mu**2*(1/Pkz0(kk))*(Pkz0(kk+dk)-Pkz0(kk-dk))/(2*dk))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def dF(kk,mu):\n",
    "    return (1./(8*pi*pi))*pow(kk,2)*deriv_i(kk,mu,zc)*deriv_j(kk,mu,zc)*Veff(kk,mu,zc)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#2D integration function\n",
    "def integrate2D(dfun, kgrid, mugrid):\n",
    "    \n",
    "    muint = [scipy.integrate.simps(dfun.T[i], mugrid) for i in range(kgrid.size)]\n",
    "    return scipy.integrate.simps(muint, kgrid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mugrid = np.linspace(-1., 1., 200) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%%time\n",
    "#   Fisher matrix   # \n",
    "\n",
    "Npar = 4\n",
    "#create array of zeros\n",
    "s = (Npar,Npar)\n",
    "\n",
    "for zi in range(0,Nzbins):\n",
    "    zc = zlist[zi]\n",
    "    Tsys = Trcv+Tsky(zc)\n",
    "    nbar = factor*nbarlist[zi]\n",
    "    kgrid = np.linspace(kmin(zc), kmax(zc), 400)\n",
    "    K, MU = np.meshgrid(kgrid, mugrid)\n",
    "    Fishermat = np.zeros(s)\n",
    "    for i in range(0,Npar):  \n",
    "        def deriv_i(kk,mu,zc):\n",
    "            if i==0:  return dlnP_dlnfsig8(kk,mu,zc)\n",
    "            elif i==1:  return dlnP_dlnbHIsig8(kk,mu,zc)\n",
    "            elif i==2:  return dlnP_dlnDA(kk,mu,zc)\n",
    "            elif i==3:  return dlnP_dlnH(kk,mu,zc)\n",
    "            else: print \"index out of range\"\n",
    "        for  j in range(0,Npar):\n",
    "            if j>=i:\n",
    "                def deriv_j(kk,mu,zc):\n",
    "                    if j==0:  return dlnP_dlnfsig8(kk,mu,zc)\n",
    "                    elif j==1:  return dlnP_dlnbHIsig8(kk,mu,zc)\n",
    "                    elif j==2:  return dlnP_dlnDA(kk,mu,zc)\n",
    "                    elif j==3:  return dlnP_dlnH(kk,mu,zc)\n",
    "                    else: print \"index out of range\" \n",
    "                Fishermat[i][j] = integrate2D(dF(K,MU),kgrid,mugrid)                \n",
    "            else: Fishermat[i,j] = Fishermat[j,i]\n",
    "    print zc, np.sqrt(linalg.inv(Fishermat)[0,0]),\\\n",
    "              np.sqrt(linalg.inv(Fishermat)[2,2]),\\\n",
    "              np.sqrt(linalg.inv(Fishermat)[3,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}