PythonicGeodynamicsMorra/Ch03 Fast Computing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How fast is your Computing Machine?\n",
    "\n",
    "Our computers do basically two things, one is storing numbers, the other is doing mathematical operations on them. They are so useful because they do these things extraordinarily fast and extraordinarily well, in the sense that they always do them correctly following our instructions. For this reason the main way in which we use them is to give them repetitive instructions for up to billions of times.\n",
    "\n",
    "Since the speed of the CPU has doubled every 1.5 years until a few years ago (Moore’s Law), we may have lost track of how really fast the computer in our hands is. When starting this journey into numerical modeling, it is important to first quickly experiment this and not only realize how remarkably fast it really is, but also to perceive its limits when we reach it.\n",
    "\n",
    "In this chapter, we start using iPython because its interactive setting allows us to feel better about what the computer does. Later we will step into writing and running a program. Let us initially load NumPy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# this function simply counts from 0 to max\n",
    "def count(max):\n",
    "    i=0\n",
    "    while i<max:\n",
    "        i=i+1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The %timeit command is used to calculate the average time of performing certain operations in python. As an example let us use this command to calculate the time that it takes to count from 0 to max in the above function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.59 ms ± 400 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit count(100000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "51.2 ms ± 2.2 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit count(1000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "543 ms ± 5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit count(10000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# also this function counts, but increasingly larger numbers\n",
    "def addCount(max):\n",
    "    a=0.0\n",
    "    for i in range(max):\n",
    "        a=a+i\n",
    "    #print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When I try to launch count(1000) or count(100000) on my laptop I simply do not observe any difference. Just the time to call the function, and I get the result. It is unnoticeable that the second call required 100 times more calculation than the first one. Only when I arrive at one million or better ten millions, I finally observe a lag. In fact what I notice is the time that it takes for one processor of my laptop to perform this serial calculation. You can test yours and perceive the power of your computer.\n",
    "\n",
    "We can also be more quantitative about the time required by our little routine. iPython offers a very simple and powerful magic command called %timeit. Let us test it for testing the time required by our loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.3 ms ± 164 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit addCount(100000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "52.3 ms ± 5.73 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit addCount(1000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "484 ms ± 8.31 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit addCount(10000000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's focus on the Numerical Tools of Python. Among others, NumPy adds to standard Python the following features:   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">A. Multidimensional vectorized arrays<br />\n",
    ">B. Mathematical functions operating on an array or portions of it<br />\n",
    ">C. Linear Algebra, Fourier development, Random Functions<br />\n",
    "      D. Input/Output functions to efficiently create and read memory mapped files<br />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NumPy Types\n",
    "\n",
    "Without any further specification, Python automatically uses float type machine precision. We might not know which precision this is, in which case it is important to discover it. A very straightforward way to discover the precision of our machine is by testing for which small ε, a and a + ε become indistinguishable. We can do so by iteratively decreasing ε of a certain ratio:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Precision: 1.0\n",
      "Precision: 0.1111111111111111\n",
      "Precision: 0.012345679012345678\n",
      "Precision: 0.0013717421124828531\n",
      "Precision: 0.00015241579027587256\n",
      "Precision: 1.6935087808430286e-05\n",
      "Precision: 1.8816764231589206e-06\n",
      "Precision: 2.0907515812876894e-07\n",
      "Precision: 2.323057312541877e-08\n",
      "Precision: 2.5811747917131966e-09\n",
      "Precision: 2.867971990792441e-10\n",
      "Precision: 3.1866355453249343e-11\n",
      "Precision: 3.5407061614721493e-12\n",
      "Precision: 3.934117957191277e-13\n",
      "Precision: 4.371242174656974e-14\n",
      "Precision: 4.85693574961886e-15\n",
      "Precision: 5.396595277354289e-16\n",
      "Precision: 5.996216974838099e-17\n",
      "Precision: 6.6624633053756655e-18\n",
      "Precision: 7.402737005972962e-19\n",
      "Precision: 8.225263339969957e-20\n",
      "Precision: 9.139181488855508e-21\n",
      "Precision: 1.0154646098728343e-21\n",
      "Precision: 1.1282940109698158e-22\n",
      "Precision: 1.2536600121886843e-23\n",
      "Precision: 1.392955569098538e-24\n",
      "Precision: 1.5477284101094868e-25\n",
      "Precision: 1.7196982334549854e-26\n",
      "Precision: 1.9107758149499837e-27\n",
      "Precision: 2.1230842388333153e-28\n",
      "Precision: 2.3589824875925725e-29\n",
      "Precision: 2.621091652880636e-30\n",
      "Precision: 2.9123240587562624e-31\n",
      "Precision: 3.2359156208402915e-32\n",
      "Precision: 3.595461800933657e-33\n",
      "Precision: 3.994957556592953e-34\n",
      "Precision: 4.4388417295477253e-35\n",
      "Precision: 4.932046366164139e-36\n",
      "Precision: 5.480051517960155e-37\n",
      "Precision: 6.088946131066838e-38\n",
      "Precision: 6.765495701185375e-39\n",
      "Precision: 7.517217445761529e-40\n",
      "Precision: 8.352463828623921e-41\n",
      "Precision: 9.28051536513769e-42\n",
      "Precision: 1.0311683739041878e-42\n",
      "Precision: 1.1457426376713199e-43\n",
      "Precision: 1.2730473751903554e-44\n",
      "Precision: 1.4144970835448395e-45\n",
      "Precision: 1.5716634261609329e-46\n",
      "Precision: 1.74629269573437e-47\n",
      "Precision: 1.9403252174826334e-48\n",
      "Precision: 2.155916908314037e-49\n",
      "Precision: 2.395463231460041e-50\n",
      "Precision: 2.6616258127333788e-51\n",
      "Precision: 2.9573620141481986e-52\n",
      "Precision: 3.2859577934979986e-53\n",
      "Precision: 3.651064214997776e-54\n",
      "Precision: 4.056738016664196e-55\n",
      "Precision: 4.50748668518244e-56\n",
      "Precision: 5.0083185390916e-57\n",
      "Precision: 5.5647983767684445e-58\n",
      "Precision: 6.183109307520494e-59\n",
      "Precision: 6.870121452800549e-60\n",
      "Precision: 7.633468280889499e-61\n",
      "Precision: 8.481631423210554e-62\n",
      "Precision: 9.424034914678393e-63\n",
      "Precision: 1.0471149905198215e-63\n",
      "Precision: 1.1634611005775795e-64\n",
      "Precision: 1.2927345561973105e-65\n",
      "Precision: 1.4363717291081228e-66\n",
      "Precision: 1.5959685878979142e-67\n",
      "Precision: 1.7732984309976824e-68\n",
      "Precision: 1.970331589997425e-69\n",
      "Precision: 2.1892573222193612e-70\n",
      "Precision: 2.4325081357992902e-71\n",
      "Precision: 2.702786817554767e-72\n",
      "Precision: 3.003096463949741e-73\n",
      "Precision: 3.3367738488330452e-74\n",
      "Precision: 3.7075264987033837e-75\n",
      "Precision: 4.1194738874482043e-76\n",
      "Precision: 4.5771932082757826e-77\n",
      "Precision: 5.085770231417536e-78\n",
      "Precision: 5.650855812686151e-79\n",
      "Precision: 6.27872868076239e-80\n",
      "Precision: 6.9763652008471e-81\n",
      "Precision: 7.751516889830111e-82\n",
      "Precision: 8.612796544255678e-83\n",
      "Precision: 9.569773938061865e-84\n",
      "Precision: 1.0633082153402071e-84\n",
      "Precision: 1.1814535726002301e-85\n",
      "Precision: 1.3127261917780335e-86\n",
      "Precision: 1.4585846575311483e-87\n",
      "Precision: 1.6206496194790536e-88\n",
      "Precision: 1.8007217994211708e-89\n",
      "Precision: 2.0008019993568564e-90\n",
      "Precision: 2.2231133326187293e-91\n",
      "Precision: 2.4701259251319215e-92\n",
      "Precision: 2.7445843612576906e-93\n",
      "Precision: 3.049538179175212e-94\n",
      "Precision: 3.388375754639125e-95\n",
      "Precision: 3.764861949599027e-96\n",
      "Precision: 4.183179943998919e-97\n",
      "Precision: 4.6479777155543544e-98\n",
      "Precision: 5.164419683949283e-99\n",
      "Precision: 5.738244093276981e-100\n",
      "Precision: 6.375826770307757e-101\n",
      "Precision: 7.084251967008619e-102\n",
      "Precision: 7.871391074454021e-103\n",
      "Precision: 8.74599008272669e-104\n",
      "Precision: 9.717766758585212e-105\n",
      "Precision: 1.0797518620650235e-105\n",
      "Precision: 1.1997242911833593e-106\n",
      "Precision: 1.3330269902037326e-107\n",
      "Precision: 1.4811411002263697e-108\n",
      "Precision: 1.6457123335848551e-109\n",
      "Precision: 1.828569259538728e-110\n",
      "Precision: 2.0317436217096978e-111\n",
      "Precision: 2.2574929130107753e-112\n",
      "Precision: 2.5083254589008614e-113\n",
      "Precision: 2.787028287667624e-114\n",
      "Precision: 3.096698097408471e-115\n",
      "Precision: 3.44077566378719e-116\n",
      "Precision: 3.823084070874656e-117\n",
      "Precision: 4.247871189860729e-118\n",
      "Precision: 4.7198568776230323e-119\n",
      "Precision: 5.244285419581147e-120\n",
      "Precision: 5.8269837995346075e-121\n",
      "Precision: 6.474426443927341e-122\n",
      "Precision: 7.193807159919269e-123\n",
      "Precision: 7.993119066576966e-124\n",
      "Precision: 8.88124340730774e-125\n",
      "Precision: 9.868048230341933e-126\n",
      "Precision: 1.096449803371326e-126\n",
      "Precision: 1.2182775593014733e-127\n",
      "Precision: 1.3536417325571926e-128\n",
      "Precision: 1.5040463695079916e-129\n",
      "Precision: 1.6711626327866573e-130\n",
      "Precision: 1.8568473697629525e-131\n",
      "Precision: 2.0631637441810583e-132\n",
      "Precision: 2.292404160201176e-133\n",
      "Precision: 2.5471157335568623e-134\n",
      "Precision: 2.830128592840958e-135\n",
      "Precision: 3.1445873253788422e-136\n",
      "Precision: 3.4939859170876023e-137\n",
      "Precision: 3.88220657454178e-138\n",
      "Precision: 4.3135628606019784e-139\n",
      "Precision: 4.792847622891088e-140\n",
      "Precision: 5.325386247656764e-141\n",
      "Precision: 5.917095830729738e-142\n",
      "Precision: 6.574550923033042e-143\n",
      "Precision: 7.305056581147824e-144\n",
      "Precision: 8.116729534608694e-145\n",
      "Precision: 9.018588371787438e-146\n",
      "Precision: 1.0020653746430487e-146\n",
      "Precision: 1.1134059718256097e-147\n",
      "Precision: 1.2371177464728997e-148\n",
      "Precision: 1.3745752738587775e-149\n",
      "Precision: 1.5273058598430862e-150\n",
      "Precision: 1.6970065109367624e-151\n",
      "Precision: 1.885562789929736e-152\n",
      "Precision: 2.0950697665885954e-153\n",
      "Precision: 2.3278552962095505e-154\n",
      "Precision: 2.5865058846772782e-155\n",
      "Precision: 2.873895427419198e-156\n",
      "Precision: 3.193217141576886e-157\n",
      "Precision: 3.54801904619654e-158\n",
      "Precision: 3.9422433846628226e-159\n",
      "Precision: 4.380270427403136e-160\n",
      "Precision: 4.8669671415590405e-161\n",
      "Precision: 5.407741268398934e-162\n",
      "Precision: 6.008601409332149e-163\n",
      "Precision: 6.6762237881468325e-164\n",
      "Precision: 7.418026431274259e-165\n",
      "Precision: 8.242251590304732e-166\n",
      "Precision: 9.158057322560813e-167\n",
      "Precision: 1.0175619247289793e-167\n",
      "Precision: 1.130624360809977e-168\n",
      "Precision: 1.2562492897888634e-169\n",
      "Precision: 1.3958325442098483e-170\n",
      "Precision: 1.5509250491220536e-171\n",
      "Precision: 1.7232500545800595e-172\n",
      "Precision: 1.9147222828667327e-173\n",
      "Precision: 2.1274692031852587e-174\n",
      "Precision: 2.363854670205843e-175\n",
      "Precision: 2.6265051891176034e-176\n",
      "Precision: 2.9183390990195595e-177\n",
      "Precision: 3.242598998910622e-178\n",
      "Precision: 3.602887776567358e-179\n",
      "Precision: 4.003208640630398e-180\n",
      "Precision: 4.448009600700442e-181\n",
      "Precision: 4.942232889667158e-182\n",
      "Precision: 5.491369877407954e-183\n",
      "Precision: 6.101522086008837e-184\n",
      "Precision: 6.779468984454264e-185\n",
      "Precision: 7.532743316060293e-186\n",
      "Precision: 8.369714795622548e-187\n",
      "Precision: 9.299683106247276e-188\n",
      "Precision: 1.0332981229163639e-188\n",
      "Precision: 1.1481090254626265e-189\n",
      "Precision: 1.275676694958474e-190\n",
      "Precision: 1.41741854995386e-191\n",
      "Precision: 1.5749094999487333e-192\n",
      "Precision: 1.7498994443874815e-193\n",
      "Precision: 1.9443327159860904e-194\n",
      "Precision: 2.160369684428989e-195\n",
      "Precision: 2.4004107604766547e-196\n",
      "Precision: 2.667123067196283e-197\n",
      "Precision: 2.9634700746625367e-198\n",
      "Precision: 3.2927445274028187e-199\n",
      "Precision: 3.6586050304475763e-200\n",
      "Precision: 4.065116700497307e-201\n",
      "Precision: 4.516796333885897e-202\n",
      "Precision: 5.018662593206552e-203\n",
      "Precision: 5.576291770229502e-204\n",
      "Precision: 6.195879744699447e-205\n",
      "Precision: 6.88431082744383e-206\n",
      "Precision: 7.649234252715366e-207\n",
      "Precision: 8.49914916968374e-208\n",
      "Precision: 9.443499077426378e-209\n",
      "Precision: 1.0492776752695976e-209\n",
      "Precision: 1.1658640836328863e-210\n",
      "Precision: 1.2954045373698738e-211\n",
      "Precision: 1.4393383748554154e-212\n",
      "Precision: 1.5992648609504617e-213\n",
      "Precision: 1.776960956611624e-214\n",
      "Precision: 1.9744010629018045e-215\n",
      "Precision: 2.1937789587797828e-216\n",
      "Precision: 2.437532176421981e-217\n",
      "Precision: 2.708369084913312e-218\n",
      "Precision: 3.009298983237013e-219\n",
      "Precision: 3.3436655369300148e-220\n",
      "Precision: 3.7151839299222385e-221\n",
      "Precision: 4.127982144358043e-222\n",
      "Precision: 4.586646827064493e-223\n",
      "Precision: 5.096274252293881e-224\n",
      "Precision: 5.662526946993201e-225\n",
      "Precision: 6.291696607770223e-226\n",
      "Precision: 6.990774008633581e-227\n",
      "Precision: 7.767526676259535e-228\n",
      "Precision: 8.630585195843928e-229\n",
      "Precision: 9.589539106493253e-230\n",
      "Precision: 1.065504345165917e-230\n",
      "Precision: 1.1838937168510189e-231\n",
      "Precision: 1.3154374631677987e-232\n",
      "Precision: 1.461597181297554e-233\n",
      "Precision: 1.6239968681083932e-234\n",
      "Precision: 1.8044409645648814e-235\n",
      "Precision: 2.0049344050720902e-236\n",
      "Precision: 2.2277048945245448e-237\n",
      "Precision: 2.4752276605828277e-238\n",
      "Precision: 2.750252956203142e-239\n",
      "Precision: 3.055836618003491e-240\n",
      "Precision: 3.3953740200038785e-241\n",
      "Precision: 3.772637800004309e-242\n",
      "Precision: 4.191819777782566e-243\n",
      "Precision: 4.6575775308695176e-244\n",
      "Precision: 5.175086145410575e-245\n",
      "Precision: 5.750095717122861e-246\n",
      "Precision: 6.388995241247624e-247\n",
      "Precision: 7.098883601386249e-248\n",
      "Precision: 7.887648445984721e-249\n",
      "Precision: 8.764053828871912e-250\n",
      "Precision: 9.737837587635458e-251\n",
      "Precision: 1.0819819541817174e-251\n",
      "Precision: 1.2022021713130193e-252\n",
      "Precision: 1.3357801903477992e-253\n",
      "Precision: 1.4842002114975546e-254\n",
      "Precision: 1.649111346108394e-255\n",
      "Precision: 1.8323459401204377e-256\n",
      "Precision: 2.035939933467153e-257\n",
      "Precision: 2.26215548163017e-258\n",
      "Precision: 2.5135060907001887e-259\n",
      "Precision: 2.792784545222432e-260\n",
      "Precision: 3.1030939391360357e-261\n",
      "Precision: 3.447882154595595e-262\n",
      "Precision: 3.830980171772884e-263\n",
      "Precision: 4.256644635303204e-264\n",
      "Precision: 4.729605150336893e-265\n",
      "Precision: 5.255116833707659e-266\n",
      "Precision: 5.839018704119621e-267\n",
      "Precision: 6.487798560132913e-268\n",
      "Precision: 7.208665066814347e-269\n",
      "Precision: 8.00962785201594e-270\n",
      "Precision: 8.899586502239934e-271\n",
      "Precision: 9.888429446933259e-272\n",
      "Precision: 1.0987143829925844e-272\n",
      "Precision: 1.2207937588806493e-273\n",
      "Precision: 1.3564375098673882e-274\n",
      "Precision: 1.5071527887415426e-275\n",
      "Precision: 1.674614209712825e-276\n",
      "Precision: 1.8606824552364721e-277\n",
      "Precision: 2.0674249502627467e-278\n",
      "Precision: 2.297138833625274e-279\n",
      "Precision: 2.55237648180586e-280\n",
      "Precision: 2.835973868673178e-281\n",
      "Precision: 3.1510820763035312e-282\n",
      "Precision: 3.5012023070039234e-283\n",
      "Precision: 3.890224785559915e-284\n",
      "Precision: 4.322471983955461e-285\n",
      "Precision: 4.8027466488394015e-286\n",
      "Precision: 5.336385165377113e-287\n",
      "Precision: 5.929316850419014e-288\n",
      "Precision: 6.588129833798905e-289\n",
      "Precision: 7.320144259776561e-290\n",
      "Precision: 8.133493621973957e-291\n",
      "Precision: 9.037215135526619e-292\n",
      "Precision: 1.004135015058513e-292\n",
      "Precision: 1.1157055722872368e-293\n",
      "Precision: 1.2396728580969298e-294\n",
      "Precision: 1.3774142867743663e-295\n",
      "Precision: 1.5304603186381847e-296\n",
      "Precision: 1.7005114651535384e-297\n",
      "Precision: 1.8894571835039315e-298\n",
      "Precision: 2.0993968705599237e-299\n",
      "Precision: 2.3326631895110265e-300\n",
      "Precision: 2.591847988345585e-301\n",
      "Precision: 2.879831098161761e-302\n",
      "Precision: 3.1998123312908456e-303\n",
      "Precision: 3.555347034767606e-304\n",
      "Precision: 3.9503855941862293e-305\n",
      "Precision: 4.389317326873588e-306\n",
      "Precision: 4.877019252081764e-307\n",
      "Precision: 5.418910280090849e-308\n",
      "Precision: 6.021011422323166e-309\n",
      "Precision: 6.6900126914702e-310\n",
      "Precision: 7.433347434967e-311\n",
      "Precision: 8.25927492774e-312\n",
      "Precision: 9.17697214195e-313\n",
      "Precision: 1.0196635713e-313\n",
      "Precision: 1.1329595237e-314\n",
      "Precision: 1.258843915e-315\n",
      "Precision: 1.39871546e-316\n",
      "Precision: 1.5541284e-317\n",
      "Precision: 1.72681e-318\n",
      "Precision: 1.91865e-319\n",
      "Precision: 2.132e-320\n",
      "Precision: 2.367e-321\n",
      "Precision: 2.6e-322\n",
      "Precision: 3e-323\n",
      "Precision: 5e-324\n"
     ]
    }
   ],
   "source": [
    "epsilon, ratio = 1.0, 9.0 \n",
    "while (epsilon):\n",
    "    print('Precision:',epsilon) \n",
    "    epsilon /= ratio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last line indicates that the precision that standard Python assumes on my machine is 10^{−324}. This is indeed associated to a 64-bit float number. One can verify this by setting a smaller precision in the first line of this sequence. Setting initially:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Precision: 1.0\n",
      "Precision: 0.11111111\n",
      "Precision: 0.012345679\n",
      "Precision: 0.0013717421\n",
      "Precision: 0.00015241579\n",
      "Precision: 1.6935088e-05\n",
      "Precision: 1.8816764e-06\n",
      "Precision: 2.0907515e-07\n",
      "Precision: 2.3230573e-08\n",
      "Precision: 2.5811748e-09\n",
      "Precision: 2.867972e-10\n",
      "Precision: 3.1866357e-11\n",
      "Precision: 3.5407064e-12\n",
      "Precision: 3.934118e-13\n",
      "Precision: 4.3712423e-14\n",
      "Precision: 4.856936e-15\n",
      "Precision: 5.3965957e-16\n",
      "Precision: 5.996218e-17\n",
      "Precision: 6.662464e-18\n",
      "Precision: 7.402738e-19\n",
      "Precision: 8.225264e-20\n",
      "Precision: 9.1391825e-21\n",
      "Precision: 1.0154648e-21\n",
      "Precision: 1.1282942e-22\n",
      "Precision: 1.2536603e-23\n",
      "Precision: 1.3929559e-24\n",
      "Precision: 1.5477288e-25\n",
      "Precision: 1.7196987e-26\n",
      "Precision: 1.9107763e-27\n",
      "Precision: 2.1230847e-28\n",
      "Precision: 2.3589831e-29\n",
      "Precision: 2.6210923e-30\n",
      "Precision: 2.9123249e-31\n",
      "Precision: 3.2359164e-32\n",
      "Precision: 3.5954625e-33\n",
      "Precision: 3.9949584e-34\n",
      "Precision: 4.4388426e-35\n",
      "Precision: 4.9320472e-36\n",
      "Precision: 5.4800523e-37\n",
      "Precision: 6.088947e-38\n",
      "Precision: 6.765497e-39\n",
      "Precision: 7.51722e-40\n",
      "Precision: 8.3524e-41\n",
      "Precision: 9.281e-42\n",
      "Precision: 1.031e-42\n",
      "Precision: 1.15e-43\n",
      "Precision: 1.3e-44\n",
      "Precision: 1e-45\n"
     ]
    }
   ],
   "source": [
    "epsilon, ratio =  np.float16(1), np.float32(9)\n",
    "while (epsilon):\n",
    "    print('Precision:',epsilon) \n",
    "    epsilon /= ratio\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly to float, there exist integers with several precisions: np.int32, np.int64, np.int128. The type of every NumPy variable can be    retrieved the dtype, e.g., eps.dtype.\n",
    "\n",
    "It is very important to keep in mind the importance of precision. All the mathematical routines implemented in Python, as well as in any other language, follow the standards determined by IEEE which imply that they are all corrected to the last digit. Therefore, when we will model a physical system we can be sure that errors will not propagate quickly to lower digits with every algorithm that does not explicitly remove this information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vectorization \n",
    "Vectorization is a concept on programming in which the programmer tries to make the minimal use of loop by replacing it vector/matrix operations to make the program run faster. It is mostly used in the programs that are slow in looping like in python. \n",
    "Furthermore, the vectorization makes the program look simple and also makes the parallization easy to implement. \n",
    "\n",
    "let us take an example of matrix multiplication one using loop, and another by the technique of vectorization. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Example: matrix multiplication using loop (non-vectorized algorithm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[19 22]\n",
      " [43 50]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "a=np.array([[1,2],[3,4]])\n",
    "b=np.array([[5,6],[7,8]])\n",
    "c=np.zeros((2,2),int)\n",
    "\n",
    "for i in np.arange(2):\n",
    "    for j in np.arange(2):\n",
    "        sum=0;\n",
    "        for k in np.arange(2):\n",
    "            sum+=a[i, k] * b[k ,j];\n",
    "        c[i,j]=sum\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Example: Matrix multiplication using vectorized algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[19 22]\n",
      " [43 50]]\n"
     ]
    }
   ],
   "source": [
    "a=np.array([[1,2],[3,4]])\n",
    "b=np.array([[5,6],[7,8]])\n",
    "c=np.matmul(a,b)\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ndarrays\n",
    "\n",
    "Let us start looking at the vectorized arrays. The core of NumPy is the object called ndarray, a storage of large quantities of data that allows to operate fast and flexible operations on it. Let’s look at some basic examples using iPython again and create a matrix (a two-dimensional array) of two rows and three columns:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dtype('int8')"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "firstArr = np.array([[4,6,3,4],[2,3,6,0]], dtype=np.int8)\n",
    "firstArr.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "firstArr.ndim # number of dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 4)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "firstArr.shape # array shape as (rows,columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dtype('int8')"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "firstArr.dtype # Out[10]: dtype(’int64’)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we find is that NumPy, if given a certain set of data, automatically decides a type (generally integer or float) and defines an array of the right shape where to store them. Given a NumPy array, it is always possible to obtain its shape and data type (dtype). The main difference between ndarrays and Python data types is that the elements of ndarrays are all of a predefined and homogeneous type, which is one of the reasons for the speed of its calculations.\n",
    "\n",
    "In general, one can set the data type at the moment of creation of the array, which is a particularly useful function when one has to handle very large models and datasets and needs to be in control of the size of the occupied memory. Data in NumPy are either int: integer, uint: unsigned integer (from 0), and float, a real number. int and uint can have sizes of 8, 16, 32, and 64 bits, while float of 16, 32, and 64 bits. If not otherwise set, automatically NumPy will set data type size to 64 bits.\n",
    "\n",
    "Most commonly used generators of ndarrays are arange, zeros and ones. Let us look at some examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,\n",
       "       17, 18, 19], dtype=int32)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "secondArr = np.arange((20),dtype=np.int32)\n",
    "secondArr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dtype('int32')"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "secondArr.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0.])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "In [14]: secondArr = np.zeros(20)\n",
    "In [15]: secondArr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dtype('float64')"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "In [16]: secondArr.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1.])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "In [17]: secondArr = np.ones(20)\n",
    "In [18]: secondArr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dtype('float64')"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "In [19]: secondArr.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Arange is the analogue of range in Python, but it creates a NumPy array of 64 bits integers. Zeros and ones instead create arrays of 64 bits float. Ones can also be created by zeros using the broadcasting properties of NumPy arrays, just by the instruction arr=np.zeros(20)+1.0. If we wish to create arrays with different types, for example due to operational reasons or of memory size, dtypes can be explicitly specified, \n",
    "e.g.,arr = np.arange(20, dtype=’float32’) or arr = np.ones(20, dtype=’int8’).\n",
    "\n",
    "Let us now check how using NumPy allows to speed the calculations of the prior section. The first temptation might be to simply to replace range with np.arange."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def addCountNP(max):\n",
    "    a=0\n",
    "    for i in np.arange(max):\n",
    "        a=a+i\n",
    "    #print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.24 s ± 103 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit addCountNP(10000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "561 ms ± 56.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit addCount(10000000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happened here is that although the arrays were ndarrays, the operation was done as for a standard Python list. The trick to speed calculations on ndarrays is to use the broadcast vectorized version of each operation. Let us look at how to add two ndarrays of one million of integers either with the standard Python loop and exploiting the NumPy broadcasting capabilities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "def addArray(a,b):\n",
    "    c=np.zeros(a.size)\n",
    "    for i in np.arange(a.size):\n",
    "        c[i]=a[i]+b[i]\n",
    "    return(c)\n",
    "\n",
    "a=np.arange(1000000)\n",
    "b=np.arange(1000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "500 ms ± 21.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit c=addArray(a,b) #standard python "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.89 ms ± 18.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit c=a+b #NumPy arrays broadcasting "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A gain of two orders of magnitude! To use the broadcasting feature of NumPy makes Python’s speed comparable to compiled codes such as C, but with the obvious huge gain in terms for code development, testing, and readability. Broadcasting is clearly not limited to the addition operation, but it works as well with all the other arithmetic operations such as c=a*b; c=a/b; c=a-b,as well as scalar-array operations like \n",
    "```python \n",
    "c=1/a; \n",
    "c=a**0.5\n",
    "```\n",
    "Besides operations on two arrays, there are numerous unary (and binary universal functions that can be used with NumPy arrays. Among unary functions, much used ones are \n",
    "```python \n",
    "np.abs()\n",
    "np.sqrt()\n",
    "np.exp()\n",
    "np.log()\n",
    "np.sign()\n",
    "```\n",
    "and all the trigonometric functions. Binary functions are \n",
    "```python \n",
    "np.add()\n",
    "np.multiply()\n",
    "np.power()\n",
    "np.maximum()\n",
    "np.mod().\n",
    "```\n",
    "Let us go back now to the initial problem. We wanted to sum all the elements of a large array (100 millions numbers). This is neither an unary or binary operation, because it is a function that projects an array into one number, the sum. The most common among these operations are already efficiently implemented in NumPy. For example, the sum is immediately obtained by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "49.1 ms ± 210 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "%timeit np.arange(10000000).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.85804383 0.75349206 0.10779009 ... 0.10493196 0.52583065 0.95307547]\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "x=np.random.rand(1001)\n",
    "print(x)\n",
    "g=1\n",
    "for y in x:\n",
    "     g = g*y\n",
    "print (g)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These operation have been in fact written in Cython, which is a C compiled operation. We will see later in this chapter how such operations, when not already implemented in NumPy, can be easily created using Cython."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Indexing and Slicing\n",
    "\n",
    "Most languages today are compiled. One of the great features of compilers is that when we have to perform a set of repetitive operations on a large set of data the compiler will take care to optimize these iterations in an exceptional way, even if we wrote a very confusing code. This is not the case with an interpreted software like Python.\n",
    "\n",
    "If we want our Python code to run fast, we have to organize the sequence of operation in a smart way. In particular a very common case like calling an **if** command inside a for loop can make Python very slow. We have seen above how looping makes Python slow. One of the secrets that makes Python fast is to use smart indexing. To use indexing in practice in this case means to vectorize the command **if** and then to use these indexes to selectively operate on arrays. This might seem a terrible setback for Python programmers, however the great advantage of forcing the developer to do this operation is that the program will be written in a completely vectorized way from the start, therefore ready to be parallelized.\n",
    "\n",
    "Indexing and slicing in NumPy is a long topic, whose full coverage goes beyond the scopes of this book. To gain a full understanding of the possibilities offered by NumPy I recommend to follow one of the online free tutorials (e.g., https://docs.scipy.org/doc/numpy-dev/user/basics.indexing.html). I will cover here some main the features that we will use more often in this book, and explain in depth some important details on the memory management associated to ndarrays.\n",
    "\n",
    "Let us now dive into the main features that interest us. Given an array, e.g., **arr**, we can access the element n with square brackets **arr[n]** and we can slice the array extracting the elements between n and m with the command **arr[n:m]**, which will return m − n elements (comprising **arr[n]** and excluding **arr[m]**. It is very important to understand that even if we associate a name to the slice of arr, this is only a view, not a new array implying that the m − n elements are not copied in a new allocated chunk of memory. This means that by changing the sliced array you will change the initial data. Let us look at an example to clearly understand how this works:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,\n",
       "       17, 18, 19])"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "secondArr=np.arange(20)\n",
    "secondArr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4, 5, 6, 7, 8, 9])"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sliceArr = secondArr[4:10]\n",
    "superslicearray = sliceArr\n",
    "sliceArr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([     4,      5,      6, 100000,      8,      9])"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "superslicearray[3]=100000\n",
    "sliceArr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([     0,      1,      2,      3,      4,      5,      6, 100000,\n",
       "            8,      9,     10,     11,     12,     13,     14,     15,\n",
       "           16,     17,     18,     19])"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "secondArr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Who is familiar with other languages might be surprised by this behavior, and believe that when defining sliceArr Python should have copied the subset of second **Arr** data into a new array. The point is that NumPy has been designed to deal with very large datasets or numerical models, therefore it uses a policy of minimization of memory usage. For this reason, if one does not force NumPy to create a copy of the sliced data, it will simply generate a view of the already existing array. NumPy can be forced to create a separate new array by adding **.copy()** to the slice, e.g., in our case **sliceArr = secondArr[4:10].copy()**.\n",
    "\n",
    "Slicing is very flexible and allows omitting, for example, the first index (e.g., **arr[:10]**) or the last index (e.g., **arr[10:]**) of a slice, which implies that the slice reaches the end of the array. It is also possible to slice the array every k elements by indicating a third parameter in the slice (e.g., **arr[4:12:3]**), and also omitting the other parameters (e.g., **arr[::3]**). Negative indexes are also admitted, which means that counting starts from the last element, backward, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10, 11, 12, 13, 14])"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "secondArr[10:-5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### N-Dimensional Indexing\n",
    "Indexing and slicing in more dimensions are just a recursive repetition of one-dimensional operations. In 2D, for example, arrays can be indexed either as **arr2d[n][m]** or, with the same effect, as **arr2d[n,m]**. The first index refers to the inner array, the second index to the outer array. For example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2,  3],\n",
       "       [ 4,  5,  6,  7],\n",
       "       [ 8,  9, 10, 11],\n",
       "       [12, 13, 14, 15],\n",
       "       [16, 17, 18, 19]])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arr2d = np.arange(20)\n",
    "sliceArr = arr2d[4:10]\n",
    "arr2d=arr2d.reshape((5,4))\n",
    "arr2d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4, 5, 6, 7, 8, 9])"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sliceArr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  0,   1,   2,   3],\n",
       "       [  4,   5,   6, 100],\n",
       "       [  8,   9,  10,  11],\n",
       "       [ 12,  13,  14,  15],\n",
       "       [ 16,  17,  18,  19]])"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sliceArr[3]=100\n",
    "arr2d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "secondArr[4*2+3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is very important to remember that the reshape command does not create a copy of the data, as well as slicing. The data remain stored in the memory as a one-dimensional array, regardless of the dimensions and shape of the array. The way in which this is done and how NumPy broadcasting operations remain extremely efficient is explained later, where strides are illustrated. Also a book that analyzes in greater detail how Python can deal with data is McKinney,W.(2012) Python for data analysis. Newton: O’Reilly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Boolean Indexing\n",
    "\n",
    "Boolean indexing is the fastest and efficient way to select, access and operate on subsets of a NumPy arrays. Let us for example create a random 5 × 5 matrix (array) and select only the positive elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.72584876,  0.41539641,  0.25705093, -0.40139017, -2.28025703],\n",
       "       [-1.47471148,  0.5574254 ,  1.32576064, -0.75029982, -1.7270752 ],\n",
       "       [ 0.28816588, -0.19908194,  1.08603933, -0.76006449,  0.27014051],\n",
       "       [-0.42589037,  0.92234231,  0.0050546 , -1.53421268,  1.5221565 ],\n",
       "       [-0.0669165 , -0.17803476,  0.73886003, -0.22918035,  0.03110576]])"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arr2d = np.random.randn(5,5)\n",
    "arr2d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[False,  True,  True, False, False],\n",
       "       [False,  True,  True, False, False],\n",
       "       [ True, False,  True, False,  True],\n",
       "       [False,  True,  True, False,  True],\n",
       "       [False, False,  True, False,  True]])"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arr2d>0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.41539641, 0.25705093, 0.5574254 , 1.32576064, 0.28816588,\n",
       "       1.08603933, 0.27014051, 0.92234231, 0.0050546 , 1.5221565 ,\n",
       "       0.73886003, 0.03110576])"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arr2d[arr2d>0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we did was to create a 5×5 array whose values were True when the elements were positive, and False in the opposite case. It was then possible to select only those elements. It is important to notice that the extracted elements have lost their 5 × 5 structure. This again results from the fact that Python stores every array, regardless of its shape, as a 1D array.\n",
    "\n",
    "This technique is normally used to operate on a certain subset of an array. For example if we desire to set to 0 all the negative elements of the above array, we can do it with one instruction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.        , 0.41539641, 0.25705093, 0.        , 0.        ],\n",
       "       [0.        , 0.5574254 , 1.32576064, 0.        , 0.        ],\n",
       "       [0.28816588, 0.        , 1.08603933, 0.        , 0.27014051],\n",
       "       [0.        , 0.92234231, 0.0050546 , 0.        , 1.5221565 ],\n",
       "       [0.        , 0.        , 0.73886003, 0.        , 0.03110576]])"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arr2d[arr2d<0]=0\n",
    "arr2d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important to emphasize that besides compactness and elegance, this technique guarantees an enormous gain in speed. To quantify this advantage, we can compare the time necessary for a standard Python for loop with the Boolean indexing tool described above for a large 2000 × 2000 random array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr2d=np.random.randn(2000,2000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.32 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "def setNegativeValuesToZero(n,m,a): \n",
    "    for i in np.arange(n):\n",
    "        for j in np.arange(m):\n",
    "            if a[i,j]<0:\n",
    "                a[i,j]=0\n",
    "%timeit -n1 -r1 setNegativeValuesToZero(2000,2000,arr2d)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "37.2 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "arr2d = np.random.randn(2000,2000)\n",
    "%timeit -n1 -r1 arr2d[arr2d<0]=0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we observe a gain of 2 orders of magnitude by using the NumPy indexing compared to standard Python loops. The message is now loud and clear: never use Python loops for large datasets and large numerical models, but always employ Python indexing/slicing features to achieve compiled code performance for speed and memory management. When not possible, rely on Cython, that we will study soon."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transposing and Axis Rotation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We already encountered the function reshape, which allows to project a 1D array as a generic n-dimensional array. One operation that cannot be achieved by reshaping is however transposing, which corresponds to swapping axis. This is a very important function when operating on 2D and 3D numerical modeling, but it is also simply essential in order to calculate the inner product. For example, to calculate Xtranspose * X , one writes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.89860875, -0.9775276 ],\n",
       "       [ 0.51299635, -0.14375967],\n",
       "       [ 3.21337996, -0.14537557],\n",
       "       [ 0.35133024,  0.83534005],\n",
       "       [-0.13308273, -0.45844327]])"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n=10\n",
    "arr2d = np.random.randn(n)\n",
    "arr2d\n",
    "arr2d=arr2d.reshape(5,int(n/5))\n",
    "arr2d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[11.53761769,  0.69201084],\n",
       "       [ 0.69201084,  1.90532435]])"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arr2d.T\n",
    "np.dot(arr2d.T,arr2d)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11.537617685403443"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(arr2d[:,0])**2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transposing, however, is not limited to swapping x and y axis. One can rotate axis as well. For example in 3D:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 0,  1,  2,  3],\n",
       "        [ 4,  5,  6,  7],\n",
       "        [ 8,  9, 10, 11]],\n",
       "\n",
       "       [[12, 13, 14, 15],\n",
       "        [16, 17, 18, 19],\n",
       "        [20, 21, 22, 23]]])"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arr3d = np.arange(2*3*4).reshape(2,3,4)\n",
    "arr3d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 3, 4)\n",
      "(4, 2, 3)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[[ 0,  4,  8],\n",
       "        [12, 16, 20]],\n",
       "\n",
       "       [[ 1,  5,  9],\n",
       "        [13, 17, 21]],\n",
       "\n",
       "       [[ 2,  6, 10],\n",
       "        [14, 18, 22]],\n",
       "\n",
       "       [[ 3,  7, 11],\n",
       "        [15, 19, 23]]])"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arr3d2 = arr3d.transpose((2,0,1)).copy()\n",
    "print(arr3d.shape)\n",
    "print(arr3d2.shape)\n",
    "arr3d[1,0,0]=1000\n",
    "arr3d2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transposing, as well as reshaping, is just a specific view of the entire array, therefore transposing will not create a new set of data, and when modifying the transposed of an array, one modifies the original array as well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Strides\n",
    "\n",
    "The technical details given here help to understand how NumPy deals with very large arrays minimizing memory occupation and maximizing access speed.\n",
    "\n",
    "It is important to emphasize that although an array can have many dimensions, the memory in our computer is structured in a purely sequential way, therefore ultimately the real structure of the array is one-dimensional, and every n-dimensional array will only be a specific view on that 1D array.\n",
    "\n",
    "NumPy is extraordinary powerful in managing n-dimensional arrays, different views, complex subsets, for example\n",
    "\n",
    "```python\n",
    "arr[arr*arr<1]. ```\n",
    "\n",
    "This power is based on the ability to directly address chunks of data by striding across the\n",
    "memory and zooming into small blocks of memory. For example, let us take a random 100 × 100 × 100 array of np.float64 data, equivalent to a 1D array of one million elements, occupying 8 bytes each. A block within the array can be extracted in microseconds, almost regardless to the size of the block:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[-8.64944929e-01, -5.90533518e-01, -4.10147374e-01, ...,\n",
       "         -5.69645473e-01, -1.01094307e+00,  1.35032788e+00],\n",
       "        [-1.19863891e-01, -1.38975586e+00, -3.94702567e-01, ...,\n",
       "          7.75620090e-01, -5.25646441e-01, -8.48567348e-01],\n",
       "        [-9.75255360e-01,  1.46582911e+00,  7.83651675e-01, ...,\n",
       "          2.05455403e-01,  6.66553318e-01, -2.51603929e-01],\n",
       "        ...,\n",
       "        [ 1.21013662e+00,  7.42753705e-01,  1.27906901e-01, ...,\n",
       "          6.89053956e-01, -5.49866401e-01,  1.23940416e+00],\n",
       "        [-7.23377423e-01,  1.61014585e-01,  9.85509721e-01, ...,\n",
       "          8.28020813e-01,  1.22753843e-01, -9.62115768e-01],\n",
       "        [-6.28193551e-01, -8.78986810e-01,  2.14976804e+00, ...,\n",
       "         -3.91930555e-01,  2.24444696e-01, -7.55916935e-01]],\n",
       "\n",
       "       [[ 7.18685731e-01,  1.69048723e+00, -1.37982218e+00, ...,\n",
       "         -1.62541779e+00, -2.09312832e+00,  9.75288858e-01],\n",
       "        [-1.39085542e+00, -9.40017465e-01,  8.79226632e-02, ...,\n",
       "          2.67830570e-01, -1.37103056e+00,  1.34539410e+00],\n",
       "        [-8.05692961e-01,  9.92133096e-01,  2.53570469e-01, ...,\n",
       "          1.20540694e-03, -5.66631476e-01,  1.26752800e+00],\n",
       "        ...,\n",
       "        [-9.04328006e-01,  9.45566196e-01, -7.49008779e-01, ...,\n",
       "          7.02133814e-01,  1.25426797e-01, -1.73760643e+00],\n",
       "        [ 1.96086003e-01,  1.03482451e+00,  3.92222191e-01, ...,\n",
       "          1.40629789e-01,  7.14684152e-01, -7.77967101e-02],\n",
       "        [ 7.46873033e-01,  1.65751508e+00, -1.34037196e+00, ...,\n",
       "         -5.67947703e-01,  2.74065342e-01, -5.27900097e-01]],\n",
       "\n",
       "       [[-8.30401835e-02,  1.11893653e+00, -9.23342096e-01, ...,\n",
       "          1.35323664e-01, -4.72578940e-01, -2.15078122e+00],\n",
       "        [-2.33825943e-02, -1.96483274e-01,  6.40807440e-01, ...,\n",
       "         -7.55958731e-01,  1.50912479e-01,  1.48147851e+00],\n",
       "        [-3.62475762e-01, -1.70104357e+00,  8.29679705e-02, ...,\n",
       "         -9.84998443e-01,  4.75440215e-01, -1.11179275e+00],\n",
       "        ...,\n",
       "        [-4.07264555e-02,  2.45605955e+00,  8.59581849e-01, ...,\n",
       "         -3.06804155e-02,  1.36347538e+00, -3.13360923e-01],\n",
       "        [-1.43686781e-01, -1.35275598e+00, -1.82640419e-01, ...,\n",
       "          6.76462242e-01,  2.08187796e-01,  3.99098625e-01],\n",
       "        [ 7.18225203e-01,  6.80858682e-02,  5.95031200e-01, ...,\n",
       "          3.04956512e-01, -4.88061504e-01, -6.60300893e-01]],\n",
       "\n",
       "       ...,\n",
       "\n",
       "       [[ 5.42652667e-01, -6.71470634e-01, -5.27070515e-01, ...,\n",
       "         -1.19118983e+00,  3.69704230e-01, -1.14344367e+00],\n",
       "        [ 1.23374239e+00,  4.55598148e-01,  2.46086822e-01, ...,\n",
       "          1.30071764e+00,  1.31839308e+00, -8.53235685e-01],\n",
       "        [ 6.21613189e-01,  6.66802312e-01,  1.24635727e+00, ...,\n",
       "         -6.66320397e-02, -4.17874221e-02,  3.56035774e-01],\n",
       "        ...,\n",
       "        [-1.50795313e-01,  7.93151739e-01, -1.30832896e+00, ...,\n",
       "          4.42285409e-01,  1.87956105e+00,  9.75540037e-01],\n",
       "        [ 8.73391897e-01,  8.18182875e-02,  5.77571158e-01, ...,\n",
       "          4.07209406e-01,  1.47091338e-01, -2.59011691e-01],\n",
       "        [-4.52699916e-01, -1.00775395e-01,  6.73023785e-01, ...,\n",
       "         -6.61058219e-01, -1.51386399e+00,  3.54500968e-01]],\n",
       "\n",
       "       [[ 7.51441068e-01, -4.25403879e-01, -9.23084821e-01, ...,\n",
       "         -2.09221410e-01,  1.60103616e-01,  2.51253356e-01],\n",
       "        [-1.58126029e+00,  1.99731880e+00, -2.75444553e-01, ...,\n",
       "         -3.71096717e-01, -4.52165506e-01, -5.51891383e-01],\n",
       "        [-4.35683423e-01, -4.54131615e-01, -7.61208255e-01, ...,\n",
       "         -2.17537162e+00,  1.11470825e+00,  2.94152945e-01],\n",
       "        ...,\n",
       "        [ 1.34100392e+00,  7.09796277e-01, -8.00938896e-01, ...,\n",
       "          1.27372491e+00,  2.65348292e-01,  4.62265025e-01],\n",
       "        [ 1.20476199e+00,  1.09746473e+00,  8.31309073e-01, ...,\n",
       "          4.91318671e-02,  1.99031597e+00, -1.55488501e+00],\n",
       "        [ 1.18800780e-01, -9.56134530e-01, -1.58216784e+00, ...,\n",
       "         -1.00567365e+00, -1.41093068e+00, -3.58525620e-01]],\n",
       "\n",
       "       [[-1.04642793e+00,  3.67117769e-01,  1.25549901e+00, ...,\n",
       "         -1.36335138e+00, -3.50241545e-01,  9.89861649e-01],\n",
       "        [-2.46137730e-01, -2.46412586e-01, -2.12807798e-01, ...,\n",
       "         -9.94574871e-01,  6.45469728e-01, -3.42567753e-01],\n",
       "        [-2.39301528e-01,  1.24980547e+00, -8.34541727e-02, ...,\n",
       "         -3.98446732e-01,  7.88851294e-01,  1.03200797e+00],\n",
       "        ...,\n",
       "        [ 7.42649877e-01, -3.17207105e+00, -2.02280992e+00, ...,\n",
       "         -9.82805931e-01, -1.87647812e+00,  8.12692188e-01],\n",
       "        [-1.35803138e-01, -9.34797178e-04,  1.20025556e+00, ...,\n",
       "          1.01406829e-01,  1.68605787e-01, -2.20869804e+00],\n",
       "        [ 5.31017358e-01,  1.80476866e-02, -1.05299813e-01, ...,\n",
       "         -7.72500415e-01,  1.33971557e-01,  5.55578609e-01]]])"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.randn(100,100,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14.7 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n1 -r1 newArray=arr3d[45:55,45:55,45:55].copy() #small block"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "23.9 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n1 -r1 newArray2=arr3d[15:85,15:85,15:85].copy() #large block"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we operate on the subset of the array however the operational time increases with the size of the block, although not in a proportional way. For example, to set the values of the subset to zero we need about 20 times more time for a subset that is over 300 times greater:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.81 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n1 -r1 arr3d[45:55,45:55,45:55]=0 #small block"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n1 -r1 arr3d[15:85,15:85,15:85]=0 #large block"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NumPy accesses n-dimensional arrays and its slices as fast as one-dimensional arrays. To understand how this is achieved, let us consider, for example an array a of 32000 integers and then create a reshaped three-dimensional 20 × 40 × 40 version assigned to b. In our case a is composed by 64-bit (8 bytes) integers, which means that one needs to proceed 8 bytes forward to access the next element along the first axis. The size of every element can be also in general assessed with a.itemsize:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "a=np.arange(32000)\n",
    "b=a.reshape(20,40,40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8 8\n"
     ]
    }
   ],
   "source": [
    "print(a.itemsize,b.itemsize)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This will be stored in memory in a buffer that contains 32000 ascending integers from 0 to 31999. As we have seen before a and b are stored in the same memory block. The way in which NumPy differentiates how to operate on them is by characterizing them by their different strides. Strides are tuples of bytes to step in each dimension when traversing an array. In practice they are the offset in bytes between an element and the neighboring one in every direction. Strides are shown explicitly using the instruction arr.strides:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(20, 40, 40) (12800, 320, 8)\n",
      "320 12800\n"
     ]
    }
   ],
   "source": [
    "print(b.shape,b.strides)\n",
    "print(8*40,8*40*40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8,)"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a.strides"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last number, 8, of the strides refers to the size of each element in the array, i.e., it is always the itemsize. The other numbers refer to the number of bytes forward necessary to access the next element along the other axes. The total number of bytes is called offset and is calculated as offset = sum(indexes * a.strides). For example, for the initial array a[n] every element is accessed at the position \n",
    "```python\n",
    "offset = n*a.strides[0]\n",
    "```\n",
    ", while b[i,j,k] is accessed calculating its associated memory offset as \n",
    "```python\n",
    "offset = i*b.strides[0]+j*b.strides[1]+k*b.strides[2]\n",
    "```\n",
    "in bytes. To obtain the element location one has to divide by b.itemsize. For example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3324"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b[2,3,4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3324.0"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(2*b.strides[0]+3*b.strides[1]+4*b.strides[2])/b.itemsize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Vector Product\n",
    "\n",
    "The inner product, np.dot() in NumPy, is the most common operation between arrays. It is equivalent to the inner product between 1-D arrays or to matrix multiplication between 2D arrays. It is defined on n-dimensional arrays as the sum product over the last axis of the first array and the second-to-last of the second array. For example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2512743135593745"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(np.random.rand(1000),b=np.random.rand(1000))/1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "will give you a number close to 0.25, since np.random.rand() is a uniformly extracted random number between 0 and 1. Similarly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.25286765, 0.26347566, 0.26084442, ..., 0.23534907, 0.24750937,\n",
       "        0.23862621],\n",
       "       [0.23163975, 0.24033687, 0.22608925, ..., 0.21354944, 0.23474152,\n",
       "        0.21846383],\n",
       "       [0.23561435, 0.25726754, 0.24283752, ..., 0.24222522, 0.24618642,\n",
       "        0.23717259],\n",
       "       ...,\n",
       "       [0.2368799 , 0.25273402, 0.23068352, ..., 0.23744693, 0.21994657,\n",
       "        0.2311918 ],\n",
       "       [0.25180896, 0.25710363, 0.27990532, ..., 0.24945943, 0.23746714,\n",
       "        0.24729834],\n",
       "       [0.26793496, 0.27259932, 0.27062371, ..., 0.25787922, 0.243902  ,\n",
       "        0.25019746]])"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(np.random.rand(100,100),b=np.random.rand(100,100))/100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "will result in a 100 × 100 matrix of numbers normally distributed around 0.25. \n",
    "\n",
    "Outer products are called by np.outer() and are an easy way to create an x, y regular mesh:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10])"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a=np.arange(0,11,1)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b=np.ones(11)\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  0.,   1.,   4.,   9.,  16.,  25.,  36.,  49.,  64.,  81., 100.],\n",
       "       [  1.,   0.,   1.,   4.,   9.,  16.,  25.,  36.,  49.,  64.,  81.],\n",
       "       [  4.,   1.,   0.,   1.,   4.,   9.,  16.,  25.,  36.,  49.,  64.],\n",
       "       [  9.,   4.,   1.,   0.,   1.,   4.,   9.,  16.,  25.,  36.,  49.],\n",
       "       [ 16.,   9.,   4.,   1.,   0.,   1.,   4.,   9.,  16.,  25.,  36.],\n",
       "       [ 25.,  16.,   9.,   4.,   1.,   0.,   1.,   4.,   9.,  16.,  25.],\n",
       "       [ 36.,  25.,  16.,   9.,   4.,   1.,   0.,   1.,   4.,   9.,  16.],\n",
       "       [ 49.,  36.,  25.,  16.,   9.,   4.,   1.,   0.,   1.,   4.,   9.],\n",
       "       [ 64.,  49.,  36.,  25.,  16.,   9.,   4.,   1.,   0.,   1.,   4.],\n",
       "       [ 81.,  64.,  49.,  36.,  25.,  16.,   9.,   4.,   1.,   0.,   1.],\n",
       "       [100.,  81.,  64.,  49.,  36.,  25.,  16.,   9.,   4.,   1.,   0.]])"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x=np.outer(a,b)\n",
    "x\n",
    "x2=np.outer(b,a)\n",
    "x2\n",
    "(x-x2)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFCFJREFUeJzt3X9sXfV5x/HPE9uJTSgJNC44gdUEtWkRyxZqdaUdVRt3pBOlQduEUqlr1U3KP9sa2KAjrdamnSaqgbZSdapk0R+gIhgKEYlL1TC5SGUaiuokVSj5ASWl5IfTmKY2EDmOEz/7w3YSB1L73nvO+Z773PdLquKc2Hmf55Y+Nfdc32PuLgBA/ZuT+gQAANlgoQNAECx0AAiChQ4AQbDQASAIFjoABMFCB4AgWOgAEAQLHQCCaC4ytmjRIu/s7CwyCQB1b/v27a+6e/tMn1foQu/s7FR/f3+RSQCoe2b269l8Hk+5AEAQLHQACIKFDgBBsNABIAgWOgAEMeOrXMzsu5I+Iemou183eewySf8tqVPSy5Juc/ff5XGCT+w8pHu37tPhoREtXtimu1Yt060rluSRKkU3ZZuZmTlqu1G6NtMdi8zsw5LekPTQOQv93yUdc/evm9ndki5193+eKdbV1eWVvGzxiZ2HtH7TcxoZO33mWFtLk+75iz/M9UFJ1U3ZZubiuinbzFyfXTPb7u5dM33ejE+5uPtPJR077/BqSQ9OfvygpFsrOrtZunfrvmkPhiSNjJ3WvVv35ZFL3k3ZZubiuinbzBy7W+1z6Je7+8Dkx0ckXX6hTzSztWbWb2b9g4ODFUUOD41UdDwrqbop28xcXDdlm5ljd2u+KOoTz9lc8Hkbd+9x9y5372pvn/EnV6dZvLCtouNZSdVN2Wbm4rop28wcu1vtQv+NmXVI0uSvR7M7pbPuWrVMbS1N0461tTTprlXL8sgl76ZsM3Nx3ZRtZo7drfa9XLZI+qykr0/+ujmzMzrH1IWDoq9Op+qmbDMzM+et0WZO0Z3Nq1wekfQRSYsk/UbSVyQ9IekxSX8g6deaeNni+RdO36TSV7kAAGb/KpcZv0N3909d4I+6Kz4rAEBu+ElRAAiChQ4AQbDQASAIFjoABMFCB4AgWOgAEAQLHQCCYKEDQBAsdAAIgoUOAEGw0AEgiGrfbbEwjx85pnv2D+jQ6JiWzGvR+qUd+ssrLgvbTdlmZmaO2m6UbqkX+uNHjunOfQc0Mj7xjpAHR8d0574DkpTrg5Kqm7LNzMzMzPXfLfVTLvfsHzjzYEwZGXfds3/gAl9R392UbWYurpuyzcyxu6Ve6IdGxyo6Xu/dlG1mLq6bss3MsbulXuhL5rVUdLzeuynbzFxcN2WbmWN3S73Q1y/tUNscm3asbY5p/dKOkN2UbWYurpuyzcyxu6W+KDp14aDoq9OpuinbzMzMeWu0mVN0Z7ynaJa4pygAVG629xQt9VMuAIDZY6EDQBAsdAAIgoUOAEGw0AEgCBY6AATBQgeAIFjoABAECx0AgmChA0AQLHQACKKmhW5md5jZ82b2CzN7xMxaszoxAEBlqn63RTNbIunzkq519xEze0zSGknfz+jcJEkDRzZr/0v36cTogFrndWjpNXeq44rVWSZK1U3ZZmZmjtpulG6tb5/bLKnNzMYkXSTpcO2ndNbAkc3au/dLGh8fkSSdGD2svXu/JEm5PiipuinbzMzMzFz/3aqfcnH3Q5Luk/SKpAFJw+7+VFYnJkn7X7rvzIMxZXx8RPtfui/LTGm6KdvMXFw3ZZuZY3erXuhmdqmk1ZKulrRY0nwz+/RbfN5aM+s3s/7BwcGKGidG3/pmqhc6npVU3ZRtZi6um7LNzLG7tVwU/ZikX7n7oLuPSdok6YPnf5K797h7l7t3tbe3VxRonffWt2q60PGspOqmbDNzcd2UbWaO3a1lob8i6QNmdpGZmaRuSXuyOa0JS6+5U3PmtE07NmdOm5Zec2eWmdJ0U7aZubhuyjYzx+5WfVHU3beZ2UZJOySdkrRTUk9WJyadvXBQ9NXpVN2UbWZm5rw12swputxTFABKjnuKAkCDYaEDQBAsdAAIgoUOAEGw0AEgCBY6AATBQgeAIFjoABAECx0AgmChA0AQLHQACIKFDgBB1HoLutzt2rVLfX19Gh4e1oIFC9Td3a3ly5eH7aZsMzMzR203SrfUC33Xrl3q7e3V2NiYJGl4eFi9vb2SlOuDkqqbss3MzMzM9d8t9VMufX19Zx6MKWNjY+rr6wvZTdlm5uK6KdvMHLtb6oU+PDxc0fF676ZsM3Nx3ZRtZo7dLfVCX7BgQUXH672bss3MxXVTtpk5drfUC727u1stLS3TjrW0tKi7uztkN2WbmYvrpmwzc+xuqS+KTl04KPrqdKpuyjYzM3PeGm3mFF3uKQoAJcc9RQGgwbDQASAIFjoABMFCB4AgWOgAEAQLHQCCYKEDQBAsdAAIgoUOAEGw0AEgCBY6AARR05tzmdlCSQ9Iuk6SS/obd382ixObcnznUb229WWdHhpV08J5umRVp+aveEeWiVJ1U7aZmZmjthulW+u7Ld4v6cfu/ldmNlfSRRmc0xnHdx7V0KYX5WPjkqTTQ6Ma2vSiJOX6oKTqpmwzMzMzc/13q37KxcwWSPqwpO9IkrufdPehrE5Mkl7b+vKZB2OKj43rta0vZ5kpTTdlm5mL66ZsM3Psbi3PoV8taVDS98xsp5k9YGbzz/8kM1trZv1m1j84OFhR4PTQaEXHs5Kqm7LNzMV1U7aZOXa3loXeLOl6Sd929xWSjku6+/xPcvced+9y96729vaKAk0L51V0PCupuinbzFxcN2WbmWN3a1noByUddPdtk7/fqIkFn5lLVnXKWqaforXM0SWrOrPMlKabss3MxXVTtpk5drfqi6LufsTMDpjZMnffJ6lb0u7sTu3shYOir06n6qZsMzMz563RZk7RrekWdGb2x5p42eJcSfslfc7df3ehz+cWdABQudnegq6mly26+88lzRgBAOSPnxQFgCBY6AAQBAsdAIJgoQNAECx0AAiChQ4AQbDQASAIFjoABMFCB4AgWOgAEAQLHQCCqPUWdLnb88zTeubRh/T6b1/V296+SDeu+Yzee+NHw3ZTtpmZmaO2G6Vb6oW+55mn9VTPt3Tq5MQdPl5/dVBP9XxLknJ9UFJ1U7aZmZmZuf67pX7K5ZlHHzrzYEw5dXJUzzz6UMhuyjYzF9dN2Wbm2N1SL/TXf/tqRcfrvZuyzczFdVO2mTl2t9QL/W1vX1TR8Xrvpmwzc3HdlG1mjt0t9UK/cc1n1Dx3+g1Vm+fO041rPhOym7LNzMV1U7aZOXa31BdFpy4cFH11OlU3ZZuZmTlvjTZzim5N9xStFPcUBYDKzfaeoqV+ygUAMHssdAAIgoUOAEGw0AEgCBY6AATBQgeAIFjoABAECx0AgmChA0AQLHQACIKFDgBB1LzQzazJzHaa2Q+zOCEAQHWyeLfFdZL2SLokg7/rTV7YdkTPbn5Jbxwb1cWXzdMNq6/Ru//kijxSpeimbDMzM0dtN0q3poVuZldKulnSv0n6x0zO6BwvbDuipx/eq1MnxyVJbxwb1dMP75WkXB+UVN2UbWZmZmau/26tT7l8Q9IXJI1ncC5v8uzml848GFNOnRzXs5tfyiOXvJuyzczFdVO2mTl2t+qFbmafkHTU3bfP8HlrzazfzPoHBwcrarxxbLSi41lJ1U3ZZubiuinbzBy7W8t36B+S9Ekze1nSo5JWmtkPzv8kd+9x9y5372pvb68ocPFl8yo6npVU3ZRtZi6um7LNzLG7VS90d1/v7le6e6ekNZJ+4u6fzuzMJN2w+ho1z51+is1z5+iG1ddkmSlNN2WbmYvrpmwzc+xuqe8pOnXhoOir06m6KdvMzMx5a7SZU3S5pygAlBz3FAWABsNCB4AgWOgAEAQLHQCCYKEDQBAsdAAIgoUOAEGw0AEgCBY6AATBQgeAIFjoABAECx0Agij9Qh/u7dWLK7u1573X6sWV3Rru7Q3dTdlmZmaO2m6UbqnfbXG4t1cD//Jl+YkTZ45Za6s6/vVrWnDLLXmcYtJuyjYzF9dN2Wbm+uzO9t0WS73QX1zZrVOHD7/pePPixXrXT/qyPLVSdFO2mbm4bso2M9dnN8Tb554aGKjoeL13U7aZubhuyjYzx+6WeqE3d3RUdLzeuynbzFxcN2WbmWN3S73Q33HH7bLW1mnHrLVV77jj9pDdlG1mLq6bss3MwbvuXth/3ve+93mlhrZs8Rc+utJ3v+e9/sJHV/rQli0V/x3VSNVN2WZmZo7arveupH6fxY4t9UVRAECQi6IAgNljoQNAECx0AAiChQ4AQbDQASAIFjoABMFCB4AgWOgAEAQLHQCCYKEDQBAsdAAIouqFbmZXmdnTZrbbzJ43s3VZntiUJ/c/qZs23qTlDy7XTRtv0pP7n8wjU5puyjYzM3PUdqN0q35zLjPrkNTh7jvM7G2Stku61d13X+hrKn1zrif3P6kN/7dBJ06fvYVTa1OrNnxwg25eenNV513mbso2MxfXTdlm5vrs5v7mXO4+4O47Jj9+XdIeSUuq/fveyv077p/2YEjSidMndP+O+7PMlKabss3MxXVTtpk5djeT59DNrFPSCknb3uLP1ppZv5n1Dw4OVvT3Hjl+pKLjWUnVTdlm5uK6KdvMHLtb80I3s4slPS7pdnd/7fw/d/ced+9y96729vaK/u4r5l9R0fGspOqmbDNzcd2UbWaO3a1poZtZiyaW+cPuvimbUzpr3fXr1No0/RZOrU2tWnd9Ltdfk3dTtpm5uG7KNjPH7jZt2LChqi80M5P0PUkH3P0rs/manp6eDWvXrp11492XvltLLl6i53/7vI6PHVfH/A7d/f67c794k6qbss3MzJy3Rps5y+5Xv/rVgQ0bNvTM9Hm1vMrlTyU9I+k5SeOTh7/o7j+60NdwCzoAqNxsX+XSXG3A3f9XklX79QCAbPGTogAQBAsdAIJgoQNAECx0AAiChQ4AQbDQASAIFjoABMFCB4AgWOgAEAQLHQCCYKEDQBDlX+i7HpP+8zppw8KJX3c9Frubss3MzBy13SDdqt+cqxC7HpN6Py+NjUz8fvjAxO8laflt8bop28zMzMxc991yf4fe97WzD8aUsZGJ4xG7KdvMXFw3ZZuZQ3fLvdCHD1Z2vN67KdvMXFw3ZZuZQ3fLvdAXXFnZ8Xrvpmwzc3HdlG1mDt0t90Lv/rLU0jb9WEvbxPGI3ZRtZi6um7LNzKG75V7oy2+TbvmmtOAqSTbx6y3fzP/iTapuyjYzM3PeGm3mBN2q7ylaDe4pCgCVm+09Rcv9HToAYNZY6AAQBAsdAIJgoQNAECx0AAiChQ4AQbDQASAIFjoABMFCB4AgWOgAEAQLHQCCqGmhm9nHzWyfmf3SzO7O6qQAAJWr+hZ0ZtYk6b8k/Zmkg5J+ZmZb3H13VicnSU/sPKR7t+7T4aERLV7YprtWLdOtK5ZkmShVN2WbmZk5artRurXcU/T9kn7p7vslycwelbRaUmYL/Ymdh7R+03MaGTstSTo0NKL1m56TpFwflFTdlG1mZmZmrv9uLU+5LJF04JzfH5w8lpl7t+4782BMGRk7rXu37ssyU5puyjYzF9dN2Wbm2N3cL4qa2Voz6zez/sHBwYq+9vDQSEXHs5Kqm7LNzMV1U7aZOXa3loV+SNJV5/z+yslj07h7j7t3uXtXe3t7RYHFC9sqOp6VVN2UbWYurpuyzcyxu7Us9J9JepeZXW1mcyWtkbQlm9OacNeqZWpraZp2rK2lSXetWpZlpjTdlG1mLq6bss3MsbtVXxR191Nm9veStkpqkvRdd38+szPT2QsHRV+dTtVN2WZmZs5bo82coss9RQGg5LinKAA0GBY6AATBQgeAIFjoABAECx0Agij0VS5mNijp11V++SJJr2Z4OvWAmRsDM8dX67zvdPcZfzKz0IVeCzPrn83LdiJh5sbAzPEVNS9PuQBAECx0AAiinhZ6T+oTSICZGwMzx1fIvHXzHDoA4Perp+/QAQC/R10s9Ea6GbWZXWVmT5vZbjN73szWpT6nophZk5ntNLMfpj6XIpjZQjPbaGZ7zWyPmd2Q+pzyZmZ3TP5z/Qsze8TMWlOfU9bM7LtmdtTMfnHOscvM7H/M7MXJXy/No136hX7Ozaj/XNK1kj5lZtemPatcnZL0T+5+raQPSPq74POea52kPalPokD3S/qxu79H0h8p+OxmtkTS5yV1uft1mnjb7TVpzyoX35f08fOO3S2pz93fJalv8veZK/1C1zk3o3b3k5KmbkYdkrsPuPuOyY9f18T/yIu5JXtCZnalpJslPZD6XIpgZgskfVjSdyTJ3U+6+1DasypEs6Q2M2uWdJGkw4nPJ3Pu/lNJx847vFrSg5MfPyjp1jza9bDQc78ZdVmZWaekFZK2pT2TQnxD0hckjac+kYJcLWlQ0vcmn2Z6wMzmpz6pPLn7IUn3SXpF0oCkYXd/Ku1ZFeZydx+Y/PiIpMvziNTDQm9IZnaxpMcl3e7ur6U+nzyZ2SckHXX37anPpUDNkq6X9G13XyHpuHL61/CymHzeeLUm/s9ssaT5ZvbptGdVPJ94aWEuLy+sh4U+q5tRR2JmLZpY5g+7+6bU51OAD0n6pJm9rImn1Faa2Q/SnlLuDko66O5T//a1URMLPrKPSfqVuw+6+5ikTZI+mPicivIbM+uQpMlfj+YRqYeFnvvNqMvEzEwTz6vucff/SH0+RXD39e5+pbt3auK/35+4e+jv3Nz9iKQDZjZ1x+BuSbsTnlIRXpH0ATO7aPKf824FvxB8ji2SPjv58Wclbc4jUvVNootSxM2oS+ZDkv5a0nNm9vPJY1909x8lPCfk4x8kPTz5jcp+SZ9LfD65cvdtZrZR0g5NvJprpwL+xKiZPSLpI5IWmdlBSV+R9HVJj5nZ32riHWdvy6XNT4oCQAz18JQLAGAWWOgAEAQLHQCCYKEDQBAsdAAIgoUOAEGw0AEgCBY6AATx/8kqAaCpOGWcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "y=np.outer(b,a)\n",
    "plt.plot(x,y,'o');plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NumPy allows composing large arrays combining blocks of one array scaled by another. This operation is called Kronecker product, and we will intensively use it to build large operators for two-dimensional continuum mechanics. For example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 1., 0., 0., 0., 0.],\n",
       "       [1., 1., 0., 0., 0., 0.],\n",
       "       [0., 0., 1., 1., 0., 0.],\n",
       "       [0., 0., 1., 1., 0., 0.],\n",
       "       [0., 0., 0., 0., 1., 1.],\n",
       "       [0., 0., 0., 0., 1., 1.]])"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.kron(np.eye(3), np.ones((2,2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Algebra\n",
    "\n",
    "All the linear algebra tools that we will need are already efficiently implemented through the standard optimized ATLAS LAPACK and BLAS libraries. Ultimately all linear algebra routines expect to operate on a 2-dimensional array. There is also a matrix type in NumPy but its use is discouraged since it is possible to obtain the same result by using arrays only.\n",
    "\n",
    "Linear algebra has many subprograms that run routines for decompositions such as Cholesky, QR and Singular Value. They can be very important for many problems, however they do not apply for the problems that we will address in this volume. The interested reader can refer to the regularly updated SciPy manual at http://docs.scipy.org/doc/numpy/reference/routines.linalg.html. \n",
    "\n",
    "Linear algebra (LA) routines are accessed through the module numpy.linalg:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy.linalg as linalg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An often employed tool is linalg.norm() both for calculating the size of a 1D (vector) and 2D (matrix) arrays. Norm, that is just the square root of the sum of the square of all the elements of a n-dimensional array, could be also calculated by the definition, but the numpy.linalg implementation can be one order of magnitude as more efficient. Let us for example benchmark the calculation of a 100 × 100 array norm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(201, 201)\n"
     ]
    }
   ],
   "source": [
    "a=np.arange(-100,101,1)\n",
    "b=np.arange(201)\n",
    "x=np.outer(a,b)\n",
    "\n",
    "x\n",
    "#normVectorized=np.sqrt(sum(sum(x**2)))\n",
    "#normVectorizedctorized\n",
    "print(x.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "292 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n 1 -r1 linalg.norm(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determinants are often employed in linear algebra because they are related to the invertibility of a matrix, and on whether the associated linear system of equations is solvable (a smaller determinant can indicate that some eigenvalues are close to zero, which makes inversion harder). However, to calculate the determinant of extremely large 2D arrays as we will often do in this volume is computationally so demanding that we will do it only for small problems. Let us calculate the determinant of few random matrices, 10 × 10 and 100 × 100, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.00878936,  0.00513914, -0.05539856,  0.00741107, -0.0006004 ,\n",
       "        0.03334587,  0.01023624, -0.02665326,  0.05143996, -0.02351945])"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dets=np.zeros(10)\n",
    "for i in np.arange(10): dets[i]=linalg.det(np.random.rand(10,10))\n",
    "dets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-2.09995750e+25, -1.16939088e+25, -1.55883240e+25, -5.88749974e+24,\n",
       "       -1.02730280e+26, -4.37716956e+25,  1.28431610e+24, -1.67252479e+25,\n",
       "        1.43810948e+25,  1.45009691e+24])"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for i in np.arange(10): dets[i]=linalg.det(np.random.rand(100,100))\n",
    "dets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cython\n",
    "\n",
    "There are situations in which a simple and straightforward NumPy solution does not exist or can be very difficult to implement. We will see some examples in the last, advanced part. When we are in this situation, it is first of all important to verify first whether this lack of vectorization makes our code sensibly slower. Since this is often the case, a powerful solution has been created by the authors of the parts that compose the Python Scientific Environment, which is called Cython.\n",
    "\n",
    "The most recent updates and documentation for Cython are available at the address cython.org). Cython is based on Pyrex, library that allowed to extend Python with C functions. Recently, however Cython developers went beyond and created an extended language where codes written in Python are optimized and compiled just at the moment of running the Python program. In other words, while the Python language is interpreted, the computer will on-the-fly compile the parts of your code written in Python in special Cython files. Still they will be written in Python, which means that the user does not need to learn a new, more challenging language, such as C. The result is a superior language, only slightly more complex than Python, that combines the best of the two worlds, C and Python.\n",
    "\n",
    "The main addition of Python to Cython are static variables. These are variables for which type and space in memory is predefined by the programmer. The advantage of using them is that the Python parser does not need to understand how to interpret them on the fly (one of the most original but computationally most inefficient characteristic of Python) while will straightaway compile and execute the fastest possible way to run the code. In practice, Cython is simply Python code with variables defined using C data types. Typically it consists of one or few functions that require being compiled with explicitly defined types and called in our code. Since the performance of Cython code are nearly equal to an equivalent C-code, but are much easier written, Cython is the definitive tool that makes prototyping in Python the easy and powerful programming tool that every scientist dreams."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cython in iPython\n",
    "\n",
    "Let us look at some examples of how it works, starting with iPython. I will use IPython 4.0.1. Different versions of IPython can make use of other magic commands.\n",
    "\n",
    "In iPython first you have to load the Cython library. This can be done with the magic command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Immediately after, one can use both the magic commands %cython, and %%cython, the first to run external Cython files, and the second to run a Cython cell. Let us see how we could have implemented the function setNegativeValuesToZero that we intro- duced earlier in the chapter. A first attempt would be to simply rewrite that function, just adding that n, m, i and l are integer types:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "def setNegativeValuesToZero(int n,int m, a):\n",
    "    cdef int i, j\n",
    "    for i in range(n):\n",
    "        for j in range(m): \n",
    "            if a[i,j]<0:\n",
    "                a[i,j]=0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As one reads on line 2, the type declaration is done simply by using the instruction cdef, then followed by the standard static C type declaration. Notice also that I replace np.arange with range, since this is the compiled version where loops are compiled. Basically, besides the addition of the int, i.e., integer, inside the call, this is exactly the same function introduced earlier. One can test its performance with the same random 2000 × 2000 array finding an improvement of a factor two compared to the pure Python version:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2000, 2000)\n",
      "1.39 s ± 103 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "arr2d=np.random.randn(2000,2000)\n",
    "print(arr2d.shape)\n",
    "%timeit setNegativeValuesToZero(2000,2000,arr2d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can one do better? Yes, but it requires some more work. I show here a version that performs practically exactly the same as the NumPy one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "cimport numpy as np\n",
    "def setNegativeValuesToZero(int n, int m, np.ndarray[double, ndim=2] a):\n",
    "    cdef int i, j\n",
    "    for i in range(n):\n",
    "        for j in range(m): \n",
    "            if a[i,j]<0:\n",
    "                a[i,j]=0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As one observes here, it is necessary to add one more line and change one type. The first is the import of numerical python within the Cython framework. This allows to set the type of the NumPy array a, which in turn allows Cython to accelerate its process. Once these two little modifications are done, one obtains the striking result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 5.28 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
      "8.75 ms ± 7.83 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "arr2d=np.random.randn(2000,2000)\n",
    "%timeit -n 1 setNegativeValuesToZero(2000,2000,arr2d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which is very close to NumPy, and in this case slightly better.\n",
    "A good Python programmer, although not Cython expert, will need only ten minutes to figure out how to modify a routine to make it ready for Cython, however for the novice programmer to go through the debugging effort that compiled programs require can be quite challenging. For this reason I recommend always using the straightforward NumPy libraries whenever necessary, and rely on Cython only when really necessary. The many more details that one has to know in order to program Cython are available in the excellent recent book Cython - A Guide for Python Programmers.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Going Parallel: mpi4py and PETSc4py\n",
    "\n",
    "mpi4py is a simple but powerful interface that allows to write MPI software natively in Python. This is however not the only option, and other projects exist such as pyMPI, that is also very simple but does not support numerical python arrays, and pypar.\n",
    "\n",
    "mpi4py is presently a translation of standard MPI-2 bindings for C++ to Python. There are two categories of variables that can be passed from one processor to the other: (i) the standard Python objects, and (ii) the NumPy arrays. Because the memory configuration of the first kind is not compact nor standard, mpi4py achieves maximum performance by building binary representations of these objects that are then sent to other processors. Clearly building this binary representation limits its use and can create a serious overhead for very large systems or data. Furthermore this process, called pickling, can also require extra memory and it needs to be tested on every machine.\n",
    "\n",
    "In order to use mpi4py, one has to install first MPI independently from the Python installation. Presently, on the windows platform, one has to install OpenMPI through Cygwin. Once this has been done, one can install mpi4py in Cygwin, or any other Linux system as well as the OSX Terminal, simply with the command pip install mpi4py. If using Anaconda the command conda install -c anaconda mpi4py=2.0.0 will add mpi4py to the Anaconda installation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpi4py import MPI \n",
    "\n",
    "comm = MPI.COMM_WORLD\n",
    "rank = comm.Get_rank()\n",
    "\n",
    "if rank == 0:\n",
    "    data = {'a': 0.1, 'b': 1, 'c':'house', 'd':1.0e100} \n",
    "    comm.send(data, dest=1, tag=111)\n",
    "    print('data sent from',rank,' :',data)\n",
    "elif rank == 1:\n",
    "    data = comm.recv(source=0, tag=111) \n",
    "    print('data received from',rank,' :',data)"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}