\documentclass[../book.tex]{subfiles}
\graphicspath{{\subfix{.}}}
\begin{document}

\chapter{Interpolating Between Data Points}
\label{12Interp}

\setcounter{section}{-1}  % Start subsections from 0

\input{prereqs/12Interp_prereqs}

\newpage
\section{Problem}
\setcounter{subsection}{-1}  % Start subsections from 0

One evening, you are chatting with a friend in the lab next to yours while you're waiting for a reaction to finish.  Your friend complains that they have a problem with data analysis. A collaborator sent them a sample of nanoparticles for your friend to added a new ligand to.  Your friend explains that they characterized the resulting particles using Raman spectroscopy and really needs to plot the difference spectrum for a paper, so that they can highlight the changes.


%has a Raman spectrum that he had recently taken containing an apparent mixture of the new nanoparticles he had synthesized and some byproducts that were previously characterized in the group.  He wanted to get a clean Raman spectrum of the new nanoparticles for a paper, but because the spectra partially overlapped, he was stuck.

``Normally, subtracting two spectra would be simple,'' your friend says, ``but the spectrum our collaborator took used a 532 nm laser and our instrument uses a 785 nm laser and different gratings. So all of the wavenumber values are different!  Also, I used up all the original nanoparticles in my synthesis and because nanoparticle syntheses are so variable, I'll never get the exact same particles as the starting ones! Like, maybe we can get ones that are close, but it would be super nice to just subtract the spectra I have already---it would be more correct and take less time and money!''

You think that this seems like a problem that could easily be solved with Python.  ``I bet we could write a quick program that uses the old data to calculate what the data points should be for each wavenumber value in the new data, and then subtract that from the new data.''

``That would definitely make my month!''

\section{Background}
\setcounter{subsection}{-1}  % Start subsections from 0

But where to start?  The simplest thing to do would be just shifting the $x$-value of each data point in the old spectrum to the nearest $x$-value of the new data, but this doesn't seem right to you.  For instance, it could easily be the case that the locations \textit{and} spacing between points is different, so that there won't be an even shift in the points, and some points may have two options of where to shift! Sketching out the problem, you get something like Figure \ref{fig:interpProblem}.

\begin{figure}[h]
    \centering
    \includegraphics[width=4in]{ChapNarrative/12Interp/Figures/TwoGaussians.png}
    \caption{A plot of two sets of data that have different locations of the data (circles) and different spacing between data points. Even though they both appear to be Gaussian lineshapes, there is no way to subtract one from the other, since the data points never line up. Note: the dotted lines are not data, they simply connect the data points.}
    \label{fig:interpProblem}
\end{figure}


You realize the core problem: \textit{All the data you have is at the points you collected...}  What happens in between them?  If you could figure this out (or take reasonable guesses) maybe you could subtract that behavior.

Thinking more about it, you consider that a spectrum is a continuously-varying thing with values that change continuously over a continuous set of wavelengths or energies.  How can you deal with such a thing in a computer?  We could keep adding more and more points, but we would never reach the continuous limit, our computer would run out of storage space, and our programs would take longer and longer to run.  So, the data that you are used to working with is a discrete set of points that will \textit{represent} these continuously changing phenomena.  You go down a rabbit hole on signals and find that this process is known as sampling; each $(x,y)$ pair (or $(t,y)$, or whatever) is like a measurement of the value of the continuous phenomenon, much like a time series that you would record in an experiment.  When you sample it, you create a discretized function.  Thinking about it a little more, you realize that this situation is the opposite: you are starting with a discretized function and trying to figure out the continuous phenomenon that it samples.  If you could propose a functional form, then you could fit it, but what could be the functional form be for a Raman spectrum?  Many of the peaks look asymmetric, what function could you use?  This seems like a problem, so maybe trying to reproduce the data is the wrong approach.  Maybe you can do what you would just do visually - look at how the data is changing point by point and try to fill in the gaps.  For instance, in Figure \ref{fig:interpProblem}, the points are drawn connected by straight lines.  Could you use that? Or perhaps another function?

A quick search for ``estimate value of data between data points'' quickly makes it clear that what you want to do is called interpolation.  Digging into interpolation, you unsurprisingly find that there are a number of different approaches from which to choose: cubic splines, linear interpolation, constant interpolation, and others.  In fact, this is reminding you a lot of the situation you found with numerical integration.  There, you could link each method with the corresponding order of complexity in a Taylor series, with each order (in principle) giving you a more accurate estimate.  You remember that the whole point of the Taylor series, when you first about learned it, was that you could use the values and the rates of change of a function at a known point to estimate the value of that function at an unknown point, which is exactly what you want to do here.  So, linear interpolation is just using the first order of the Taylor series, a line with constant slope, to estimate the value of the function between the data points.  A cubic spline uses a cubic equation, the third order, to do the same, which allows it to take into account some of the curvature of the function between the points.

You think, ``this seems like the kind of thing that Numpy or Scipy would be able to do.''  Sure enough, you find that Numpy has a 1d interpolation \cite{function}, \lstinline{interp()}, and Scipy has a whole \cite{module}, \lstinline{scipy.interpolate}.  

You decide to write a quick \cite{script} to see how they work!''

\section{Solution}
\setcounter{subsection}{-1}  % Start subsections from 0

You sketch out a plan to test the interpolation approaches that you want to try:
\begin{enumerate}
    \item Pick an analytical function that will give you exact values at any point, so you can test the accuracy of different interpolation methods.  Peaks in spectra are often well described by Gaussian functions, so that seems like a good place to start.
    \item Pick a set of $x$ values to operate as the ``data.''  These should be relatively sparsely spaced.
    \item Pick a different set of $x$ values to operate as the points at which you would like interpolated data.
    \item Do the interpolation.
    \item Plot the analytical function, the data, the interpolated data, and the difference.  Compute the Root-Mean-Squared Error (RMSE) between the interpolated data and the analytical curve.
\end{enumerate}

\subsection{Implementing linear interpolation using \lstinline{numpy.interp}}
\label{Implementing_linear_interpolation_using_numpy.interp}

The only part of this that is new to you is actually doing the interpolation.  Inspecting the documentation for \lstinline{numpy.interp()}, you see that this is simply a \cite{function} that takes as \cite{positional_argument}s the $x$-values at which you want interpolated points and the $x$- and $y$-values of the data from which you want to interpolate.\footnote{For more information on interpolation methods, check out books like Hildebrand, F. B. \textit{Introduction to Numerical Analysis}; Dover Press, 1987, or Eric Sullivan's \textit{Numerical Methods: An Inquiry-Based Approach With Python} (\href{https://numericalmethodssullivan.github.io/}{https://numericalmethodssullivan.github.io}).}  This function simply returns an array that holds the $y$-values for any $x$-value where you want an interpolation to occur.  

\begin{pyin}
    import numpy as np
\end{pyin}
\begin{pyin}
    y_interp = np.interp(x_interp, x_data, y_data)
\end{pyin}
assuming that you already have the arrays \lstinline{x_data} and \lstinline{y_data}.

\subsection{Implementing cubic spline interpolation using \lstinline{scipy.interpolate.CubicSpline}}
\label{Implementing_cubic_spline_interpolation_using_scipy.interpolate.CubicSpline}

Looking at the documentation for \lstinline{scipy.interpolate.CubicSpline}, you find that the situation is more complex.  \lstinline{CubicSpline} is a \cite{class} and the constructor function only takes the $x$- and $y$-values of the data from which you will interpolate---it does not need to know the $x$-values at which you would like interpolated $y$-values.  All that it is doing in the background is computing the coefficients of all of the cubic splines.  Once an instance of the \lstinline{CubicSpline} class is created, the resulting \cite{object} has a default \cite{method} that takes a \cite{positional_argument} that is the $x$-value(s) at which you would like interpolated $y$-values.  So, you figure the usage is two steps instead of the single step for \lstinline{numpy.interp}:

\begin{pyin}
    import scipy.interpolate as spi
\end{pyin}
\begin{pyin}
    cs = spi.CubicSpline(x_data, y_data)
\end{pyin}
\begin{pyin}
    y_interp = cs(x_interp)
\end{pyin}

Now you are ready to make a test \cite{script} to see how well these approaches do.  It turns out to be pretty simple.  You want to:
\begin{enumerate}
    \item Define an analytical function so that you know the exact $y$-value of the function at any $x$-value.  Since the Gaussian is so prevalent in science, you will just use it.
    \item Create some sparsely-separated data points to serve as the ``fake" data that you will use interpolate from.
    \item Create the exact $x$- and $y$-values for the points that you would like to interpolate to.  This will be a good way to compare the accuracy of the interpolated points.
    \item Use linear and cubic spline interpolation to calculate interoplated data points at each exact $x$-value.
    \calculate the root-mean-squared error between the interpolated and exact data.
\end{enumerate}
With this plan in hand, you quickly bang out a script to play with.

\lstinputlisting[language=Python]{../../Distribution/Code by Chapter/Ch11 - Interpolating Between Data Points/InterpExample.py}

\noindent which gives you the output:

\begin{multiout}
    Linear Interpolation RMSE: 0.73
     Cubic Interpolation RMSE: 0.13
\end{multiout}

Clearly the cubic spline interpolation is more accurate, but you want to really see what is going on, so you create some plots:
\begin{enumerate}
    \item A plot of the exact analytical curve, the data points, and the interpolated curve.
    \item The error in the interpolated y-values, defined as $y_{error} = y_{interp} - y_{exact}$.
    \item The relative error in the interpolated y-values, defined as $y_{relerror} = y_{interp} - y_{exact} / y_{exact}$.
\end{enumerate}

The plots for linear and cubic spline interpolation are given in Figure \ref{fig:Interp:InterpExamples}.

\begin{figure}[H]
    \centering
    \begin{tabular}{ccc}
        \includegraphics[width=6in]{ChapNarrative/12Interp/Figures/InterpExamples.png}
    \end{tabular}
    \caption{An example of the performance of (a) linear and (b) cubic spline interpolation for a generated Gaussian peak without noise, compared to the effects of noise on the same methods (c and d, respectively).}
    \label{fig:Interp:InterpExamples}
\end{figure}

Again, the cubic spline is more accurate than the linear interpolation.  So you wonder why anyone would use linear interpolation, but then think, maybe it doesn't handle noise as well.  So, you add noise to your simulation and get the following results:

\begin{multiout}
    Linear Interpolation RMSE: 5.48
     Cubic Interpolation RMSE: 6.46
\end{multiout}

In this case, the linear interpolation approach is slightly less affected by noise, and would probably be preferable.  Plotting the linear and cubic spline interpolation results similarly for the same Gaussian peak with noise, shown on the right in Figure \ref{fig:Interp:InterpExamples}, shows a the same story, with the cubic spline interpolation tending to be off by more in the areas between data points despite being off by the same amount at the data points.  Based on this, you figure that cubic splines will generally do better if you don't have noise and if your data is reasonably closely spaced---basically if you have a few data points on each side of your peak.

With the different approaches validated, writing the script for your friend will be quick.  All you need to do is read both data files, interpolate one spectrum using the $x$-values of the other, and save the interpolated spectrum.  You reason that it would be better to use the data file with more closely spaced points as the data, and the $x$-values of the less closely spaced data as the points at which to interpolate.  

%The only wrinkle is which wavenumber values to interpolate to.  Looking at the spectra, you see that both start at 2000 cm$^{-1}$ but that the 532 nm spectrum does not go as low in wavenumbers as the 785 nm spectrum does.  Since you can't interpolate the 532 nm spectrum at all 785 nm spectrum values, but you can do the opposite, you decide to interplate the 785 nm spectrum to the $x$-values of the 532 nm spectrum.  



\section{Code}
\setcounter{subsection}{-1}  % Start subsections from 0

\begin{figure}[H]
    \centering
    \includegraphics[]{ChapNarrative/12Interp/Figures/FlowChartCh11.pdf}
    \caption{Flowchart for code that will take two spectral data files with different x-values, interpolate one, and then subtract one spectrum from the other.}
    \label{fig:Interp:ActivityDiagram}
\end{figure}

\lstinputlisting[language=Python]{../../Distribution/Code by Chapter/Ch11 - Interpolating Between Data Points/RamanSubtraction.py}

You run this script, select the folder ``Nanoparticles" in your directory, and it produces following output:

\begin{figure}[H]
    \centering
    \includegraphics[width=3in]{Figures/raman.png}
    \caption{The results of interpolation and subtraction of the spectra.}
    \label{fig:Interp:Raman}
\end{figure}

You show this result to your friend, who says, ``That was so easy for you!  I couldn't figure out any way to do this with any of the programs that our lab uses.  Now we can be really quantitative about the effects of the ligand exchange reaction we did.  Can I use this in my paper?  We will totally plug you in the acknowledgments!"

\subsection{Line-by-line explanation}

\begin{footnotesize}
\begin{itemize}
    \item [8--11:] We need \lstinline{numpy}, \lstinline{make_subplots()}, \lstinline{quickSelectFolder()} and \lstinline{quickPopupMessage()} from \lstinline{codechembook.quickTools}, and the wavenumber symbol from the \lstinline{codechembook.symbols.chem} module.
    \item [14:] The two spectra have different $y$-scaling, so we need a normalization factor before we subtract.  Put it here so that it is easy to find and change.
    \item [17\&18:] Get the path to the folder that holds the two Raman spectrum csv files.
    \item [21--24:] Read both files and store them in separate variables.
    \item [27:] Do the actual interpolation.  Here \lstinline{x532} is the $x$-values at which we want interpolated $y$-values, and \lstinline{x785} and \lstinline{y785} are the $x$- and $y$-values of the data we want to interpolate.  We need to interpolate the 785 nm spectrum onto the $x$-values of the 532 nm spectrum because the 785 nm spectrum has $x$-values that fall outside of the range of the 532 nm spectrum and thus can not be interpolated!
    \item [30]: Here we subtract the 532 nm spectrum from the interpolated 785 nm spectrum.  We need a normalization factor for the 532 nm spectrum because it is lower in intensity than the 785 nm spectrum.  You can try to make this fancy or you can just guess and check a few times to get very close.
    \item [33--45:] Nothing fancy here, we're just plotting the two spectra in one subplot and the difference between them in another.
\end{itemize}
\end{footnotesize}

\hrule
\input{Exercises/Chapter 11 Interpolation}
\hrule

\end{document}
