{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An introduction to the ```numpy``` array library\n",
    "\n",
    "*This notebook is taken from a course of Numerical Physics at Ecole Polytechnique, many thanks to his author Michel Ferrero (CPHT, Ecole Polytechnique) for his authorisation to distribute it.*\n",
    "\n",
    "As you have read in the notebook introducing the ```python``` language, a convenient way\n",
    "to gather objects together are lists. As a physicist you will very often want to work with\n",
    "vectors and matrices. It could be tempting to use lists to do this. For example, you could\n",
    "very well define a vector ```v = [1,2,3]``` and then access its element with, e.g. ```v[1]```\n",
    "etc.\n",
    "\n",
    "However, while lists are very generic objects, they have some limitations. For example,\n",
    "they do not behave like what we would expect for a vector. To see this, try to execute the code\n",
    "below and see if you can anticipate the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# what is the result of this?\n",
    "# do you think it will be [3,2,4]?\n",
    "print([1, 2, 3] + [2, 0, 1])\n",
    "\n",
    "# what about this?\n",
    "print(2 * [1, 2, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you see, lists do not behave like vectors. The reason for this is that lists can contain anything. You can\n",
    "create a list ```[1,'a',3.2]``` that mixes different types. What would a multiplication by 2 mean in this case?\n",
    "Finally, because of this flexibility, lists are also rather slow.\n",
    "\n",
    "To overcome these limitations, there is a very useful library called ```numpy``` (for ```num```eric ```py```thon)\n",
    "that provides multidimensional arrays. Arrays can be used to create vectors and matrices. As we will see, the\n",
    "```numpy``` library also has tools to do mathematical operations, linear algebra, Fourier transforms, etc.\n",
    "In order to use `numpy` you first have to import it. Often you do this and give `numpy` the alias name `np`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create your first arrays\n",
    "\n",
    "The example below shows how to create arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1, 2, 3])               # create from a list\n",
    "b = np.zeros(3)                       # an array full of zeros (1-dimensional)\n",
    "c = np.zeros((2, 3))                  # a matrix full of zeros (2-dimensional)\n",
    "d = np.zeros((2, 3), dtype=complex)   # by default an array contains floats\n",
    "\n",
    "print(\"a =\", a)\n",
    "print(\"b =\", b)\n",
    "print(\"c =\\n\", c)\n",
    "print(\"d =\\n\", d)\n",
    "print(\"shape of d:\", d.shape)         # in order to know the dimensions of an array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Other common functions to create arrays\n",
    "\n",
    "Often we do not want to create arrays manually. These functions can be useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a range of numbers (like the python range command)\n",
    "a = np.arange(5, 25, 2)         # like range this will stop at 23\n",
    "print(a)\n",
    "\n",
    "# Create a discretized interval\n",
    "b = np.linspace(1, 2, 5)       # start, stop, number of points\n",
    "print(b)\n",
    "\n",
    "# You can get the step by adding an argument\n",
    "c, step = np.linspace(1, 2, 5, retstep=True)\n",
    "print(\"step = \", step)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Indexing and slicing arrays\n",
    "\n",
    "You access an element of an array with the bracket operator. You have to be careful,\n",
    "just like for a `list`, the index starts at 0 for arrays. The bracket also allows you\n",
    "to very easily take partial views of an array. Note that if you modify a slice of an\n",
    "array, the original array will be modified too! This is one of the powers of the `numpy`\n",
    "arrays. Slicing does not introduce any useless copy of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.random.randint(10, size=(3,3))     # this is a random 3x3 matrix\n",
    "print(a)\n",
    "\n",
    "# Indexing\n",
    "print(a[0, 0])                            # this is the first matrix element\n",
    "print(a[1, 2])                            # second row, third column\n",
    "\n",
    "# Slicing (partial views)\n",
    "print(a[0, :])                            # this is the first line\n",
    "print(a[:, 1])                            # this is the second column\n",
    "print(a[0:2, 0:2])                        # this is the upper left 2x2 corner\n",
    "\n",
    "b = a[0:2, :]                             # this is the upper 2x3 slice\n",
    "b[0, 1] = 2                               # b is a slice of a, if modified a will be modified too\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Standard operations with arrays\n",
    "\n",
    "The rule is that when you write operations on arrays it applies the\n",
    "operation on all elements, one after the other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(7)\n",
    "print(a)\n",
    "print(a**2)\n",
    "print(2 * a)\n",
    "\n",
    "a = np.array([1, 2, 3])\n",
    "b = np.array([2, 1, 4])\n",
    "\n",
    "print(a + b)\n",
    "print(a * b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applying a function on an array\n",
    "\n",
    "The standard operations above act on all elements of an array. Of course, one would also\n",
    "want to be able to do this with functions, like e.g. `cos`. If you use the cosine function\n",
    "that is provided in the `math` library of the standard library of `python`, it will\n",
    "not be able to act on an array. For this reason, all mathematical functions have been\n",
    "included in the `numpy` module and generalized so that they can also act on arrays.\n",
    "Therefore I strongly recommend that you always use mathematical functions from the\n",
    "`numpy` module. Here is an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.linspace(0, 2*np.pi, 10)    # note that np.pi gives the pi constant\n",
    "b = np.cos(a)                      # this works and returns an array of cosines\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It may happen that you want to act on the elements of an array with a function that is\n",
    "not part of the `numpy` library. For example, a function that you have defined yourself.\n",
    "In that case, you can *vectorize* it: it turns the function into a new function that\n",
    "can also act on arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is my own function (acting on numbers only)\n",
    "def my_function(x):\n",
    "    if x > 0: return 1\n",
    "    else: return -1\n",
    "\n",
    "a = np.array([-10, -2, 4, 24, 5])\n",
    "\n",
    "# This would not work: my_function(a)\n",
    "# so I \"vectorize\" the function: it can now act on arrays\n",
    "vfunc = np.vectorize(my_function)\n",
    "\n",
    "# Now this works as expected\n",
    "print(vfunc(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Copying arrays: a common mistake\n",
    "\n",
    "When you want to copy an array, be sure to use `.copy()` as in this example:\n",
    "\n",
    "```python\n",
    "a = np.array([1,2,3])\n",
    "b = a.copy()\n",
    "```\n",
    "\n",
    "Indeed, you might have a surprise otherwise. Check the output of this code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1, 2, 3])\n",
    "print(\"original a: \", a)\n",
    "\n",
    "b = a                              # this *identifies* b and a\n",
    "b[0] = 4                           # b is the same as a so this line modifies a too\n",
    "print(\"a has been modified: \", a)\n",
    "\n",
    "b = a.copy()\n",
    "b[0] = 12\n",
    "print(\"a is unchanged: \", a)\n",
    "print(\"b has changed: \", b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scalar product and matrix product\n",
    "\n",
    "As you saw above, if you consider two vectors or two matrices with the\n",
    "same size and write `a * b` it will do an elementwise multiplication.\n",
    "If you instead want to do a scalar product or a matrix multiplication,\n",
    "you should write it like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1, 2, 3])\n",
    "b = np.array([2, 1, 1])\n",
    "\n",
    "print(\"a scalar product b = \", np.dot(a,b)) # or a.dot(b)\n",
    "\n",
    "a = np.array([[1, 2], [2, 3]])\n",
    "b = np.array([[1, 1], [2, 2]])\n",
    "\n",
    "print(\"matrix product = \\n\", np.dot(a,b))\n",
    "\n",
    "# If you're using Python version 3.5 or higher, then the '@' symbol can be used for matrix multiplication\n",
    "# a @ b == np.dot(a, b)\n",
    "# c = np.array([[4, 7], [-10, 15]])\n",
    "# a @ b @ c == a.dot(b).dot(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Saving and loading arrays to text files\n",
    "\n",
    "It can often be useful to save the data of an array to a file. Typically, you could\n",
    "have a first production code that generates and saves some data to a file and\n",
    "a second code that reads the data and analyzes and plots it. These two functions\n",
    "can then be useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save an array in a text file\n",
    "a = np.array([[1,2],[3,4]])\n",
    "np.savetxt(\"my_data.dat\", a)\n",
    "\n",
    "# read the array from the text file\n",
    "b = np.loadtxt(\"my_data.dat\")\n",
    "print(b)\n",
    "\n",
    "# if you want each column separately, then pass unpack=True\n",
    "x, y = np.loadtxt(\"my_data.dat\", unpack=True)\n",
    "print(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More information\n",
    "\n",
    "Of course arrays can do a lot more than what has been shown here.\n",
    "In particular, there are many functions in `numpy` that allow to operate\n",
    "on arrays. You can do Fourier transforms, polynomials fits, linear\n",
    "algebra operations and many more things. Here are some useful links\n",
    "to learn more:\n",
    "\n",
    "- If you want to learn a bit more you can follow this (still quite short)\n",
    "  [tutorial](http://docs.scipy.org/doc/numpy/user/quickstart.html).\n",
    "  \n",
    "- This is the [tutorial](http://www.scipy-lectures.org/intro/numpy/index.html) from\n",
    "  the scipy lectures.\n",
    "  \n",
    "- You can of course read the more complete\n",
    "  [user's guide](http://docs.scipy.org/doc/numpy/user/index.html)\n",
    "\n",
    "- If you want details about a function you know the name of, you can type it\n",
    "  in this [search page](http://docs.scipy.org/doc/numpy/search.html). You can\n",
    "  also just google the numpy function.\n",
    "\n",
    "- If you don't know the name of the function and look for something that might\n",
    "  be in numpy, you can browse through the reference and see if you\n",
    "  can find it: [reference](http://docs.scipy.org/doc/numpy/reference/index.html)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}