{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import timeit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 1 : Python lists vs Numpy Lists\n",
    "Aim : Learn the difference between python lists and numpy arrays. Practice creating lists of given lengths and studying the data.\n",
    "\n",
    "When using python, you will be faced with two different types of lists.\n",
    "- **Python Lists** are simple to use and modify. Empty lists can be initialised using  `list = []`. To add elements to them, simply use `list.append(value)` to add your value to the list This value can be a number, a string, or even another list. Two lists can be added together by using `new_list = list1 + list2`.  \n",
    "The advantage of python lists is that they can freely grow and shrink in size, making them realy usefull for quick codes where you do not have a precise idea of the length of your list.\n",
    "- **Numpy Arrays** come from the numpy package. These arrays have a fixed size, which is an advantage and a disadvantage. They cannot grow or shrink in size but come with a multitude of in-built methods for rapid data manipulation. They can be initialised from a python list `list = np.array(python_list)`, filled with zeros (`np.zeros()`) or ones (`np.ones()`). Due to their nature, everytime you initialize a numpy array from nothing, you are required to enter the shape you want the array to have. This is done by giving it the `shape` parameter.   \n",
    "`shape = 10` or `shape = (10,)` creates a normal \"1D\" list of size 10. `shape = (3,5)` creates a 3x5 matrix (3 rows and 5 columns). By doing so, you can create a list of any dimension you want.  \n",
    "\n",
    "In both cases, to access elements of the list, precise the indice of the element you wish to access. In a 1D list, `list[0]` will give you the first element, `list[4]` will give you the 5th element and `list[-1]` the last one. `list[1:4]` will give you a list with all the elements from indice 1 to 3, `list[:5]` from the first to the 5th element and `list[1:5:2]` every 2 elements between indice 1 and 4.   \n",
    "If your list is a 2D matrix, then doing `list[4]` will give you a list : the 5th row. So `list[4][2]` gives the 3rd element of the 5th row. Accessing columns is often esier with numpy arrays, where `list[:,4]` return the 5th column (all elements who have 4 on their second indice)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's comparate the performance of python lists and numpy arrays. To do this, we need to create a matrix of random numbers of a size NxN, using python lists and numpy arrays,\n",
    "and then measure the time it takes. In order to measure the time of the code evaluation, we will use the following wrapper of our code:\n",
    "\n",
    "```python\n",
    "import timeit\n",
    "start_time = timeit.default_timer()\n",
    "# Your code here\n",
    "end_time = timeit.default_timer()\n",
    "print(\"Time taken : \", end_time - start_time)\n",
    "```\n",
    "\n",
    "Alternatively, you can use the `%%timeit` magic command if you are using a Jupyter notebook. This will run the code multiple times and give you an average time.\n",
    "\n",
    "```python\n",
    "%%timeit\n",
    "# Your code here\n",
    "```\n",
    "\n",
    "Create the NxN matrix of a size N = 100, 1000, and 3000, filled with random numbers from a starndard normal distribution (mean = 0, std = 1) using python lists and numpy arrays.\n",
    "You can generate a single the random number using `np.random.normal()` command. In the case of python lists you can use the nested `for` loops to create the NxN matrix. In the case of numpy - you can just put the `size` parameter into a `np.random.normal()` function call. Compare the time it takes to generate the matrix with both methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100 # Change N here for larger arrays\n",
    "\n",
    "# Use Python lists\n",
    "start = ...\n",
    "# Your code\n",
    "end = ...\n",
    "print(\"Time taken with python lists : \", end-start)\n",
    "\n",
    "# Use Numpy Array to create\n",
    "start = ...\n",
    "# Your code\n",
    "end = ...\n",
    "print(\"Time taken with numpy arrays : \", end-start)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try other values of N and plot the time it takes to generate the matrix with both methods. You can use `matplotlib` package to plot the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_values = [100, 1000, 3000, 5000, 10000]\n",
    "python_timings = []\n",
    "numpy_timings = []\n",
    "for N in N_values:\n",
    "    # Your code\n",
    "    python_time = end - start  # Measure time for python lists\n",
    "    # Your code\n",
    "    numpy_time = end - start  # Measure time for numpy arrays\n",
    "\n",
    "    python_timings.append(python_time)\n",
    "    numpy_timings.append(numpy_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting the results\n",
    "plt.plot(N_values, python_timings, label='Python Lists')\n",
    "plt.plot(N_values, numpy_timings, label='Numpy Arrays')\n",
    "plt.xlabel('N (size of NxN matrix)')\n",
    "plt.ylabel('Time (seconds)')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's calculate the mean and the standard deviation of each row in the generated matrix. To do this, you can use the `np.mean()` and `np.std()` functions (or `array.mean()` and `array.std()` methods) for numpy arrays. The direction over which the mean or std is calculated can be specified using the `axis` parameter. For example, `np.mean(array, axis=1)` will calculate the mean of each row, while `np.mean(array, axis=0)` will calculate the mean of each column.\n",
    "\n",
    "Additionally, let's try to find the mininum and maximum values in the generated matrix, as well as their positions (indices) in the array. You can use the `np.min()`, `np.max()`, `np.argmin()`, and `np.argmax()` functions (or `array.min()`, `array.max()`, `array.argmin()`, and `array.argmax()` methods) for numpy arrays. The `np.argmin()` and `np.argmax()` functions return the indices of the minimum and maximum values in the flattened array. To get the indices in the original array shape, you can use the `np.unravel_index()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100\n",
    "data = np.random.normal(size=(N, N))\n",
    "\n",
    "mean = ...\n",
    "std = ...\n",
    "\n",
    "print(\"Mean of each row: \")\n",
    "print(mean)\n",
    "print(\"Standard deviation of each row: \")\n",
    "print(std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_value = ...\n",
    "max_index_flattened = ...\n",
    "max_index = ...\n",
    "\n",
    "min_value = ...\n",
    "min_index_flattened = ...\n",
    "min_index = ...\n",
    "\n",
    "print(f\"Flattened index {max_index_flattened} vs. 2D index: {max_index}\")\n",
    "print(f\"Maximum value: {max_value} at index {max_index}\")\n",
    "print(f\"Minimum value: {min_value} at index {min_index}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 2. Gaussian distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this part, we practice with creating gaussian distributions, plotting them and comparing them to the gaussian function. First, let's create a 1D array of N=10000 points\n",
    "sampled from the normal distribution with a certain mean and standard deviation. You can use the `np.random.normal(mean, std)` function to generate these points. Then, plot a histogram of these points using the `plt.hist()` function from the `matplotlib` package. Use the parameter `density=True` to normalize the distribution. Control the number of bins in the histogram using the `bins` parameter. Draw the vertical line for the mean and the standard deviation on the histogram using `plt.axvline()` function and the place of the expected mean value. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10000\n",
    "mean = 4\n",
    "std = 2\n",
    "num_bins = 101\n",
    "\n",
    "# Your code here\n",
    "data = ...\n",
    "plt.hist(...)\n",
    "plt.axvline(...)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's compare the histogram to the theoretical gaussian function. The gaussian function is given by:\n",
    "\n",
    "$$f(x) = \\frac{1}{\\sigma \\sqrt{2\\pi}} e^{-\\frac{(x - \\mu)^2}{2\\sigma^2}}$$\n",
    "\n",
    "where $\\mu$ is the mean and $\\sigma$ is the standard deviation. Let's try to analyze the difference between the \"empirical\" distribution (the histogram) and the theoretical one (the gaussian function). To do this, let's get the bin centers from the histogram and calculate the gaussian function at these points. Then, plot the histogram and the gaussian function on the same graph for comparison. Both bin centers and the histogram heights can be obtained from the output of the `plt.hist()` function as follows:\n",
    "\n",
    "```python\n",
    "freqs, bin_edges, _ = plt.hist(...)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "freqs, bins, _ = plt.hist(...)\n",
    "gaussian = ... # Evaluate the gaussian function at bins using the formula\n",
    "\n",
    "plt.plot(bins, gaussian, color='red', label='Gaussian Function')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's plot the difference between the empirical and theoretical distributions. To do this, calculate the difference between the histogram heights and the gaussian function values at the bin centers. Then, plot this difference as a function of the bin centers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "difference = ... # Calculate the difference between freqs and gaussian\n",
    "plt.plot(bins, difference, label='Difference')\n",
    "plt.xlabel('Value')\n",
    "plt.ylabel('Difference')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
