space.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Header files:\n",
    "#     mytypes.h - defines constant multiples of PI and the absolute value function, all available within Python\n",
    "\n",
    "#     int_de.h - defines a function:\n",
    "#         int int_de()\n",
    "\n",
    "#     calcpath.h - declares functions and defines the following constants\n",
    "# NEQS = 8 \n",
    "# kmax = 20000 \n",
    "\n",
    "#     randoms.h - declares the following functions:\n",
    "#         double rrandom()\n",
    "    \n",
    "#     spectrum.h - declares one struct:\n",
    "#         struct params\n",
    "#     and the following functions:\n",
    "#         int spectrum\n",
    "#         int derivs1\n",
    "#         int scalarsys\n",
    "#         int tensorsys\n",
    "#     and the following constants:\n",
    "# Y = 50                  # Y = value of k/aH at which to initialize mode fcns \n",
    "# Amp = 2.0803249e-9      # scalar amplitude at knorm \n",
    "# knorm = 0.05            # normalization scale\n",
    "# ki_file = \"ks.dat\"      # file containing k-values for integration\n",
    "# k_file = \"ks_eval.dat\"  # file containing k-values at which to evaluate spectrum\n",
    "# kinos = 214             # total number of k-values to use for integration \n",
    "# knos = 1575             # total number of k-values to evaluate\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy\n",
    "import random\n",
    "from scipy.integrate import solve_ivp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 218,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "OUTFILE1_NAME = \"nr.dat\"  \n",
    "OUTFILE2_NAME = \"esigma.dat\"\n",
    "\n",
    "NEQS = 8\n",
    "kmax = 20000\n",
    "\n",
    "NUMPOINTS = 20000\n",
    "\n",
    "NUMEFOLDSMAX = 60.0\n",
    "NUMEFOLDSMIN = 40.0\n",
    "\n",
    "SMALLNUM = 0.0001\n",
    "VERYSMALLNUM = 1e-18\n",
    "LOTSOFEFOLDS = 1000.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Calc:\n",
    "    def __init__(self):\n",
    "        self.Y = numpy.zeros (NEQS, dtype=float, order='C')\n",
    "        self.initY = numpy.zeros (NEQS, dtype=float, order='C')\n",
    "        self.ret = \"\"\n",
    "        self.npoints = 0\n",
    "        self.Nefolds = 0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Fills vector * yinit with randomly chosen initial values of the flow\n",
    "# parameters, as well as a randomly chosen amount of inflation, Nefolds.\n",
    "\n",
    "# POSSIBLE ISSUE WITH EXACT DECIMAL NUMBERS\n",
    "\n",
    "def pick_init_vals ():\n",
    "    init_vals = numpy.zeros (NEQS, dtype=float, order='C')\n",
    "    \n",
    "    init_vals[0] = 0.0\n",
    "    init_vals[1] = 1.0\n",
    "    init_vals[2] = random.uniform(0.0 , 0.8)\n",
    "    init_vals[3] = random.uniform(-0.5 , 0.5)\n",
    "    init_vals[4] = random.uniform(-0.05 , 0.05)\n",
    "    \n",
    "    for i in range (5 , NEQS):\n",
    "        upper_val = 1 * 0.05 * 0.1 ** (i-5)  - 0.5 * 0.05 * 0.1 ** (i-5)\n",
    "        init_vals[i] = random.uniform(-upper_val , upper_val)\n",
    "        \n",
    "    init_Nefolds = random.uniform(NUMEFOLDSMIN , NUMEFOLDSMAX)\n",
    "    \n",
    "    return init_vals , init_Nefolds\n",
    "\n",
    "\n",
    "def derivs (t, y):\n",
    "    dydN = numpy.zeros(NEQS , dtype=float , order='C')\n",
    "    \n",
    "    if y[2] >= 1.0:\n",
    "        dydN = numpy.zeros(NEQS , dtype=float , order='C')\n",
    "    else:\n",
    "        if y[2] > VERYSMALLNUM:\n",
    "            dydN[0] = - numpy.sqrt ( y[2] / (4 * numpy.pi ) )\n",
    "        else:\n",
    "            dydN[0] = 0.0\n",
    "        \n",
    "        dydN[1] = y[1] * y[2]\n",
    "        dydN[2] = y[2] * ( y[3] + 2.0 * y[2] )\n",
    "        dydN[3] = 2.0 * y[4] - 5.0 * y[2] * y[3] - 12.0 * y[2] * y[2]\n",
    "        \n",
    "        for i in range (4 , NEQS-1):\n",
    "            dydN[i] = ( 0.5 * (i-3) * y[3] + (i-4) * y[2] ) * y[i] + y[i+1]\n",
    "            \n",
    "        dydN[NEQS-1] = ( 0.5 * (NEQS-4) * y[3] + (NEQS-5) * y[2] ) * y[NEQS-1]\n",
    "    \n",
    "    return dydN\n",
    "\n",
    "def int_de (y , N , Nend , derivs):\n",
    "    h = 1e-6\n",
    "    ydoub = y.copy()\n",
    "        \n",
    "    Nsol = solve_ivp(derivs, (N , Nend), ydoub, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False , first_step=h)\n",
    "    \n",
    "    sol_status = Nsol['status']\n",
    "    sol_length = len ( Nsol['t'] )\n",
    "    sol_x = Nsol['t']\n",
    "    sol_Y = Nsol['y']\n",
    "    \n",
    "    if sol_length > kmax:\n",
    "        sol_length = kmax\n",
    "        sol_x = sol_x [0:kmax]\n",
    "        sol_Y = sol_Y [: , 0:kmax]\n",
    "        \n",
    "    return sol_status , sol_length , sol_x , sol_Y\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 229,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "nontrivial\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "nontrivial\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "nontrivial\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "nontrivial\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "nontrivial\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "nontrivial\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      "asymptote\n",
      ".\n",
      "asymptote\n",
      "asymptote\n"
     ]
    }
   ],
   "source": [
    "calc = Calc()\n",
    "\n",
    "path = numpy.zeros (( NEQS , kmax ))\n",
    "\n",
    "# Open output files\n",
    "\n",
    "try:\n",
    "    outfile1 = open (OUTFILE1_NAME, \"w\")\n",
    "except IOError as e:\n",
    "    print (\"Could not open file\" , str(OUTFILE1_NAME) , \", errno =\" , str(e) , \".\")\n",
    "    \n",
    "try:\n",
    "    outfile2 = open (OUTFILE2_NAME, \"w\")\n",
    "except IOError as e:\n",
    "    print (\"Could not open file\" , str(OUTFILE2_NAME) , \", errno =\" , str(e) , \".\")\n",
    "    \n",
    "# Allocate buffers\n",
    "\n",
    "y = numpy.zeros(NEQS, dtype=float, order='C')\n",
    "yinit = numpy.zeros(NEQS, dtype=float, order='C')\n",
    "\n",
    "# iters = total number of iterations\n",
    "# points = points saved with n < NMAX\n",
    "# asymcount = points with 0 < n < NMAX , r = 0\n",
    "# nontrivcount = nontrivial points\n",
    "# insuffcount = points where slow roll breaks down before N efolds\n",
    "# noconvcount = points that do not converge to either a late\n",
    "# time attractor or end of inflation.\n",
    "    \n",
    "iters = 0\n",
    "points = 0\n",
    "errcount = 0\n",
    "outcount = 0\n",
    "asymcount = 0\n",
    "nontrivcount = 0\n",
    "insuffcount = 0\n",
    "noconvcount = 0\n",
    "badncount = 0\n",
    "savedone = 0\n",
    "\n",
    "# Currently this loop is set to depend on the number of\n",
    "# nontrivial points.  This can be changed to depend on total\n",
    "# number of models, or other criteria.\n",
    "\n",
    "deyancount = 0\n",
    "\n",
    "while nontrivcount < NUMPOINTS:\n",
    "    deyancount += 1\n",
    "    \n",
    "    iters += 1\n",
    "    \n",
    "    if iters % 100 == 0:\n",
    "        if iters % 1000 == 0:\n",
    "            print (\"\\n asymcount =\" , asymcount , \",nontrivcount =\" , nontrivcount , \",insuffcount =\" , insuffcount , \", noconvcount =\" , noconvcount , \", badncount =\" , badncount , \", errcount =\" , errcount , \"\\n\")\n",
    "            print (\"\\n\", iters)\n",
    "        else:\n",
    "            print (\".\")\n",
    "            \n",
    "# Select a random initial condition such that the slow\n",
    "# roll hierarchy converges.\n",
    "\n",
    "    yinit , calc.Nefolds = pick_init_vals()\n",
    "    \n",
    "    y = yinit.copy()\n",
    "    \n",
    "    calc.ret = calcpath ( calc.Nefolds , y , path , calc.npoints)\n",
    "    \n",
    "    if deyancount > 100:\n",
    "        break\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 223,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def calcpath (Nefolds , y , path , count):\n",
    "# \tdouble Nstart,Nend;\n",
    "# \tgsl_vector * xp;\n",
    "# \tgsl_matrix * yp;\n",
    "# \tenum calcpath_retcode retval=internal_error;\n",
    "# \tint i,j,k,z,kount;\n",
    "# \tdouble Hnorm;\n",
    "\n",
    "    retval = \"internal_error\"\n",
    "    \n",
    "    if NEQS < 6:\n",
    "        raise Exception(\"calcpath(): NEQS must be at least 6\\n\")\n",
    "    \n",
    "    yp = numpy.zeros((NEQS , kmax) , dtype=float , order='C')\n",
    "    xp = numpy.zeros(kmax , dtype=float , order='C')\n",
    "    \n",
    "    Nstart = LOTSOFEFOLDS\n",
    "    Nend = 0.0\n",
    "    \n",
    "    z , kount , xp , yp = int_de ( y , Nstart , Nend , derivs )\n",
    "    \n",
    "    y = yp [ : , -1 ]\n",
    "    \n",
    "    if z < 0:\n",
    "        retval = \"internal_error\"\n",
    "        z = 0\n",
    "#         goto end;\n",
    "    else:\n",
    "#         Find when epsilon passes through unity.\n",
    "        i = check_convergence ( yp , kount )\n",
    "        \n",
    "        if not i:\n",
    "#             We never found an end to inflation, so we must be at a late-time attractor.\n",
    "            if y[2] > SMALLNUM or y[3] < 0.0:\n",
    "#                 The system did not evolve to a known asymptote.\n",
    "                retval = \"noconverge\"\n",
    "            else:\n",
    "                retval = \"asymptote\"\n",
    "        else:\n",
    "#             We found an end to inflation: integrate backwards Nefolds e-folds from that point.\n",
    "            Nstart = xp[i-2] - xp[i-1]\n",
    "            Nend = Nefolds\n",
    "            \n",
    "            y = yp [ : , i-2 ].copy()\n",
    "            \n",
    "            z , kount , xp , yp = int_de ( y , Nstart , Nend , derivs )\n",
    "            \n",
    "            if z < 0:\n",
    "                retval = \"internal_error\"\n",
    "                z = 0\n",
    "            else:\n",
    "                if check_convergence ( yp , kount ):\n",
    "                    retval = \"insuff\"\n",
    "                else:\n",
    "                    retval = \"nontrivial\"\n",
    "                    \n",
    "    print (retval)\n",
    "                \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def check_convergence (yy , kount):\n",
    "    for i in range (kount):\n",
    "        if numpy.absolute ( yy [ 2 , i ] ) >= 1.0:\n",
    "            return i\n",
    "            break\n",
    "        \n",
    "    return 0\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.081451381839780623"
      ]
     },
     "execution_count": 140,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "4 * (numpy.log(2) + 0.5772156649) - 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a = numpy.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10])"
      ]
     },
     "execution_count": 146,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}