{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Acceleration with numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys\n",
    "sys.path.insert(0, os.path.abspath(\"..\"))\n",
    "\n",
    "import geoclide as gc\n",
    "from time import process_time\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## BBox -Ray intersection test, multiples bboxes and 1 ray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here we create 100000 bounding boxes and 1 ray\n",
    "nx = 100\n",
    "ny = 100\n",
    "nz = 10\n",
    "x = np.linspace(0., nx-1, nx, np.float64)\n",
    "y = np.linspace(0., ny-1, ny, np.float64)\n",
    "z = np.linspace(0., nz-1, nz, np.float64)\n",
    "x_, y_, z_ = np.meshgrid(x,y,z, indexing='ij')\n",
    "pmin_arr = np.vstack((x_.ravel(), y_.ravel(), z_.ravel())).T\n",
    "x = np.linspace(1., nx, nx, np.float64)\n",
    "y = np.linspace(1., ny, ny, np.float64)\n",
    "z = np.linspace(1., nz, nz, np.float64)\n",
    "x_, y_, z_ = np.meshgrid(x,y,z, indexing='ij')\n",
    "pmax_arr = np.vstack((x_.ravel(), y_.ravel(), z_.ravel())).T\n",
    "r0 = gc.Ray(gc.Point(-2., 0., 0.25), gc.normalize(gc.Vector(1, 0., 0.5)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Intersection tests using a loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "elapsed time (s) using loop:  0.8166620819999997\n"
     ]
    }
   ],
   "source": [
    "start = process_time()\n",
    "nboxes = pmin_arr.shape[0]\n",
    "t0_ = np.zeros(nboxes, dtype=np.float64)\n",
    "t1_ = np.zeros_like(t0_)\n",
    "is_int_ = np.full(nboxes, False, dtype=bool)\n",
    "for ib in range (0, nboxes):\n",
    "    bi = gc.BBox(gc.Point(pmin_arr[ib,:]), gc.Point(pmax_arr[ib,:]))\n",
    "    t0_[ib], t1_[ib], is_int_[ib] = bi.intersect(r0, ds_output=False)\n",
    "end = process_time()\n",
    "print(\"elapsed time (s) using loop: \", end - start)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Intersection tests using numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "elapsed time (s) using numpy:  0.011210285999999847\n"
     ]
    }
   ],
   "source": [
    "start = process_time()\n",
    "pmin = gc.Point(pmin_arr)\n",
    "pmax = gc.Point(pmax_arr)\n",
    "b_set = gc.BBox(pmin, pmax)\n",
    "t0, t1, is_int1 = b_set.intersect(r0, ds_output=False)\n",
    "end = process_time()\n",
    "print(\"elapsed time (s) using numpy: \", end - start)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## BBox - Ray intersection test, multiples bboxes and multiple rays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Case 1: for each ray, perform intersection test with all the bounding boxes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We create 400 bounding boxes and 400 rays\n",
    "nx = 20\n",
    "ny = 20\n",
    "nz = 1\n",
    "x = np.linspace(0., nx-1, nx, np.float64)\n",
    "y = np.linspace(0., ny-1, ny, np.float64)\n",
    "z = np.linspace(0., nz-1, nz, np.float64)\n",
    "x_, y_, z_ = np.meshgrid(x,y,z, indexing='ij')\n",
    "pmin_arr = np.vstack((x_.ravel(), y_.ravel(), z_.ravel())).T\n",
    "x = np.linspace(1., nx, nx, np.float64)\n",
    "y = np.linspace(1., ny, ny, np.float64)\n",
    "z = np.linspace(1., nz, nz, np.float64)\n",
    "x_, y_, z_ = np.meshgrid(x,y,z, indexing='ij')\n",
    "pmax_arr = np.vstack((x_.ravel(), y_.ravel(), z_.ravel())).T\n",
    "nboxes = pmin_arr.shape[0]\n",
    "x_, y_, z_ = np.meshgrid(np.linspace(0.5, nx-0.5, nx, np.float64),\n",
    "                        np.linspace(0.5, ny-0.5, ny, np.float64),\n",
    "                        nz+1, indexing='ij')\n",
    "\n",
    "o_set_arr = np.vstack((x_.ravel(), y_.ravel(), z_.ravel())).T\n",
    "nrays = o_set_arr.shape[0]\n",
    "d_set_arr = np.zeros_like(o_set_arr)\n",
    "d_set_arr[:,0] = 0.\n",
    "d_set_arr[:,1] = 0.\n",
    "d_set_arr[:,2] = -1.\n",
    "o_set = gc.Point(o_set_arr)\n",
    "d_set = gc.Vector(d_set_arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The tests using loops"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "case 1 - elapsed time (s) using loops: 0.3339018809999996\n"
     ]
    }
   ],
   "source": [
    "start = process_time()\n",
    "t0_ = np.zeros((nboxes, nrays), dtype=np.float64)\n",
    "t1_ = np.zeros_like(t0_)\n",
    "is_int_ = np.full((nboxes,nrays), False, dtype=bool)\n",
    "list_rays = []\n",
    "for ir in range(0, nrays):\n",
    "  list_rays.append(gc.Ray(gc.Point(o_set_arr[ir,:]),\n",
    "                          gc.normalize(gc.Vector(d_set_arr[ir,:]))))\n",
    "for ib in range (0, nboxes):\n",
    "  bi = gc.BBox(gc.Point(pmin_arr[ib,:]), gc.Point(pmax_arr[ib,:]))\n",
    "  for ir in range(0, nrays):\n",
    "      t0_[ib,ir], t1_[ib,ir], is_int_[ib,ir] = bi.intersect(list_rays[ir], ds_output=False)\n",
    "end = process_time()\n",
    "print(\"case 1 - elapsed time (s) using loops:\", end-start)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The tests using numpy calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "case 1 - elapsed time (s) using numpy: 0.0032663020000001097\n"
     ]
    }
   ],
   "source": [
    "start = process_time()\n",
    "r_set = gc.Ray(o_set, d_set)\n",
    "pmin = gc.Point(pmin_arr)\n",
    "pmax = gc.Point(pmax_arr)\n",
    "b_set = gc.BBox(pmin, pmax)\n",
    "t0, t1, is_int1 = b_set.intersect(r_set, ds_output=False)\n",
    "end = process_time()\n",
    "time_fast = end-start\n",
    "print(\"case 1 - elapsed time (s) using numpy:\", end-start)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Case 2: diagonal calculations i.e. test only between ray(i) and bbox(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We create 40000 bounding boxes and 40000 rays\n",
    "nx = 200\n",
    "ny = 200\n",
    "nz = 1\n",
    "x = np.linspace(0., nx-1, nx, np.float64)\n",
    "y = np.linspace(0., ny-1, ny, np.float64)\n",
    "z = np.linspace(0., nz-1, nz, np.float64)\n",
    "x_, y_, z_ = np.meshgrid(x,y,z, indexing='ij')\n",
    "pmin_arr = np.vstack((x_.ravel(), y_.ravel(), z_.ravel())).T\n",
    "x = np.linspace(1., nx, nx, np.float64)\n",
    "y = np.linspace(1., ny, ny, np.float64)\n",
    "z = np.linspace(1., nz, nz, np.float64)\n",
    "x_, y_, z_ = np.meshgrid(x,y,z, indexing='ij')\n",
    "pmax_arr = np.vstack((x_.ravel(), y_.ravel(), z_.ravel())).T\n",
    "nboxes = pmin_arr.shape[0]\n",
    "x_, y_, z_ = np.meshgrid(np.linspace(0.5, nx-0.5, nx, np.float64),\n",
    "                        np.linspace(0.5, ny-0.5, ny, np.float64),\n",
    "                        nz+1, indexing='ij')\n",
    "\n",
    "o_set_arr = np.vstack((x_.ravel(), y_.ravel(), z_.ravel())).T\n",
    "nrays = o_set_arr.shape[0]\n",
    "d_set_arr = np.zeros_like(o_set_arr)\n",
    "d_set_arr[:,0] = 0.\n",
    "d_set_arr[:,1] = 0.\n",
    "d_set_arr[:,2] = -1.\n",
    "o_set = gc.Point(o_set_arr)\n",
    "d_set = gc.Vector(d_set_arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The tests using loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "case 2 - elapsed time (s) using loops: 0.421498755\n"
     ]
    }
   ],
   "source": [
    "start = process_time()\n",
    "t0_ = np.zeros((nboxes), dtype=np.float64)\n",
    "t1_ = np.zeros_like(t0_)\n",
    "is_int_ = np.full((nboxes), False, dtype=bool)\n",
    "list_rays = []\n",
    "for ib in range(0, nboxes):\n",
    "    bi = gc.BBox(gc.Point(pmin_arr[ib,:]), gc.Point(pmax_arr[ib,:]))\n",
    "    ri = gc.Ray(gc.Point(o_set_arr[ib,:]), gc.Vector(d_set_arr[ib,:]))\n",
    "    t0_[ib], t1_[ib], is_int_[ib] = bi.intersect(ri, ds_output=False)\n",
    "end = process_time()\n",
    "print(\"case 2 - elapsed time (s) using loops:\", end-start)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The tests using numpy calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "case 2 - elapsed time (s) using numpy: 0.001745509000000034\n"
     ]
    }
   ],
   "source": [
    "start = process_time()\n",
    "r_set = gc.Ray(o_set, d_set)\n",
    "pmin = gc.Point(pmin_arr)\n",
    "pmax = gc.Point(pmax_arr)\n",
    "b_set = gc.BBox(pmin, pmax)\n",
    "t0, t1, is_int1 = b_set.intersect(r_set, diag_calc=True, ds_output=False)\n",
    "end = process_time()\n",
    "print(\"case 2 - elapsed time (s) using numpy:\", end-start)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "geoclide-py313",
   "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.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}