earth.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import meshio\n",
    "import numpy as np\n",
    "from multiprocessing import Pool\n",
    "import time\n",
    "\n",
    "from ugs import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "infile = \"earth/data/globe.vtk\"\n",
    "outfile = \"costed_earth.vtk\"\n",
    "mesh = meshio.read(infile)\n",
    "\n",
    "max_distance = 500000  #\n",
    "                      # A smaller value means visiting far fewer nodes, so it speeds things up a lot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We need to define a way to calculate cost\n",
    "@functools.lru_cache(maxsize=32768)\n",
    "def distance_with_elevation_scaled(mesh, current, _next):\n",
    "    # The travel_cost can be any function, including just the distance.\n",
    "    # Here, we exagerate the elevation difference, to make changing elevation more costly\n",
    "    if current == _next:\n",
    "        return 0\n",
    "    \n",
    "    z_scaling = 100.  # 100. is a random number to pick, but has quite a big impact on the resulting paths.\n",
    "    \n",
    "    new_current = np.append(mesh.points[current][:2], mesh.point_data['Z'][current] * z_scaling)\n",
    "    new_next    = np.append(mesh.points[_next][:2],   mesh.point_data['Z'][_next]   * z_scaling)\n",
    "    \n",
    "    return int(np.linalg.norm(new_current - new_next))  # return as Int, just for niceness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We need to define a way to calculate cost\n",
    "@functools.lru_cache(maxsize=32768)\n",
    "def elevation_only(mesh, current, _next):\n",
    "    # Only take into account elevation changes for costs\n",
    "    if current == _next:\n",
    "        return 0\n",
    "    return int(abs(mesh.point_data['Z'][current] - mesh.point_data['Z'][_next]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose which cost function you want\n",
    "travel_cost = elevation_only"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare for parallel execution\n",
    "\n",
    "We want to use a parallel map function to calculate the cost for all points. We can make this easier by 'baking in' the parameters of a function we know, and leave only the `starting point` as a variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "\n",
    "# We bake in the mesh and travel_cost function into the get_from_point function\n",
    "get_from_point_in_mesh = partial(get_from_point, mesh = mesh, travel_cost_function = travel_cost, max_distance = max_distance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2087, 0)\n"
     ]
    }
   ],
   "source": [
    "# Now we can use the get_from_point function with only a point ID as a parameter, by going via the new get_from_point_in_mesh function\n",
    "# Here we do a test, and see the output format: (point, total cost of all paths to that point)\n",
    "print(get_from_point_in_mesh(2087))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare input data\n",
    "\n",
    "We don't want to calculate the LEC of a point below sea-level, so here we find all the points above sea-level, so they can be used as starting points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "points_above_sealevel = np.nonzero(mesh.point_data['Z'] >= 0)[0][::-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total starting points available:  183646\n"
     ]
    }
   ],
   "source": [
    "print(\"Total starting points available: \", points_above_sealevel.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Process ForkPoolWorker-5:\n",
      "Process ForkPoolWorker-3:\n",
      "Process ForkPoolWorker-4:\n",
      "Process ForkPoolWorker-2:\n",
      "Process ForkPoolWorker-6:\n",
      "Process ForkPoolWorker-1:\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "Traceback (most recent call last):\n",
      "  File \"/usr/lib/python3.5/multiprocessing/process.py\", line 249, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/usr/lib/python3.5/multiprocessing/process.py\", line 249, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/usr/lib/python3.5/multiprocessing/process.py\", line 249, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/usr/lib/python3.5/multiprocessing/process.py\", line 249, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/usr/lib/python3.5/multiprocessing/process.py\", line 249, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/usr/lib/python3.5/multiprocessing/process.py\", line 249, in _bootstrap\n",
      "    self.run()\n",
      "  File \"/usr/lib/python3.5/multiprocessing/process.py\", line 93, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/usr/lib/python3.5/multiprocessing/process.py\", line 93, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/usr/lib/python3.5/multiprocessing/process.py\", line 93, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/usr/lib/python3.5/multiprocessing/process.py\", line 93, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/usr/lib/python3.5/multiprocessing/process.py\", line 93, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/usr/lib/python3.5/multiprocessing/process.py\", line 93, in run\n",
      "    self._target(*self._args, **self._kwargs)\n",
      "  File \"/usr/lib/python3.5/multiprocessing/pool.py\", line 119, in worker\n",
      "    result = (True, func(*args, **kwds))\n",
      "  File \"/usr/lib/python3.5/multiprocessing/pool.py\", line 119, in worker\n",
      "    result = (True, func(*args, **kwds))\n",
      "  File \"/usr/lib/python3.5/multiprocessing/pool.py\", line 119, in worker\n",
      "    result = (True, func(*args, **kwds))\n",
      "  File \"/usr/lib/python3.5/multiprocessing/pool.py\", line 119, in worker\n",
      "    result = (True, func(*args, **kwds))\n",
      "  File \"/usr/lib/python3.5/multiprocessing/pool.py\", line 119, in worker\n",
      "    result = (True, func(*args, **kwds))\n",
      "  File \"/usr/lib/python3.5/multiprocessing/pool.py\", line 44, in mapstar\n",
      "    return list(map(*args))\n",
      "  File \"/usr/lib/python3.5/multiprocessing/pool.py\", line 44, in mapstar\n",
      "    return list(map(*args))\n",
      "  File \"/usr/lib/python3.5/multiprocessing/pool.py\", line 44, in mapstar\n",
      "    return list(map(*args))\n",
      "  File \"/usr/lib/python3.5/multiprocessing/pool.py\", line 44, in mapstar\n",
      "    return list(map(*args))\n",
      "  File \"/live/share/ugs.py\", line 160, in get_from_point\n",
      "    total_cost = get_total_cost_for_point(point, mesh, travel_cost_function, max_distance)\n",
      "  File \"/usr/lib/python3.5/multiprocessing/pool.py\", line 44, in mapstar\n",
      "    return list(map(*args))\n",
      "  File \"/live/share/ugs.py\", line 160, in get_from_point\n",
      "    total_cost = get_total_cost_for_point(point, mesh, travel_cost_function, max_distance)\n",
      "  File \"/live/share/ugs.py\", line 160, in get_from_point\n",
      "    total_cost = get_total_cost_for_point(point, mesh, travel_cost_function, max_distance)\n",
      "  File \"/live/share/ugs.py\", line 127, in get_total_cost_for_point\n",
      "    came_from, cost_so_far, dist_so_far = cost_search(mesh, start, travel_cost_function, max_distance)\n",
      "  File \"/live/share/ugs.py\", line 160, in get_from_point\n",
      "    total_cost = get_total_cost_for_point(point, mesh, travel_cost_function, max_distance)\n",
      "  File \"/live/share/ugs.py\", line 160, in get_from_point\n",
      "    total_cost = get_total_cost_for_point(point, mesh, travel_cost_function, max_distance)\n",
      "  File \"/live/share/ugs.py\", line 127, in get_total_cost_for_point\n",
      "    came_from, cost_so_far, dist_so_far = cost_search(mesh, start, travel_cost_function, max_distance)\n",
      "  File \"/live/share/ugs.py\", line 127, in get_total_cost_for_point\n",
      "    came_from, cost_so_far, dist_so_far = cost_search(mesh, start, travel_cost_function, max_distance)\n",
      "  File \"/live/share/ugs.py\", line 75, in cost_search\n",
      "    for _next in graph_neighbours(mesh, current):\n",
      "  File \"/live/share/ugs.py\", line 127, in get_total_cost_for_point\n",
      "    came_from, cost_so_far, dist_so_far = cost_search(mesh, start, travel_cost_function, max_distance)\n",
      "  File \"/live/share/ugs.py\", line 127, in get_total_cost_for_point\n",
      "    came_from, cost_so_far, dist_so_far = cost_search(mesh, start, travel_cost_function, max_distance)\n",
      "  File \"/live/share/ugs.py\", line 75, in cost_search\n",
      "    for _next in graph_neighbours(mesh, current):\n",
      "  File \"/live/share/ugs.py\", line 75, in cost_search\n",
      "    for _next in graph_neighbours(mesh, current):\n",
      "  File \"/live/share/ugs.py\", line 77, in cost_search\n",
      "    new_cost = cost_so_far[current] + travel_cost_function(mesh, current, _next)\n",
      "  File \"/live/share/ugs.py\", line 75, in cost_search\n",
      "    for _next in graph_neighbours(mesh, current):\n",
      "  File \"/live/share/ugs.py\", line 19, in graph_neighbours\n",
      "    aa= mesh.cells['triangle']==current\n",
      "KeyboardInterrupt\n",
      "  File \"/live/share/ugs.py\", line 20, in graph_neighbours\n",
      "    a = np.where(aa)[0]\n",
      "KeyboardInterrupt\n",
      "  File \"/live/share/ugs.py\", line 19, in graph_neighbours\n",
      "    aa= mesh.cells['triangle']==current\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "KeyboardInterrupt\n",
      "  File \"/usr/lib/python3.5/multiprocessing/pool.py\", line 119, in worker\n",
      "    result = (True, func(*args, **kwds))\n",
      "  File \"/usr/lib/python3.5/multiprocessing/pool.py\", line 44, in mapstar\n",
      "    return list(map(*args))\n",
      "  File \"/live/share/ugs.py\", line 160, in get_from_point\n",
      "    total_cost = get_total_cost_for_point(point, mesh, travel_cost_function, max_distance)\n",
      "  File \"/live/share/ugs.py\", line 127, in get_total_cost_for_point\n",
      "    came_from, cost_so_far, dist_so_far = cost_search(mesh, start, travel_cost_function, max_distance)\n",
      "  File \"/live/share/ugs.py\", line 79, in cost_search\n",
      "    new_dist = dist_so_far[current] + distance(mesh, current, _next)\n",
      "  File \"/live/share/ugs.py\", line 7, in distance\n",
      "    @functools.lru_cache(maxsize=32768)\n",
      "KeyboardInterrupt\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-10-ed2c60ace9ec>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     16\u001b[0m     \u001b[0mstart_time\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     17\u001b[0m     \u001b[0;31m# Here we use a parallel map function to send out the chunk across the CPUs in the Pool\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 18\u001b[0;31m     \u001b[0mcosts\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mget_from_point_in_mesh\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpoints_above_sealevel\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mstart\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mstop\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     19\u001b[0m     \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"From \"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstart\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\" to \"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstop\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\" took \"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mstart_time\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"seconds, starting with point\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpoints_above_sealevel\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mstart\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\". Percent complete: \"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m100\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfloat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstop\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mpoints_above_sealevel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     20\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.5/multiprocessing/pool.py\u001b[0m in \u001b[0;36mmap\u001b[0;34m(self, func, iterable, chunksize)\u001b[0m\n\u001b[1;32m    258\u001b[0m         \u001b[0;32min\u001b[0m \u001b[0ma\u001b[0m \u001b[0mlist\u001b[0m \u001b[0mthat\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0mreturned\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    259\u001b[0m         '''\n\u001b[0;32m--> 260\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_map_async\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0miterable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmapstar\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mchunksize\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    261\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    262\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mstarmap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0miterable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mchunksize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.5/multiprocessing/pool.py\u001b[0m in \u001b[0;36mget\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    600\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    601\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 602\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwait\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    603\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mready\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    604\u001b[0m             \u001b[0;32mraise\u001b[0m \u001b[0mTimeoutError\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.5/multiprocessing/pool.py\u001b[0m in \u001b[0;36mwait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    597\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    598\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mwait\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 599\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_event\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwait\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    600\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    601\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.5/threading.py\u001b[0m in \u001b[0;36mwait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    547\u001b[0m             \u001b[0msignaled\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_flag\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    548\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0msignaled\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 549\u001b[0;31m                 \u001b[0msignaled\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cond\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwait\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    550\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0msignaled\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    551\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.5/threading.py\u001b[0m in \u001b[0;36mwait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    291\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m    \u001b[0;31m# restore state no matter what (e.g., KeyboardInterrupt)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    292\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mtimeout\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 293\u001b[0;31m                 \u001b[0mwaiter\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0macquire\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    294\u001b[0m                 \u001b[0mgotit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    295\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# Output data variables\n",
    "all_costs = []\n",
    "mesh.point_data['cost'] = np.zeros_like(mesh.point_data['Z'])\n",
    "\n",
    "# For parallel execution, we chunk the data up into increments (inc), so we can monitor the progress\n",
    "# Larger chunks may be better able to take advantage of the LRU caching, but can lead to some CPUs finishing their batch \n",
    "# much earlier, and they just wait around until the next chunk is loaded.\n",
    "start = 0\n",
    "inc = 100\n",
    "stop = inc\n",
    "\n",
    "# Run on 11 CPUs\n",
    "p = Pool(6)\n",
    "\n",
    "while start <  points_above_sealevel.shape[0]-1:\n",
    "    start_time = time.time()\n",
    "    # Here we use a parallel map function to send out the chunk across the CPUs in the Pool\n",
    "    costs = p.map(get_from_point_in_mesh, points_above_sealevel[start:stop])\n",
    "    print(\"From \", start, \" to \", stop, \" took \", time.time() - start_time, \"seconds, starting with point\", points_above_sealevel[start], \". Percent complete: \", 100*(float(stop)/points_above_sealevel.shape[0]))\n",
    "    \n",
    "    # Save the data\n",
    "    all_costs.extend(costs)\n",
    "        \n",
    "    # Write out data progressively, so we can see progress in Paraview\n",
    "    for i in all_costs:\n",
    "        mesh.point_data['cost'][i[0]] = i[1]\n",
    "    meshio.write(outfile, mesh)\n",
    "    \n",
    "    # move to the next chunk of data\n",
    "    start += inc\n",
    "    stop += inc\n",
    "    if stop >= points_above_sealevel.shape[0]:\n",
    "        stop = points_above_sealevel.shape[0]-1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "## Output\n",
    "\n",
    "Show some of the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print all_costs[:100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}