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

\chapter{Computing Numerical Derivatives}
\label{14NumDer}

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

\input{prereqs/14NumDer_prereqs}

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

You just wrote a \cite{script} that found the peaks of the reduction and oxidation waves of a CV experiment and used them to calculate the redox potential, $E_{1/2}$.  You expected your advisor to be impressed, but instead they only had a series of questions, the main one being how you can know if the redox event is reversible without \textit{looking} at the waves. 

After some discussion, you both come to the conclusion that you need to automatically calculate the peak oxidative and reductive currents. So now you need to figure out how to automatically extract these values and compare them to see if they're `close enough' to establish reversibility.  You know how to find the peak positions (Chapter \ref{13PeakFind}), but you cannot simply read out the $y$-values at these points. This is because CV scans have a capacitive charging component that makes a significant contribution to the observed currents.  You will need to figure out what this capacitive component is, and then subtract it from the peak currents to find the current due to the chemical redox event.\footnote{If you care to know more, we recommend the open-access tutorial article \textit{A Practical Beginner’s Guide to Cyclic Voltammetry} by Elgishi et al. in \textit{J. Chem. Ed.}}
%You are working up a promising electrocatalyst and need to do a series of CV measurements.  You want to create a script that automates the workup.  The first thing that you need to figure out is the capacitive charging component of the signal, which is the derivative of charge with respect to time.  Since all you have is data, you can do this by fitting, but you can also just do it by taking the derivative of the data.

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

The capacitive charging portion of the CV wave is typically treated as a linear extrapolation of the straight part of the scan at higher potential than the peak for oxidation and at lower potential than the peak for reduction.  Normally, people will fit a line to the linear regions and then use the fit to compute the capacitive contribution to the total current at $E_{pa}$ or $E_{pc}$.  Subtracting the capacitive contribution from the total current gives the electrochemical current.  For a reversible process, $I_{pa}$ and $I_{pc}$ should be equal.

\begin{figure}[H]
    \centering
    \includegraphics[width=3in]{ChapNarrative/14NumDer/Figures/CVpeaks_annotated.png}
    \caption{A diagram of your task---you need to figure out the slope of the capacitive parts of the CV and extrapolate them to the peak to determine $I_{pa}$ and $I_{pc}$.}
    \label{fig:NumDer:CVexample}
\end{figure}

It is not immediately apparent how to use Python to do this.  You can certainly use the lmfit \cite{library} to do the fitting, but how do you automatically figure out the right \textit{part} of the wave to fit, so that you can extrapolate for the capacitive components?

Since the redox potential of each sample will likely be somewhat different, and your sweep rates are likely to also be different, then you can't just always fit the same range of potentials for each file.  ``I need to figure out where in the wave the slope is basically constant,'' you mutter to no one in particular.  So, you need to calculate the slope and look for regions where it is essentially the same.  ``Better yet, I can look for when slope of the slope---the second derivative---is zero, then I know the slope is not changing!''  Now you need to figure out how to calculate derivatives of data.

Remembering back to your calculus class, you recall a number of rules for derivatives, like
\begin{align*}
    \diff{}{x}x^n &= nx^{n-1}\\
&\textrm{or}\\
    \diff{}{x}\cos{x} &= -\sin{x}.\\
\end{align*}
How do you apply these rules to numerical data?  Thinking back a little farther, you remember the definition of a derivative at any point $x$:
\begin{equation}
    \diff{}{x}f(x) = \lim_{h\to0}\frac{f(x+h) - f(x)}{h}
\end{equation}
In analytical math, you can take this limit, but with discrete data, there is a minimum spacing between points, and such, $h$ can only get as small as $x_{i+1}-x_{i}$, where $i$ is the index of an array of data points.  So, if you define $x$ as $x_i$, and replace $h$ with $x_{i+1}-x_{i}$, then you get
\begin{equation}
    \diff{}{x}f(x)|_{x_i} \approx \frac{f(x_{i+1}) - f(x_i)}{x_{i+1}-x_{i}}.
\end{equation}
This is just the old formula for the slope that you learned in algebra!
\begin{equation}
    m = \frac{y_2-y_1}{x_2-x_1}
\end{equation}
Is it really true that numerical differentiation is as easy as finding the slope?

By this definition, you can sketch out what your derivative process should look like, as shown in Figure \ref{fig:NumDer:ForwardDerivativeExample}.  While this makes sense initially, as you start to think about it, you notice some issues.  First, and most importantly, your sketch has seven data points but only six derivatives!  This makes sense, because the derivative connects every pair of points, but it seems pretty inconvenient.  Secondly, the derivative you calculate is the slope between two adjacent points, but this doesn't actually tell you the value of the derivative at any given point.  Perhaps this tells you the value of the derivative at some unknown point between the two data points, like in interpolation,\footnote{See Chapter \ref{12Interp}} but that isn't necessarily what you want.  
\begin{figure}
    \centering
    \includegraphics[width=2in]{ChapNarrative/14NumDer/Figures/ForwardDerivativeExamples.png}
    \caption{A basic sketch of numerical derivatives, in this case, the ``forward derivative.''}
    \label{fig:NumDer:ForwardDerivativeExample}
\end{figure}

``Maybe I can average the derivative before the point with the derivative after the point to get a better estimate of the derivative at the point?''  Working out the math, you reason:
\begin{equation}
    \diff{}{x}f(x)|_{x_i} \approx  \frac{\frac{f(x_{i+1}) - f(x_i)}{x_{i+1}-x_{i}} + \frac{f(x_{i}) - f(x_{i-1})}{x_{i}-x_{i-1}}}{2}.
\end{equation}
At least when your $x$-values are equally spaced, this simplifies to:
\begin{equation}
    \diff{}{x}f(x)|_{x_i} \approx  \frac{f(x_{i+1}) - f(x_{i-1})}{x_{i+1}-x_{i-1}}
\end{equation}
which is just the slope of the line connecting the point just before each point to the point just after it.

You realize that this can provide you with a solution for the slope for all points but the first and last point, so now you are two points short for the number of derivative points compared to the number of data points.  However, you figure that you can probably get the first point's derivative by computing the slope to the next point, and the last point's derivative by computing the slope to the previous point, and that will give you a value for the derivative at each data $x$-value.

Finally, since the second derivative is just the derivative of the derivative, you can perform this same procedure over and over to achieve whatever order you need. Of course, while you know you could code this, you also think that taking numerical derivatives is common enough that there must be a library for it.  A quick Google search for `numerical derivatives Python' shows that Numpy---that old familiar standby---has functions for this.

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

Your task is fairly simple when looking at it broadly.  Starting with your peak finding code,\footnote{See Chapter \ref{13PeakFind}} you can break the task up into six major activities: reading the data from the file, splitting it into reduction and oxidation waves, finding the peaks, subtracting the capacitive charging contribution, printing the results, and plotting a figure:
\begin{figure}[H]
    \centering
    \includegraphics[]{ChapNarrative/14NumDer/Figures/FlowChartCh13code.pdf}
    \caption{A flowchart showing the big picture of what you are trying to accomplish when expanding your CV analysis code to include correcting the peak currents for capacitive charging effects The boxes with the dashed outline and grey coloring correspond to regions in Figure \ref{fig:FlowChart:NumDerRefactor}, while the box with and black coloring corresponds to the flowchart in Figure \ref{fig:NumDer:NumDerFunc}.}
    \label{fig:FlowChart:NumDerOverview}
\end{figure}
You could just add the capacitive charging part to the code that you already have, but you think that would not be a great option because you do these experiments a lot and would likely end up with many slightly different copies of the same code on your computer.  This could lead to headaches and errors down the road.

Since you plan to keep doing these experiments and analyses, you would like to break down the steps into a set of tools that not only are easy to maintain, but allow you to do the analysis interactively if needed.  You realize that your flowchart provides a clear way to break up the task into logically-related subtasks!  You could write a function that splits the data, another function that finds the peaks, and a third function that determines the capacitive charging contribution.  You already have the code that would need to be in the first two functions, they're just in a script instead of their own functions. You could re-write them to be function, which will make them more convenient to use here.  Reading a little about programming, you realize that what you want to do is called \cite{refactoring}, so you dig in to it more.

\subsection{Refactoring code to turn it into functions you can reuse}
\label{Refactoring_code_to_turn_it_into_functions_you_can_reuse}

Looking at your flowchart for your CV peak finding code, you see that it is easy to find which steps go with splitting the data and which steps go with finding the peaks:
\begin{figure}[H]
    \centering
    \includegraphics[]{ChapNarrative/14NumDer/Figures/FlowChartCh13refactor.pdf}
    \caption{Your flowchart from Chapter \ref{13PeakFind}, annotated to identify the logical blocks to split the oxidation and reduction waves (dotted box) and the peak finding (gray shaded box).}
    \label{fig:FlowChart:NumDerRefactor}
\end{figure}

% We can take as an example the FTIR smoothing we did in Chapter \ref{11SmoothCrim}.  Perhaps we have many tasks that involve using at least some of this capacity.  Rather than copying and pasting this code that we already have into other code that we're developing, thus wasting time, potentially inducing new errors, and losing any future changes we make to the original code, we can simply refactor it so that parts of it can be called from external code.  Starting from the code (header removed to save space):
% \begin{lstlisting}
% import numpy as np
% import scipy.signal as sps
% from pathlib import Path
% from codechembook.quickplots import quickScatter
% from codechembook.chemtools import chemSymbols as cs

% # Get the spectrum file names
% data_file = Path("../../Data/10Integrals/ftir_0.CSV")

% # Read the file
% wavenumber, absorption = np.genfromtxt(data_file, delimiter = ",", unpack = True)

% # Apply Savitzky-Golay smoothing
% smooth_absorption = sps.savgol_filter(absorption, 100, 3)

% # Set the baseline limits and determine the values from slope and y-intercept
% base_min = 2460
% base_max = 2580
% abs_min = smooth_absorption[np.argmin(np.abs(wavenumber - base_min))]
% abs_max = smooth_absorption[np.argmin(np.abs(wavenumber - base_max))]
% m = (abs_max - abs_min) / (base_max - base_min)
% b = abs_min - m * base_min
% baseline = m * wavenumber + b

% # Plot all of the spectra to allow the limits of integration to be determined
% fig = quickScatter(x = [wavenumber]*3, y = [absorption, smooth_absorption, baseline], 
%                    name = ["raw data", "smoothed data", "baseline"], output = None)

% # go in and modify each trace to have the appearance we want.
% fig.update_traces(selector = 0, line = dict(color = "black"))
% fig.update_traces(selector = 1, line = dict(color = "gray"))
% fig.update_traces(selector = 2, line = dict(color = "black", dash = "dash"))

% # Find the y range for just the data we want to see
% ymax = np.max(absorption[np.logical_and(wavenumber > 2400, wavenumber < 2700)])
% ymin = np.min(absorption[np.logical_and(wavenumber > 2400, wavenumber < 2700)])
% fig.update_yaxes(title = "absorption", range = [ymin - 0.1 * (ymax-ymin), 1.1 * ymax]) # update y axis
% fig.update_xaxes(title = f"wavenumber /{cs('cm-1')}", range = [base_min, base_max]) # update x axis

% fig.update_layout(template = "simple_white")
% fig.show("png")
% fig.write_image(format = "png", file = "SmoothIR.png")
% \end{lstlisting}
% we can see that the code is really divided into three logical blocks:
% \begin{enumerate}
%     \item We read the file and immediately smooth it.
%     \item We determine the baseline based on two points that are provided for the start and the end of the peak.
%     \item We plot the data, the smoothed data, and the baseline.
% \end{enumerate}
% This can thus be divided into three functions:
% \begin{enumerate}
%     \item A function that reads and smooths the data.  It can take as an argument only the file name (assuming we want to apply the same smoothing each time) and then return the arrays for \lstinline{wavenumber}, \lstinline{absorption}, and \lstinline{smooth_absorption} back to the main code.
%     \item A function that computes the baseline.  It can take as a arguments \lstinline{wavenumber}, \lstinline{smooth_absorption}, \lstinline{base_min}, and \lstinline{base_max}, and return \lstinline{baseline}.
%     \item A function that plots the data.  It can take as arguments     \lstinline{wavenumber}, \lstinline{absorption}, \lstinline{smooth_absorption}, the range that we want to plot on the x axis, and the file name.
% \end{enumerate}
% We can quite easily modify this code to accommodate this structure.

The process for refactoring this code will be to identify the logically-related code segments, move them to a new function that needs to be defined, and determine the parameters that need to be passed and the variables that need to be returned.  Looking at the code in Chapter \ref{13PeakFind}, you see that lines 22--36 deal with finding the turnaround point of the potential and creating the oxidation and reduction waves, so you create a new function called \lstinline{SplitWaves()} and move those lines into its body.  This code only needs the potential and current arrays to function, so you include them in the function definition as \cite{positional_argument}s.  It produces four new arrays, the current and potential for reduction and oxidation, so you include those in a \lstinline{return} statement at the end.  Finally, you write a helpful \cite{docstring}.  Now you have a function that does the work, and you can replace all of the code that you just moved with a simple function \cite{call} to it in the main code!  Following the same logic, you create a function called \lstinline{AnalyzeCV()} to do the peak finding, calculate $E_{1/2}$, and determine the indices of the peaks.  Replacing this code with a call to the new function gets your code working again, but now in a refactored form.  You can import these functions at the command line, if you want to work with data interactively, and save yourself a ton of time.  You can also import the functions into a new code that performs the capacitive charging correction.  The rest of the code, which handles reading files and output of results, you leave in a block at the end of the script.  This block is denoted by an \lstinline{if __name__ == '__main__':} statement that makes sure it only runs if you run the script directly instead of importing the functions.\footnote{See Section \ref{Importing_any_Python_code_from_anywhere_on_your_computer}.}  The code you come up with is:

\lstinputlisting[language=python]{../../Distribution/Code by Chapter/Ch13 - Computing Numerical Derivatives/CVPeakFind_refactored.py}

% You can see here that the actual execution of the code comes down to just three function calls: \lstinline{SmoothSpectrum()}, \lstinline{GetBaseline()}, and \lstinline{PlotSmoothed()}.  The rest is handled in the functions, which are only slightly modified from the ones in your original script.  If you ran this refactored code directly, it would run no differently than the original code.  However, if you imported it, say because you wanted to use the \lstinline{GetBaseline()} function in a different code, then none of the rest of the code would execute.

% \begin{lstlisting}
% import numpy as np
% import scipy.signal as sps
% from pathlib import Path
% from codechembook.quickplots import quickScatter
% from codechembook.chemtools import chemSymbols as cs

% def SmoothSpectrum(filename):
%     """<helpful docstring>"""
%     # Get the spectrum file names
%     data_file = Path(filename)

%     # Read the file
%     wavenumber, absorption = np.genfromtxt(data_file, delimiter = ",", unpack = True)

%     # Apply Savitzky-Golay smoothing
%     smooth_absorption = sps.savgol_filter(absorption, 100, 3)
%     return wavenumber, absorption, smooth absorption

% def GetBaseline(wavenumber, smooth_absorption, base_min, base_max):
%     """<helpful docstring>"""
%     # Set the baseline limits and determine the values from slope and y-intercept
%     base_min = 2460
%     base_max = 2580
%     abs_min = smooth_absorption[np.argmin(np.abs(wavenumber - base_min))]
%     abs_max = smooth_absorption[np.argmin(np.abs(wavenumber - base_max))]
%     m = (abs_max - abs_min) / (base_max - base_min)
%     b = abs_min - m * base_min
%     baseline = m * wavenumber + b
%     return baseline

% def PlotSmoothed(wavenumber, absorption, smooth_absorption, xrange_min, xrange_max, filename):
%     "<helpful docstring>"
%     # Plot all of the spectra to allow the limits of integration to be determined
%     fig = quickScatter(x = [wavenumber]*3, y = [absorption, smooth_absorption, baseline], 
%                        name = ["raw data", "smoothed data", "baseline"], output = None)

%     # go in and modify each trace to have the appearance we want.
%     fig.update_traces(selector = 0, line = dict(color = "black"))
%     fig.update_traces(selector = 1, line = dict(color = "gray"))
%     fig.update_traces(selector = 2, line = dict(color = "black", dash = "dash"))

%     # Find the y range for just the data we want to see
%     ymax = np.max(absorption[np.logical_and(wavenumber > xrange_min, wavenumber < xrange_max)])
%     ymin = np.min(absorption[np.logical_and(wavenumber > xrange_min, wavenumber < xrange_max)])
%     fig.update_yaxes(title = "absorption", range = [ymin - 0.1 * (ymax-ymin), 1.1 * ymax]) # update y axis
%     fig.update_xaxes(title = f"wavenumber /{cs('cm-1')}", range = [base_min, base_max]) # update x axis

%     fig.update_layout(template = "simple_white")
%     fig.show("png")
%     fig.write_image(format = "png", file = filename)

% if __name__ == "__main__":
%     wave, abs, sm_abs = SmoothSpectrum("../../Data/10Integrals/ftir_0.CSV")
%     basel = GetBaseline(wave, sm_abs)
%     PlotSmoothed(wave, abs, sm_abs, 2400, 2700, "SmoothIR.png")
% \end{lstlisting}

\subsection{Importing objects from Python files in arbitrary locations using \texttt{codechembook.quickTools.importFromPy()}}
\label{Importing_objects_from_Python_files_in_arbitrary_locations}

It is often the case that you want to use functions or variables defined in Python scripts located in some folder of your file system other than the current working directory (CWD).  We discussed this already in Section \ref{Importing_any_Python_code_from_anywhere_on_your_computer}, but briefly, our options are typically:
\begin{enumerate}
    \item To copy the file to the CWD. This is not ideal since you now have multiple copies to keep track of, making it more likely that you have multiple different versions of the code to deal with.
    \item To add all the locations to the PythonPATH.  The code to do this is clunky to use in every script that needs this capability.
\end{enumerate}
The \lstinline{codechembook.quickTools} module features a function to simplify this process,\footnote{Technically, what this function is doing is adding the requested object's variable names to the `namespace' of the Python kernel. This allows you to use the variable names, just as you normally would. Exploring this more is well beyond the scope of this book, but we mention it here, so that you have terms to search, if you want to learn more. } \lstinline{importFromPy()}.  This function has a required parameter to specify the name of the file that has the code we want to import, and then a series of arbitrary positional parameters that are the names of the objects (functions, variables) in that code that you want to import.  When you run the function, the specified Python code will run, and the specified objects will then be available to the rest of the code that calls \lstinline{importFromPy()}.


\subsection{Computing derivatives using \texttt{numpy.gradient()}}
\label{Computing_derivatives_using_numpy.gradient}

Turning your attention back to derivatives, you realize that computing these slopes is very simple, particularly if your points are evenly spaced.  You could easily code this up yourself, but why do that when you could use \lstinline{numpy.gradient()}.  Reading the documentation, you see that this \cite{function} performs derivatives in multiple dimensions (including one) and has some useful features built in.\footnote{For more information on numerical derivatives, 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}).}  You find that this function is very simple---it has one required \cite{argument}, the $y$-values, though the second \cite{positional_argument} is the $x$-values which will generally be useful if your points are not spaced by the default of 1 or are unequally spaced.  The second argument could also be a single value per dimension representing the spacing between points, but only if that spacing is constant in each dimension.  It also has one \cite{keyword_argument}, \lstinline{edge_order}, that specifies how the first and last points are handled.  The default, 1, treats it as described above, while 2 represents a higher order approximation that is often more accurate.  You decide to do some quick tests by creating some data representing a function whose derivative is known analytically, $\cos(x)$ whose derivative is $-\sin(x)$.  You can then test the approach to see how it works.  You whip up a quick script, which produces Figure \ref{fig:NumDer:CosDerivativeExample}:
\lstinputlisting[language=Python]{../../Distribution/Code by Chapter/Ch13 - Computing Numerical Derivatives/CosDerivativeExample.py}
You note that the derivatives are fairly accurate and that, in this case, the 2nd order edge calculation works much better.  As would be expected, increasing the number of points that you use to sample substantially decreases the error.  

\begin{figure}[H]
    \centering
    \includegraphics[width=3in]{ChapNarrative/14NumDer/Figures/CosDerivativeExample.png}
    \caption{The results of a test of the accuracy of numerical derivatives, comparing the known analytical solution to numerical solutions using first and second order edge corrections.}
    \label{fig:NumDer:CosDerivativeExample}
\end{figure}

With this understanding, you are now ready to write your code.  

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

You start with the refactored peak finding code you wrote earlier as a template.  You can write another function that separately calculates the currents so that your code is equally modular for current as it is for potentials.  Your plan of attack is:
\begin{enumerate}
    \item Calculate the first and second derivatives (easy now).
    \item Find the part of each wave where the second derivative is nearly zero.
    \item Fit the wave where the second derivative is nearly zero to a linear model.
    \item Use the linear fit parameters to calculate the capacitive current at the peak potentials.
    \item Subtract those currents from the peak currents.
\end{enumerate}
Part 2 is really the only part that you don't know how to do.  Naively, you would think that you can just find the parts of the second derivative that are zero, but even clean CV data has some noise that will make it so the second derivative is usually not \textit{exactly} zero.  Your first idea is to look for where the second derivative is less than some threshold value.  It would be good to actually visualize what the second derivative is, so you quickly whip up a code to plot it and get the following:
\begin{figure}[H]
    \centering
    \includegraphics[width=4in]{Figures/SecondDerivativeCV.png}
    \caption{Plots of the CV, the first derivative, and the second derivative, as well as the peak locations.}
    \label{fig:NumDer:SecondDerivativeCV}
\end{figure}
You see the problem: the second derivative is also nearly zero by accident in many different parts of the wave because the noise is amplified by taking the derivative.  So, if you just look for where the second derivative is near zero, you will get points all over the place and trying to fit them will give you wrong slopes.  So, you need to narrow down the points somehow.  One thing that you know is that, for the reduction wave, the linear part should be somewhere between the beginning of the scan and the peak, and for the oxidation wave, it should be somewhere between the turnaround point and the peak.  This is helpful, but still doesn't solve the problem.  One thing you could do is find the largest contiguous range of points that are all nearly zero.  This should filter out all of the points where random noise made the value accidentally zero.  You could also smooth the data, or smooth the derivatives, to try to reduce this noise some.  All of this would require a lot of user intervention, though, for different files with different noise levels.  You would like something that is relatively immune to noise.

\begin{mdframed}
    \textbf{Derivatives from Local Regression Smoothing}
    
    \noindent For particularly noisy data, it can be useful to obtain derivatives directly from the smoothing process rather than from the raw or smoothed data.  Local regression smoothing techniques use analytical functions that have analytical formulae for their derivatives, and \lstinline{scipy.signal.savgol_filter()} can return derivatives \textit{via} the \lstinline{deriv} keyword.
\end{mdframed}

After thinking about it for a while, you come to a realization.  You can account for the noise automatically by using the standard deviation.  When the second derivative is really approximately zero, essentially all of the variance from zero will be noise.  When the second derivative is changing, then you will get an additional contribution to the standard deviation from the fact that the underlying signal is also changing.  This will necessarily increase the standard deviation.  So, you could do the following:
\begin{enumerate}
    \item Approximate the standard deviation at each point by computing the standard deviation in a moving window of points.  The standard deviation of each point will then be the standard deviation in the range of points around it, say $\pm$10 points.
    \item Find the minimum standard deviation of each wave.
    \item Start at the beginning of each wave, and find the first point where the standard deviation drops below some multiple of the minimum standard deviation.
    \item Keep going until you find the last contiguous point at which the standard deviation remains below that multiple of the minimum standard deviation.
    \item This range will be the range that you fit.
\end{enumerate}

You modify your code to also plot the standard deviation:
\begin{figure}[H]
    \centering
    \includegraphics[width=3in]{Figures/stdCV.png}
    \caption{A plot of the standard deviations of the reduction and oxidation waves with respect to potential.}
    \label{fig:NumDer:stdCV}
\end{figure}
You can now clearly see the parts of the curves that you want to fit!  


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

\begin{figure}[H]
    \centering
    \includegraphics[]{ChapNarrative/14NumDer/Figures/FlowChartCh13func.pdf}
    \caption{Flowchart for code that will determine the capacitive charging component of the current at any potential point.}
    \label{fig:NumDer:NumDerFunc}
\end{figure}


Finally, you put it all together and get the code working on a single file:
\lstinputlisting[language=python]{../../Distribution/Code by Chapter/Ch13 - Computing Numerical Derivatives/CVPeakFindandCurrent.py}

You can run this code the same way as the previous version,\footnote{See Chapter \ref{13PeakFind}.} and it outputs the following plot and text:
\begin{figure}[H]
    \centering
    \includegraphics[width=3in]{ChapNarrative/14NumDer/Figures/CV.png}

    \caption{Results of applying the above code. )top) The data visualization that is produced and (bottom) the text output in the \cite{console}.}
    \label{fig:NumDer:AnalyzedCV}
\end{figure}
\begin{multiout}
Epa = -1.236 V, ipa -2.371e-04 A
Epc = -0.885 V, ipc 2.081e-04 A
delta E = 0.351 V, E1/2 = -1.061 V
\end{multiout}
The currents are close enough to call this reversible.  You'll have to test this code a bit more on some additional data files to make sure that the process is robust, but you are even more sure than before that you have saved yourself a ton of time!

\subsection{Line-by-line explanation}

\begin{footnotesize}
\begin{enumerate}
    \item [8--13:] We will need \lstinline{numpy} and \lstinline{make_subplots()} as usual.  We also need \lstinline{LinearModel} from \lstinline{lmfit} for the fitting and extrapolation.  We'll use \lstinline{quickOpenFilename()}, \lstinline{quickPopupMessage()}, and \lstinline{importFromPy()} from the \lstinline{codechembook.quickTools} module.  Finally, we'll import the capital delta symbol from the \lstinline{codechembook.symbols.greek} module and, in addition, 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 [14:] We use \lstinline{importFromPy} to get the \lstinline{SplitWaves()} and \lstinline{AnalyzeCV} functions from the \cite{refactor}ed version of the code in Chapter \ref{13PeakFind}.
    \item [17--61:] Compute the capacitive charging component of the current.  Since we need to do this for the reduction and oxidation wave, we just create a function that we can call twice rather than reproduce the same code again.
    \begin{enumerate}
        \item [35\&36:] Take the first and second derivatives using \lstinline{numpy.gradient} with default parameters.
        \item [40\&41:] Find the standard deviation of the second derivative at each relevant point between the start of the wave and the peak.  We will ignore anything after the peak because it definitely should not be fit to determine the charging.  At each point, we define a window of points before and after it to use to compute the standard deviation.
        \item [44:] Now we want to get the indices where the second derivative is nearly zero.  We will compute this by using list comprehensions to find the points at which the standard deviation of the second derivative is less than twice its minimum value.
        \item [49:] We want to find the first contiguous range over which the standard deviation is nearly zero.  There are many ways to do this, but a simple way here is to compare the array of indices where the standard deviation is nearly zero to an array of indices that starts at the same first index and just increases by one.  At some data point, the standard deviation will be more than twice the minimum and that index will not be included in the first array.  Subtracting the second array from the first one thus gives you the difference in indices.  This should be zero for the contiguous part and non-zero at the first index where the standard deviation is more than twice the minimum.  This line will return an array that holds these differences.
        \item [53:] Next, we need to know the indices that correspond to the contiguous range where the standard deviation of the second derivative is nearly zero.  We do that here using a list comprehension.  We zip up the array of valid indices and the array that tells us whether the range is contiguous and return a new list that holds only the indices for the contiguous range.
        \item [56\&57:] Now we can simply fit the contiguous range we identified.
        \item [60:] Use the best fit slope and intercept to compute an array holding the capacitive current at each value of the potential.
        \item [61:] Return the capacitive component of the current to the main code.
    \end{enumerate}
    \item [64--74:] This is all the same as the refactored version of the code in Chapter \ref{13PeakFind}, which itself is based on the code discussed in that chapter, so see the discussion of that code for the details.
    \item [87\&88:] Using our function, compute the capacitive currents for the reduction and oxidation waves.
    \item [91--93:] Print the potentials and currents that we found.
    \item [96--125:] Make a plot that shows the peaks found and the capacitive extrapolations so that we can do a sanity check.
\end{enumerate}
\end{footnotesize}



\hrule
\input{Exercises/Chapter 13 Numerical Derivatives}
\hrule

\end{document}
