Lecture2conv.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Programming in Python, Mathematics, Statistics, and Modeling\n",
    "## Lecture 2:\n",
    "This lecture contains the contents: \n",
    "* Linear Algebra\n",
    "    * Solving linear system with Python\n",
    "    * Simple Linear Ordinary Differential Equation\n",
    "    * Least Squares solution \n",
    "\n",
    "This lecture starts with a brief recap of the most important things mentioned in the self-study slides on Linear Algebra. Some concepts will be illustrated further by implementing them in Python. A elobarote in-class exercise will be carried out simulating gas-diffusion in a layer reservoir (represented by a linear system, i.e. Linear Ordinary Differential Equation). Finally the Least Squares method will be explained using the concept of orthogonality. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recap of lecture on programming and self-study "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_03/Slide2.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise 1: Given a Python list, remove all occurrence of 20 from the list\n",
    "\n",
    "list1 = [5, 20, 15, 20, 25, 50, 20]\n",
    "\n",
    "# expected output [5, 15, 25, 50]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise 2: Unpack the following tuple into 4 variables\n",
    "aTuple = (10, 20, 30, 40)\n",
    "\n",
    "# Your code here\n",
    "\n",
    "print(a) # should print 10\n",
    "print(b) # should print 20\n",
    "print(c) # should print 30\n",
    "print(d) # should print 40"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise 3: Reverse a given string\n",
    "str1 = \"PYnative\"\n",
    "print(\"Original String is:\", str1)\n",
    "\n",
    "# expected output: evitanYP\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Excersise 4: Given a Python dictionary, Change Brad’s salary to 8500\n",
    "sampleDict = {\n",
    "     'emp1': {'name': 'Jhon', 'salary': 7500},\n",
    "     'emp2': {'name': 'Emma', 'salary': 8000},\n",
    "     'emp3': {'name': 'Brad', 'salary': 6500}\n",
    "}\n",
    "\n",
    "# Your code here\n",
    "\n",
    "print(sampleDict['emp3'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_02/Slide2.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_02/Slide4.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_02/Slide5.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_02/Slide6.PNG\" width=\"50%\"> \n",
    "<img style=\"float: right;\" src=\"lecture_02/Slide7.PNG\" width=\"50%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_02/Slide8.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_02/Slide9.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_02/Slide12.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve the linear system using Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "# Define linear system and right-hande-side (rhs) for following system:\n",
    "# x_1   - 3x_2  + 3x_3  =  2\n",
    "# 0x_1  + 6x_2  - 1x_3  =  0\n",
    "# 7x_1   - 0x_2  + 9x_3  =  7\n",
    "# Allocate memory:\n",
    "# rows x cols = m x n = ???\n",
    "m = ?\n",
    "n = ?\n",
    "A = np.zeros((?, ?))\n",
    "rhs = np.zeros((?,))\n",
    "\n",
    "# Fill matrix and rhs vector:\n",
    "# your code here!\n",
    "\n",
    "# Solve system:\n",
    "# hint: use np.linalg.??? (when you write \"np.linalg.\" you can \n",
    "# press \"Tab\" to do code complete, it will show which functions are available!)\n",
    "solution = np.linalg.???(A, rhs)\n",
    "\n",
    "print(\"Solution: [x1, x2, x3] = [{:3.2f}, {:3.2f}, {:3.2f}]\".format(solution[0], \n",
    "                                                                    solution[1],\n",
    "                                                                    solution[2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_02/Slide18.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_02/Slide19.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_02/Slide21.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What does Python when there is no solution?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "# Define linear system and right-hande-side (rhs):\n",
    "A = np.zeros((2, 2))\n",
    "rhs = np.zeros((2,))\n",
    "\n",
    "A[0, :] = [2, -1]\n",
    "A[1, :] = [-4, 2]\n",
    "rhs[:] = [0, -2]\n",
    "\n",
    "# Solve system:\n",
    "solution = np.linalg.solve(A, rhs)\n",
    "\n",
    "print(\"Solution: [x1, x2] = [{:3.2f}, {:3.2f}]\".format(solution[0], solution[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_02/Slide22.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What does Python when there are infinite number of solutions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "# Define linear system and right-hande-side (rhs):\n",
    "A = np.zeros((2, 2))\n",
    "rhs = np.zeros((2,))\n",
    "\n",
    "A[0, :] = [2, -1]\n",
    "A[1, :] = [-4, 2]\n",
    "rhs[:] = [1, -2]\n",
    "\n",
    "# Solve system:\n",
    "solution = np.linalg.solve(A, rhs)\n",
    "\n",
    "print(\"Solution: [x1, x2] = [{:3.2f}, {:3.2f}]\".format(solution[0], solution[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_02/Slide24.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_02/Slide27.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_02/Slide28.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_02/Slide29.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import pi, sin, cos\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# Some transformation (use the angle of 60 degrees):\n",
    "angle = ?\n",
    "\n",
    "L = np.zeros((?, ?))\n",
    "L[0, :] = [?, ?]  # first row\n",
    "L[1, :] = [?, ?]  # second row\n",
    "\n",
    "# Original vector:\n",
    "vec1 = [1, 1]\n",
    "\n",
    "# Result (which method to use for matrix-vector product)\n",
    "resvec = np.???(L, vec1)\n",
    "\n",
    "# Plot result:\n",
    "# Create dictionary that holds information on font:\n",
    "font = {'family': 'serif',\n",
    "        'color':  'darkred',\n",
    "        'weight': 'normal',\n",
    "        'size': 16,\n",
    "        }\n",
    "\n",
    "# Plot solution using the imported function matplotlib.pyplot referenced by plt:\n",
    "plt.figure(num=None, figsize=(15, 5), dpi=80, facecolor='w', edgecolor='k')\n",
    "\n",
    "plt.plot([0, vec1[0]], [0, vec1[1]], 'red', linewidth=5)\n",
    "plt.plot([0, resvec[0]], [0, resvec[1]], 'blue', linewidth=5)\n",
    "plt.title('Visualization of Transformation', fontdict=font)\n",
    "plt.xlabel('x_1', fontdict=font)\n",
    "plt.ylabel('x_2', fontdict=font)\n",
    "plt.axis('equal')\n",
    "plt.gca().legend(('org vec','res vec'))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_02/Slide30.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_02/Slide31.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Projections and Least Squares\n",
    "### Definitions:\n",
    "* Column space of **A**, denoted as C(**A**), is the space span by the columns of **A**.\n",
    "* For any linear system **Ax**=**b**, if **b** is not in C(**A**), the system doesn't have a solution. \n",
    "* Any two vectors **u** and **v** that are orthogonal to eachother have a dot (or inner) product of zero.\n",
    "* Shortest distance between a plane and a point above a plane is the length of the orthogonal line between the particular point and plane. \n",
    "\n",
    "### Simple example: Data points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"least_squares/Slide1.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"least_squares/Slide2.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"least_squares/Slide3.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example in Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# Some parameters:\n",
    "num_points = 100\n",
    "length = 100\n",
    "segm_length = length / num_points\n",
    "indep_var = np.linspace(segm_length/2, length - segm_length/2, num_points)\n",
    "\n",
    "# Noise related parameters:\n",
    "np.random.seed(0)\n",
    "use_trend = True\n",
    "if use_trend:\n",
    "    noise_fac = 0.5\n",
    "else: \n",
    "    noise_fac = 25\n",
    "noise = np.random.uniform(low=-noise_fac, high=noise_fac, size=(num_points,))\n",
    "\n",
    "# Original line parameters:\n",
    "slope = 5\n",
    "y_intersect = 10\n",
    "\n",
    "if use_trend:\n",
    "    # Define noisy data (with trend):\n",
    "    data = y_intersect + indep_var * (slope + noise)\n",
    "    extra_txt = \" with Trend\"\n",
    "else:\n",
    "    # Define noisy data (without trend):\n",
    "    data = y_intersect + indep_var * slope + noise\n",
    "    extra_txt = \" without Trend\"\n",
    "\n",
    "# Plot data:\n",
    "plt.figure(num=None, figsize=(15, 5), dpi=80, facecolor='w', edgecolor='k')\n",
    "\n",
    "font = {'family': 'serif',\n",
    "        'color':  'darkred',\n",
    "        'weight': 'normal',\n",
    "        'size': 16,\n",
    "        }\n",
    "\n",
    "plt.plot(indep_var, data, 'or', markersize=3)\n",
    "plt.title('Noisy data' + extra_txt, fontdict=font)\n",
    "plt.xlabel('x', fontdict=font)\n",
    "plt.ylabel('data', fontdict=font)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform least squares:\n",
    "A = np.zeros((num_points, 2))\n",
    "A[:, 0] = np.ones((num_points, ))\n",
    "A[:, 1] = indep_var\n",
    "\n",
    "# Solving the normal equations:\n",
    "solution = np.linalg.solve(np.dot(A.T, A), np.dot(A.T, data))\n",
    "\n",
    "# Plot solution:\n",
    "plt.figure(num=None, figsize=(15, 5), dpi=80, facecolor='w', edgecolor='k')\n",
    "\n",
    "plt.plot(indep_var, data, 'or', markersize=3)\n",
    "plt.plot(indep_var, solution[0] + solution[1]*indep_var, 'black', linewidth=2)\n",
    "plt.title('Linear Regression' + extra_txt, fontdict=font)\n",
    "plt.xlabel('x', fontdict=font)\n",
    "plt.ylabel('data', fontdict=font)\n",
    "plt.gca().legend(('data','regres'))\n",
    "plt.show()\n",
    "\n",
    "print(\"Original:\\t\\tY-intersect is: {:4.2f}, and Slope is: {:4.2f}\".format(y_intersect, slope))\n",
    "print(\"Linear regression:\\tY-intersect is: {:4.2f}, and Slope is: {:4.2f}\".format(solution[0], solution[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some programming:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_03/Slide5.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "# Create matrix of coefficients:\n",
    "A = np.zeros((3, 3))\n",
    "A[0, :] = [0.8, 0.1, 0.04]\n",
    "A[1, :] = [0.15, 0.7, 0.13]\n",
    "A[2, :] = [0.05, 0.2, 0.83]\n",
    "\n",
    "# Create initial gas distribution:\n",
    "gas_distr = np.array([1, 0, 0])\n",
    "\n",
    "# Time step size:\n",
    "dt = 1  # [years]\n",
    "    \n",
    "# Calculate distribution at several timesteps:\n",
    "time_after_injection = [6, 15, 50]  # [years]\n",
    "nr_time_steps = [int(ii / dt) for ii in time_after_injection]  # matrix powers (and time-steps) need to be integers :)! \n",
    "\n",
    "for n in nr_time_steps:\n",
    "# your code here .....\n",
    "\n",
    "    print('Gas distribution after {:d} years is :\\t'.format(n * dt), new_distr)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Last example where we know the condition after X years and \n",
    "# can calculate back the concentration (solving the linear system for init_gas_distr)\n",
    "\n",
    "# Create matrix of coefficients:\n",
    "A = np.zeros((3, 3))\n",
    "A[0, :] = [0.8, 0.1, 0.04]\n",
    "A[1, :] = [0.15, 0.7, 0.13]\n",
    "A[2, :] = [0.05, 0.2, 0.83]\n",
    "\n",
    "# Create initial gas distribution\n",
    "gas_distr = np.array([1, 0, 0])\n",
    "\n",
    "# Calculate distribution at new time-step:\n",
    "time_step = 6\n",
    "new_distr = np.dot(np.linalg.???(A, time_step), gas_distr)\n",
    "print('Gas distribution after {:d} years is :\\t'.format(time_step), new_distr)\n",
    "\n",
    "# Use new distribution to calculate back the initial distribution (which method???):\n",
    "init_distr = np.linalg.???(matrix_prod_pow(A, time_step), new_distr)\n",
    "init_distr[init_distr<1e-14] = 0\n",
    "\n",
    "print('Initial gas distribution:\\t\\t', init_distr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# Plot gas distribution in each layer over 25 years time-period:\n",
    "num_steps = 25\n",
    "nr_time_steps = np.linspace(1, num_steps, num_steps, dtype=int)\n",
    "num_layers = 3\n",
    "gas_distr_mat = np.zeros((num_steps, num_layers))\n",
    "\n",
    "# your code here .....\n",
    "\n",
    "font = {'family': 'serif',\n",
    "        'color':  'darkred',\n",
    "        'weight': 'normal',\n",
    "        'size': 16,\n",
    "        }\n",
    "\n",
    "# Plot solution:\n",
    "plt.figure(num=None, figsize=(15, 5), dpi=80, facecolor='w', edgecolor='k')\n",
    "\n",
    "plt.plot(nr_time_steps, gas_distr_mat[:, 0], 'r', linewidth=2)\n",
    "plt.plot(nr_time_steps, gas_distr_mat[:, 1], 'blue', linewidth=2)\n",
    "plt.plot(nr_time_steps, gas_distr_mat[:, 2], 'black', linewidth=2)\n",
    "plt.title('Gas distribution over time', fontdict=font)\n",
    "plt.xlabel('time [years]', fontdict=font)\n",
    "plt.ylabel('gas concentration [-]', fontdict=font)\n",
    "plt.gca().legend(('layer 1','layer 2', 'layer 3'))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}