post qlens parameters - cnfw pmode 3 (mvir, c, rc_kpc) - get  rs, rc, rcore, c_hat_200 and rm2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gravitational Lensing Research\n",
    "# Program to take a qlens-generated sample chain and derive alternative parameters\n",
    "# d3: This version takes the cnfw pmode=1 format (mvir, c, beta) and derives scale radius (r_s),\n",
    "# cnfw param r_c (r_c), 'core radius' (r_core, defined as where log slope=-1) and revised concentration (c_hat)\n",
    "# c_hat_200 is the defined as r_vir / r_{-2}, and r_{-2} is where the log slope is -2.\n",
    "# Uses python 3.6\n",
    "# input from cNFW needed is mvir, c, rc_kpc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "import numpy as np\n",
    "import os\n",
    "import pandas as pd\n",
    "import timeit\n",
    "from IPython.core.debugger import set_trace\n",
    "import datetime\n",
    "from shutil import copyfile\n",
    "from scipy.optimize import minimize_scalar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/kevin/CloudStation/KEVIN/UCI/Research/gravlensing/A611_temp/chains_Abell_611_v30.ip.cNFW_d5\n"
     ]
    }
   ],
   "source": [
    "# The current working directory should CONTAIN the \"chains_...\" subdirectory for the run of interest.\n",
    "print(os.getcwd())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/kevin/CloudStation/KEVIN/UCI/Research/gravlensing/RXCJ2248\n"
     ]
    }
   ],
   "source": [
    "label = 'RXCJ2248.cnfw.v8'  # Don't use this on NFW, there is no 'beta' parameter.\n",
    "path = '/Users/kevin/CloudStation/KEVIN/UCI/Research/gravlensing/RXCJ2248' \n",
    "os.chdir(path)\n",
    "print(os.getcwd())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "# verify that the key files exist\n",
    "assert os.path.exists('chains_' + label)\n",
    "for ext in ['', '.nparam', '.paramnames', '.latex_paramnames', '.ranges']:\n",
    "    assert os.path.exists('chains_' + label + '/' + label + ext)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create the new directory, if it doesn't already exist\n",
    "app = '_d5'\n",
    "newlabel = label + app\n",
    "if not os.path.exists(os.getcwd() + '/chains_' + newlabel):\n",
    "        os.makedirs(os.getcwd() + '/chains_' + newlabel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "# copy key files\n",
    "for ext in ['', '.nparam', '.paramnames', '.latex_paramnames', '.ranges']:\n",
    "    copyfile('chains_' + label + '/' + label + ext, 'chains_' + newlabel + '/' + newlabel + ext) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create 'newpath', for convenience\n",
    "newpath = path + '/chains_' + newlabel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we are deriving scale radius and core radius\n",
    "addl_pnames = ['r_s', 'r_c', 'r_core', 'c_hat_200', 'rm2']\n",
    "addl_latex_pnames = ['r_{s}', 'r_{c}', 'r_{core}', '\\hat{c}_{200}', 'r_{-2}']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "FlatLambdaCDM(H0=70 km / (Mpc s), Om0=0.3, Tcmb0=0 K, Neff=3.04, m_nu=None, Ob0=None)"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from astropy.cosmology import FlatLambdaCDM\n",
    "import astropy.units as u\n",
    "cosmo = FlatLambdaCDM(H0=70 * u.km / u.s /u.Mpc, Om0=0.3)\n",
    "cosmo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "zlens = 0.288"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "182.36871306056253 solMass / kpc3\n",
      "4.329314910896306 kpc / arcsec\n",
      "rho_200 =  2959632.997002065 solMass / arcsec3\n"
     ]
    }
   ],
   "source": [
    "rho_crit = cosmo.critical_density(zlens).to(u.M_sun / u.kpc**3)\n",
    "print(rho_crit)\n",
    "kpc_per_arcsec = cosmo.kpc_proper_per_arcmin(zlens) * u.arcmin / (60. * u.arcsec)\n",
    "print (kpc_per_arcsec)\n",
    "rho_200 = rho_crit * 200 * kpc_per_arcsec**3  # units of M_sun arcsec^-3\n",
    "print('rho_200 = ', rho_200)\n",
    "# get rid of units now\n",
    "rho_200 = rho_200.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in the chain and append the new parameter names\n",
    "samples = np.loadtxt(newpath + '/' + newlabel, comments=\"#\", delimiter=None, unpack=False)\n",
    "with open(newpath + '/' + newlabel + '.paramnames') as afile:\n",
    "    paramnames = afile.readlines()\n",
    "with open(newpath + '/' + newlabel + '.latex_paramnames') as afile:\n",
    "    latex_paramnames = afile.readlines()\n",
    "# append the new parameter names\n",
    "for nm, latex_nm in zip(addl_pnames, addl_latex_pnames):\n",
    "    paramnames.append(nm + '\\n')\n",
    "    latex_paramnames.append(nm + '\\t' + latex_nm + '\\n')\n",
    "    with open(newpath + '/' + newlabel + '.paramnames', 'a') as afile:\n",
    "        afile.write(nm + '\\n')\n",
    "    with open(newpath + '/' + newlabel + '.latex_paramnames', 'a') as afile:\n",
    "        afile.write(nm + '\\t' + latex_nm + '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/kevin/CloudStation/KEVIN/UCI/Research/gravlensing/RXCJ2248/chains_RXCJ2248.cnfw.v8_d5\n"
     ]
    }
   ],
   "source": [
    "# change the working directory to the new path\n",
    "os.chdir(newpath)\n",
    "print(os.getcwd())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0  derived parameters, adding  5\n"
     ]
    }
   ],
   "source": [
    "# increase the number of derived parameters in the .nparam file\n",
    "pnum, dpnum = np.loadtxt(newlabel + '.nparam', comments=\"#\", delimiter=None, unpack=False)\n",
    "addl_np = len(addl_pnames)\n",
    "print(dpnum, ' derived parameters, adding ', addl_np)\n",
    "dpnum += addl_np\n",
    "np.savetxt(newlabel + '.nparam', [np.int(pnum), np.int(dpnum)], fmt='%i')\n",
    "\n",
    "# add new ranges to the .ranges file. Not sure if these should be calculated, just using 0 and 1e30.\n",
    "with open(newlabel + '.ranges', 'a') as afile:\n",
    "    for nm in addl_pnames:\n",
    "        afile.write('0  1e30\\n')\n",
    "\n",
    "# Strip off the \\n character.\n",
    "paramnames = [item[:-1] for item in paramnames] \n",
    "latex_paramnames = [item[:-1] for item in latex_paramnames]  \n",
    "samples = np.array(samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(95405, 14)\n",
      "['mvir', 'c', 'rc_kpc', 'q', 'theta', 'xc', 'yc', 'shear1', 'shear2', 'mtot1', 'mtot2', 'raw_chisq', 'r_s', 'r_c', 'r_core', 'c_hat_200', 'rm2']\n"
     ]
    }
   ],
   "source": [
    "ns, p = samples.shape # note that 'p' here is two more than the number of parameters.\n",
    "# as it counts the first column (weight) and the last (chisquare).\n",
    "print(samples.shape)\n",
    "print(paramnames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check to make sure the names are there. Also finds log(paramname).\n",
    "assert all([any([name.find(item)!=-1 for item in ['mvir', 'c', 'rc_kpc']]) for name in paramnames[0:3]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extend the samples to have addl_np additional columns\n",
    "samples2 = np.append(samples, np.zeros((ns, addl_np)), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "# move the chisquare to be the last column\n",
    "samples2[:,-1] = samples2[:, p-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the new parameters\n",
    "target_cols = [1, 2, 3]\n",
    "targets = [col for col in samples2[:, target_cols].T] \n",
    "# for the cc model with n constant, we need cols [1,2, 3, 4, 5], which are k0, gamma, a1, s and q\n",
    "# for the pmode=3 cNFW model we need cols [1, 2, 3], which are mvir, c, and rc_kpc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert log parameters to non-log. But we are not saving them. Let mkdist do the tranform for plots, as\n",
    "# you need to change the latex names.\n",
    "for i,pm in enumerate(targets):\n",
    "    if paramnames[i][0:4] == 'log(':\n",
    "        targets[i] = 10** pm\n",
    "#         samples[:,i] = 10** samples[:,i]\n",
    "#         paramnames[i-1] = paramnames[i-1][4:len(paramnames[i-1])-1]\n",
    "        \n",
    "# print(paramnames)        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "mvir, c, rc_kpc = targets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([9.22921e+14, 1.33811e+15, 9.66055e+14, ..., 1.68365e+15,\n",
       "        1.76934e+15, 1.63059e+15]),\n",
       " array([12.4721 , 13.5587 , 12.6208 , ...,  7.01562,  6.36054,  7.47692]),\n",
       " array([ 85.6659, 171.578 ,  84.6951, ...,  54.4686,  46.7505,  60.8804]))"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mvir, c, rc_kpc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.00448616"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.amin(rc_kpc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the equation for virial mass and solve for the virial radius:\n",
    "$m_{200} = \\frac{4}{3} \\pi r_{200}^3 \\rho_{200}$, and $r_s = r_{200} / c$,\n",
    "Therefore\n",
    "$$ r_s = (\\frac{3}{4\\pi})^{1/3} (m_{200} / (c  \\rho_{200}))^{1/3}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6203504908994001"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(3 / (np.pi * 4.))**(1./3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.0"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Now need a function to find c_hat. First, the cnfw log slope.\n",
    "def cnfwlogslope(r, rs, rc): \n",
    "    return - (r * (3 * r + 2* rc + rs)) / ((r + rc) * (r + rs))\n",
    "# test:\n",
    "cnfwlogslope(15, 15, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "64.67915297409839\n",
      "8.357816691600547\n"
     ]
    }
   ],
   "source": [
    "# A function to find r_{beta}, the radius at which the log slope is 'beta'.\n",
    "def rbeta_cnfw(r_s, r_c, beta):\n",
    "    if beta == -1:\n",
    "        return r_c * (np.sqrt(1 + 8. * r_s / r_c) - 1) / 4.\n",
    "    elif beta == -2:\n",
    "        return r_s * (np.sqrt(1 + 8. * r_c / r_s) + 1) / 2.\n",
    "    else:\n",
    "        rb = minimize_scalar(lambda r: np.abs(beta - cnfwlogslope(r, r_s, r_c)), method='bounded', bounds=(0,5000))\n",
    "        if rb.success !=True:\n",
    "            print('Optimize Failure')\n",
    "            print(rb)\n",
    "            print(r_s, r_c, beta)\n",
    "        return rb.x\n",
    "# test\n",
    "print(rbeta_cnfw(36.8, 24.5, -2)) # should be 64.679 per Mathematica\n",
    "print(rbeta_cnfw(20., 12., -1))# should be 8.358"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the new parameters\n",
    "rvir = (3 / (np.pi * 4.))**(1./3) * mvir**(1./3) / (rho_200**(1./3))\n",
    "r_s = rvir / c\n",
    "r_c = rc_kpc / kpc_per_arcsec.value\n",
    "beta = r_s / r_c\n",
    "rm2 = [rbeta_cnfw(r_sp, r_cp, -2) for r_sp, r_cp in zip(r_s, r_c)]\n",
    "c_hat = rvir / rm2 \n",
    "r_core = [rbeta_cnfw(r_sp, r_cp, -1) for r_sp, r_cp in zip(r_s, r_c)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples2[:, p-1:p-1+addl_np] = np.array([r_s, r_c, r_core, c_hat, rm2]).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(95405, 19)"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "samples2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['mvir', 'c', 'rc_kpc', 'q', 'theta', 'xc', 'yc', 'shear1', 'shear2', 'mtot1', 'mtot2', 'raw_chisq', 'r_s', 'r_c', 'r_core', 'c_hat_200', 'rm2']\n",
      "17\n"
     ]
    }
   ],
   "source": [
    "print(paramnames)\n",
    "print(len(paramnames))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write the revised samples back to the file\n",
    "np.savetxt(newlabel, samples2, fmt=\"%.6g\", delimiter='   ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best fit chi squared is  70.1092\n"
     ]
    }
   ],
   "source": [
    "for i, name in enumerate(paramnames):\n",
    "    if name == 'raw_chisq':\n",
    "        print('Best fit chi squared is ', np.amin(samples2[:,i+1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'RXCJ2248.cnfw.v8_d5'"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newlabel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "ename": "SystemExit",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "An exception has occurred, use %tb to see the full traceback.\n",
      "\u001b[0;31mSystemExit\u001b[0m\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/kevin/anaconda/envs/py36/lib/python3.6/site-packages/IPython/core/interactiveshell.py:2918: UserWarning: To exit: use 'exit', 'quit', or Ctrl-D.\n",
      "  warn(\"To exit: use 'exit', 'quit', or Ctrl-D.\", stacklevel=1)\n"
     ]
    }
   ],
   "source": [
    "sys.exit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If there is a transform that you want to make, and a dat file for conversion, specify it here\n",
    "transform = True\n",
    "t_file = 'rcx2_transform.dat'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Success\n"
     ]
    }
   ],
   "source": [
    "# make the error file using mkdist\n",
    "if transform:\n",
    "    code = os.system('mkdist ' + newlabel + ' -E -T:../' + t_file + ' | tee ' + newlabel + '.err')\n",
    "else:\n",
    "    code = os.system('mkdist ' + newlabel + ' | tee ' + newlabel + '.err')    \n",
    "if code != 0:\n",
    "    print('Error! Code returned ', code)\n",
    "else:\n",
    "    print('Success')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Success\n"
     ]
    }
   ],
   "source": [
    "# make the various distributions and plots\n",
    "# having trouble with 'ranges much smaller than actual ranges'\n",
    "# so try adding the option '-B5e-4'\n",
    "if transform:\n",
    "    code = os.system('mkdist ' + newlabel + ' -n40 -N40  -B5e-4 -T:../' + t_file + ' |tee tmp.txt')\n",
    "else:\n",
    "    code = os.system('mkdist ' + newlabel + ' -n40 -N40  -B5e-4'+ ' |tee tmp.txt')\n",
    "if code != 0:\n",
    "    print('Error! Code returned ', code)\n",
    "else:\n",
    "    with open('tmp.txt', 'r') as file:\n",
    "        tmp = file.readlines()\n",
    "    if 'WARNING' in tmp or 'Error:' in tmp:\n",
    "        print('Completed with WARNINGS or Error')\n",
    "    elif len(tmp) == 1:\n",
    "        print('Warning: output was only one line:')\n",
    "        print(tmp)\n",
    "    else:\n",
    "        print('Success')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.exit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Success\n"
     ]
    }
   ],
   "source": [
    "# to run the python triangle plot pdf maker\n",
    "code = os.system('python ' + newlabel + '_tri.py')\n",
    "if code != 0:\n",
    "    print('Error! Code returned ', code)\n",
    "else:\n",
    "    print('Success')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:py36]",
   "language": "python",
   "name": "conda-env-py36-py"
  },
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}