nanograv11yr_makehdf5.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-10-17T04:57:09.627662Z",
     "start_time": "2017-10-17T04:57:08.655422Z"
    }
   },
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "can't multiply sequence by non-int of type 'float'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-78e8810e44fa>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpyplot\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 10\u001b[0;31m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrcParams\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'savefig.dpi'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m1.5\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrcParams\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'savefig.dpi'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mnumpy\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: can't multiply sequence by non-int of type 'float'"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "%load_ext Cython\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "from __future__ import division\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "matplotlib.rcParams['savefig.dpi'] = 1.5 * matplotlib.rcParams['savefig.dpi']\n",
    "\n",
    "import numpy as np\n",
    "import sys,os,glob,h5py\n",
    "import libstempo as T2\n",
    "import libstempo.plot as LP\n",
    "\n",
    "import NX01_psr\n",
    "import NX01_datafile\n",
    "\n",
    "try:\n",
    "    from IPython.core.display import clear_output\n",
    "    have_ipython = True\n",
    "except ImportError:\n",
    "    have_ipython = False\n",
    "\n",
    "Msol = 1.98855*10.0**30.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parentpath = '/Users/staylor/Research/NANOGrav/stochastic_11yr_analysis/data/'\n",
    "\n",
    "parpath = parentpath + 'partim_no_noise'\n",
    "timpath = parentpath + 'partim_no_noise'\n",
    "noisepath = parentpath + 'nano_11_noisefiles_30_marg'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the parfiles and timfiles\n",
    "parfiles = sorted(glob.glob(parpath+'/*.par'))\n",
    "timfiles = sorted(glob.glob(timpath+'/*.tim'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the noise files\n",
    "noisefiles = sorted(glob.glob(noisepath+'/*.txt'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(parfiles), len(timfiles), len(noisefiles)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## This code block is a one-time deal to make par files stripped of EFACS, EQUADS, and ECORRs (don't do this if you have already got stripped par files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "stripped_pars = list(parfiles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in range(len(stripped_pars)):\n",
    "    stripped_pars[ii] = stripped_pars[ii].replace('9yv1.gls.par', '9yv1.gls.strip.par')\n",
    "    stripped_pars[ii] = stripped_pars[ii].replace('9yv1.t2.gls.par', '9yv1.t2.gls.strip.par')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for ii in range(len(stripped_pars)):\n",
    "    os.system('awk \\'($1 !~ /T2EFAC/ && $1 !~ /T2EQUAD/ && $1 !~ /ECORR/ && $1 !~ /RNAMP/ && $1 !~ /RNIDX/ ) {{print $0}} \\' {0} > {1}'.format(parfiles[ii],stripped_pars[ii]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "parfiles = sorted(glob.glob(parpath+'/*.strip.par'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now moving on to processing these pulsars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#######################################\n",
    "# PASSING THROUGH TEMPO2 VIA libstempo\n",
    "#######################################\n",
    "\n",
    "t2psr = []\n",
    "for ii in range(len(parfiles)):\n",
    "    \n",
    "    t2psr.append( T2.tempopulsar(parfile = parfiles[ii], timfile = timfiles[ii], \n",
    "                                 maxobs=30000, ephem='DE421') )\n",
    "    \n",
    "    #if np.any(np.isfinite(t2psr[ii].residuals())==False)==True:\n",
    "    #    t2psr[ii] = T2.tempopulsar(parfile = parfiles[ii], timfile = timfiles[ii])\n",
    "                 \n",
    "    clear_output()\n",
    "    print '\\r', '{0} of {1}'.format(ii+1,len(parfiles))\n",
    "    sys.stdout.flush()                 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check out some plots if you want\n",
    "#LP.plotres(t2psr[0])\n",
    "plt.errorbar(t2psr[0].toas(),t2psr[0].residuals()/1e-6,t2psr[0].toaerrs,fmt='.',alpha=0.2)\n",
    "plt.xlabel(r'MJD')\n",
    "plt.ylabel(r'Residuals [$\\mu$s]')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Cell for time-slicing of data\n",
    "\n",
    "Nyears = 'tot'\n",
    "\n",
    "Tmin = np.min([np.min(p.toas()) for p in t2psr])\n",
    "if Nyears == 'tot':\n",
    "    Tcutoff = np.max([np.max(p.toas()) for p in t2psr])\n",
    "else:\n",
    "    Tcutoff = Tmin + Nyears*365.25\n",
    "\n",
    "ind_slice = []\n",
    "for ii,p in enumerate(t2psr):\n",
    "    \n",
    "    mask = np.logical_and(p.toas() >= Tmin, p.toas() <= Tcutoff)\n",
    "    \n",
    "    if np.sum(mask) > 0:\n",
    "        Tobs = (p.toas()[mask].max() - p.toas()[mask].min()) / 365.25\n",
    "        \n",
    "        if Tobs > 1.0:\n",
    "            ind_slice.append(ii)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Pass all tempopulsar objects to NX01 pulsar class\n",
    "psr = [NX01_psr.PsrObj(t2psr[ii]) for ii in ind_slice]\n",
    "[p.grab_all_vars(jitterbin=1.0,makeGmat=False,\n",
    "                 fastDesign=False,planetssb=True,\n",
    "                 startMJD=Tmin, endMJD=Tcutoff) for p in psr] # according to the 9 year paper,\n",
    "                                                              # the jitterbin used is 1s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Pass all tempopulsar objects to NX01 pulsar class\n",
    "#psr = [NX01_psr.PsrObj(p) for p in t2psr]\n",
    "#[p.grab_all_vars(jitterbin=1.0,makeGmat=False,\n",
    "#                 fastDesign=False,planetssb=True,\n",
    "#                 startMJD=Tmin, endMJD=Tcutoff) for p in psr] # according to the 9 year paper,\n",
    "#                                                              # the jitterbin used is 1s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fill in the locations of par, tim, and noise files\n",
    "for ii,jj in enumerate(ind_slice):\n",
    "    psr[ii].parfile = parfiles[jj] \n",
    "    psr[ii].timfile = timfiles[jj]\n",
    "    psr[ii].noisefile = noisefiles[jj]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Only need to execute if you want roemer delays \n",
    "# from other ephemerides.\n",
    "# Dummy libstempo passes to grab delays.\n",
    "\n",
    "ephems = ['DE421', 'DE430', 'DE435', 'DE436']\n",
    "\n",
    "for ii,jj in enumerate(ind_slice):\n",
    "    \n",
    "    for eph in ephems:\n",
    "        if eph != psr[ii].ephemname:\n",
    "    \n",
    "            dummy_t2psr = T2.tempopulsar(parfile = parfiles[jj], timfile = timfiles[jj], \n",
    "                                         maxobs=30000, ephem=eph)\n",
    "        \n",
    "            psr[ii].roemer[eph] = np.double(dummy_t2psr.roemer)\n",
    "            if psr[ii].tmask is not None:\n",
    "                psr[ii].roemer[eph] = psr[ii].roemer[eph][psr[ii].tmask]\n",
    "            if psr[ii].isort is not None:\n",
    "                psr[ii].roemer[eph] = psr[ii].roemer[eph][psr[ii].isort]\n",
    "            \n",
    "                 \n",
    "    clear_output()\n",
    "    print '\\r', '{0} of {1}'.format(ii+1,len(parfiles))\n",
    "    sys.stdout.flush()                 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Make HDF5 file directory\n",
    "dirname = os.getcwd() + '/11yr_psr_hdf5_files/DE421/tm_svdstabilized/'\n",
    "if not os.path.exists(dirname):\n",
    "    os.makedirs(dirname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dump all pulsars into HDF5 files \n",
    "for ii,p in enumerate(psr):\n",
    "    df = NX01_datafile.DataFile(dirname + '/' + p.name + '.hdf5')\n",
    "    df.addTempoPulsar(p)\n",
    "    \n",
    "    clear_output()\n",
    "    print '\\r', '{0} of {1}'.format(ii+1,len(psr))\n",
    "    sys.stdout.flush()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create information text file for pulsar hdf5 file locations\n",
    "fil = open(dirname + '/psrList_nano11yr.txt','w')\n",
    "print >>fil, \"NAME\", \"HDF5-PATH\", \"PARFILE-PATH\", \"TIMFILE-PATH\"\n",
    "print >>fil, \"#############################################\"\n",
    "for p in psr:\n",
    "    print >>fil, p.name, dirname+'/'+p.name+'.hdf5'\n",
    "fil.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rankings = np.genfromtxt(parentpath+'psrlist.txt',dtype=str)\n",
    "fil = open(dirname + '/psrList_nano11_ranked.txt','w')\n",
    "print >>fil, \"NAME\", \"HDF5-PATH\", \"PARFILE-PATH\", \"TIMFILE-PATH\"\n",
    "print >>fil, \"#############################################\"\n",
    "for p in rankings:\n",
    "    print >>fil, p, dirname+'/'+p+'.hdf5'\n",
    "fil.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rankings = np.genfromtxt(parentpath+'psrlist_Tg3yr.txt',dtype=str)\n",
    "fil = open(dirname + '/psrList_nano11_ranked_Tg3yr.txt','w')\n",
    "print >>fil, \"NAME\", \"HDF5-PATH\", \"PARFILE-PATH\", \"TIMFILE-PATH\"\n",
    "print >>fil, \"#############################################\"\n",
    "for p in rankings:\n",
    "    print >>fil, p, dirname+'/'+p+'.hdf5'\n",
    "fil.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "gmu_uls = np.array([3.87595541e-13, 4.71210227e-13, 5.71803686e-13,\n",
    "                    6.95797061e-13, 8.32522097e-13, 1.02163421e-12,\n",
    "                    1.26938114e-12, 1.54532333e-12, 1.86832268e-12,\n",
    "                    2.25810134e-12, 2.83286775e-12, 3.51803788e-12,\n",
    "                    4.52328994e-12, 6.02546201e-12, 7.67416070e-12,\n",
    "                    1.04064636e-11, 1.39940969e-11, 2.06332855e-11,\n",
    "                    3.08028537e-11, 5.48868339e-11])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.loglog(gmu_uls,10.0**prob)\n",
    "plt.fill_betweenx(10.0**prob,gmu_uls,1e-9,alpha=0.4)\n",
    "plt.xlim(1e-13,1e-9)\n",
    "plt.ylim(1e-3,1.0)\n",
    "plt.xlabel(r'$G\\mu / c^2$')\n",
    "plt.ylabel(r'$p$')\n",
    "plt.show()"
   ]
  }
 ],
 "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.13"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}