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

\chapter{Fitting Simulations to Data}
\label{06AnalFit}


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

\input{prereqs/06AnalFit_prereqs}

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

The UV-vis spectra you have taken during the titration of your compound with \ce{HCl} demonstrate how clean the titration is, but you also want to know the p$K_a$ of your compound. Since you also recorded the p$H$ values after each acid addition, you can generate the titration curve shown in Figure \ref{fig:AnalFit:TitrationCurve}.

Normally you would simply look at the titration curve, find the equivalence point, and then determine the p$K_a$ by finding the half-equivalence point, just as you had been taught in your undergraduate laboratory class. However, you have three issues with this approach.  First, it is somewhat slow... you have to manually go in and analyze the data yourself and you are starting to get a taste for having Python do things for you. Second, it always feels a bit like guessing. Looking at the data you collected, you wonder how accurately you can really find the half-equivalence point using the procedure you learned as an undergraduate. This is especially problematic since there is a bit of noise in the data, making knowing the `true' value difficult. Finally, you are getting used to knowing the uncertainties in parameters using the lmfit \cite{library}.  You think it would be great to have the uncertainties in your compound's p$K_a$ value, but there is no built-in model for titration behavior in \lstinline{lmfit}. 

\begin{figure}[h]
    \centering
    \includegraphics[width = 4in, trim={0 0 0 100},clip]{ChapNarrative/06AnalFit/Figures/ExpDataOnly.png}
    \caption{You would like to fit this titration curve, so that you can extract the p$K_a$ associated with the complex you titrated.}
    \label{fig:AnalFit:TitrationCurve}
\end{figure}


Looking at the lmfit documentation, you this is doable, as lmfit allows you to define custom models that you can fit to data.  So you simply have to code a model that will reproduce the titration response you wish to fit. 

You know that your titration curve is determined by the Henderson-Hasselbalch equation:
\vspace{-0em}
\begin{align}
    pH = pK_a + log\left(\frac{[A^-]}{[HA]}\right)
    \label{nar:analfit:eq:HH}
\end{align}
\vspace{-0em}
where $[A^-]$ and $[HA]$ are the concentration of the deprotonated and protonated forms of your compound, respectively.

While you don't know the relative concentrations of the deprotonated and protonated for of your compound at every point in the titration, you know the initial concentration of the compound and that it was in the deprotonated form. In other words, you know an initial concentration, $[A^-]_i$. You also know the amount of a strong acid (i.e., HCl) added at each point, so you also know the concentration of \ce{H+} you added. In other words, you know an initial concentration, [\ce{H+}]$_i$.  Thus, for any given point on the titration curve, you can solve for the concentration of protonated and deprotonated forms present using an ICE table and the value of the equilibrium constant ($K_a = 10^{-pK_a}$). You don't yet know the p$K_a$, but you think you can use it as a fit parameter and fit your custom model to your titration data. 

You work out the familiar ICE table in your lab notebook:

\begin{tiny}
\begin{center}
\begin{tabular}{cccccc}
            & HA   & \ce{<-->}  & $A^-$   & + & \ce{H+}  \\
initial     & 0 &            & $[A^-]_i$ &   & [\ce{H+}]$_i$ \\
change      & x      &            & -x     &   & -x \\
equilibrium & x &        & $[A^-]_i$ - x & & [\ce{H+}]$_i$ - x
    
\end{tabular}
\end{center}
\end{tiny}
Just as you did when you took general chemistry, you solve for the value of $x$:
\begin{tiny}
\begin{align}
    10^{-pK_a} &= \frac{([A^-]_i - x)([H^+]_i - x)}{x}\nonumber\\
    10^{-pK_a}\cdot x &= [A^-]_i \cdot [H^+]_i - [A^-]
_i\cdot x - [H^+]_i\cdot x + x^2\nonumber\\
0 & = x^2 - ([A^-]_i + [H^+]_i + 10^{-pK_a})\cdot x + [A^-]_i\cdot[H^+]_i\nonumber\\
%\end{align}
%\begin{align}
&\textrm{using the quadratic equation to solve for x...
}\nonumber\\
x & =\frac{-b \pm \sqrt{b^2 - 4ac}}{2a} \\
&\textrm{where:}\nonumber\\
    &\ a = 1\nonumber\\
    &\ b = -1\cdot ([A^-]_i + [H^+]_i + 10^{-pK_a})\nonumber\\
    &\ c = [A^-]_i\cdot[H^+]_i\nonumber\\
&\textrm{and then you finally find the concentrations of acid and base using:}\nonumber\\
[A^-] &= [A^-]_i - x \\
&\textrm{ and } \nonumber\\
[HA] &= x \label{nar:analfit:eq:finalX}
\end{align}
\end{tiny}
You also know that, while there will be two roots for $x$, only one is going to be physically possible. In this particular case, correct value of $x$:
\begin{enumerate}
    \item must be positive
    \item must not be greater than either the starting value of the deprotonated species ($[A^-]_i$) or the amount of acid added ($[H^+]_i$)
    \item must not be a complex number (i.e., no imaginary numbers)
\end{enumerate}
Once you determine the correct value for $x$, you can use it to solve for the concentrations of the protonated and deprotonated forms, and then use Equation \ref{nar:analfit:eq:HH} to determine the p$H$.  In this manner, you can predict the titration curve for any species with a \textit{known} p$K_a$. 

You figure that, if you can predict the value of p$H$ as a function of acid concentration for a \textit{known} p$K_a$, then you can use this same approach to find which value of p$K_a$ gives the best fit to experimental data for a compound with an \textit{unknown} p$K_a$.  To do so, you fit the above behavior to your data to extract the p$K_a$ using \lstinline{lmfit}. This might sound complicated, but you have almost all the knowledge you need to do it! 


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

\subsection{Planning}
Starting a rough outline of what you want to do, you arrive at the following:
\begin{enumerate}
    \item \textit{Write code to predict the titration behavior, using the mathematical process outlined above.} An easy way to validate this would be to simulate the titration behavior for a series of equally spaced $x$-values (i.e., acid additions). If you write a function to do this, you can even treat this it its own \cite{library} and \citett{import} this behavior into any code you like, such as your code to fit the experimental data.

    \item \textit{Import the data you will fit your predicted behavior to.} You are already fairly conversant with procedures for this.

    \item \textit{Fit the prediction of the titration behavior to your data} in order to extract an estimate of the p$K_a$, and the uncertainties in this estimate. You already know how to do this for models that are built-in to \lstinline{lmfit}, and so you simply need to learn how to substitute in models you create yourself. 

    \item \textit{Plot the results}, including the estimate of the p$K_a$ and the uncertainties in this estimate. 
\end{enumerate}
Looking at the list, you decide that you already know how to do the second and final steps. You know how to do most of step 3, though there are a few new things to introduce. So the bulk of what you need to learn concerns how to code up the expected titration behavior, validate it, and import it for fitting. They also say that learning a few new coding ideas can make this a bit easier. 



%When fitting models to data, we make initial guesses at the model's parameters, evaluate the model for these guesses, compare the results to experiment, and then change our guess so that the model more closely matches the experiment and try again.  We will do this until our agreement with experiment is `good enough.'\footnote{How good is good enough? This depends on exactly how the fitting is implemented, but will not really be addressed in this book, since the default approach is good enough for us.} The practical implication of this approach is that we will need to evaluate Equation \ref{nar:analfit:eq:HH} and find the roots of a quadratic \textit{many} times. The best way to accomplish performing the same task over and over in programming is to collect the task inside a user-defined \cite{function} and then just call the function every time you want to do the task. 


% IF, elif, else 
\subsection{Automatically running different code under different conditions using \texttt{if-else} statements}
\label{Automatically_running_different_code_under_different_conditions_using_if-then-else_statements}

Evaluating the roots of the quadratic means testing if they meet certain conditions (i.e., is the root a real number). We can do this using the \lstinline{if} \cite{control_structure}.  The structure is simple: we check if something is \lstinline{True}; if it is, then we run the code \cite{block} that immediately follows. However, if the value is \lstinline{False}, we skip the immediate code and move on to the next block.  We can illustrate this below:
\begin{lstlisting}
mylist = ["True", 1, "Blue"]
if True:
    print ("This was True")
if False:
    print("This was False")
if 1 in mylist:
    print("1 is in mylist")
if 3 in mylist:
    print("3 is in mylist")
\end{lstlisting}
Running this code will output:
\begin{multiout}
This was True
1 is in mylist
\end{multiout}
You see `This was True' because the block after an \lstinline{if} statement runs when the code to the right of \lstinline{if} evaluates to \lstinline{True}. Since you directly provided the \cite{literal} \lstinline{True}, Python runs the block immediately afterwards. However, when you supply the literal \lstinline{False}, then it does not run, and we don't see 'This was false,' printed out. You do see `1 is in mylist' because \lstinline{1 in mylist} evalulated to \lstinline{True}, but since the list does not have 3, \lstinline{3 in list} evaluates to \lstinline{False}, and you don't see the results of the following block. 

It is nice to be able to check for truth, but we can make the \lstinline{if} control structure more powerful using \lstinline{else} and \lstinline{elif} (else-if)---both of which must be used in conjunction (and following) an \lstinline{if} block. To see the utility of these two new additions, consider the following code that uses 3 \lstinline{if} blocks.
\begin{lstlisting}
a = 0
if a > 0:
    print("The variable is greater than 0.")

if a < 0:
    print("The variable is less than 0.")

if a == 0:
    print("The variable is 0.")
\end{lstlisting}
\vspace{-1em}
\begin{pyout}
'The variable is 0.'
\end{pyout}

Looking at this code, you can see that you didn't really have to perform the last test: if the first and second are not true, the third must be true. Also, if the first one was \lstinline{True}, then we didn't really need to check the other two.  This is where  \lstinline{elif} and \lstinline{else} help us out. The first stands for 'else if' and basically is only run if the tests above it were False, while the last is run if every statement before it was false. 
\begin{lstlisting}
a = 0
if a > 0:
    print("The variable is greater than 0.")
elif a < 0:
    print("The variable is less than 0.")
else:
    print("The variable is 0.")

print("Done checking values.")
\end{lstlisting}
\vspace{-1em}
\begin{multiout}
The variable is 0.
Done checking values.
\end{multiout}

In this example, the first if statement was \lstinline{False}, so the program checked the conditional for the \lstinline{elif} statement, which was \lstinline{False}, so it automatically ran what was in the \lstinline{else} block, then the final print statement. If the variable had a value greater than 0, the program would have run the code associated with the if statement, then skipped the \lstinline{elif} and \lstinline{else} blocks and to go straight to the final print statement. You can play around with other examples to get a sense of how these structures work together.
\footnote{There are other reasons to use \lstinline{elif} and \lstinline{else}, and you can consult Chapter \ref{FlowControl} to learn more.}


% % linspace
% \subsection{Generating an array of equally spaced numbers using \lstinline{numpy.linspace()}}
% \label{Sec:N:06AnalFit:linspace}

% To validate our function that predicts titration behavior, we will simulate a titration experiment.  Since a titration experiment would often use a series of equivalent acid additions, this is equivalent to a series of equally spaced independent variables (a.k.a., the $x$-values). Though we could type these out by hand, this is irritating to do, and would be burdensome, or impossible, for large numbers of independent variables.  Fortunately, Numpy includes a function for generating such points: \lstinline{linspace()}.

% The basic usage is that we provide three \ul{positional arguments}, the first is the starting value, the second is the stopping value, and the third is the number of equally spaced points desired. This will return a \textbf{Numpy array} that meets these specifications. 

% \begin{pyin}
% import numpy as np
% \end{pyin}
% \begin{pyin}
% np.linspace(3.14, 100, 11)
% \end{pyin}
% \begin{pyout}
% array([ 3.14 , 12.826, 22.512, 32.198, 41.884, 51.57, 61.256, 70.942,  80.628, 90.314, 100.])
% \end{pyout}
% Here you can see we have 11 equally spaced values, starting at the first (3.14) and ending at the last (100).  So, both the upper and lower limit are \ul{inclusive}. Finally, while these \ul{arguments} can be either integer or float \textbf{data types}, \lstinline{linespace()} always returns a \lstinline{numpy array} of floats.\footnote{This is true even if only integers are used and the requested values could be represented at integers.  If you want to see, try \lstinline{np.linspace(1, 11, 11).}}

\subsection{Setting the working directory using \texttt{os.chdir()}}
\label{Setting_the_working_directory_using_os.chdir}

Any time you run Python code, the Python interpreter program must be active in \textit{some} location in your \cite{file_system} (even if you are not reading or creating files).  In other words, there must be some folder on your computer in which the files can be directly accessed by the Python program without having to include their \cite{path}.  That folder is known as the \cite{working_directory}. You can think of this as the specific folder on your computer where Python is currently active. You can find your working directory by typing \lstinline{pwd} (\ul{p}rint \ul{w}orking \ul{d}irectory) into the console in the Spyder IDE, which will give you the path associated with the working directory. If you have imported \lstinline{Path} from the pathlib \cite{library}, then you can also use \lstinline{Path.cwd()} (\ul{c}urrent \ul{w}orking \ul{d}irectory).  The latter will return a Path \cite{object} that is the current working directory. 

To change the working directory, we need to use a different library, \lstinline{os}, which provides the \lstinline{chdir()} function.\footnote{You might find code demonstrating this capability with pathlib, but at the time that this book was written, this only worked for Macs and not Windows PCs. The approach using \lstinline{os} works on both.}  This function can accept a path object or a string object, and will set the current working directory to the directory the argument specifies, as follows:

\begin{pyin}
from pathlib import Path
\end{pyin}
\begin{pyin}
Path.cwd()
\end{pyin}
This will output your current working directory. We can then change it to a different path

\begin{pyin}
import os
\end{pyin}
\begin{pyin}
new_path = Path('/path/to/new/folder')
\end{pyin}
\begin{pyin}
os.chdir(new_path)
\end{pyin}
\begin{pyin}
Path.cwd()
\end{pyin}
As long as the \lstinline{new_path} value is different than the original working directory you will have moved to a new working directory, as seen in the second output.


\subsection{Importing any Python code from anywhere on your computer}
\label{Importing_any_Python_code_from_anywhere_on_your_computer}

It is very convenient to use \cite{function}s from \acite{libraries}, because we can use the \lstinline{from <library> import <component> as <alias>} construction.  But, you realize that, all of these libraries are just Python code!  That is, they are held in a file that ends in .py!  The \cite{Numpy} library you have been using so much is just in a file called numpy.py. This must mean that we could import our own code from its own .py file.  In the context of this chapter, it might be nice to store a function that simulates the titration response in a file called titration.py and then import it into our code and use it.  However, Python cannot just import any .py file---largely for security reasons. There are two cases in which you can import a Python file into your script. 
\begin{enumerate}
    \item The .py file is in a folder that is part of the so-called \cite{pythonpath}.  This is a list in which Python keeps track of where it can read .py files from.  Part of installing a new library using \lstinline{pip/pip3/conda install} involves insuring that library is added to the PythonPATH.  When you try to import a library, Python checks the paths in it to see if the file you are trying to import are there.  If so, you can import it.  If not, then the next condition needs to be met. 

    \item The .py file you want to import is in the current \cite{working_directory}. 
\end{enumerate}

If neither of the above is true, then you will get an error when trying to import the Python file because Python doesn't know where to find it. 

\begin{wrapfigure}{r}{.5\textwidth}
    \begin{mdframed}
        \textbf{Why not Jupyter notebooks: revisited.}
        
        \noindent In Chapter \ref{Chap:N:00bArrays} we gave a few reasons as to why we don't use Jupyter notebooks for this book.  One of the reasons given was that it is hard to import code from Jupyter notebooks into other code. The reason is that a Jupyter notebook is not a .py file, and so we cannot easily import it like we will do in multiple cases in this book.
    \end{mdframed}
\end{wrapfigure}

In summary, you have two options for importing your own Python code as a library: you can add the .py file to one of the \acite{directories} that is in the PythonPATH, or you can ensure that the .py file you want to import is in the current working directory. The former is the `normal' behavior we have been using so far.  The latter is the behavior we introduce in this chapter. 

Another important thing to note: When importing components from an external Python file, the code in the external file will only be run the very first time you run the code. After that, you have access to the components that you specify in the \lstinline{import} \cite{statement}. For instance, if the function for simulating p$H$ titration is called \lstinline{pH_response}, we can use \lstinline{from titration_sim import pH_response}.  This statement will result in running the entire titration\_sim.py file and then give you access to the pH\_response function. 

Finally, when you import from a Python file, the whole code will run.  This may not be a desirable situation, since your function could be embedded in a code that works on different data, or that takes some time to run, or that needs an input file that you don't have or don't want to give it.  To avoid this, you can use a ``main entry point,'' which is some code that defines how the script should run if it is called directly as we have normally done so far.  If the same Python script is imported, then the code in the code block denoted by the main entry point will be ignored.  It is denoted with a line with an \lstinline{if} statement like:
\begin{pyin}
if __name__ == "__main__":
\end{pyin}
with the code that you want in the code block executed if the expression evaluates to \lstinline{True}. This will be \lstinline{True} when the file is run directly, and \lstinline{False} otherwise (so, \lstinline{False} when run in an \lstinline{import} statement). In this way, you can ensure that if you run the script directly, as we have been doing so far, it will run as expected, but when you import elements from it, the code in the main entry point will not be run.



%
% try except
%
\subsection{Handling errors without crashing using \lstinline{try}-\lstinline{except}-\lstinline{finally} }
\label{Handling_errors_without_crashing_using_try-except-finally_}

There are going to be times that you want to perform an operation that you know may fail sometimes.  In real life, we do this often.  We might attempt a new chemical synthesis, knowing it could fail.  However, even though the synthesis didn't work, you still want to keep doing chemistry. In Python, there are times that a step in your \cite{script} might fail, but you don't want the failure to stop your entire script from running.  Python provides a means to do this via the \lstinline{try-except-finally} \cite{control_structure}. 

The basic idea behind this control structure is to arrange code into 2-3 \cite{block}s. The first (required) is the \lstinline{try} block. The second (required) is the \lstinline{except} block. The third (optional) is the \lstinline{finally} block. 

When encountering this control structure, Python first attempts code in the \lstinline{try} block.  If this code runs without error, then the \lstinline{except} block is skipped. However, if an error occurs in the \lstinline{try} block, then the block is excited and the code in the \lstinline{except} block is attempted. Because of this behavior, the \lstinline{try} and \lstinline{except} blocks cannot exist in isolation.  You must have both.

The \lstinline{finally} block is optional, and if it is included, then it runs after completion of either the \lstinline{try} or \lstinline{except} blocks. 

\begin{lstlisting}
     var = "hi"
     try: 
         var  = 2 + var
     except:
         print("it appears that your variable is not a number.")
     finally:
         print(f"the value of var is {var}")
\end{lstlisting}

If you change this code so that \lstinline{var} starts as a number, you will get a different output.\footnote{In this usage, the \lstinline{try, except} blocks look like another approach to error handling, but we will see that they can be used in a more flexible manner as we move through this book. You can read more details about this flow control structure in Chapter \ref{FlowControl}. For now, we will use it to handle importing from a Python file that may or may not be in our current working directory.}


\medskip
Ok! With these concepts out of the way, you can now produce a script that will fit your titration curve and give an estimate (and uncertainty) for the p$K_a$. 


\subsection{Plotting in a single line of code using \texttt{codechembook.quickPlots.quickScatter()}}
\label{Plotting_in_a_single_line_of_code_using_codechembook.quickPlots.quickScatter}

The shortest code that you can produce to plot a single set of $(x,y)$-data in \cite{Plotly} is three lines:
\begin{lstlisting}
fig = make_subplots()
fig.add_scatter(x = xdata, y = ydata)
fig.show()
\end{lstlisting}
Notably, this code produces a figure with default formatting, which is not great for presenting scientific data, so you will at least want an additional \lstinline{update_layout()} function, and if you want to include axis titles, you will also need \lstinline{update_xaxes()} and \lstinline{update_yaxes()}.  This is starting to be a lot of work just to make a plot quickly.

The \lstinline{codechembook} library contains a module, \lstinline{quickPlots}, that streamlines the creation of many plots so that the same plot can be created in a single statement.  As an example, consider the function \lstinline{quickScatter()}: it largely mirrors the usage of \lstinline{fig.add_scatter()}, but includes many keyword arguments with reasonable defaults so that
\begin{pyin}
quickScatter(x = xdata, y = ydata) 
\end{pyin}
will create and show a nice looking plot.  The plot can be customized using a number of keyword arguments:
\begin{itemize}
    \item \lstinline{xlabel}, \lstinline{ylabel}: accept strings containing the text for the axis labels.  If these keywords are not used, the variable names passed to the \lstinline{x} and \lstinline{y} keyword arguments will be used as the labels.
    \item \lstinline{name}: the name of the trace, which appears in interactive mode when the mouse hovers over a data point.  This name also will appear in the legend if it is enabled, which it is by default.
    \item \lstinline{mode}: the drawing mode, either 'lines' or 'markers' corresponding to the same keyword in \lstinline{add_scatter()}.
    \item \lstinline{template}: the template to use, where the default is the built-in `simple\_white' template that you have been using.
    \item \lstinline{output}: how to show the plot.  By default, a .png file that will show in the Spyder Plots pane, but it will take any output option that \lstinline{fig.show()} can take (e.g., `browser'), as well as \lstinline{None} if you do not want the plot to show.
\end{itemize}
Importantly, \lstinline{quickScatter()} returns a Plotly Figure object that you can store as a variable and update using the usual Plotly Figure object methods for options that \lstinline{quickScatter()} does not natively support.  Thus, you can use \lstinline{quickScatter()} as a quick starting point for making more complex single-subplot figures.  In this case, passing \lstinline{None} to the \lstinline{output} keyword is convenient, since the plot will not show before it is finished being formatted.

If you want to plot multiple traces simultaneously, you have three options:
\begin{enumerate}
    \item You can use \lstinline{quickScatter()} as described above, but with \lstinline{output = None}, and store the Plotly Figure object as a variable.  You can then use \lstinline{add_scatter()} to add traces to the plot before you show it.
    \item If all of your traces have the same $x$-data, you can pass a list of $y$-data arrays to the keyword \lstinline{y} and each item in the list will be plotted as a separate trace.  Note that all $y$-data arrays must be the same length as the $x$-data array.  Also, in this case, you can pass lists of the same length containing the modes and names of the traces to the respective keywords.
    \item If your traces have different lengths, you can pass lists of $x$-data arrays and $y$-data arrays to the corresponding keywords.  As explained for the previous option, you can also use lists for the arguments to the \lstinline{mode} and \lstinline{name} parameters.
\end{enumerate}
The \lstinline{quickScatter()} function does not support subplots.  For this, see the \lstinline{quickGrid()} function in the same module.  You will also find shortcuts for making bar charts, histograms, and other useful figures.

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

Even though you know enough now to just write a complete fitting \cite{script}, it will save a lot of time and headache to first write and test the \cite{function} that models the titration behavior.  That way, you can make sure it works correctly before attempting to fit to it.  You can do this by using it to simulate some data and check to make sure that the result is what you expect.

Thinking through the this simulation, you can see that there are 4 steps:
\begin{enumerate}
    \item Set up the parameters of the experiment: the starting concentration of the deprotonated species ($[A^-]_i$) and the concentrations of the added acid for each point ([\ce{H+}]$_i$). In titrations, the latter is often expressed in terms of equivalents (let's call these \lstinline{eqs}), but you can convert the equivalents to concentration using \lstinline{eqs}$\cdot [A^-]_i$.

    \item Solve for the change in concentration at each point.  This is done using the quadratic equation, and you will find two possible values of $x$ (call them \lstinline{x1} and \lstinline{x2})

    \item Determine which value of $x$ is correct for each point.

    \item Use the value of $x$ to find the values of pH for each point and return these values. 
\end{enumerate}

Working from the above, you sketch out the plan, as the following flowchart:
\begin{figure}[h]
    \centering
    \includegraphics[]{ChapNarrative/06AnalFit/Figures/FlowChartCh7func.pdf}
    \caption{Flowchart for the function that will simulate titration behavior.}
    \label{fig:AnalFit:ActivityDiagramForSimulation}
\end{figure}

Your code that contains a function that performs these steps and uses the code to simulate a titration experiment and plot the results is:
%code
\lstinputlisting[language=Python]{../../Distribution/Code by Chapter/Ch07 - Fitting Simulations to Data/TitrationModel.py}
Running this code yields the following Figure: \\
\begin{figure}[H]
    \centering
    \includegraphics[width = 4in]{ChapNarrative/06AnalFit/Figures/Predicted.png}\\

    \caption{Results of a simulation of titration behavior. This was done using a function you wrote to simulate this behavior. This function will serve as the model for fitting experimental titration data.}
    \label{fig:AnalFit:SimTitration}
\end{figure}
You can see you have a titration curve with 20 points, for a species with a p$K_a$ of 7, but we can explain this in more detail.

%\newpage
\subsection{Line-by-line explanation}

\begin{footnotesize}
\begin{itemize}
    \item [6 \& 7:] We import \acite{libraries} we need to use.  We are introducing the \lstinline{quickScatter()} function of the \lstinline{codechembook} library. This provides a quicker way to generate scatter plots.

    \item [9:] Define the function for simulating the titration curve. We specify the \ul{parameters} that are expected. \lstinline{eqs} is a positional parameter, while \lstinline{pKa} and \lstinline{base_i} are keyword parameters. 

    \item[10-21:] This is the \ul{docstring} for the function. It explains the purpose of the function, the parameters that are in the function, and what is returned by the function.

    \item [23:] The a-term in the quadratic formula. As long as we have a monoprotic acid, a will be equal to 1 in our quadratic equation.

    \item [24:] the b-term in the quadratic formula

    \item [25:] the c-term in the quadratic formula. Because \lstinline{eqs} is a \textbf{Numpy array}, we will get another \textbf{Numpy array} here, with a unique value for each value of \lstinline{eqs}. 

    \item [28 \& 29:] We use the quadratic formula to calculate the two roots for $x$.  One of these is the correct root for our problem. For both of these, because \lstinline{c} is a \textbf{Numpy array}, we are also going to get a \textbf{Numpy array} back.  Between these two \textbf{Numpy arrays}, we will have all the roots for all the values of \lstinline{eqs}. 

    \item [32:] We make a list that has no items. It will eventually hold the correct roots.

    \item [33:] We start a \lstinline{for} \ul{loop}, using \lstinline{zip()} to move though the equivalents (\lstinline{eqs}) and roots (\lstinline{x1} and \lstinline{x2}) of the quadratic equation. In this for loop, we will run a \lstinline{if-elif-else} series of statements to identify the correct root, using the three conditions out outlined near the start of the chapter, and add it to the list. 

    \begin{itemize}

        \item [34 \& 35:] We first check to see if the root in \lstinline{x1} satisfies our first condition: that it must be positive.  If it does not, then we know that the second root (\lstinline{e2}) is the correct root, so we can add it to our list of correct roots. 
    
        \item [36 \& 37:] If the root in \lstinline{x1} passed the first condition, then we need to check the second condition: that it is smaller than both the initial concentration of base and the added acid. The later is done by noting the number of equivalents added, and using the \lstinline{base_i} as a scaling factor to get units of concentration for the acid. Also, note that we are checking two conditions, and also using \lstinline{or} here.  Again, if the first root fails, then the second root must the the correct root and we can add it our list of correct roots. 
    
        \item [38 \& 39:] If the root in \lstinline{x1} passed the second condition, then we need to check the third condition: that it not be complex. If the first root fails, then the second root must the the correct root and we can add it to our list of correct roots. 
    
        \item [40 \& 41:] If the first root did not fail any of our criteria, then it must be the correct root, and we can add to to our list of correct roots.

    \end{itemize}

    \item [43:] We convert the list of correct roots to a \textbf{Numpy array}. 

    \item [45:] We return the \textbf{Numpy array}, converted to pH.  This is the end of the function.

    \item [48:] We want to include a test of this function, to ensure it is working.  But we also know that we are going to want to call this function from other scripts. It will likely be the case that you do not wish the test to run each time you call the function, so we include a \lstinline{if __name__ == '__main__':} block here.  This block will only run if you directly run the script---which is different from reading the script during an \cite{import} statement (which is what you would do to include this function in another script). 

    \item [50:] We need a series of independent variables to pass to our function. In the context in which we have developed this function, the independent variable is the equivalents of acid we have added.  It makes the most sense to check equivalence just past 0 (0.05) and just before 1 (0.95), so we use \lstinline{numpy} \lstinline{linspace} to generate 20 equally spaced points between these two limits. 

    \item [53 \& 54:] We wish to plot the data we created.  We use \lstinline{quickScatter()} to do so.  We supply the x and y data as \cite{keyword_argument}s, using the function we just wrote to create the y-values within the call to the \lstinline{quickScatter()} function. We also can specify the mode we want, the x-axis and y-axis title, using the \lstinline{xlabel} and \lstinline{ylabel} keyword arguments. While this function does return a Plotly Figure object, the function also (by default) shows the figure.  So as long as we just want to show the plot, we need only call the function, which we have done. 
\end{itemize}
\end{footnotesize}

\medskip
With a function that can simulate a titration curve, all that is left is to import that function into a fitting script that will will turn this function into a \texttt{lmfit} model and use it to fit your data.

\begin{figure}[h]
    \centering
    \includegraphics[]{ChapNarrative/06AnalFit/Figures/FlowChartCh7code.pdf}
    \caption{Flowchart for code that uses the function you wrote to simulate titration data as a model to fit the experimental titration data.}
    \label{fig:AnalFit:ActivityDiagramForFitting}
\end{figure}

%code
\lstinputlisting[language=Python]{../../Distribution/Code by Chapter/Ch07 - Fitting Simulations to Data/TitrationFit.py}

%Running this code results in two outputs.  The first is output to the console window in Spyder and is the contents of the fit report (generated at line 36).  The second is a 

Running this code brings up a dialog box to select the titration data to fit.  You navigate to the Catalyst/Titration/pH folder and select titration.csv.  The code outputs the following plot, where the data is the orange points, the fit is the blue line, and the black line indicates the estimated value of the p$K_a$. Below this is shown the residuals.

%\includegraphics[width = 3in]{ChapNarrative/03AnalFit/Figures/FitReport.png}
\begin{figure}[H]
    \centering
    \includegraphics[width = 4in]{ChapNarrative/06AnalFit/Figures/FitTitration.png}
    \caption{Results of using the function you wrote to simulate titration behavior as a model to fit the experimental data. Shown is (top) the fit to the data and (bottom) the residuals of the fit. }
    \label{fig:AnalFit:FitResult}
\end{figure}

%You can refer to the \lstinline{lmfit} documentation for details in the fit report, though we note here that it returns the reduced chi-square value (from which you can judge the goodness of fit) and the final parameter estimate and uncertainty (bottom of the report). 

\subsection{Line-by-line explanation}
\begin{footnotesize}
\begin{itemize}
    \item [8--14:] Import the \acite{libraries} and portions of \acite{libraries} we will need.  Use aliases, as we like. 

    \item [17:] Use the \lstinline{quickOpenFilename()} function to point to the data we want to fit, and then pass the Path object it returns to \lstinline{genfromtxt()} to extract and unpack the data into two arrays. 

    \item [20--28:] We wish to import the function we wrote above so that we can use it as a fit model.  The process is:
    \begin{itemize}
        \item [20--21:] We start a \lstinline{try} block, in which we attempt to import the pH\_response function we wrote from our TitrationModel.py. This will work as long as the file (TitrationModel.py) is in the current working directory. Otherwise, this will fail, which will push us to the \lstinline{except} block.

        \item [22--28:] This \lstinline{except} block will handle the error that occurs when the current working directory does not have (TitrationModel.py).

        \begin{itemize}
            \item [23:] We open up a popup box so that the user is aware that the import attempt failed, and that they need identify which folder contains this TitrationModel.py. 

            \item [24:] We use the \lstinline{quickOpenFilename()} function to allow the user to identify the desired folder, using the system's native file explorer. 

            \item [25:] We first create a Path object pointing to the directory that is the current working directly (\lstinline{Path('.')}), use the use the \lstinline{resolve()} method of Path objects to identify where this directory is, and then store this as a variable.

            \item [26:] We then change the directory to the location identified by the user on line 24.

            \item [27:] We then again attempt to import the pH\_response function from the TitrationModel.py, which should work now because we set the current working directory to be in the same place as this file. 

            \item [28:] We then change the working directory back to the place we originally were (which we stored on line 25).  Thus, we imported the function we wanted, and we are still in the same working directly we started in. 
        \end{itemize}
    \end{itemize}
    
    \item [31:] We make an \lstinline{lmfit} model object, specifying the function that provides the desired behavior, and also the fact that we will have two independent parameters. The function is the \lstinline{pH_response} function we imported within our \lstinline{try-except} control structure.  The independent parameter is the equivalence values.

    \item [34:] We create an \texttt{lmfit} parameter object.  Initially it is empty. 

    \item [35:] We add the fit parameter, p$K_a$, to the parameter object.  Here we give the parameter a name and an initial value to start the fitting process.  We set the initial value to just the average of the pH values as a rough guess. Depending on the model and data you use, you might need to find other ways to estimate the initial guess for your model parameters. More information on parameters, including constraints, can be found in the online documentation for lmfit.

    \item [36:] We set the value of the initial concentration, in units of molar. 

    \item [39:] Fit our model to the data, and pass the fit object to a variable. We need to specify the data we use, the parameters we will use, and any independent variables.  

    \item [42:] Using the \lstinline{fitplot()} function from the \lstinline{codechembook} library, we get a Figure object of the fit.  Using the \ul{keyword argument} \lstinline{output}, we can specify that we don't want the function to show the plot.  But we still return a Figure object, and assign this to a variable. 

    \item [45--48:] It would be nice to have a line that indicates the best-fitting p$K_a$ value, and we construct this here and add it to the Figure object.
    \begin{itemize}
        \item [42:] For the $x$-values, we provide just two: the minimum and maximum values of the experimental \lstinline{eqs}.  This ensures the line we draw spans the full range of experimental values on the $x$-axis.

        \item [43:] For the $y$-values, we also provide a pair of values.  In this case, both are the final value of our parameter. This value is an attribute of a dictionary entry, with the key equal to our parameter name, and this dictionary is, itself, an attribute of the fit result object (in this case \lstinline{pH_fit}). We need to provide this twice, because we have two $x$-values, and we need the same number of $x$ and $y$ values for the plot to work. 

        \item [44:] We want this to be a line (also, no markers), and we do not need to include it in a legend. 

        \item [45:] Finally, we want to change the color of our line to be black, and so we supply this information as a dictionary to the line keyword---this is how we can format the appearance of our line. The second right parentheses ")" is what ends the function call for \lstinline{add_scatter()}
    \end{itemize}

    \item [51--54:] It would also be convenient to label the plot with the p$K_a$ value and uncertainty that was extracted from the fit. So, we create an annotation and specify its content, format, and placement.
    \begin{itemize}
        \item [48:] we specify the $x$ position of the annotation, and place it at the maximum $x$-value. We also state that the anchor for the annotation is on the right (i.e., the $x$-position species the far right position of the annotation).

        \item [49:] We specify the $y$-position of the annotation, placing it at the best fit value for the p$K_a$. We access this value just as was done on line 43. We also specify that the anchor is the bottom, meaning that the $y$-position specifies the lowest potion of the annotation. 

        \item [53:] We provide the string we wish to write to the \lstinline{text} keyword argument. We use an f-string, adding the final fit values via the corresponding \lstinline{ModelResult} object attributes, but formatted to 2 decimal places.  We also do the same for the uncertainty. 

        \item [54:] Annotations default to showing a line (called: arrow) pointing to the point being annotated. We don't want this, so we set this value to False, and then end the \lstinline{add_annotation()} function call with the closing parentheses. 
    \end{itemize}
     
    \item [55:] We show a static version of the figure. 
    \end{itemize}
\end{footnotesize}

\hrule
\input{Exercises/Chapter 07 Fitting data when your model cannot be represented by a single equation}
\hrule

\end{document}
