Euclid-like.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": 1,
   "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": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import scipy\n",
    "from scipy.interpolate import interp1d\n",
    "from __future__ import division"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "font = {'size'   : 16, 'family':'STIXGeneral'}\n",
    "axislabelfontsize='x-large'\n",
    "plt.rc('font', **font)\n",
    "plt.rcParams['legend.fontsize']='medium'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy import integrate\n",
    "from scipy import linalg\n",
    "\n",
    "pi=np.pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Fiducial cosmological parameters\n",
    "c=3e5\n",
    "hubble=0.675\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",
    "\n",
    "nss = 0.9655\n",
    "alphass = 0.0\n",
    "betass = 0.0\n",
    "\n",
    "gamma=0.545\n",
    "\n",
    "kpivot = 0.05 \n",
    "#kpivot = 0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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, nrun=alphass, nrunrun=betass, pivot_scalar=kpivot)\n",
    "pars.set_for_lmax(2500, lens_potential_accuracy=0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#calculate results for these parameters\n",
    "results = camb.get_results(pars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Get matter power spectrum at z=0: P(k,z=0)\n",
    "#Not non-linear corrections couples to smaller scales than you want\n",
    "pars.set_matter_power(redshifts=[0.], kmax=1.0)\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-5, maxkh=0.5, npoints = 800)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Derivatives\n",
    "\n",
    "def dlogPd_dns(kk):\n",
    "    return np.log(kk/kpivot)\n",
    "\n",
    "def dlogPd_dalphas(kk):\n",
    "    return (1./2)*pow(np.log(kk/kpivot),2)\n",
    "\n",
    "def dlogPd_dbetas(kk):\n",
    "    return (1./6)*pow(np.log(kk/kpivot),3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.7  0.8  0.9  1.   1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9  2. ]\n",
      "ztest = 1.2\n",
      "Number of redshift bins = 14\n"
     ]
    }
   ],
   "source": [
    "#Redshift bins\n",
    "zlist = np.arange(0.7,2.05,0.1)\n",
    "ztest = zlist[5]\n",
    "Nzbins = len(zlist)\n",
    "\n",
    "Area = 15000.0 #deg^2\n",
    "omegatot = Area*pow(pi/180,2)\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, 8, 6, 4, 3.5, 2, 1]\n",
    "\n",
    "biaslist = [1.083, 1.125, 1.104, 1.126, 1.208, 1.243, 1.282, 1.292, 1.363, 1.497, 1.486, \\\n",
    "            1.491, 1.573, 1.568]\n",
    "\n",
    "Dzbin = 0.1\n",
    "\n",
    "\n",
    "print zlist\n",
    "print \"ztest =\", ztest\n",
    "print \"Number of redshift bins =\", Nzbins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3881.21332089\n"
     ]
    }
   ],
   "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)\n",
    "\n",
    "print rcom(ztest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.901505284449\n"
     ]
    }
   ],
   "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)\n",
    "\n",
    "print get_growth(ztest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.558232251833\n"
     ]
    }
   ],
   "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)\n",
    "\n",
    "print Dgz(ztest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Construct Pgg(k,μ,z) \n",
    "def Pgg(kk,mu,zc):    \n",
    "    return pow(bg+get_growth(zc)*mu**2,2)*pow(Dgz(zc),2)*Pkz0(kk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def kmin(zc):\n",
    "    return 2*pi*pow(Vsur(zc),-1/3)\n",
    "\n",
    "kNL = 0.1 #Mpc^{-1}    \n",
    "    \n",
    "def kmax(zc):\n",
    "    return kNL*pow(1+zc,2/(2+nss)) #non-linear cutoff "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7 0.00300391239432 0.143027630415\n",
      "0.8 0.00285567190509 0.14864884108\n",
      "0.9 0.00274250846588 0.154169231521\n",
      "1.0 0.0026545118315 0.159595782335\n",
      "1.1 0.00258515759285 0.164934667994\n",
      "1.2 0.00252997099126 0.170191383714\n",
      "1.3 0.00248578002456 0.17537084751\n",
      "1.4 0.00245027417685 0.180477483171\n",
      "1.5 0.00242173160884 0.185515288356\n",
      "1.6 0.0023988439947 0.190487890971\n",
      "1.7 0.00238060039977 0.195398596213\n",
      "1.8 0.00236620816162 0.200250426112\n",
      "1.9 0.00235503768985 0.205046153\n",
      "2.0 0.00234658314854 0.209788328003\n"
     ]
    }
   ],
   "source": [
    "for zi in range(0,Nzbins):\n",
    "    zc = zlist[zi]\n",
    "    print zc, kmin(zc), kmax(zc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def Pshot(nbar):\n",
    "    return 1/nbar\n",
    "\n",
    "def Veff(kk,mu,zc):\n",
    "    return Vsur(zc)*(Pgg(kk,mu,zc)/(Pgg(kk,mu,zc)+Pshot(nbar)))**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "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": 21,
   "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": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mugrid = np.linspace(-1., 1., 513) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Fisher matrix parameters \n",
    "params = [\"0:ns\",\"1:alphas\",\"2:betas\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "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": 25,
   "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": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%%time\n",
    "#   Fisher matrix   # \n",
    "\n",
    "Npar = 2 #Npar = 3\n",
    "#create array of zeros\n",
    "s = (Npar,Npar)\n",
    "\n",
    "Fishermat = np.zeros(s)\n",
    "\n",
    "for zi in range(0,Nzbins):\n",
    "    zc = zlist[zi]\n",
    "    nbar = factor*nbarlist[zi]\n",
    "    bg = biaslist[zi]\n",
    "    kgrid = np.linspace(kmin(zc), kmax(zc), 400)\n",
    "    K, MU = np.meshgrid(kgrid, mugrid)\n",
    "    for i in range(0,Npar):  \n",
    "        def deriv_i(kk,mu,zc):\n",
    "            if i==0:  return dlogPd_dns(kk)\n",
    "            elif i==1:  return dlogPd_dalphas(kk)\n",
    "            elif i==2:  return dlogPd_dbetas(kk)\n",
    "            else: print \"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 dlogPd_dns(kk)\n",
    "                    elif j==1:  return dlogPd_dalphas(kk)\n",
    "                    elif j==2:  return dlogPd_dbetas(kk)\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]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2664131.9402093   1340721.35385563]\n",
      " [ 1340721.35385563   743896.81698823]]\n"
     ]
    }
   ],
   "source": [
    "print Fishermat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  4.03621557e-06  -7.27445028e-06]\n",
      " [ -7.27445028e-06   1.44549763e-05]]\n"
     ]
    }
   ],
   "source": [
    "print linalg.inv(Fishermat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.00200903349255\n",
      "0.00380197004973\n"
     ]
    }
   ],
   "source": [
    "#print marginalised uncertainties on (ns,alphas,betas)\n",
    "for i in range(0,Npar):\n",
    "    print np.sqrt(linalg.inv(Fishermat)[i,i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "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
}