Lesson_5_20200410.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving Equations and Curve Fitting.\n",
    "\n",
    "This week we will be solving equations and fitting data to curves.  Both of these procedures rely on linear algebra techniques to efficiently solve problems numerically.\n",
    "\n",
    "## Solving systems of linear equations\n",
    "\n",
    "Sometimes physics problems can involve solving many, many equations simultaneously.  This can be done analytically by hand using linear algebra, but we will be using the `numpy.linalg` library to solve a problem quickly.  Consider this system of coupled equations:\n",
    "\n",
    "$2 x_1 + x_2 = 0$\n",
    "\n",
    "$-5 x_1 + 3 x_2 = 10$\n",
    "\n",
    "We can write this as a single matrix equation and then use matrix methods to calculate a solution for both $x_1$ and $x_2$:\n",
    "\n",
    "$\\begin{pmatrix} 2 & 1 \\\\ -5 & 3 \\end{pmatrix}\\begin{pmatrix} x_1 \\\\ x_2 \\end{pmatrix} = \\begin{pmatrix} 0 \\\\ 10 \\end{pmatrix}$\n",
    "\n",
    "The example below demonstrates the use of numpy.linalg to solve for the two $x$ values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x1 = 1.0\n",
      "x2 = 4.0\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def solve_two_equations(a1, a2, b1, b2, z1, z2):\n",
    "    '''Solves two equations of the form a1 * x1 + a2 * x2 = z1 and b1 * x1 + b2 * x2 = z2'''\n",
    "    left = [[a1, a2], [b1, b2]] #matrix\n",
    "    right = [z1, z2] #vector\n",
    "    return np.linalg.solve(left, right)\n",
    "\n",
    "[x1, x2] = solve_two_equations(2, 1, -5, 3, 6, 7)\n",
    "print('x1 =', x1)\n",
    "print('x2 =', x2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 1\n",
    "From the circuit below, we can write down a set of 6 linear equations for the currents $i1$ through $i6$ in the Wheatstone bridge circuit using Ohm's Law, $V = IR$, and use Kirchoff's Laws to write the equations.  Kirchoff's Laws state that all currents entering a point in a circuit must add up to zero, and all voltages in a loop must also add up to zero.\n",
    "\n",
    "<img src=\"Bridge Circuit 6.png\" style=\"width: 200px\">\n",
    "\n",
    "Here is an example of six independent linear equations:\n",
    "\n",
    "$i6 - i1 - i2 = 0$ (sum of currents at point at top of circuit)\n",
    "\n",
    "$i2 \\times R2 - i3 \\times R3 - i1 \\times R1 = 0$ (sum of voltages around top loop of the bridge)\n",
    "\n",
    "$i1 - i3 - i4 = 0$\n",
    "\n",
    "$i4 + i5 - i6 = 0$\n",
    "\n",
    "$i3 \\times R3 - i4 \\times R4 + i5 \\times R5 = 0$\n",
    "\n",
    "$i1 \\times R1 + i4 \\times R4 = V$\n",
    "\n",
    "Using these equations you will create a function that takes the values $R1$–$R6$ and $V$ as arguments and generates a 6x6 `numpy` array (**M**) from your equations and an 6x1 array with the constants (**r**), such that your matrix equation is **M i = r**. \n",
    "\n",
    "To represent a system of equations to `numpy`, remember that it is the *coefficients* of the variables that distinguish one equation from another.  So our 6x6 matrix of equations is really a 6x6 matrix of coefficients, with each row representing one equation, and each position in a row implicitly associated with one of the independent variables.  In our case, our variables are $i1$–$i6$, so as an example, the second equation above would be passed to `numpy` as `[-R1, R2, -R3, 0, 0, 0]`.  Make sure you understand why that is before moving on, and write out your matrix of coefficients _on paper_ before writing your function\n",
    "\n",
    "\n",
    "Have your function use the `numpy.linalg.solve()` function solve your equations for the i's.  Call your function with R1 = 6 ohm, R2 = 4 ohm, R3 = 18 ohm, R4 = 7 ohm, R5 = 3 ohm and V = 12 V.  Print the full solution, as well as specifically printing the current and voltage across R3.  (Any printing of output like this should also print some explanatory text along with it, so that it's not necessary to look at the details of the code in order to understand the output.) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitting functions to data\n",
    "\n",
    "Often in physics we will have some experimental data, and we would like to model the data with some particular mathematical function.  For example, we know that a ball thrown in the air should follow a parabolic path as a function of time: $y = y_0 + v_0 t + \\frac{1}{2}g t^2$.  Thus, if we have a bunch of $(t_i, y_i)$ data points, we could fit the aforementioned function to the data, with $y_0$ and $v_0$ as free parameters, and thus use the data to determine how fast the ball was thrown and from what initial height. \n",
    "\n",
    "In other cases, we might have a dataset, but we might not know *a priori* what functional form is the best one to fit the data, so we might try a few different options.  Or we might have two competing theories that predict different functional forms for the data, and we might want to see which one in fact ends up giving a better fit to the data.\n",
    "\n",
    "In many computer languages (python included), doing this sort of fitting proceeds in two steps:\n",
    "\n",
    "* Write a function that takes the independent variable as input (e.g. as a `numpy` array), along with any values for the function parameters, and returns an array of corresponding values of the fitting function at each of the input data points. \n",
    "\n",
    "* Pass that function, along with initial guesses for the parameters to be fit, to a routine that finds the \"best fit\" values of the fitting parameters. \n",
    "\n",
    "For now we will leave what we mean by \"best fit\" vague, and work through a simple example of this sort of fitting, using the `scipy` fucntion `scipy.optimize.curve_fit()`. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 2\n",
    "\n",
    "Write a function that takes an array of time values as input, along with the parameters $y_0$ and $v_0$, and returns an array of the $y$ position values, following the equation above for $y(t)$.  So that it will be usable with `scipy.optimize.curve_fit()`, make sure that the time values (independent variable) is the first parameter in your function declaration.\n",
    "\n",
    "Now generate a set of time values (using, e.g., `np.linspace`), and pass those to your function, returning the corresponding $y$ values.  Graph the result.  Make the time values closely spaced enough so that the curve looks relatively smooth. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 3\n",
    "\n",
    "Using the same array of time values and same function from above, generate a \"noisy\" dataset using `numpy.random.normal` like last week.  That is, generate values from your function (or just use the ones created in the last step), and then add random noise to the values and store the results in a new array.  Plot these noisy values (as points, not connected by a line).  If necessary, adjust the noise level so that the overall shape of the curve is apparent, but that there is clearly some scatter in the data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 4\n",
    "\n",
    "Now use `scipy.optimize.curve_fit()` to find the best fit to your \"noisy\" datapoints generated in the previous exercise.  Looking at the documentation for that will show you an example of its usage, but briefly, you pass it the name of your function, the x array (independent variable), and the y array (dependent variable, i.e. function values), and it returns an array of the parameters that provide the best fit, along with an estimate of the *covariance matrix*.  The latter is an estimate of the uncertainty in the fitted parameters.  The diagonal elements can be thought of as uncertainties on the individual parameters in isolation, and the off-diagonal elements provide an indication of how linked each pair of parameters is, i.e. how changing the value of one would influence the value of the other. \n",
    "\n",
    "Repeat your plot from the previous exercise, showing the noisy data as points.  Then overplot a line showing the best-fit curve based on the parameters returned by `curve_fit`.  Note that you can calculate this best-fit line easily using your original function, just by passing it the same independent-variable array, and the parameters returned from `curve_fit`, which will be in the same order as needed for your function.  Also print out the values of the best-fit parameters, and comment on whether they are close to the values you used to generate the noisy data in the first place.  (They should be!) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitting lines from a spectrum, and comparing fits\n",
    "\n",
    "Now let's try a somewhat more useful example - fitting the lines in a spectrum.  We'll load in data for a real spectrum of the mineral magnesite, and fit some of the lines.  We'll compare fits to see which of two possible functional forms is a better fit to the line profile. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First let's just get the data and see what we're working with.  The `numpy` routine `genfromtxt` makes it easy to read data from a text file, which is something you'll often want to do when working with experimental data.  When you read in data in this way, the routine returns an n_columns x n_rows 2D `numpy` array.  It's often simplest to go ahead and assign columns from this to separate variable names right away, e.g. if you read in a two-column datafile into array $A$, you may then want to do `x = A[:, 0]` and `y = A[:, 1]`.  (Note that in `numpy` arrays, the *second* dimension refers to columns.  So that line `x = A[:, 0]` means \"take every row from A, but only column 0\".)\n",
    "\n",
    "#### Exercise 5\n",
    "\n",
    "Import the Raman spectrum for the mineral magnesite from the file Magnesite.txt using the `np.genfromtxt` function as below and then plot the entire spectrum.  The x axis is in units of cm$^{-1}$.  (Inverse centimeters are called \"wavenumbers\", and are a unit of frequency of light - they are more common in chemistry than physics, but you'll sometimes see them crop up.)\n",
    "\n",
    "Once you've plotted the whole spectrum, select a subset of data around one of the strong spectral lines, and plot only that line - that's a line we'll use for fitting in the next step.  Save the x and y values for this line in separate arrays for use next.\n",
    "\n",
    "A nice way to do this is to use logical conditions to index your data: \n",
    "\n",
    "`indices_in_range = np.logical_and(x > xmin, x < xmax)`\n",
    "\n",
    "`data_i_want = data[indices_in_range]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#get data from file and parse into physical variables\n",
    "def get_spectrum_data(filename):\n",
    "    \"\"\"Returns intensity versus wavenumber data from text file in current directory\n",
    "    \n",
    "    INPUTS: Filename as a string.\n",
    "    \n",
    "    OUTPUTS: \n",
    "    WAVENUMBER is an array of the wavenumber values at which the intensity is measured.\n",
    "    INTENSITY is an array of the intensity values at each wavenumber.\n",
    "    \"\"\"\n",
    "    \n",
    "    data = np.genfromtxt('Magnesite.txt')\n",
    "\n",
    "    wavenumber = data[:,0]\n",
    "    intensity = data[:, 1]\n",
    "    \n",
    "    return wavenumber, intensity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  Exercise 6\n",
    "\n",
    "Now let's try fitting one of the lines in this spectrum.  \n",
    "\n",
    "First we'll try fitting the line with a Gaussian function.  In the cell below is an example function that takes as arguments the central frequency, peak amplitude, and width of the Gaussian, as well as a constant additive offset (since we can see above that the background level doesn't go all the way down to zero).  It also needs to take as its first argument the independent variable, the frequency.  Use `scipy.optimize.curve_fit()` to fit the Gaussian.  You may need to make an initial guess using the `p0=` argument (look at the documentation for `curve_fit` as necessary).\n",
    "\n",
    "Before running the fitting it's a good idea to just calculate and plot the function, to make sure that it is behaving as you expect.  Once it looks good with you manually specifying the parameters, you can go ahead and run the fit to find the best-fitting parameters.  This manual running of your function is also a good way to determine some good initial guesses for the parameters, which can be helpful to allow the fit to converge quickly to the optimal solution. \n",
    "\n",
    "Print the fit parameters and plot the data and its fit, following the earlier example. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Gaussian function with four fit parameters: center frequency, overall amplitude, width parameter, constant vertical offset.\n",
    "def gaussian(freq, center_freq, amp, width, const):\n",
    "    return amp*np.exp(-(freq-center_freq)**2/(2*width**2)) + const\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 7\n",
    "\n",
    "Now do the same exercise as above but define a function to calculate a Lorentzian \n",
    "line profile, which has the functional form\n",
    "\n",
    "$L = \\frac{1}{2\\pi} \\frac{w}{(\\nu - \\nu_0)^2 + (w/2)^2}$\n",
    "\n",
    "where $w$ is the width of the line and $\\nu_0$ is its central frequency.   As before, allow a multiplicative\n",
    "scaling amplitude and a constant additive offset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Goodness of fit\n",
    "\n",
    "Now, which of these functions is more appropriate for our dataset?  A Lorentzian profile is expected theoretically if the main line broadening mechanism is the energy uncertainty from the Heisenberg uncertainty principle.   A Gaussian profile is expected if the main thing broadending the line is Doppler shift due to the random motions of individual atoms.  So determining the best-fit profile can tell us something about the physical situation in the gas. \n",
    "\n",
    "Your eye can probably tell you already which is the better fit, but it's worth using a quantitative tool. The most common way to do this is by calculating a statistic called $\\chi^2$ (pronounced \"chi squared\").  There are many different definitions of this statistic, beginning with just summing up the squared deviations between the data and fit:\n",
    "\n",
    "$\\chi_{basic}^2 = \\Sigma (y_{data,i} - y_{fit,i})^2$\n",
    "\n",
    "where the sum is over all the data points.  A smaller $\\chi^2$ indicates a better fit.\n",
    "\n",
    "More commonly, the statistic is weighted: each $y_i$ data point has an uncertainty $\\sigma_i$ and we use these to weight the sum:\n",
    "\n",
    "$\\chi^2 = \\Sigma \\frac{(y_{data,i} - y_{fit,i})^2}{\\sigma_i^2}.$\n",
    "\n",
    "The effect of the weighting is that the sum is affected less by points far from the expected value if they have large uncertainties.\n",
    "\n",
    "In my anecdotal experience, the reduced $\\chi^2$ defined in Taylor is the most commonly used.  This statistic is:\n",
    "\n",
    "$\\tilde{\\chi}^2 = \\chi^2/d,$\n",
    "\n",
    "where $d$ represents the number of degrees of freedom in the fit or comparison.\n",
    "\n",
    "#### Exercise 8\n",
    "\n",
    "Write a function to calculate the reduced $\\chi^2$ and call it to print the $\\tilde{\\chi}^2$ value for each of your two fits above.  Use an uncertainty of 0.02 for the measured intensity values (arbitrary units).  Which is a better fit, a Lorentzian or a Gaussian?  Which physical mechanism seems dominant for generating the shape of this spectral line?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 9 (optional)\n",
    "\n",
    "Now that we have a good fit for our data, what would the signal (which is currently in the frequency domain) look like in the time domain?  Take the inverse Fourier transform of your fit to the Raman data and plot it. (Don't worry about getting the x axis units correct; to do that, we'd first have to convert from cm$^{-1}$ to SI units.)  In the case of Magnesite, the 1095 cm$^{-1}$ and 735 cm$^{-1}$ peaks are associated with stretching and bending of the carbonate group in Magnesite, and the other two main peaks are associated with vibrational modes of the lattice.  The time domain shows you how these vibrational modes behave, and you can measure this behavior using time-domain techniques like pump-probe spectroscopy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}