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

\chapter{Finding Peaks}
\label{13PeakFind}

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

\input{prereqs/13PeakFind_prereqs}

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

You need to do a series of cyclic voltammetry experiments on your catalysts under a number of different conditions.  You will need to evaluate different scan rates, different catalyst derivatives, and different substrates, adding up to more than a \textit{hundred} individual experiments.  You want to figure out the $E_{1/2}$ for each of them, but with so many experiments, that will take a lot of slow, hand-on data analysis. You are not looking forward to quantifying all these results by hand, so you start thinking about how to automate it using your growing coding skills.

\begin{figure}[H]
    \centering
    \includegraphics[width=4in]{ChapNarrative/13PeakFind/Figures/CVexample.png}
    \caption{A diagram of your task---you need to figure out the value of $E_{1/2}$, the potential halfway between the andodic and cathodic peak potentials.}
    \label{fig:PeakFind:CVexample}
\end{figure}

% You need to figure out three parameters from the CV:
% \begin{enumerate}
%     \item The formal potential $E_{1/2}$, which is the average of the potentials of the cathodic and anodic peaks ($E_{pc}$ and $E_{pa}$, respectively).
%     \item The anodic peak current $i_{pa}$, which is the difference between the current at the anodic peak and the projection of the slope of the current prior to the reduction wave.
%     \tiem The cathodic peak current $i_{pc}$, which is the difference between the current at the cathodic peak and the projection of the slope of the current prior to the oxidation wave.
% \end{enumerate}

% Certainly the first thing you need to do is identify the peaks programmatically.  Once you have their potential values, you can easily calculate $E_{1/2}$.  Projecting the slopes will be a little trickier but you can either fit them or use the numerical derivatives we learned in the last chapter to extrapolate them.  So how to find the peaks?

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

You think about how you would find a peak just by looking at a cyclic voltammegram (CV)---you do it without even really thinking about it.  It's just the highest data point in a range of data, or in other words, the data point that has a higher value than any other point (at least in its vicinity).  So, it seems like the easiest thing to do is just finding the index of the \textit{global} maximum potential value in the reduction wave, and the minimum potential value in the oxidation wave, and using it to index the potential and current arrays.  Thinking about this in very general terms, you draw something with peaks and valleys (Figure \ref{fig:PeakFind:Prominence}) and immediately realize this would only find the leftmost peak in your drawing.

\begin{figure}[H]
    \centering
    \includegraphics[width=4in]{ChapNarrative/13PeakFind/Figures/Prominence.png}
    \caption{(Left) An example of how prominence is determined for a given peak.  (Right) The prominence of a series of peaks.}
    \label{fig:PeakFind:Prominence}
\end{figure}

How to handle multiple peaks?  Then finding the maximum or minimum is will only find one of them.  Time for a more sophisticated approach that finds \textit{local} maxima and minima.  Again, as a human you can usually easily spot them, but how do you get a computer to do it?

You realize that there are a few more possible problems. Your catalyst is not very soluble in all the solvents you will use in your experiment. This means that the resulting data could be fairly noisy in those solvents.  Thus, the local maximum is likely not going to be the actual maximum of the underlying signal, but rather a point close to the maximum of the underlying signal that just happens to have a particularly large contribution from noise. This is like what is shown in Figure \ref{fig:PeakFind:ProminenceNoise}.  What can you do to mitigate these problems?

\begin{figure}[H]
    \centering
    \includegraphics[width=6in]{ChapNarrative/13PeakFind/Figures/ProminenceNoise.png}
    \caption{Examples of peak filtering with difference prominence requirements.}
    \label{fig:PeakFind:ProminenceNoise}
\end{figure}


You could fit the peak to a function that has the peak center as a parameter, if your peak follows a well-behaved analytical function or if you just fit part of the peak (say, to a parabola).  However, this would probably require a lot of user intervention and it also seems like it might take a bit of tuning of parameters, so maybe there is a more generic method.

Another approach could be to simply find every point that is higher than the points adjacent to it.  This would pick way \textit{more} peaks than you want, but you could filter out the ones you don't want with a number of different filters until you have narrowed your list down to only the ones you want. Talking with \first, you learn there are a number of different filters that you could imagine applying that would cover a range of peaks and a range of different signal-to-noise ratios.  They point you to the \lstinline{scipy.signal} signal module.

%By now you have learned to think back to calculus class.  You remembered learning that you can find maxima and minima of a function by finding the zeroes of its derivative.  You can then use the sign of the second derivative at each zero to determine if it is a maximum or a minimum.  This seems like a general approach that could work well, too.  This could either be done by smoothing first to get rid of noise, then taking the derivative, or by using the derivatives of the smoothing functions themselves (for instance, you can get them for free from \lstinline{scipy.signal.savgol_filter}).  

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

\subsection{Finding local maxima using \lstinline{scipy.signal.find_peaks()}}
\label{Finding_local_maxima_using_scipy.signal.find_peaks}
In \lstinline{scipy.signal} you find a function \lstinline{find_peaks()}, which is a complex function that appears to do what you want.  It seems to just go through your data and look for points that meet several criteria, the most basic of which is that its value is higher than those of its nearest neighbors.  It will return to you an array of \cite{indices} of the peaks.  However, since this could find way too many peaks if your data is noisy, there are other more complex criteria by which you can filter out most or all of the spurious peaks that aren't real or that you don't actually want. You can filter by the $y$-value a number of different ways using a few options: \lstinline{height} is the minimum $y$-value (or a range of $y$-values if a two-element array is specified) of allowable peak heights, \lstinline{threshold} is similar but specifically refers how much higher than its nearest neighbors it must be, and \lstinline{prominence} is the difference in $y$-values between the peak height and the highest valley between it and the next nearest peak.  


So, if you know that the peak you want is much higher than any neighboring peaks, you can use \lstinline{threshold} to filter out baseline noise or \lstinline{prominence} to filter out noise that might be causing multiple peaks to appear near the top of the peak that you want.  \lstinline{prominence} seems a particularly useful parameter in many chemistry contexts.


You can also filter on properties related to the $x$-axis, with \lstinline{distance} setting the minimum number of points between neighboring peaks and \lstinline{width} setting the minimum (if a single value is given) or a range (if two values are given) of peak widths that are acceptable.  \lstinline{width} is particularly useful if you have an expectation of the widths of your peaks---often noise is much narrower than signal features.



\subsection{Looping until some conditions are met with \lstinline{while} loops}
\label{Looping_until_some_conditions_are_met_with_while_loops}

The above seems like a solid plan, but you wonder how to handle the fact that the cathodic trace is going to have a negative peak? You realize that you can just take the cathodic scan, and multiply the $y$-values by -1, so that you can treat them as positive, but this raises another question... how do you know when you are dealing with points on the cathodic scan? You would need a way to loop through all the values of potential in the scan, and then determine when you moved from the anodic to the cathodic scan.

You already know how to use \citett{for} \cite{loop}s, but the problem with them is that \lstinline{for} loops will run for exactly as many iterations as the length of the object that is being iterated over. This means you need to have an object that is exactly as long as the number of time you want the for loop to run.  Sometimes, there is no way to know how many times you need the loop to run ahead of time: either you are depending on input from a user, data from a source outside of the computer, or simply executing some analysis in which it is not clear how many iterations you will need. In the case of the CV experiment we are considering, you don't know how many iterations there will be until you move from the anodic to the cathodic scan. 

A \lstinline{while} loop addresses this case.  A \lstinline{while} loop will keep iterating until some specific logical expression evaluates to \lstinline{False}.  In Python, \lstinline{while} loops are written like this:

\begin{lstlisting}
while <conditional expression>:
    <code to execute repetitively>
<rest of script>
\end{lstlisting}

As long as \lstinline{<conditional expression>} evaluates to \lstinline{False}, \lstinline{<code to execute repetitively>} will continue to loop.  The while loop will terminate once \lstinline{<conditional expression>} evaluates to False.  This evaluation is performed only at the end of the loop.  If something changes in the middle of an iteration, the iteration will complete, and then the loop will exit to move on with the rest of the main code.  If the \lstinline{<conditional expression>} evaluates to False even before the first iteration, the code in the loop will never run at all.\footnote{See Chapter \ref{fig:FlowControl:While} for additional information.}

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

\begin{figure}[H]
    \centering
    \includegraphics[width=\textwidth]{ChapNarrative/13PeakFind/Figures/FlowChartCh12.pdf}
    \caption{Flowchart for code that reads in cyclic voltammetry data, splits the data into anodic and cathodic scans, and then identify the locations of peaks in the anodic and cathodic scans.}
    \label{fig:PeakFind:ActivityDiagram}
\end{figure}

\lstinputlisting[language=python]{../../Distribution/Code by Chapter/Ch12 - Finding Peaks/CVPeakFind.py}

You run this code, navigate to the eChem subfolder of your Catalyst folder, and select the file `CVdata.csv' to use as test data for this script.  In a flash, the following command line output appears:
\begin{multiout}
    Epa = -1.236 V
    Epc = -0.885 V
    delta E = 0.351 V, E1/2 = -1.061 V
\end{multiout}
along with a plot:
\begin{figure}
    \centering
    \includegraphics[width=3in]{ChapNarrative/13PeakFind/Figures/CVpeaks.png}
    \caption{The CV data you analyzed along with the peaks marked.}
    \label{fig:PeakFind:CVplot}
\end{figure}

This looks great!  You just saved yourself a lot of time opening and analyzing CV data by hand.


\subsection{Line-by-line explanation}

\begin{footnotesize}
\begin{enumerate}
    \item [8--13:] We need \lstinline{numpy} and \lstinline{make_subplots()} like usual, as well as \lstinline{scipy.signal} to get \lstinline{find_peaks()}.  We'll also import the capital delta symbol from the \lstinline{codechembook.symbols.greek} module and separately import the HTML formatting functions \lstinline{textsub()} and \lstinline{textit}.  We could just import the whole symbols module, but the code looks cleaner this way compared to typing out \lstinline{cs.typesettingHTML.textsub()} each time we want to use it.
    \item [16\&17:] Get the filename for the CV data csv file by notifying the user what to do and then opening the file dialog.
    \item [20:] Read the raw CV data.  This data has both the oxidation and reduction waves in it, so we will have to separate them next.
    \item [23--29:] We want separate arrays for the reduction and oxidation waves. First, we need to find the index at which the potential starts going positive (after going negative for the first part of the scan).  We do this by looking for where the difference between the last and the current potential switches signs:
    \begin{enumerate}
        \item [23:] We want to count data points so that we know the index.  First we need to initialize the count.  We set it at 1 instead of 0 because we are computing the direction of the change of the voltage by subtracting the current point from the point behind it, and there is no point behind the point at index 0, so we have to ignore it.
        \item [24:] We also initialize the stop condition to \lstinline{False}, since the loop will terminate when it becomes \lstinline{True}.
        \item [25--29:] Now we loop until we reach the end of the voltage array or the stop condition is no longer \lstinline{False}.  If we don't include the first condition, the counter could become greater than the length of the voltage array, causing an error.  If the voltage at this index is less than the voltage at the previous index, then we are still scanning to the negative direction, and we increment the index and continue the loop.  If this is not true, then we have stopped scanning in the negative direction and have found the turnaround point, so we set \lstinline{stop} to be True, which terminates the loop.
    \end{enumerate}
    \item [32:] Now we can split the CV data into reduction and oxidation waves.  We will slice from the beginning to the turnaround point to get the reduction wave.
    \item [33--36:] The oxidation wave is slightly more difficult.  It is possible that there are two points that are exactly the same at the turnaround: the last voltage for the negative scan and the first voltage for the positive scan.  We use an \lstinline{if-else} statement to test for this case and slice accordingly.
    \item [39\&40:] Now, we will use \lstinline{find_peaks()} on the reduction and oxidation waves, looking for any peak that has at least 10\% prominence (i.e., 10\% greater than any surrounding base).  Since \lstinline{find_peaks()} only takes the prominence in terms of absolute difference between the base and the peak, we can calculate this in terms of a fraction of the maximum value of the wave.  \lstinline{find_peaks()} returns the indices of the peaks as an \lstinline{ndarray} and a dictionary of the peak properties, so we need to accept both of them on the LHS of the assignment operator.
    \item [43\&44:] We don't want to know the indices of the peaks, just the E values, so we use the peak indices to index the data and return the E value of each peak.
    \item [47--49:] Print out the results to the command line.
    \item [52--71:] Plot the data with the found peaks marked on top.
\end{enumerate}
\end{footnotesize}

\hrule
\input{Exercises/Chapter 12 Peak Finding}
\hrule

\end{document}
