test_nuggets.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python\n",
    "\n",
    "\n",
    "\"\"\" McNuggets: the post-processing tool for brewster\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import forwardmodel\n",
    "import ciamod\n",
    "import pickle\n",
    "from scipy.io.idl import readsav\n",
    "from scipy import interpolate\n",
    "from scipy.interpolate import interp1d\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def teffRM(theta,runargs):\n",
    "    \n",
    "    invmr = theta[0:7]\n",
    "    logf = 0.0 #theta[5]\n",
    "    #logbeta = theta[6]\n",
    "    logg = theta[7]\n",
    "    r2d2 = theta[8]\n",
    "    dlam = theta[9]\n",
    "    gam = theta[10]\n",
    "    intemp = theta[11:]\n",
    "    pcover, cloudparams,do_clouds,gasnum,cloudnum,inlinetemps,\\\n",
    "        coarsePress,press,inwavenum,linelist,cia,ciatemps,\\\n",
    "        use_disort = runargs\n",
    "\n",
    "\n",
    "    # interp temp onto finer grid coarsePress => press\n",
    "    # spline fit with max smoothing\n",
    "    tfit = sp.interpolate.splrep(np.log10(coarsePress),np.log10(intemp),s=0)\n",
    "    temp = 10.**(np.asfortranarray(sp.interpolate.splev(np.log10(press),tfit,der=0),dtype='d'))\n",
    "\n",
    "    # get the ngas\n",
    "    ngas = invmr.shape[0] + 1\n",
    "    # Hard code nlayers\n",
    "    nlayers = press.shape[0]\n",
    "    # interp temp onto finer grid coarsePress => press\n",
    "    # spline fit with max smoothing\n",
    "    tfit = sp.interpolate.splrep(np.log10(coarsePress),intemp,s=0)\n",
    "    temp = np.asfortranarray(sp.interpolate.splev(np.log10(press),tfit,der=0),dtype='d')\n",
    "    # now loop through gases and get VMR for model\n",
    "    # check if its a fixed VMR or a profile\n",
    "    # VMR is log10(VMR) !!!\n",
    "    logVMR = np.empty((ngas,nlayers),dtype='d')\n",
    "    alkratio = 16.2 #  from Asplund et al (2009)\n",
    "    if invmr.size > invmr.shape[0]:\n",
    "        # now sort Na and K\n",
    "        tmpvmr = np.empty((ngas,nlayers),dtype='d')\n",
    "        tmpvmr[0:(ngas-2),:] = invmr[0:(ngas-2),:]\n",
    "        tmpvmr[ngas-2,:] = np.log10(10.**invmr[ngas-2,:] / (alkratio+1.))\n",
    "        tmpvmr[ngas-1,:] = np.log10(10.**invmr[ngas-2,:] * (alkratio / (alkratio+1.)))                                \n",
    "        for i in range(0,ngas):\n",
    "            vfit = sp.interpolate.splrep(np.log10(coarsepress),tmpvmr[i,:],s=0)\n",
    "            logVMR[i,:] = sp.interpolate.splev(np.log10(press),vfit,der=0)\n",
    "    else:\n",
    "        # now sort Na and K\n",
    "        tmpvmr = np.empty(ngas,dtype='d')\n",
    "        tmpvmr[0:(ngas-2)] = invmr[0:(ngas-2)]\n",
    "        tmpvmr[ngas-2] = np.log10(10.**invmr[ngas-2] / (alkratio+1.))\n",
    "        tmpvmr[ngas-1] = np.log10(10.**invmr[ngas-2] * (alkratio / (alkratio+1.)))\n",
    "        for i in range(0,ngas):                              \n",
    "            logVMR[i,:] = tmpvmr[i]\n",
    "\n",
    "    # now need to translate cloudparams in to cloud profile even\n",
    "    # if do_clouds is zero..\n",
    "    # 5 entries for cloudparams for simple slab model are:\n",
    "    # 0) log10(number density)\n",
    "    # 1) top layer id (or pressure)\n",
    "    # 2) base ID (these are both in 61 layers)\n",
    "    # 3) rg\n",
    "    # 4) rsig\n",
    "    if (do_clouds == 1):\n",
    "        npatch = cloudparams.shape[0]\n",
    "        ncloud = cloudparams.shape[1]\n",
    "        cloudrad = np.empty((npatch,nlayers,ncloud),dtype='d')\n",
    "        cloudsig = np.empty_like(cloudrad)\n",
    "        cloudprof = np.zeros_like(cloudrad)\n",
    "        ndens= np.reshape(cloudparams['f0'],(npatch,ncloud))\n",
    "        c1 = np.reshape(cloudparams['f1'],(npatch,ncloud))\n",
    "        c2 = np.reshape(cloudparams['f2'],(npatch,ncloud))\n",
    "        rad = np.reshape(cloudparams['f3'],(npatch,ncloud))\n",
    "        sig = np.reshape(cloudparams['f4'],(npatch,ncloud))\n",
    "        for i in range(0, npatch):\n",
    "            for j in range(0, ncloud):\n",
    "                b1 = c1[i,j] - 1\n",
    "                b2 = c2[i,j] -1 \n",
    "                cloudprof[i,b1:b2+1,j] = ndens[i,j]\n",
    "                cloudrad[i,:,j] = rad[i,j]\n",
    "                cloudsig[i,:,j] = sig[i,j]        \n",
    "    else:\n",
    "        npatch = 1\n",
    "        ncloud = 1\n",
    "        cloudrad = np.ones((npatch,nlayers,ncloud),dtype='d')\n",
    "        cloudsig = np.ones_like(cloudrad)\n",
    "        cloudprof = np.ones_like(cloudrad)\n",
    "\n",
    "    modcall = [temp,logg,r2d2,gasnum,logVMR,pcover,do_clouds,cloudnum,cloudrad,cloudsig,cloudprof,\\\n",
    "               inlinetemps,press,inwavenum,linelist,cia,ciatemps,use_disort]\n",
    "    # now we can call the forward model\n",
    "    outspec = forwardmodel.marv(*modcall)\n",
    "\n",
    "    wave = np.array(outspec[0,::-1])\n",
    "    flux = np.array(outspec[1,::-1])\n",
    "\n",
    "    # now calculate Fbol by summing the spectrum across its wave bins\n",
    "    fbol = 0.0\n",
    "    \n",
    "    for j in range(1, (wave.size - 1)):\n",
    "        sbin = ((wave[j] - wave[j-1]) + (wave[j+1] - wave[j])) / 2. \n",
    "    \n",
    "        fbol = (sbin * flux[j]) + fbol\n",
    "\n",
    "    # now get T_eff\n",
    "    t_ff = ((fbol/(r2d2 * 5.670367e-8))**(1./4.))\n",
    "\n",
    "    # and Radius\n",
    "    parallax = 5.84\n",
    "    sigpi = 0.03\n",
    "    sigphot = 0.02\n",
    "    \n",
    "    sigR2D2 = sigphot * r2d2 * (-1./2.5)* np.log(10.)\n",
    "\n",
    "    sigD = sigpi * 3.086e16\n",
    "    D = parallax * 3.086e16\n",
    "\n",
    "    R = np.sqrt(((np.random.randn() * sigR2D2)+ r2d2)) \\\n",
    "        * ((np.random.randn()* sigD) + D)\n",
    "\n",
    "    g = (10.**logg)/100.\n",
    "\n",
    "    # and mass\n",
    "\n",
    "    M = (R**2 * g/(6.67E-11))/1.898E27\n",
    "    R = R / 71492e3\n",
    "    result = np.concatenate((theta,np.array([t_ff, R, M])),axis=0)\n",
    "    \n",
    "    return result\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: '../run_analysis/570D_retrieval_result.pk1'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-3-3a58fde07192>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mwith\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'../run_analysis/570D_retrieval_result.pk1'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'rb'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0minput\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m     \u001b[0msampler\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpickle\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mload\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minput\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mndim\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msampler\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mchain\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0msamples\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msampler\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mchain\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m15000\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mndim\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: '../run_analysis/570D_retrieval_result.pk1'"
     ]
    }
   ],
   "source": [
    "with open('./Gl229B_nocloud.pk1', 'rb') as input:\n",
    "    sampler = pickle.load(input) \n",
    "\n",
    "ndim = sampler.chain.shape[2]\n",
    "samples = sampler.chain[:,15000:,:].reshape((-1, ndim))\n",
    "\n",
    "slen = samples.shape[0]\n",
    "#samplus = np.zeros([slen,ndim+3])\n",
    "\n",
    "#samplus[:,0:ndim] = samples\n",
    "#print samplus[0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up run arguments\n",
    "\n",
    "w1 = 1.0\n",
    "w2 = 20.0\n",
    "pcover = 1.0\n",
    "do_clouds = 0\n",
    "use_disort = 0\n",
    "\n",
    "\n",
    "# set up pressure grids\n",
    "logcoarsePress = np.arange(-4.0, 2.5, 0.53)\n",
    "#logcoarsePress = np.arange(-4.0, 3.0, 0.5)\n",
    "coarsePress = 1000.* pow(10,logcoarsePress)\n",
    "logfinePress = np.arange(-4.0, 2.4, 0.1)\n",
    "finePress = 1000.* pow(10,logfinePress)\n",
    "press = finePress\n",
    "\n",
    "# now the linelist\n",
    "# Set up number of gases, and point at the lists. see gaslist.dat\n",
    "ngas = 8\n",
    "gasnum = np.asfortranarray(np.array([1,2,4,5,6,3,20,21],dtype='i'))\n",
    "lists = [\"../Linelists/xsecarrH2O_1wno_500_10000.save\",\"../Linelists/xsecarrCH4_1wno_500_10000.save\",\"../Linelists/xsecarrCO_1wno_500_10000_02.save\",\"../Linelists/xsecarrCO2_1wno_500_10000_02.save\" ,\"../Linelists/xsecarrNH3_1wno_500_10000.save\",\"../Linelists/xsecarrH2S_1wno_500_10000.save\",\"../Linelists/xsecarrK_new_1wno_500_10000_02.save\",\"../Linelists/xsecarrNa_new_1wno_500_10000_02.save\"]\n",
    "# get the basic framework from water list\n",
    "x=readsav('../Linelists/xsecarrH2O_1wno_500_10000.save')\n",
    "inlinelist=x.xsecarr  #3D array with Nwavenubmers x Ntemps x Npressure\n",
    "inlinetemps=np.asfortranarray(x.t,dtype='float64')\n",
    "inpress=1000.*x.p\n",
    "rawwavenum=x.wno\n",
    "wn1 = 10000./w2\n",
    "wn2 = 10000. / w1\n",
    "inwavenum = np.asfortranarray(rawwavenum[np.where(np.logical_not(np.logical_or(rawwavenum[:] > wn2, rawwavenum[:] < wn1)))],dtype='float64')\n",
    "ntemps = inlinetemps.size\n",
    "npress= finePress.size\n",
    "nwave = inwavenum.size\n",
    "r1 = np.amin(np.where(np.logical_not(np.logical_or(rawwavenum[:] > wn2, rawwavenum[:] < wn1))))\n",
    "r2 = np.amax(np.where(np.logical_not(np.logical_or(rawwavenum[:] > wn2, rawwavenum[:] < wn1))))\n",
    "\n",
    "# Here we are interpolating the linelist onto our fine pressure scale. \n",
    "linelist = (np.ones([ngas,npress,ntemps,nwave],order='F')).astype('float64', order='F')\n",
    "for gas in range (0,ngas):\n",
    "    inlinelist=readsav(lists[gas]).xsecarr\n",
    "    for i in range (0,ntemps):\n",
    "        for j in range (r1,r2+1):\n",
    "            pfit = interp1d(np.log10(inpress),np.log10(inlinelist[:,i,j]))\n",
    "            linelist[gas,:,i,(j-r1)] = np.asfortranarray(pfit(np.log10(finePress)))\n",
    "\n",
    "linelist[np.isnan(linelist)] = -50.0\n",
    "\n",
    "# cloudparams is structured array with 5 entries\n",
    "# each one has a patch*cloud entries\n",
    "cloudparams = np.ones(5)\n",
    "# 5 entries in cloudparams for simple slab model are:\n",
    "# 0) log10(number density)\n",
    "# 1) top layer id (or pressure)\n",
    "# 2) base ID (these are both in 61 layers)\n",
    "# 3) rg\n",
    "# 4) rsig\n",
    "cloudparams[0] = -20.\n",
    "cloudparams[1] = 10\n",
    "cloudparams[2] = 12\n",
    "cloudparams[3] = 1e-4\n",
    "cloudparams[4] = 1e-5\n",
    "# hardwired gas and cloud IDs\n",
    "cloudnum = np.array([1],dtype='i')\n",
    "\n",
    "# Get the cia bits\n",
    "tmpcia, ciatemps = ciamod.read_cia(\"CIA_DS_aug_2015.dat\",inwavenum)\n",
    "cia = np.asfortranarray(np.empty((4,ciatemps.size,nwave)),dtype='float32')\n",
    "cia[:,:,:] = tmpcia[:,:,:nwave] \n",
    "ciatemps = np.asfortranarray(ciatemps, dtype='float32')\n",
    "\n",
    "runargs = pcover, cloudparams,do_clouds,gasnum,cloudnum,inlinetemps,coarsePress,press,inwavenum,linelist,cia,ciatemps,use_disort\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -3.28799047e+00,  -3.19601596e+00,  -1.05987170e+01,\n",
       "        -6.58080820e+00,  -4.42955954e+00,  -6.14382385e+00,\n",
       "        -5.48272079e+00,   5.22060480e+00,   1.83316777e-19,\n",
       "         1.41338934e-03,   1.70582381e+05,   9.84474294e+02,\n",
       "         8.44922181e+02,   5.93459204e+02,   4.88045942e+02,\n",
       "         4.88125470e+02,   6.82786148e+02,   4.60758758e+02,\n",
       "         6.09931653e+02,   6.71726569e+02,   9.53092268e+02,\n",
       "         1.25451657e+03,   1.58464588e+03,   1.47568056e+03,\n",
       "         7.63414910e+02,   1.09211839e+00,   8.00271026e+01])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "teffRM(samples[1,:],runargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1000000, 24)\n"
     ]
    }
   ],
   "source": [
    "print samples.shape\n",
    "results = []\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  -3.28799047   -3.19601596  -10.59871702 ...,  763.41490958\n",
      "     1.09280919   80.12837437]\n",
      " [  -3.28799047   -3.19601596  -10.59871702 ...,  763.41490958\n",
      "     1.09003059   79.72142068]\n",
      " [  -3.28799047   -3.19601596  -10.59871702 ...,  763.41490958\n",
      "     1.08400631   78.84266109]\n",
      " ..., \n",
      " [  -3.26480176   -3.17703479   -7.50721097 ...,  777.25621832\n",
      "     1.04712037   78.38416175]\n",
      " [  -3.26586804   -3.17750718   -7.46675549 ...,  777.01852628\n",
      "     1.03500999   76.42277237]\n",
      " [  -3.26586804   -3.17750718   -7.46675549 ...,  777.01852628\n",
      "     1.02983568   75.66056369]]\n"
     ]
    }
   ],
   "source": [
    "for i in range (0,200):\n",
    "    results.append(teffRM(samples[i,0:ndim],runargs))     \n",
    "samplus = np.array(results)\n",
    "print samplus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[   -3.28799047    -3.19601596   -10.59871702 ...,  1254.51656828\n",
      "   1584.64588311  1475.6805571 ]\n",
      " [   -3.28799047    -3.19601596   -10.59871702 ...,  1254.51656828\n",
      "   1584.64588311  1475.6805571 ]\n",
      " [   -3.28799047    -3.19601596   -10.59871702 ...,  1254.51656828\n",
      "   1584.64588311  1475.6805571 ]\n",
      " ..., \n",
      " [   -3.270054      -3.21041688    -9.09278101 ...,  1278.51371455\n",
      "   1588.88456932  1770.25438089]\n",
      " [   -3.270054      -3.21041688    -9.09278101 ...,  1278.51371455\n",
      "   1588.88456932  1770.25438089]\n",
      " [   -3.270054      -3.21041688    -9.09278101 ...,  1278.51371455\n",
      "   1588.88456932  1770.25438089]]\n"
     ]
    }
   ],
   "source": [
    "print samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1000000, 24)\n"
     ]
    }
   ],
   "source": [
    "print samples.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}