qlens kappa distribution calc v5.1 - RXCJ2248.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gravitational Lensing Research\n",
    "# Program to take a qlens-generated sample chain, derive kappa vs. radius values. No plotting.\n",
    "# v5.1: Only DM. Bug fix from v5.\n",
    "# Switching to kpc rather than arc seconds.\n",
    "# Also calculate the mass enclosed versus radius.\n",
    "# Also calculate the radius containing 1% of the M200 mass.\n",
    "# Saves the equal weighted chain.\n",
    "# This version is for RXCJ2248"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "import numpy as np\n",
    "import os\n",
    "import timeit\n",
    "from IPython.core.debugger import set_trace\n",
    "import datetime\n",
    "import concurrent.futures as cf\n",
    "import multiprocessing as mp\n",
    "import copy\n",
    "from scipy import interpolate\n",
    "from scipy.optimize import minimize_scalar\n",
    "import pandas as pd\n",
    "from astropy.cosmology import FlatLambdaCDM\n",
    "import astropy.units as u\n",
    "from astropy import constants as const"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for parallel processing\n",
    "nproc = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Changing the output radii values to kpc. Make sure redshift is correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/kevin/CloudStation/KEVIN/UCI/Research/gravlensing/RXCJ2248/chains_RXCJ2248.cnfw.v10c\n"
     ]
    }
   ],
   "source": [
    "name = 'RXCJ2248'\n",
    "dirname = 'RXCJ2248'\n",
    "suffix = '.cnfw.v10c'\n",
    "label = name + suffix\n",
    "# Set your path here, to the \"chains_...\" directory of interest\n",
    "#path = '/home/grant/Documents/qlens-beta/' + name + '/chains_' + label\n",
    "# path = '/home/kea/KEVIN/UCI/Research/gravlensing/' + dirname + '/chains_' + label\n",
    "path = '/Users/kevin/CloudStation/KEVIN/UCI/Research/gravlensing/' + dirname + '/chains_' + label\n",
    "os.chdir(path)\n",
    "print(os.getcwd())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "zlens =  0.348\n"
     ]
    }
   ],
   "source": [
    "rs_dic = {'A611': 0.288, 'A2537': 0.294, 'RXCJ2248': 0.348, 'MS2137': 0.314, 'A383': 0.189,\n",
    "          'A2261': 0.225, 'M2129': 0.589, 'A1703': 0.280,  'M1720': 0.387}\n",
    "\n",
    "zlens = rs_dic[name]\n",
    "print('zlens = ', zlens)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FlatLambdaCDM(H0=70 km / (Mpc s), Om0=0.3, Tcmb0=0 K, Neff=3.04, m_nu=None, Ob0=None)\n",
      "195.12769323056327 solMass / kpc3\n",
      "4.92139759498602 kpc / arcsec\n",
      "rho_200 =  4651727.483682498 solMass / arcsec3\n"
     ]
    }
   ],
   "source": [
    "\n",
    "cosmo = FlatLambdaCDM(H0=70 * u.km / u.s /u.Mpc, Om0=0.3)\n",
    "print(cosmo)\n",
    "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",
    "kpc_per_arcsec = kpc_per_arcsec.value\n",
    "rho_200 = rho_200.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2226281196.5328746 solMass / kpc2\n"
     ]
    }
   ],
   "source": [
    "# Calculate Sigma_crit, assuming a source redshift of 2.0.\n",
    "sigma_crit = (const.c**2 * cosmo.angular_diameter_distance(2.0) / (4. * np.pi * const.G *\\\n",
    "                           cosmo.angular_diameter_distance(zlens) * cosmo.angular_diameter_distance_z1z2(zlens, 2.0))\\\n",
    "             ).to(u.M_sun / u.kpc**2)\n",
    "print(sigma_crit)   \n",
    "# get rid of units now\n",
    "sigma_crit = sigma_crit.value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Review the starting and stopping values!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Covers radii of   8.56 kpc to    494.22 kpc.\n",
      "Covers radii of   1.74 arcsec to 100.42 arcsec.\n",
      "QLens stopping value is 145.20 arc seconds\n"
     ]
    }
   ],
   "source": [
    "# radii values, in KPC, log spaced. Based on 0.5 * innermost and 2.0 * outermost data points.\n",
    "rstart = 17.12 * 0.5     \n",
    "rstop  = 247.11 * 2.0\n",
    "nr = 12  \n",
    "radii_kpc = np.logspace(np.log10(rstart), np.log10(rstop), nr, endpoint=True)\n",
    "radii_as = radii_kpc / kpc_per_arcsec\n",
    "# QLens uses no endpoint, so it's stopping value needs to be calculated differently\n",
    "ql_stop_as = radii_kpc[-1] * radii_kpc[-1] / radii_kpc[-2] / kpc_per_arcsec\n",
    "radii_kpc = np.round(radii_kpc, 2)\n",
    "print(f'Covers radii of {radii_kpc[0]:6.2f} kpc to    {radii_kpc[-1]:6.2f} kpc.')\n",
    "print(f'Covers radii of {radii_as[0]:6.2f} arcsec to {radii_as[-1]:6.2f} arcsec.')\n",
    "print(f'QLens stopping value is {ql_stop_as:6.2f} arc seconds')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.allclose(radii_as, np.logspace(np.log10(radii_kpc[0] / kpc_per_arcsec), np.log10(ql_stop_as),\n",
    "                                         nr, endpoint=False), atol=0.01, rtol=0.001)\n",
    "radii_as = np.round(radii_as, 2)\n",
    "radii_kpc = np.round(radii_kpc, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def weighted_quantile(values, quantiles, sample_weight=None, \n",
    "                      values_sorted=False, old_style=False):\n",
    "    \"\"\" Very close to numpy.percentile, but supports weights.\n",
    "    NOTE: quantiles should be in [0, 1]!\n",
    "    :param values: numpy.array with data\n",
    "    :param quantiles: array-like with many quantiles needed\n",
    "    :param sample_weight: array-like of the same length as `array`\n",
    "    :param values_sorted: bool, if True, then will avoid sorting of\n",
    "        initial array\n",
    "    :param old_style: if True, will correct output to be consistent\n",
    "        with numpy.percentile.\n",
    "    :return: numpy.array with computed quantiles.\n",
    "    \"\"\"\n",
    "    values = np.array(values)\n",
    "    quantiles = np.array(quantiles)\n",
    "    if sample_weight is None:\n",
    "        sample_weight = np.ones(len(values))\n",
    "    sample_weight = np.array(sample_weight)\n",
    "    assert np.all(quantiles >= 0) and np.all(quantiles <= 1), \\\n",
    "        'quantiles should be in [0, 1]'\n",
    "\n",
    "    if not values_sorted:\n",
    "        sorter = np.argsort(values)\n",
    "        values = values[sorter]\n",
    "        sample_weight = sample_weight[sorter]\n",
    "\n",
    "    weighted_quantiles = np.cumsum(sample_weight) - 0.5 * sample_weight\n",
    "    if old_style:\n",
    "        # To be convenient with numpy.percentile\n",
    "        weighted_quantiles -= weighted_quantiles[0]\n",
    "        weighted_quantiles /= weighted_quantiles[-1]\n",
    "    else:\n",
    "        weighted_quantiles /= np.sum(sample_weight)\n",
    "    return np.interp(quantiles, weighted_quantiles, values)\n",
    "\n",
    "# now calculated the median and +/- 1-sigma values, accounting for sample weights.\n",
    "def find_bands(p, samples, n_rad, weighted=True):\n",
    "    # 'p' is the offset of the first column to sort by.\n",
    "    # 'n_rad' is the number of columns to sort.\n",
    "    # The first column is the weight, unless weighted=False.\n",
    "    if not weighted:\n",
    "        p += 1\n",
    "        samples = np.hstack((np.ones((samples.shape[0], 1)), samples))\n",
    "    total = [np.sum(samples[:, i + p] * samples[:, 0]) for i in range(n_rad)]\n",
    "    weights = samples[:, 0]\n",
    "    sigma_minus  = np.zeros(n_rad)\n",
    "    sigma_plus = np.zeros(n_rad)\n",
    "    median = np.zeros(n_rad)\n",
    "    for i in range(n_rad):\n",
    "        sigma_minus[i], median[i], sigma_plus[i] = weighted_quantile(samples[:, i+p], [.16, .50, .84], weights)\n",
    "    return sigma_minus, sigma_plus, median"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['mvir', 'c', 'rc_kpc', 'q', 'theta', 'xc', 'yc', 'shear1', 'shear2', 'mtot1', 'mtot2']\n"
     ]
    }
   ],
   "source": [
    "# Open the paramnames file and display\n",
    "with open(path + '/' + label + '.paramnames') as afile:\n",
    "    paramnames = afile.readlines()\n",
    "paramnames = [item[:-1] for item in paramnames]  # Strip off the \\n character.\n",
    "print(paramnames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ql_script5(i, p, params, tstart, ns):\n",
    "    # assembles and runs the QLens script, returning ONE vector of kappa values (dmo)\n",
    "    # Note that the first and pth columns are not parameters.\n",
    "    # the pth column is the chi square.\n",
    "    # Note that \"like\" is no longer on this list; equal weighted samples. \n",
    "    # Note that since we are now just looking at dark matter, this greatly simplifies the\n",
    "    # script, eliminating parts b and c.\n",
    "    # This could actually be done analytically.\n",
    "    \n",
    "    mvir, c, rc_kpc, q, theta, xc, yc, shear1, shear2, mtot1, mtot2, loglike = params[0:p]\n",
    "    # NOTE: Check the above line; it may be subtly different for each halo, depending on the number of\n",
    "    # perturbers, derived parameters, etc.\n",
    "    \n",
    "    pars0 = ' '.join([str(item) for item in [mvir, c, rc_kpc, q, theta, xc, yc]]) # Halo parameters\n",
    "    \n",
    "    pars1 = ' '.join([str(item) for item in [shear1, shear2]]) \n",
    "\n",
    "    # Note that since we are now just doing dark matter of the main halo, we don't need to enter \n",
    "    # the code for the BCG or perturbers!\n",
    "\n",
    "    kappaname = 'plotdata/kappadist_temp_' + str(mp.current_process().name) + '.tmp'\n",
    "    kappaname_dmo = 'plotdata/kappadist_temp_' + str(mp.current_process().name) + '.dmo.tmp'\n",
    "\n",
    "    script_part_a1 = '''\n",
    "\n",
    "    zlens=0.348    # REVIEW\n",
    "    zsrc=2.0\n",
    "    zsrc_ref=2.0\n",
    "    shear_components on\n",
    "    major_axis_along_y off    \n",
    "\n",
    "    sci_notation off\n",
    "    gridtype cartesian  # REVIEW\n",
    "    grid -70 65 -45 55  # REVIEW\n",
    "    imgdata read ../RXCJ2248_SL_Data_v3.dat   #REVIEW\n",
    "    '''\n",
    "\n",
    "    script_part_a2 = '''\n",
    "    lens clear\n",
    "    lens cnfw pmode=3 '''\n",
    "    # here we insert the parameters for the first lens (halo)\n",
    "\n",
    "    script_part_b1 =''' shear='''\n",
    "    # here we insert the paramters for the shear\n",
    "    \n",
    "    script_part_b2 ='''\n",
    "    lens pjaffe '''\n",
    "    # here we insert the parameters for the second lens (BCG)\n",
    "\n",
    "    script_part_b3 ='''\n",
    "    lens pjaffe '''\n",
    "    # pert 1\n",
    "\n",
    "    script_part_b4 ='''\n",
    "    lens pjaffe  '''\n",
    "    # pert 2\n",
    "    \n",
    "    script_part_b5 ='''\n",
    "    lens pjaffe  '''\n",
    "    # pert 3\n",
    "    \n",
    "#     script_part_c = '''\n",
    "# lens perturber pjaffe\t  0.275/anchor=5,0\t0.275/anchor=5,1\t0.0096/anchor=5,2    0.67\t     80.6\t-10.88\t10.22\t# 4\n",
    "\n",
    "# lens perturber pjaffe\t  0.209/anchor=5,0    \t0.209/anchor=5,1\t0.0073/anchor=5,2    0.84\t     128.4\t -16.79\t 0.60\t# 5\n",
    "\n",
    "# lens perturber pjaffe\t  0.158/anchor=5,0    \t0.158/anchor=5,1\t0.0055/anchor=5,2\t 0.9\t\t 131.6\t 1.13\t -2.78\t# 6\n",
    "\n",
    "# lens perturber pjaffe\t  0.240/anchor=5,0    \t0.240/anchor=5,1\t0.0084/anchor=5,2\t 0.79\t      61.7\t -13.68\t 12.87\t# 7\n",
    "\n",
    "#     '''\n",
    "    # Other perturbers, if any, can go in part c above.\n",
    "    \n",
    "    script_part_d = '''\n",
    "    plotkappa ''' + str(rstart / kpc_per_arcsec) + ' ' + str(ql_stop_as)\\\n",
    "    + ' ' + str(nr) + ' '\n",
    "\n",
    "    script_part_e = '''\n",
    "    cosmology\n",
    "    '''\n",
    "    \n",
    "    script = script_part_a1 + script_part_a2 + pars0 \\\n",
    "     + script_part_d + kappaname_dmo + ' lens=0\\n' \\\n",
    "     + script_part_e\n",
    "    \n",
    "    #         + script_part_b1 + pars1\\\n",
    "#         + script_part_b2 + pars2 + script_part_b3 + pars3  \\\n",
    "#         + script_part_b4 + pars4 + script_part_b5 + pars5  \\\n",
    "#         + script_part_c \\\n",
    "#         + script_part_d + kappaname + ' \\n' \\\n",
    "    \n",
    "    sc_name = 'plotdata/kappascript_' + str(mp.current_process().name) + '.in'\n",
    "    \n",
    "    # double check redshift\n",
    "    rs_pos = script.find('zlens=')\n",
    "    if rs_pos == -1:\n",
    "        print(\"Can't locate 'zlens=' in script.\")\n",
    "        sys.exit()    \n",
    "    eol = script.find('# REVIEW', rs_pos)\n",
    "    assert float(script[rs_pos + 6: eol]) == zlens\n",
    "    \n",
    "    with open(sc_name, 'w') as outfile:\n",
    "        outfile.write(script)\n",
    "\n",
    "    # Now, run the script, saving the output in a temporary file\n",
    "    os.system('qlens ' + sc_name + ' -q | tee plotdata/scriptout' + str(mp.current_process().name) + '.tmp')\n",
    "    \n",
    "    # read back the cosmology and determine the r_200 of lens 0\n",
    "    with open('plotdata/scriptout' + str(mp.current_process().name) + '.tmp', 'r') as outf:\n",
    "        output = ''.join(outf.readlines())\n",
    "    cos_pos = output.find('cosmology')\n",
    "    lens0_pos = output.find('Lens 0:', cos_pos)\n",
    "    r200_pos = output.find('r_200 = ', lens0_pos)\n",
    "    kpc_pos = output.find('kpc', r200_pos)\n",
    "    r200 = np.float(output[r200_pos + 7 : kpc_pos])\n",
    "    m200 = 4. * np.pi / 3. * r200**3 * rho_200 / kpc_per_arcsec**3 \n",
    "    \n",
    "    # read in kappa file back in, and store results\n",
    "#     kdat = np.loadtxt(kappaname)\n",
    "#     kappa = kdat[:, 1]\n",
    "    kdat_dmo = np.loadtxt(kappaname_dmo)\n",
    "    kappa_dmo = kdat_dmo[:, 1]\n",
    "    encl_mass = kdat_dmo[:, 4]\n",
    "    # The above assumes that the QLens data is in M_sun.\n",
    "#     kappa_bary = kappa - kappa_dmo\n",
    "    \n",
    "    # Calculate the radius containing 1% of m200:\n",
    "    f_encl = interpolate.interp1d(radii_kpc, encl_mass, fill_value='extrapolate')\n",
    "    r_1_res = minimize_scalar(lambda r: np.abs(f_encl(r) - 0.01 * m200))\n",
    "    if r_1_res.success != True:\n",
    "        set_trace()\n",
    "    r_1pct = r_1_res.x\n",
    "    if r_1pct <0. or r_1pct > 4000.:\n",
    "        set_trace()\n",
    "    \n",
    "    if i==0:\n",
    "        assert np.all(np.abs((kdat_dmo[:, 0] * kpc_per_arcsec - radii_kpc) / radii_kpc) \\\n",
    "                      < .01)\n",
    "    if i%50 == 1:\n",
    "        tavg = (timeit.default_timer() - tstart) / (i+1.)\n",
    "        proj_fin =  tavg * (ns -i -1)\n",
    "        sys.stdout.write(\"\\r{0:5.2%} complete. Avg time per iteration: {1:8.6f} sec. Projected finish in {2:4.2f} min. \"\\\n",
    "                   .format((float(i) / ns), tavg,  proj_fin/60))\n",
    "    sys.stdout.flush() # required in multiprocessing or else it buffers this output\n",
    "    return kappa_dmo, encl_mass, r_1pct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now doing this calculation not in a function.\n",
    "# Import the chain. This is the EQUAL weighted chain.\n",
    "samples = np.loadtxt(path + '/' + label + 'post_equal_weights.dat', comments=\"#\", delimiter=None, unpack=False)\n",
    "samples = np.array(samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the ranges, and scale the parameters.\n",
    "ranges = np.loadtxt(path + '/' + label + '.ranges', comments=\"#\", delimiter=None, unpack=False)\n",
    "for i, (lwr, upr) in enumerate(ranges):  # All columns but the last one in samples.\n",
    "    if upr != 1.e30:\n",
    "        samples[:, i] = samples[:, i] * (upr - lwr) + lwr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(40878, 12)\n",
      "Median parameter values: \n",
      "['mvir: 1.839e+15', 'c: 5.39', 'rc_kpc: 19.82', 'q: 0.5044', 'theta: -35.9', 'xc: -0.5391', 'yc: 0.3222', 'shear1: 0.005871', 'shear2: 0.02715', 'mtot1: 8.159e+11', 'mtot2: 7.412e+09']\n"
     ]
    }
   ],
   "source": [
    "ns, p = samples.shape\n",
    "print(samples.shape)\n",
    "print('Median parameter values: ')\n",
    "print([f'{pn}: {mv:.4g}' for (pn, mv) in zip(paramnames, [np.median(i) for i in samples.T])])\n",
    "\n",
    "chain_columns = paramnames + ['loglike']\n",
    "\n",
    "# for testing only, you can uncomment the two lines below, and it will only use the last 500 samples,\n",
    "# which is useful to quickly see if things are working right.\n",
    "\n",
    "# ns = 500\n",
    "# samples = samples[-ns:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                                                                      \n",
      "Kappa calculation completed in  20.08  minutes.\n"
     ]
    }
   ],
   "source": [
    "# make sure there is a directory called plotdata. If not, create it.\n",
    "if not os.path.exists(os.getcwd()+'/plotdata'):\n",
    "    os.makedirs(os.getcwd()+'/plotdata')\n",
    "\n",
    "results = []\n",
    "tstart = timeit.default_timer()\n",
    "with cf.ProcessPoolExecutor(nproc) as pp:\n",
    "    for i, params in enumerate(samples[:, 0:p]):\n",
    "        results.append((i, pp.submit(ql_script5, i, p, params, tstart, ns)))\n",
    "\n",
    "# code waits here until all processes are finished\n",
    "sys.stdout.write('\\r                                                                                      ')\n",
    "print(\"\\nKappa calculation completed in \", np.round((timeit.default_timer() - tstart) / 60, 2), \" minutes.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Updating data array...\n",
      "Data updated.\n"
     ]
    }
   ],
   "source": [
    "all_names = []  #['all_%07.2f' % rad for rad in radii_kpc]\n",
    "dm_names = ['dm_%07.2f' % rad for rad in radii_kpc]\n",
    "bary_names = []  #['bary_%07.2f' % rad for rad in radii_kpc]\n",
    "m_names = ['m_encl_%07.2f' % rad for rad in radii_kpc]\n",
    "col_names = all_names + dm_names + bary_names + m_names + ['r_1pct']\n",
    "\n",
    "# Use a numpy array for holding everything. DataFrames are too slow.\n",
    "samples2 = np.append(samples, np.full((ns, nr * 2 + 1), np.nan), axis=1)  # Note: changed *4 to *3 here\n",
    "\n",
    "old_i = -1\n",
    "print('Updating data array...')\n",
    "for (i,k) in results:\n",
    "    k_dm, m_encl, r_1= k.result()\n",
    "    reslist = list(k_dm) + list(m_encl)+ [r_1]\n",
    "    samples2[i, p:] = reslist\n",
    "    if old_i +1 != i:\n",
    "        print('WARNING: Out of order at i=', i)\n",
    "    old_i = i\n",
    "print('Data updated.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Making bands...\n",
      "Bands made.\n"
     ]
    }
   ],
   "source": [
    "print('Making bands...')\n",
    "# sm, sp, med = find_bands(0, samples2[:, p: p + nr], nr, weighted=False)\n",
    "sm_dmo, sp_dmo, med_dmo = find_bands(0, samples2[:, p + nr * 0: p + nr * 1], nr, weighted=False)\n",
    "# sm_bary, sp_bary, med_bary = find_bands(0, samples2[:, p + nr * 2: p + nr *3], nr, weighted=False)\n",
    "sm_mencl, sp_mencl, med_mencl = find_bands(0, samples2[:, p + nr * 1: p + nr * 2], nr, weighted=False)\n",
    "sm_r1, sp_r1, med_r1 = find_bands(0, samples2[:, p + nr * 2: p + nr * 2 + 1], 1, weighted=False)\n",
    "print(\"Bands made.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mvir</th>\n",
       "      <th>c</th>\n",
       "      <th>rc_kpc</th>\n",
       "      <th>q</th>\n",
       "      <th>theta</th>\n",
       "      <th>xc</th>\n",
       "      <th>yc</th>\n",
       "      <th>shear1</th>\n",
       "      <th>shear2</th>\n",
       "      <th>mtot1</th>\n",
       "      <th>...</th>\n",
       "      <th>m_encl_0025.87</th>\n",
       "      <th>m_encl_0037.41</th>\n",
       "      <th>m_encl_0054.09</th>\n",
       "      <th>m_encl_0078.21</th>\n",
       "      <th>m_encl_0113.08</th>\n",
       "      <th>m_encl_0163.50</th>\n",
       "      <th>m_encl_0236.41</th>\n",
       "      <th>m_encl_0341.81</th>\n",
       "      <th>m_encl_0494.22</th>\n",
       "      <th>r_1pct</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>4.087800e+04</td>\n",
       "      <td>40878.000000</td>\n",
       "      <td>40878.000000</td>\n",
       "      <td>40878.000000</td>\n",
       "      <td>40878.000000</td>\n",
       "      <td>40878.000000</td>\n",
       "      <td>40878.000000</td>\n",
       "      <td>40878.000000</td>\n",
       "      <td>40878.000000</td>\n",
       "      <td>4.087800e+04</td>\n",
       "      <td>...</td>\n",
       "      <td>4.087800e+04</td>\n",
       "      <td>4.087800e+04</td>\n",
       "      <td>4.087800e+04</td>\n",
       "      <td>4.087800e+04</td>\n",
       "      <td>4.087800e+04</td>\n",
       "      <td>4.087800e+04</td>\n",
       "      <td>4.087800e+04</td>\n",
       "      <td>4.087800e+04</td>\n",
       "      <td>4.087800e+04</td>\n",
       "      <td>40878.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>1.848519e+15</td>\n",
       "      <td>5.449458</td>\n",
       "      <td>21.537840</td>\n",
       "      <td>0.504281</td>\n",
       "      <td>-35.889833</td>\n",
       "      <td>-0.545280</td>\n",
       "      <td>0.324391</td>\n",
       "      <td>0.005858</td>\n",
       "      <td>0.027077</td>\n",
       "      <td>8.950429e+11</td>\n",
       "      <td>...</td>\n",
       "      <td>7.416650e+12</td>\n",
       "      <td>1.449281e+13</td>\n",
       "      <td>2.768320e+13</td>\n",
       "      <td>5.145163e+13</td>\n",
       "      <td>9.260865e+13</td>\n",
       "      <td>1.606802e+14</td>\n",
       "      <td>2.676117e+14</td>\n",
       "      <td>4.264283e+14</td>\n",
       "      <td>6.488860e+14</td>\n",
       "      <td>42.464579</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>1.408857e+14</td>\n",
       "      <td>0.654340</td>\n",
       "      <td>11.393482</td>\n",
       "      <td>0.013269</td>\n",
       "      <td>0.410105</td>\n",
       "      <td>0.208035</td>\n",
       "      <td>0.164816</td>\n",
       "      <td>0.004859</td>\n",
       "      <td>0.012385</td>\n",
       "      <td>7.122293e+11</td>\n",
       "      <td>...</td>\n",
       "      <td>3.775497e+11</td>\n",
       "      <td>5.520848e+11</td>\n",
       "      <td>7.281558e+11</td>\n",
       "      <td>8.549323e+11</td>\n",
       "      <td>9.585740e+11</td>\n",
       "      <td>1.479274e+12</td>\n",
       "      <td>3.382544e+12</td>\n",
       "      <td>8.020427e+12</td>\n",
       "      <td>1.755874e+13</td>\n",
       "      <td>1.592729</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>1.418113e+15</td>\n",
       "      <td>3.707292</td>\n",
       "      <td>0.001222</td>\n",
       "      <td>0.447728</td>\n",
       "      <td>-37.510151</td>\n",
       "      <td>-1.444772</td>\n",
       "      <td>-0.329024</td>\n",
       "      <td>-0.019284</td>\n",
       "      <td>-0.022612</td>\n",
       "      <td>1.000578e+10</td>\n",
       "      <td>...</td>\n",
       "      <td>6.385700e+12</td>\n",
       "      <td>1.283326e+13</td>\n",
       "      <td>2.525078e+13</td>\n",
       "      <td>4.828014e+13</td>\n",
       "      <td>8.899978e+13</td>\n",
       "      <td>1.547543e+14</td>\n",
       "      <td>2.539628e+14</td>\n",
       "      <td>3.952640e+14</td>\n",
       "      <td>5.837414e+14</td>\n",
       "      <td>37.067278</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>1.748689e+15</td>\n",
       "      <td>4.976196</td>\n",
       "      <td>13.104794</td>\n",
       "      <td>0.495481</td>\n",
       "      <td>-36.166841</td>\n",
       "      <td>-0.683223</td>\n",
       "      <td>0.213828</td>\n",
       "      <td>0.002598</td>\n",
       "      <td>0.018746</td>\n",
       "      <td>2.023390e+11</td>\n",
       "      <td>...</td>\n",
       "      <td>7.130122e+12</td>\n",
       "      <td>1.407233e+13</td>\n",
       "      <td>2.713318e+13</td>\n",
       "      <td>5.081082e+13</td>\n",
       "      <td>9.191570e+13</td>\n",
       "      <td>1.596806e+14</td>\n",
       "      <td>2.653463e+14</td>\n",
       "      <td>4.210450e+14</td>\n",
       "      <td>6.369294e+14</td>\n",
       "      <td>41.364692</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>1.839097e+15</td>\n",
       "      <td>5.390324</td>\n",
       "      <td>19.824124</td>\n",
       "      <td>0.504414</td>\n",
       "      <td>-35.899245</td>\n",
       "      <td>-0.539149</td>\n",
       "      <td>0.322193</td>\n",
       "      <td>0.005871</td>\n",
       "      <td>0.027146</td>\n",
       "      <td>8.158775e+11</td>\n",
       "      <td>...</td>\n",
       "      <td>7.393340e+12</td>\n",
       "      <td>1.447783e+13</td>\n",
       "      <td>2.768957e+13</td>\n",
       "      <td>5.148628e+13</td>\n",
       "      <td>9.264584e+13</td>\n",
       "      <td>1.607010e+14</td>\n",
       "      <td>2.676292e+14</td>\n",
       "      <td>4.264279e+14</td>\n",
       "      <td>6.486312e+14</td>\n",
       "      <td>42.390886</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>1.937086e+15</td>\n",
       "      <td>5.865426</td>\n",
       "      <td>28.383153</td>\n",
       "      <td>0.513305</td>\n",
       "      <td>-35.615958</td>\n",
       "      <td>-0.402344</td>\n",
       "      <td>0.434084</td>\n",
       "      <td>0.009110</td>\n",
       "      <td>0.035447</td>\n",
       "      <td>1.453420e+12</td>\n",
       "      <td>...</td>\n",
       "      <td>7.673426e+12</td>\n",
       "      <td>1.488529e+13</td>\n",
       "      <td>2.822511e+13</td>\n",
       "      <td>5.211160e+13</td>\n",
       "      <td>9.333716e+13</td>\n",
       "      <td>1.617046e+14</td>\n",
       "      <td>2.699096e+14</td>\n",
       "      <td>4.317851e+14</td>\n",
       "      <td>6.604379e+14</td>\n",
       "      <td>43.486476</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>2.464483e+15</td>\n",
       "      <td>9.418072</td>\n",
       "      <td>91.256955</td>\n",
       "      <td>0.559011</td>\n",
       "      <td>-34.032077</td>\n",
       "      <td>0.216247</td>\n",
       "      <td>1.002830</td>\n",
       "      <td>0.030492</td>\n",
       "      <td>0.077162</td>\n",
       "      <td>3.386619e+12</td>\n",
       "      <td>...</td>\n",
       "      <td>8.913909e+12</td>\n",
       "      <td>1.651625e+13</td>\n",
       "      <td>3.013513e+13</td>\n",
       "      <td>5.417483e+13</td>\n",
       "      <td>9.554669e+13</td>\n",
       "      <td>1.658738e+14</td>\n",
       "      <td>2.811612e+14</td>\n",
       "      <td>4.605071e+14</td>\n",
       "      <td>7.231346e+14</td>\n",
       "      <td>50.183709</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>8 rows × 37 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "               mvir             c        rc_kpc             q         theta  \\\n",
       "count  4.087800e+04  40878.000000  40878.000000  40878.000000  40878.000000   \n",
       "mean   1.848519e+15      5.449458     21.537840      0.504281    -35.889833   \n",
       "std    1.408857e+14      0.654340     11.393482      0.013269      0.410105   \n",
       "min    1.418113e+15      3.707292      0.001222      0.447728    -37.510151   \n",
       "25%    1.748689e+15      4.976196     13.104794      0.495481    -36.166841   \n",
       "50%    1.839097e+15      5.390324     19.824124      0.504414    -35.899245   \n",
       "75%    1.937086e+15      5.865426     28.383153      0.513305    -35.615958   \n",
       "max    2.464483e+15      9.418072     91.256955      0.559011    -34.032077   \n",
       "\n",
       "                 xc            yc        shear1        shear2         mtot1  \\\n",
       "count  40878.000000  40878.000000  40878.000000  40878.000000  4.087800e+04   \n",
       "mean      -0.545280      0.324391      0.005858      0.027077  8.950429e+11   \n",
       "std        0.208035      0.164816      0.004859      0.012385  7.122293e+11   \n",
       "min       -1.444772     -0.329024     -0.019284     -0.022612  1.000578e+10   \n",
       "25%       -0.683223      0.213828      0.002598      0.018746  2.023390e+11   \n",
       "50%       -0.539149      0.322193      0.005871      0.027146  8.158775e+11   \n",
       "75%       -0.402344      0.434084      0.009110      0.035447  1.453420e+12   \n",
       "max        0.216247      1.002830      0.030492      0.077162  3.386619e+12   \n",
       "\n",
       "       ...  m_encl_0025.87  m_encl_0037.41  m_encl_0054.09  m_encl_0078.21  \\\n",
       "count  ...    4.087800e+04    4.087800e+04    4.087800e+04    4.087800e+04   \n",
       "mean   ...    7.416650e+12    1.449281e+13    2.768320e+13    5.145163e+13   \n",
       "std    ...    3.775497e+11    5.520848e+11    7.281558e+11    8.549323e+11   \n",
       "min    ...    6.385700e+12    1.283326e+13    2.525078e+13    4.828014e+13   \n",
       "25%    ...    7.130122e+12    1.407233e+13    2.713318e+13    5.081082e+13   \n",
       "50%    ...    7.393340e+12    1.447783e+13    2.768957e+13    5.148628e+13   \n",
       "75%    ...    7.673426e+12    1.488529e+13    2.822511e+13    5.211160e+13   \n",
       "max    ...    8.913909e+12    1.651625e+13    3.013513e+13    5.417483e+13   \n",
       "\n",
       "       m_encl_0113.08  m_encl_0163.50  m_encl_0236.41  m_encl_0341.81  \\\n",
       "count    4.087800e+04    4.087800e+04    4.087800e+04    4.087800e+04   \n",
       "mean     9.260865e+13    1.606802e+14    2.676117e+14    4.264283e+14   \n",
       "std      9.585740e+11    1.479274e+12    3.382544e+12    8.020427e+12   \n",
       "min      8.899978e+13    1.547543e+14    2.539628e+14    3.952640e+14   \n",
       "25%      9.191570e+13    1.596806e+14    2.653463e+14    4.210450e+14   \n",
       "50%      9.264584e+13    1.607010e+14    2.676292e+14    4.264279e+14   \n",
       "75%      9.333716e+13    1.617046e+14    2.699096e+14    4.317851e+14   \n",
       "max      9.554669e+13    1.658738e+14    2.811612e+14    4.605071e+14   \n",
       "\n",
       "       m_encl_0494.22        r_1pct  \n",
       "count    4.087800e+04  40878.000000  \n",
       "mean     6.488860e+14     42.464579  \n",
       "std      1.755874e+13      1.592729  \n",
       "min      5.837414e+14     37.067278  \n",
       "25%      6.369294e+14     41.364692  \n",
       "50%      6.486312e+14     42.390886  \n",
       "75%      6.604379e+14     43.486476  \n",
       "max      7.231346e+14     50.183709  \n",
       "\n",
       "[8 rows x 37 columns]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Multiply the kappas in the samples2 array by Sigma_crit, to yield surface density in M_sun kpc^-2\n",
    "samples3 = samples2\n",
    "samples3[:, p: p + nr] = samples2[:, p: p + nr] * sigma_crit  # Note, this is altered for dmo, otherwise need nr*3\n",
    "df = pd.DataFrame(samples3, columns = chain_columns + col_names)\n",
    "df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.7015575 , 1.645916  , 1.566243  , 1.4590855 , 1.3238275 ,\n",
       "       1.1630325 , 0.9848686 , 0.7997282 , 0.6184622 , 0.45300535,\n",
       "       0.312952  , 0.2037076 ])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "med_dmo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The values above should start at approx. 1.0 to 2.0, then decrease gradually to approx. 0.1 to 0.3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Results saved to bands file, and dataframe saved to .csv file with label  RXCJ2248.cnfw.v10c\n",
      "Path:  /Users/kevin/CloudStation/KEVIN/UCI/Research/gravlensing/RXCJ2248/chains_RXCJ2248.cnfw.v10c\n"
     ]
    }
   ],
   "source": [
    "# np.savez(label + '.all.bands5', radii_kpc, sm, sp, med)\n",
    "np.savez(label + '.dmo.bands5', radii_kpc, sm_dmo, sp_dmo, med_dmo)\n",
    "# np.savez(label + '.bary.bands5', radii_kpc, sm_bary, sp_bary, med_bary)\n",
    "np.savez(label + '.mencl.bands5', radii_kpc, sm_mencl, sp_mencl, med_mencl)\n",
    "np.savez(label + '.r1pct.bands5', radii_kpc, sm_r1, sp_r1, med_r1)\n",
    "with open(label + '.column_names', 'w') as colfile:\n",
    "    colfile.writelines('\\n'.join(list(df.columns.values)))\n",
    "df.to_csv(label + '.csv')\n",
    "print('Results saved to bands file, and dataframe saved to .csv file with label ', label)\n",
    "print('Path: ', path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use Notebook \"Kappa band plotter v5\" to plot the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['mvir', 'c', 'rc_kpc', 'q', 'theta', 'xc', 'yc', 'shear1', 'shear2',\n",
       "       'mtot1', 'mtot2', 'loglike', 'dm_0008.56', 'dm_0012.38', 'dm_0017.90',\n",
       "       'dm_0025.87', 'dm_0037.41', 'dm_0054.09', 'dm_0078.21', 'dm_0113.08',\n",
       "       'dm_0163.50', 'dm_0236.41', 'dm_0341.81', 'dm_0494.22',\n",
       "       'm_encl_0008.56', 'm_encl_0012.38', 'm_encl_0017.90', 'm_encl_0025.87',\n",
       "       'm_encl_0037.41', 'm_encl_0054.09', 'm_encl_0078.21', 'm_encl_0113.08',\n",
       "       'm_encl_0163.50', 'm_encl_0236.41', 'm_encl_0341.81', 'm_encl_0494.22',\n",
       "       'r_1pct'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/kevin/opt/anaconda3/envs/py37/lib/python3.7/site-packages/pandas/core/frame.py:4238: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  return super().rename(**kwargs)\n"
     ]
    }
   ],
   "source": [
    "# Make a dataframe for Manoj to use in SIDM calculations.\n",
    "dfsidm = df[['mvir', 'c', 'rc_kpc', 'mtot1', 'loglike',  *dm_names]]\n",
    "dfsidm.rename({'mtot1': 'mBCG'}, axis='columns', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original dataframe is 40878 rows. Attempting to thin to 5109\n"
     ]
    }
   ],
   "source": [
    "means = dfsidm.mean(axis=0, numeric_only=True)\n",
    "stdev = dfsidm.std(axis=0, numeric_only=True)\n",
    "r, c = dfsidm.shape\n",
    "r_thin_target = 5000.\n",
    "stride = max(1, int(r / r_thin_target))\n",
    "r_thin = int(r / stride)\n",
    "print(f'Original dataframe is {r} rows. Attempting to thin to {r_thin}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5110, 17)\n",
      "RMS mean change is 0.3%\n",
      "RMS stdev change is 1.0%\n"
     ]
    }
   ],
   "source": [
    "dfthin = dfsidm[::stride]\n",
    "print(dfthin.shape)\n",
    "thin_means = dfthin.mean(axis=0, numeric_only=True)\n",
    "thin_stdev = dfthin.std(axis=0, numeric_only=True)\n",
    "mean_chg = np.sqrt(np.mean(((means - thin_means) / means)**2))\n",
    "sd_chg = np.sqrt(np.mean(((stdev - thin_stdev) / stdev)**2))\n",
    "print(f'RMS mean change is {mean_chg:.1%}')\n",
    "print(f'RMS stdev change is {sd_chg:.1%}')\n",
    "assert mean_chg < 0.01\n",
    "assert sd_chg < 0.02"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saved thinned SIDM dataframe to RXCJ2248.cnfw.v10c.SIDM.csv\n"
     ]
    }
   ],
   "source": [
    "# save if error tolerance met\n",
    "dfthin.to_csv(label + '.SIDM.csv')\n",
    "print(f'Saved thinned SIDM dataframe to {label}.SIDM.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mvir</th>\n",
       "      <th>c</th>\n",
       "      <th>rc_kpc</th>\n",
       "      <th>mBCG</th>\n",
       "      <th>loglike</th>\n",
       "      <th>dm_0008.56</th>\n",
       "      <th>dm_0012.38</th>\n",
       "      <th>dm_0017.90</th>\n",
       "      <th>dm_0025.87</th>\n",
       "      <th>dm_0037.41</th>\n",
       "      <th>dm_0054.09</th>\n",
       "      <th>dm_0078.21</th>\n",
       "      <th>dm_0113.08</th>\n",
       "      <th>dm_0163.50</th>\n",
       "      <th>dm_0236.41</th>\n",
       "      <th>dm_0341.81</th>\n",
       "      <th>dm_0494.22</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>5.110000e+03</td>\n",
       "      <td>5110.000000</td>\n",
       "      <td>5110.000000</td>\n",
       "      <td>5.110000e+03</td>\n",
       "      <td>5110.000000</td>\n",
       "      <td>5.110000e+03</td>\n",
       "      <td>5.110000e+03</td>\n",
       "      <td>5.110000e+03</td>\n",
       "      <td>5.110000e+03</td>\n",
       "      <td>5.110000e+03</td>\n",
       "      <td>5.110000e+03</td>\n",
       "      <td>5.110000e+03</td>\n",
       "      <td>5.110000e+03</td>\n",
       "      <td>5.110000e+03</td>\n",
       "      <td>5.110000e+03</td>\n",
       "      <td>5.110000e+03</td>\n",
       "      <td>5.110000e+03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>1.849007e+15</td>\n",
       "      <td>5.445033</td>\n",
       "      <td>21.496627</td>\n",
       "      <td>9.072186e+11</td>\n",
       "      <td>-42.760813</td>\n",
       "      <td>3.817064e+09</td>\n",
       "      <td>3.679471e+09</td>\n",
       "      <td>3.491438e+09</td>\n",
       "      <td>3.246454e+09</td>\n",
       "      <td>2.943129e+09</td>\n",
       "      <td>2.587198e+09</td>\n",
       "      <td>2.192480e+09</td>\n",
       "      <td>1.780261e+09</td>\n",
       "      <td>1.376748e+09</td>\n",
       "      <td>1.008605e+09</td>\n",
       "      <td>6.973203e+08</td>\n",
       "      <td>4.543723e+08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>1.385936e+14</td>\n",
       "      <td>0.642713</td>\n",
       "      <td>11.234691</td>\n",
       "      <td>7.100991e+11</td>\n",
       "      <td>2.501500</td>\n",
       "      <td>2.741335e+08</td>\n",
       "      <td>2.175188e+08</td>\n",
       "      <td>1.589643e+08</td>\n",
       "      <td>1.039179e+08</td>\n",
       "      <td>5.842773e+07</td>\n",
       "      <td>2.812180e+07</td>\n",
       "      <td>1.750103e+07</td>\n",
       "      <td>1.887757e+07</td>\n",
       "      <td>2.159716e+07</td>\n",
       "      <td>2.421675e+07</td>\n",
       "      <td>2.503857e+07</td>\n",
       "      <td>2.284468e+07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>1.468804e+15</td>\n",
       "      <td>3.796233</td>\n",
       "      <td>0.002746</td>\n",
       "      <td>1.001039e+10</td>\n",
       "      <td>-56.173201</td>\n",
       "      <td>3.217709e+09</td>\n",
       "      <td>3.169329e+09</td>\n",
       "      <td>3.089798e+09</td>\n",
       "      <td>2.949224e+09</td>\n",
       "      <td>2.747563e+09</td>\n",
       "      <td>2.478254e+09</td>\n",
       "      <td>2.129020e+09</td>\n",
       "      <td>1.706331e+09</td>\n",
       "      <td>1.297385e+09</td>\n",
       "      <td>9.277867e+08</td>\n",
       "      <td>6.184912e+08</td>\n",
       "      <td>3.855518e+08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>1.748949e+15</td>\n",
       "      <td>4.982282</td>\n",
       "      <td>13.319280</td>\n",
       "      <td>2.165970e+11</td>\n",
       "      <td>-44.235759</td>\n",
       "      <td>3.610724e+09</td>\n",
       "      <td>3.515567e+09</td>\n",
       "      <td>3.371268e+09</td>\n",
       "      <td>3.167213e+09</td>\n",
       "      <td>2.899294e+09</td>\n",
       "      <td>2.567398e+09</td>\n",
       "      <td>2.180883e+09</td>\n",
       "      <td>1.767772e+09</td>\n",
       "      <td>1.362375e+09</td>\n",
       "      <td>9.922644e+08</td>\n",
       "      <td>6.797724e+08</td>\n",
       "      <td>4.383164e+08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>1.842735e+15</td>\n",
       "      <td>5.382210</td>\n",
       "      <td>19.650624</td>\n",
       "      <td>8.429651e+11</td>\n",
       "      <td>-42.518356</td>\n",
       "      <td>3.788704e+09</td>\n",
       "      <td>3.665377e+09</td>\n",
       "      <td>3.488203e+09</td>\n",
       "      <td>3.248359e+09</td>\n",
       "      <td>2.946643e+09</td>\n",
       "      <td>2.588565e+09</td>\n",
       "      <td>2.192469e+09</td>\n",
       "      <td>1.780434e+09</td>\n",
       "      <td>1.377092e+09</td>\n",
       "      <td>1.009136e+09</td>\n",
       "      <td>6.972862e+08</td>\n",
       "      <td>4.542048e+08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>1.935888e+15</td>\n",
       "      <td>5.853338</td>\n",
       "      <td>28.033491</td>\n",
       "      <td>1.462913e+12</td>\n",
       "      <td>-41.040539</td>\n",
       "      <td>3.989747e+09</td>\n",
       "      <td>3.824820e+09</td>\n",
       "      <td>3.603987e+09</td>\n",
       "      <td>3.323439e+09</td>\n",
       "      <td>2.988215e+09</td>\n",
       "      <td>2.608130e+09</td>\n",
       "      <td>2.204327e+09</td>\n",
       "      <td>1.792831e+09</td>\n",
       "      <td>1.391188e+09</td>\n",
       "      <td>1.024988e+09</td>\n",
       "      <td>7.139524e+08</td>\n",
       "      <td>4.694942e+08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>2.446971e+15</td>\n",
       "      <td>8.061289</td>\n",
       "      <td>74.583217</td>\n",
       "      <td>3.147034e+12</td>\n",
       "      <td>-36.130069</td>\n",
       "      <td>5.001336e+09</td>\n",
       "      <td>4.517091e+09</td>\n",
       "      <td>4.033487e+09</td>\n",
       "      <td>3.552219e+09</td>\n",
       "      <td>3.083537e+09</td>\n",
       "      <td>2.668554e+09</td>\n",
       "      <td>2.253402e+09</td>\n",
       "      <td>1.846411e+09</td>\n",
       "      <td>1.465729e+09</td>\n",
       "      <td>1.113716e+09</td>\n",
       "      <td>8.034876e+08</td>\n",
       "      <td>5.481138e+08</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               mvir            c       rc_kpc          mBCG      loglike  \\\n",
       "count  5.110000e+03  5110.000000  5110.000000  5.110000e+03  5110.000000   \n",
       "mean   1.849007e+15     5.445033    21.496627  9.072186e+11   -42.760813   \n",
       "std    1.385936e+14     0.642713    11.234691  7.100991e+11     2.501500   \n",
       "min    1.468804e+15     3.796233     0.002746  1.001039e+10   -56.173201   \n",
       "25%    1.748949e+15     4.982282    13.319280  2.165970e+11   -44.235759   \n",
       "50%    1.842735e+15     5.382210    19.650624  8.429651e+11   -42.518356   \n",
       "75%    1.935888e+15     5.853338    28.033491  1.462913e+12   -41.040539   \n",
       "max    2.446971e+15     8.061289    74.583217  3.147034e+12   -36.130069   \n",
       "\n",
       "         dm_0008.56    dm_0012.38    dm_0017.90    dm_0025.87    dm_0037.41  \\\n",
       "count  5.110000e+03  5.110000e+03  5.110000e+03  5.110000e+03  5.110000e+03   \n",
       "mean   3.817064e+09  3.679471e+09  3.491438e+09  3.246454e+09  2.943129e+09   \n",
       "std    2.741335e+08  2.175188e+08  1.589643e+08  1.039179e+08  5.842773e+07   \n",
       "min    3.217709e+09  3.169329e+09  3.089798e+09  2.949224e+09  2.747563e+09   \n",
       "25%    3.610724e+09  3.515567e+09  3.371268e+09  3.167213e+09  2.899294e+09   \n",
       "50%    3.788704e+09  3.665377e+09  3.488203e+09  3.248359e+09  2.946643e+09   \n",
       "75%    3.989747e+09  3.824820e+09  3.603987e+09  3.323439e+09  2.988215e+09   \n",
       "max    5.001336e+09  4.517091e+09  4.033487e+09  3.552219e+09  3.083537e+09   \n",
       "\n",
       "         dm_0054.09    dm_0078.21    dm_0113.08    dm_0163.50    dm_0236.41  \\\n",
       "count  5.110000e+03  5.110000e+03  5.110000e+03  5.110000e+03  5.110000e+03   \n",
       "mean   2.587198e+09  2.192480e+09  1.780261e+09  1.376748e+09  1.008605e+09   \n",
       "std    2.812180e+07  1.750103e+07  1.887757e+07  2.159716e+07  2.421675e+07   \n",
       "min    2.478254e+09  2.129020e+09  1.706331e+09  1.297385e+09  9.277867e+08   \n",
       "25%    2.567398e+09  2.180883e+09  1.767772e+09  1.362375e+09  9.922644e+08   \n",
       "50%    2.588565e+09  2.192469e+09  1.780434e+09  1.377092e+09  1.009136e+09   \n",
       "75%    2.608130e+09  2.204327e+09  1.792831e+09  1.391188e+09  1.024988e+09   \n",
       "max    2.668554e+09  2.253402e+09  1.846411e+09  1.465729e+09  1.113716e+09   \n",
       "\n",
       "         dm_0341.81    dm_0494.22  \n",
       "count  5.110000e+03  5.110000e+03  \n",
       "mean   6.973203e+08  4.543723e+08  \n",
       "std    2.503857e+07  2.284468e+07  \n",
       "min    6.184912e+08  3.855518e+08  \n",
       "25%    6.797724e+08  4.383164e+08  \n",
       "50%    6.972862e+08  4.542048e+08  \n",
       "75%    7.139524e+08  4.694942e+08  \n",
       "max    8.034876e+08  5.481138e+08  "
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dfthin.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Other statistics that should accompany the csv file:\n",
    "#  Offset\n",
    "# BCG scale radius and core radius"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "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/opt/anaconda3/envs/py37/lib/python3.7/site-packages/IPython/core/interactiveshell.py:3327: 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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# copy to the 'kappas' folder\n",
    "os.chdir(path)\n",
    "os.system('cp ' + label + '*.bands5.npz  \"' + path + '/../../kappas/\"')\n",
    "os.system('cp ' + label + '.column_names  \"' + path + '/../../kappas/\"')\n",
    "os.system('cp ' + label + '.csv  \"' + path + '/../../kappas/\"')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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
}