notebook/computing_subhalo_clusters.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SubhaloClusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyspark.sql.functions as F\n",
    "from pyspark.sql.functions import col, expr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6.72 ms, sys: 2.47 ms, total: 9.19 ms\n",
      "Wall time: 43.5 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "halos = spark.read.parquet('/datasets/uchuu/RockstarExtendedParquet/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Correspondence between the snapshot number and redshift for all 50 directories: http://www.skiesanduniverses.org/resources/Uchuu_snapshot_redshift_scalefactor.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parameters:\n",
    "- Redshift\n",
    "- Minimum vpeak\n",
    "- Minimum and maximum halo mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "redshift = 0.00\n",
    "min_vpeak = 70\n",
    "\n",
    "# Number of hosts to evaluate with this range: 4\n",
    "cmass_min = 2.0e15\n",
    "cmass_max = 2.03e15\n",
    "\n",
    "# Number of hosts to evaluate with this range: 1140\n",
    "#cmass_min = 1.0e15\n",
    "#cmass_max = 5.0e15\n",
    "\n",
    "# Number of hosts to evaluate with this range: 211022\n",
    "#cmass_min = 1.0e14\n",
    "#cmass_max = 5.0e15\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select HOST halos for `Z=0` (redshift=0) in the Mvir range `cmass_min` - `cmass_max`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "host_halos = halos.where((col('redshift') == redshift)\n",
    "                    & (col('pid') == -1)\n",
    "                    & (col('Mvir') > cmass_min)\n",
    "                    & (col('Mvir') < cmass_max))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We extend the selected hosts with the new 'CvirKly' column and we rename columns to avoid naming collisions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "hosts = (host_halos\n",
    "           .withColumn('CvirKly', expr('Rvir/Rs_Klypin'))\n",
    "           .withColumnRenamed('x', 'host_x')\n",
    "           .withColumnRenamed('y', 'host_y')\n",
    "           .withColumnRenamed('z', 'host_z')\n",
    "           .withColumnRenamed('Rvir', 'host_Rvir')\n",
    "           .select('id', 'host_x', 'host_y', 'host_z', 'host_Rvir'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We select the halos that correspond to satellites (pid != -1) and with a given value of 'Vpeak':"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "satellites = (halos\n",
    "           .where((col('redshift') == redshift) & (col('Vpeak') > min_vpeak) & (col('pid') != -1))\n",
    "           .select('pid', 'x', 'y', 'z', 'Vpeak'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want you can take a reduce random sample of the halos, eg. to use 1% of the data:\n",
    "\n",
    "    myhalos = (halos\n",
    "           .where((col('redshift') == redshift) & (col('Vpeak') > min_vpeak) & (col('pid') != -1))\n",
    "           .select('pid', 'x', 'y', 'z', 'Vpeak')\n",
    "           .sample(0.01))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since it is small, we will cache in memory the hosts dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DataFrame[id: bigint, host_x: double, host_y: double, host_z: double, host_Rvir: double]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hosts.cache()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Number of hosts:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hosts.count()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since hosts is a very reduced subset of data, to calculate the median the simplest way is using `Alternative 4: pandas` to calculate the median (see `analyzing_rockstar_extended_data_with_pyspark` notebook to all the alternatives): "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 194 ms, sys: 58.9 ms, total: 253 ms\n",
      "Wall time: 4.67 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "hosts_pdf = host_halos.withColumn('CvirKly', expr('Rvir/Rs_Klypin')).select('Mvir', 'Rvir', 'CvirKly').toPandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 690 µs, sys: 0 ns, total: 690 µs\n",
      "Wall time: 680 µs\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "Mvir_median_pandas = hosts_pdf['Mvir'].median()\n",
    "Rvir_median_pandas = hosts_pdf['Rvir'].median()\n",
    "CvirKly_median_pandas = hosts_pdf['CvirKly'].median()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Properties Mean host: Mvir median = 2015500000000000.0 Rvir = 2568.325 CvirKly = 5.595388327119377\n"
     ]
    }
   ],
   "source": [
    "print('Properties Mean host:','Mvir median =', Mvir_median_pandas, 'Rvir =', Rvir_median_pandas, 'CvirKly =', CvirKly_median_pandas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Results obtained with VAEX:\n",
    "- Number of hosts = 4\n",
    "- Properties Mean host: Mvir median = 2015529411764706.0 Rvir = 2568.3263725490197 CvirKly = 5.600783982528668"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "subhalos = hosts.join(satellites, hosts.id == satellites.pid, 'left')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = subhalos.withColumn('radius', \n",
    "                             F.sqrt(( F.pow(col('host_x') - col('x'), 2)\n",
    "                                    + F.pow(col('host_y') - col('y'), 2)\n",
    "                                    + F.pow(col('host_z') - col('z'), 2)))\n",
    "                             /(col('host_Rvir')/1000.))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+---------------+--------------------+-----+\n",
      "|             id|(host_Rvir / 1000.0)|count|\n",
      "+---------------+--------------------+-----+\n",
      "|  3814028194986|             2.57321|  857|\n",
      "|214679739377433|             2.56777|  734|\n",
      "|139500640393732|  2.5677600000000003|  728|\n",
      "|101876720636316|             2.56888|  805|\n",
      "+---------------+--------------------+-----+\n",
      "\n",
      "CPU times: user 4.94 ms, sys: 2.37 ms, total: 7.3 ms\n",
      "Wall time: 20.6 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "result.groupBy('id', 'host_Rvir').count().selectExpr('id', 'host_Rvir/1000.0', 'count').show()\n",
    "# Recordar dividir host_Rvir entre 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 5.61 ms, sys: 4.81 ms, total: 10.4 ms\n",
      "Wall time: 16.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "radius_quantiles = result.approxQuantile('radius', [0.025, 0.5, 0.975], 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6.69 ms, sys: 1.47 ms, total: 8.16 ms\n",
      "Wall time: 14.5 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "rows = result.agg(F.stddev('radius').alias('std')).collect()\n",
    "radius_std = rows[0].std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Statistics radius:\n",
      "\tsigma = 0.22744851346072806\n",
      "\tquantile_0.025 = 0.20744410621018508\n",
      "\tmedian = 0.6424474540350691\n",
      "\tquantile_0.975 = 0.9838158737789053\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(f'Statistics radius:\\n'\n",
    "      f'\\tsigma = {radius_std}\\n'\n",
    "      f'\\tquantile_0.025 = {radius_quantiles[0]}\\n'\n",
    "      f'\\tmedian = {radius_quantiles[1]}\\n'\n",
    "      f'\\tquantile_0.975 = {radius_quantiles[2]}\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing executing times: VAEX vs Spark"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "VAEX running on FT vs Spark running on Hadoop3. \n",
    "\n",
    "The times taken at the Hadoop3 cluster will vary depending on how many resources are assigned to the application. In this case I run them when compiting with other resource intensive application from other user. When running alone the times are even smaller.\n",
    "\n",
    "The benchmark was run using the smallest example that computes satellites around only 4 hosts:\n",
    "\n",
    "Results for vaex stored on `snap-50-subhalo-clusters.o4218594` in FT, in 24 hours the calculation could only finish calculating the first host, so to have an estimation of the total time we assume the the other 3 will take a similar time. Addionally the VAEX calculation will have an additional time to compute the global statistics at the end that is not taken into account but that will increase even further the total time.\n",
    "\n",
    "|-                         | VAEX                     | Spark                     |\n",
    "|--------------------------|--------------------------|---------------------------|\n",
    "| Preliminary computations |  109 min                 | less than 1 min           |\n",
    "| Finding halos            |  9 hours only first host | less than 2 min ALL hosts |\n",
    "| Total time               |  more than 38 hours (est)| less than 3 min           |\n",
    "\n",
    "\n",
    "For the medium case example with 1140 hosts the time needed by Spark is less than 6 minutes.\n"
   ]
  }
 ],
 "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": 4
}