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

\chapter{Integrating Data Numerically}
\label{10NumInt}

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

\input{prereqs/10NumInt_prereqs}

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

The more code you write, the more you start to see the range of problems you can now solve and the benefits of automating data analysis workflows over manually working up large amounts of data in vendor-specific software.  On your way into work one day, you overhear a conversation between two people in a lab down the hall from you:

``...I really need to measure the rate of this click reaction, but there's no clean peak for the reactant or product and $^{15}$N-azides are too expensive for me to try NMR.  I could track the disappearance of the azide by FTIR, but the software doesn't do kinetics.''

You pipe up, ``You measure the kinetics pretty easily by writing your own script to analyze the FTIR data.  It is actually easier than you'd think, and probably faster than doing it by hand in some software.''

You all get some coffee and talk through the problem.  Tracking the azide concentration is straightforward by FTIR, because its spectral signature is well-removed from the complicated fingerprint region.  The peak shape is complex, so fitting it to extract the intensity is probably out of the question, but you can definitely integrate the band shape it after subtracting out the background.  So, all you need is a series of spectra at different time points, and you can produce a plot of the concentration of azide versus time!

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

Thinking back to your calculus class, you recall that integration was introduced to you with the ``rectangle'' rule, where you draw a series of rectangles under the curve of the function and add up the areas of the rectangles to approximate the area under the curve.  You spend a little time reading up on integration in Python and find that several \cite{method}s for integration are provided in the SciPy \cite{library}: the trapezoid rule, Simpson's rule, Gaussian quadrature...  Which one should you use?  

\begin{figure}[H]
    \centering
    \includegraphics[width=4in]{Figures/IntegrationExamples.png}
    \caption{The different approaches to numerical integration that are commonly used.}
\label{fig:NumInt:Approaches}
\end{figure}

You come across Gaussian quadrature and read that can be exactly correct for integrals of polynomials, but you are integrating experimental data with noise.  It seems that it is not as frequently used for such data because noise can cause its results to be inaccurate.  The rectangle rule, trapezoid rule, and Simpson's rule all basically come down to what you assume is happening to your $y$-values between the data points.  For the rectangle rule, all the $y$-values are the same between any two data points.  For the trapezoid, rule, you assume you can connect each $x$,$y$-pair with a straight line.  For Simpson's rule, things initially seem much more complicated, but you realize that you are basically connecting every three $x$,$y$ pairs with a parabola.\footnote{You can derive all three rules by integrating a Taylor series, a general way to estimate the value of a mathematical function at an unknown point based on knowing the value of the function and its derivatives at a known point!  Many numerical methods are applications of the Taylor series up to first or second order, with higher orders being generally more accurate but slower.}

So, the accuracy of your result (which can rarely be perfect when working with numerical data like this, compared to integrating analytical functions) will depend on how your data changes between the point and the method that you use, and often, there is not one method that is always better.  Luckily, it is easy to try different methods to see what is most accurate, in which situations.


% SOLUTION %

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

\subsection{Integrating data using \texttt{scipy.integrate}}
\label{Integrating_data_using_scipy.integrate}

You already saw in the \lstinline{scipy.integrate} module a function for integration of differential equations (Chapter \ref{07SimKin}), but looking at the module documentation you noticed that you can also integrate analytical functions and numerical data with it as well.  \lstinline{scipy.integrate} contains built-in functions to integrate numerical data using the trapezoid rule or Simpson's rule.  Apparently, the rectangle rule is too simple to bother with.  So, how do you actually do it?\footnote{For more information on integration methods, check out books like Hildebrand, F. B. \textit{Introduction to Numerical Analysis}; Dover Press, 1987, or Eric Sullivan's \textit{Numerical Methods: An Inquiry-Based Approach With Python} (\href{https://numericalmethodssullivan.github.io/}{https://numericalmethodssullivan.github.io}).}

Trapezoid rule integration is done by \lstinline{scipy.integrate.trapezoid()}.  This function has only one required \cite{argument}, \lstinline{y}, which is the $y$-values of the data to be integrated.  If you provide only this argument, the \cite{function} assumes that your data points are spaced by 1 on the $x$-axis.  Usually, this is not the case, so you will be supplying the $x$-values using the \lstinline{x} keyword.  If your points are equally spaced on the $x$-axis and you would prefer it, you can also provide the spacing as a single number \lstinline{dx} instead.  If you have multidimensional data and only want to integrate one dimension, then you can specify that with the \lstinline{axis} \cite{keyword}.  The function will return the integral(s) as a single value, or \lstinline{ndarray} if the data is multidimensional.

Simpson's rule integration is handled very similarly by \lstinline{scipy.integrate.simpson()}.  In fact, you can switch from trapezoid to Simpson's just by replacing \lstinline{trapezoid} with \lstinline{simpson}.  You decide to do a quick test, integrating $y = x$ from 0 to 1.  You know the answer should be 0.5, so you try:
\begin{pyin}
import numpy as np
\end{pyin}
\begin{pyin}
import scipy.integrate as spi
\end{pyin}
\begin{pyin}
x = np.linspace(0, 1, 10)
\end{pyin}
\begin{pyin}
spi.trapezoid(x, x = x)
\end{pyin}
\begin{pyout}
0.5
\end{pyout}
\begin{pyin}
spi.simpson(x, x = x)
\end{pyin}
\begin{pyout}
0.49999999999999994
\end{pyout}
So far so good!

The rectangle rule is not implemented in \lstinline{scipy.integrate}.  However, it can be performed in a single line.  For the forward rectangle rule, where the height of the rectangle is set by the $y$-value of the left side of the rectangle, the integral is simply \lstinline{np.sum(y[:-1] * (x[1:] - x[:-1]))}
Here, if we have $N$ data points, then we have $N-1$ rectangles.  We can compute the widths of each rectangle simply by subtracting all of the right hand side $x$-values (\lstinline{x[1:]}, skipping the first point) from all of the left hand side $x$-values (\lstinline{x[:-1]}, skipping the last point) to produce an $N-1$ array that holds the difference between each point.  We can then produce an $N-1$-length slice of \lstinline{y} by omitting the last data point, multiply that slice by the widths, and sum each product to get the integral.  If we instead wanted to use the right hand side $y$-value for the height, then we replace \lstinline{y[:-1]} with \lstinline{y[1:]}.
\begin{pyin}
np.sum(x[:-1] * (x[1:] - x[:-1]))
\end{pyin}

\subsection{Generating test data to evaluate the accuracy of numerical analysis}
\label{Generating_test_data_to_evaluate_the_accuracy_of_numerical_analysis}

Even though your azide peaks are not well reproduced by Gaussian lineshapes, the Gaussian is sort of similar and so may still serve as a good test case to evaluate the different approaches to numerical integration. If your Gaussian function is 
\begin{equation}
f\left(x\right) = A\exp\left(-\left(x-x_0\right)/2\sigma^2\right),
\end{equation}
\noindent where $A$ is the amplitude, $x_0$ is the center, and $\sigma$ is the standard deviation, then the area can be calculated simply as:

\begin{equation}
area = A\sigma\sqrt{2\pi}.
\end{equation}

So, you again have an analytical solution to compare to. You just need some text data. You whip up a quick test script to make this data and confirm that the integration approaches will work:

\lstinputlisting[language=Python, lastline = 49]{../../Distribution/Code by Chapter/Ch09 - Integrating Data Numerically/GaussianIntegral.py}

\subsection{Line-by-Line Explanation}

\begin{footnotesize}
\textbf{Lines}: \textbf{Explanation}
\begin{itemize}
    \item [6--7:] We need to import \lstinline{numpy} for arrays and the \lstinline{scipy.integrate} module for integration.
    \item [10--12:] Here we define variables holding the parameters for a Gaussian curve.
    \item [15--17:] We need to compute the $x$,$y$-data points for the Gaussian curve and its area using the analytical solution.
    \item [20:] We compute the rectangle rule integral.  If there are $N$ points in our curve, then there are $N-1$ rectangles, since each rectangle has to span two points.  Of the two $y$-values for each rectangle, we will use the one with the lower $x$-value (left rectangle rule).  We multiply this height by the width, which we calculate for each pair of x-values in case they are not equally spaced.  We can do all of this in one line using a \cite{vectorized} approach.
    \item [21:] Now we compute the trapezoid rule integral by calling the \lstinline{scipy.integrate} function \lstinline{trapezoid()}.  The first parameter is positional and takes an array holding the $y$-values, while the $x$-values are optional and are passed to a keyword parameter \lstinline{x}.
    \item [21:] We compute the Simpson's rule integral by calling the \lstinline{scipy.integrate} function \lstinline{simpson()}.  Conveniently, it is used the same way as \lstinline{trapezoid()}.
    \item [25-29:] Print the analytically-determined area and the error of the numerical methods for comparison.
\end{itemize}
\end{footnotesize}

The rectangle rule calculation is so simple that it apparently wasn't worth including in Scipy! Running this script gives the following output:

\begin{multiout}
Area = 2.506628
Errors:     Absolute --> % Diff
Rectangle: -0.036055 --> -1.4%
Trapezoid: -0.036055 --> -1.4%
Simpson's: -0.522722 --> -20.9%
\end{multiout}

You are surprised to find that the Simpson's rule calculation was less accurate than the other two, despite being presumably a better approximation.  It seems that, in this case, the parabola approximation is wose than the linear one!  Since you have the script, you try several other combinations of Gaussian parameters and $x$-coordinates, finding similar results.  Apparently, it will be fine to use the trapezoid rule for bands like this. . .

\begin{mdframed}
    \textbf{How Many Data Points Do You Need?}
    
    \noindent In most cases, your codes will run so fast that you can create test data with a very large number of data points and not worry about it taking long to calculate.  In some cases, this is not true, and you will want to use the minimum number of points necessary to describe the phenomenon you are modeling.  There are good guidelines for choosing the number, or frequency, of data points to use that come from signal analysis.  A good place to start is the Nyquist-Shannon Sampling Theorem, which links the spacing between points to the maximum rate of change of the data.  In essence, you do not want your data points to be so far apart that there could be an important shape in the data that falls completely between two data points and is, therefore, not represented at all by the data points.
\end{mdframed}

% \subsection{Combining logical expressions with \lstinline{and} and \lstinline{or} for advanced slicing}
% \label{Sec:N:10NumInt:CombiningLogicalExpressions}

% Sometimes we want to apply multiple tests to an object.  For instance, we want to know if a point $x$ is greater than 7 and less than 10.  Mathematically, we would thus write $10 > x > 7$, but we would realize that this is really two distinct logical operations: $10 > x$ and $x > 7$.  For the expression to evaluate to true, $x$ must satisfy the first \textit{and} the second.  In Python, the same thing is accomplished using the \lstinline{and} operator.  Similarly, if we want an expression to evaluate to true if one of two conditions is met, we can use the \lstinline{or} operator.  Particularly in Python, there are also shortcuts, so that \lstinline{10 > x > 7} is a shorthand for \lstinline{10 > x and x > 7}.  Even more complicated tests can be created using more \lstinline{and} and \lstinline {or} statements, like \lstinline{10 > x and x > 7 or x == 0}.  For these statements, there is an order of operations that is followed, but it is usually better to use parentheses to make it very clear what order you intend the conditionals to be evaluated.

% So far, you have sliced arrays using a single conditional statement, but in this case, you want to take a slice that corresponds to wavenumbers less than \textcolor{red}{2200} but greater than \textcolor{red}{2000}.  A quick test shows the problem:

% \begin{pyin}
% wavenumber = np.arange(400, 4000, 1)
% \end{pyin}
% \begin{pyin}
% azide_range = wavenumber[2200 > wavenumber > 2000]
% \end{pyin}
% \begin{multiout}
% Traceback (most recent call last):

%   Cell In[73], line 1
%     peak_range = 2000 < xaxis < 2400

% ValueError: The truth value of an array with more than one element is 
% ambiguous. Use a.any() or a.all()
% \end{multiout}
% Trying 

% \begin{pyin}
% azide_range = wavenumber[2200 > wavenumber and wavenumber > 2000]
% \end{pyin} yields the same result.

% A quick search shows you the problem: you gave Python a mathematically reasonable statement, but since \texttt{wavenumber} is a numpy array, Python evaluates one conditional and returns an array of equal size but containing only \texttt{True} and \texttt{False}.  When you link two conditionals like this, Python interprets it as you looking for values that satisfy the first condition AND the second one, simultaneously, and implicitly inserts the \lstinline{and} operator for you.  Unfortunately, the \lstinline{and} operator is only designed to evaluate to a single \texttt{True} or \texttt{False} value, not an array of values as you'd like.  Basically, it is asking ``are these arrays both true?'' rather than ``are the $n$th items in each of these arrays both true?''  This obviously is a problem.

% A little more digging shows you that what you really want to do is an ``element-wise'' logical operation, performing a logical operation on two corresponding items (elements) in two different arrays to form a new array holding the result of that logical operation.  The use of the `elementwise and' operator, `\lstinline{&}', provides this functionality:

% \begin{pyin}
% azide_range = (wavenumber > 2000) & (wavenumber < 2200)
% \end{pyin}
% Basically, you use \lstinline{&} rather than the Python \lstinline{and} command and everything else works as expected.  Similarly, \lstinline{or} can be replaced by `\lstinline{|}'

% Now, you can simply use \lstinline{azide_range} to slice your x and y data arrays to get just the azide peak and pass that to the \texttt{Scipy} integrator.  

\subsection{Creating customized versions of existing functions using wrappers}
\label{Creating_customized_versions_of_existing_functions_using_wrappers}

Now that you know how to integrate, you will need to define the range over which the integral should be evaluated.  You note that the Scipy integration routines just integrate the whole array, and provide no direct means to set the limits of integration.  You figure that you can just use \cite{mask}ing to select a range of data, much like you did when fitting the UV/vis spectra.\footnote{See Chapter \ref{05CompFit}.} This seems like a pain to code every time, and you could see yourself using such a function to integrate HPLC and GC chromatograms, mass spectra, and maybe even NMR data, so you decide to make your own function that calls the Scipy integration functions that you call \lstinline{IntegrateRange}.  These kinds of functions turn out to be common, and are called \cite{wrapper}s because they wrap an existing function with some helper code to make the function more convenient to use.  You could, of course, include the extra code in your \cite{script}, but if you use the \cite{function} often enough, the wrapper will likely save you considerable time.

\lstinputlisting[language=Python, lastline=80]{../../Distribution/Code by Chapter/Ch09 - Integrating Data Numerically/integrateRange.py}

\subsection{Line-by-Line Explanation}

\begin{footnotesize}
\textbf{Lines}: \textbf{Explanation}
\begin{itemize}
    \item [1--2:] We need to import \lstinline{numpy} and \lstinline{scipy.integrate} to get \lstinline{trapezoid()} and \lstinline{simpson()}.  Since \lstinline{y} and \lstinline{x} must be Numpy arrays, we will not have to import \lstinline{numpy} separately.
    \item [4:] Define the wrapper function.  It requires the $y$-values, the $x$-values, and a list of the limits of integration.  We include a keyword parameter for a string specifying which method we want to use (default: trapezoid).
    \item [21:] We sort the list so that it doesn't matter in which order the user provides it.
    \item [24:] We determine the elements of the \lstinline{x} and \lstinline{y} arrays that are between the limits of integration using \cite{mask}ing as we did in Chapter \ref{05CompFit}.
    \item [27--35:] Depending on which method the user specifies, calculate and return the appropriate integral.
\end{itemize}
\end{footnotesize}

\subsection{Simplifying the processing of iterable objects using comprehensions}
\label{Simplifying_the_processing_of_iterable_objects_using_comprehensions}

Sometimes all you want to do with a \lstinline{for} \cite{loop} is perform the same simple operation on each element of an \cite{iterable} \cite{object}.  For instance, suppose you have a list of p$H$ values, and you want to make a list of formatted strings that express those p$H$ values.  You could use a \lstinline{for} loop: 
\begin{lstlisting}
pH = [5.5, 5.7, 5.9, 6.2, 6.5, 7.1, 8.6]
pH_text = []

for trial in pH:
	pH_text.append(f'The pH is {str(trial)}.')
\end{lstlisting}
This is clear and simple, but requires writing a fair amount of code.  You'd rather do something that only takes a single line, and some searching leads you to quickly discover that Python has a way to do that called \cite{list_comprehension}s.  The basic format of the list comprehension is [<statement producing new element> for <element> in <iterable variable>], which evaluates to a new list that can be assigned to a new variable.  For the example above, you instead write:
\begin{lstlisting}
pH = [5.5, 5.7, 5.9, 6.2, 6.5, 7.1, 8.6]
pH_text = [f'The pH is {str(trial)}.' for trial in pH]
\end{lstlisting}

You further discover that list comprehensions also allow you to construct new lists that use only certain elements from an original list using a form like [<statement producing new element> for <element> in <iterable variable> if <conditional test based on element>].  The simplest example you can think of is:
\begin{pyin}
pH_basic = [trial for trial in pH if trial > 7.0]
\end{pyin}
which yields the list
\begin{pyout}
[7.1, 8.6]
\end{pyout}
This looks to you to be the list equivalent to the Numpy concept of masking.

You realize that you can combine these concepts.  Before, if you wanted to make a new list for only the trials in the above example where the solution was basic, you would write:
\begin{lstlisting}
pH = [5.5, 5.7, 5.9, 6.2, 6.5, 7.1, 8.6]
pH_text = []

for trial in pH:
	if trial > 7.0:
		pH_text.append(f'The pH is {str(trial)}.')
\end{lstlisting}
but instead you write a very terse list comprehension:
\begin{lstlisting}
pH = [5.5, 5.7, 5.9, 6.2, 6.5, 7.1, 8.6]
ph_text = [f'The pH is {str(trial)}.' for trial in pH if trial > 7.0]
\end{lstlisting}

The advantage of list comprehensions is that they do simple, boring things without taking up too much space, making it easier to see more important parts of your code.  The disadvantage is that their \cite{syntax} is not terribly intuitive, particularly to new coders or coders not familiar with Python, making your code in some ways potentially harder to read.  However, they are considered ``\cite{Pythonic}," so you are likely to encounter them even if you do not personally use them in your codes.  You can also see that they could be a great way to work with dictionaries stored in lists.

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

Now you feel like you know have you need to solve the full problem.  You come up with a plan:

\begin{enumerate}
    \item Read the data files and the associated times that they were recorded.

    \item Background subtract each spectrum.

    \item Integrate the azide peak for each spectrum and store the result.

    \item Plot the relative concentration vs. time.

    \item Fit the data to obtain a rate.
\end{enumerate}
For each spectrum, you have several pieces of data: the wavenumber array, the absorbance, the background-subtracted absorbance, the time the spectrum was recorded, and the value of the integral.  You figure that it would be best to store all of these values as a dictionary for each spectrum to keep them organized, and then make a list of these dictionaries.  When thinking about the plot, you also think it might be nice to display three subplots: one that shows the full IR spectrum (so you can remember what it looks like), one that is a zoom-in on the azide, so you can see where the integrated values are coming from, and one that is the results of the numerical integration (so you can verify that the results you obtain meet your expectations).

\begin{figure}[H]
    \centering
    \includegraphics[]{ChapNarrative/10NumInt/Figures/FlowChartCh9.pdf}
    \caption{Flowchart for calculating the integrals of the azide band in several FTIR spectra and constructing a kinetic trace showing how their concentrations vary with time.}
    \label{fig:NumInt:ActivityDiagram}
\end{figure}

Now you are finally ready to solve the problem.  You whip up the following code:

\lstinputlisting[language=Python, lastline=80]{../../Distribution/Code by Chapter/Ch09 - Integrating Data Numerically/IntegrateAzide.py}

\noindent Your new collaborator provides you with a series of csv files named ``azide\_XXXX.csv," where XXXX is the time in seconds after the start of the reaction at which the spectrum was recorded.  You place the files in a new folder in your research directory called ``Azide" and run the code.  It asks you to specify these files, and then churns away for a second or two and produces the following output:

\begin{figure}[H]
    \centering
    \includegraphics[width=6in]{./Figures/IntegrateAzide.png}
    \caption{(left) An overview of the spectra.  (middle) The feature of interest after baseline subtraction.  (right) The kinetics of the azide formation at short times.}
    \label{fig:NumInt:IntegrateAzide}
\end{figure}
You take your computer down the hall to show your colleague, who says, ``I can't believe that this wasn't produced by some professional software---I could just directly use this as a figure in my paper!  Maybe you can show me how to do this..."  You say, ``Sure, let's meet for lunch tomorrow and get started.  Bring your laptop!"

\subsection{Line-by-line explanation}

\begin{footnotesize}
\textbf{Lines}: \textbf{Explanation}
\begin{itemize}
    \item [9--14:] Import \lstinline{numpy} and \lstinline{make_subplots}.  We will also import \lstinline{quickOpenFilenames} and \lstinline{quickPopupMessage} from the \lstinline{codechembook} library for convenient file selection and user notifications.
    \item [17\&18:] We give the user a message instructing them to pick the csv files that they want to analyze, and then open the file dialog box so that they can find and select them.  This process returns a list of \lstinline{Path} variables, which we immediately sort and store as a new variable.  There is one new import here, \lstinline{customColorList}, which will be used on line 24.
    \item [21:] Set the limits over which we want to integrate.
    \item [24:] We want to make a color series that smoothly varies from one color to another.  We will use a simple function from the \lstinline{codechembook.quickPlots} module that helps us do that.  The function \lstinlnie{sampleColorScale()} will start with a Plotly color scale that is designed to continuously vary through colors and will pick out \textit{N} evenly spaced values from the scale, where \textit{N} is the number of colors you want. This value is a required parameter, and there is one optional parameter that specifies the name of the color scale (default: `bluered').  It returns a list of rgb values.
    \item [27:] We define an empty array to hold the data we are about to read.
    \item [28:] We want to make a plot that shows the original spectra, a zoom in of the peak we care about after baseline subtraction, and the kinetic trace, so we will make a 1 $\times$ 3 subplot array.
    \item [29--39:] We loop through each file one at a time to to read and plot them:
    \begin{enumerate}
        \item [29:] Zip the color scale and the filename list together so that each file is associated with a different color.
        \item [30:] Read the wavenumber and absorption data into arrays.
        \item [31:] Get the time point for each spectrum from the file name, similar to our approach in Chapter \ref{05CompFit}.  The files are named \lstinline{azide_*.csv} where \lstinline{*} is the reaction time point at which the spectrum was recorded.  We split the stem of the filename at \lstinline{_} and take the second part, which just has the text for the time.
        \item [34:] Create a dictionary that holds the wavenumber, absorption, and time point of each spectrum so that we don't have to manage several arrays later.
        \item [37--39:] Add this spectrum to the plot of the full spectra.
    \end{enumerate}
    \item [42--61:] Loop through all of the files again, this time to perform the background subtraction and integration:
    \begin{enumerate}
        \item [44\&45:] We need to figure out the array indices closest to our limits of integration.  Subtract the limit of integration from each $x$-value and look for the one closest to zero by taking the absolute value of the difference and finding the smallest one. This is done with \lstinline{numpy.argmin()}, which returns the index of the smallest value of an array.
        \item [48\&49:] Now we need to get the slope and intercept of the background.  This can be done algebraically by solving first for the slope using $\frac{rise}{run}$, then plugging that slope into $b = y - mx$ for one of the points to find the $y$-intercept.
        \item [52--55:] Here we subtract the baseline to create a new array and enter it into the dictionary.  Next, we and add this trace to our second panel of our subplot with the corresponding color.
        \item [58:] We can integrate now.  We use the \lstinline{integrateRange()} function from this chapter to simplify the process.  We store the value of the integral as a new entry in our dictionary.
        \item [61:] We will also store the maximum $y$-value for this spectrum to help us determine the $y$-range of the subplot that shows a zoom in to the peak we are integrating.
    \end{enumerate}
    \item [64]: Now we want to find the maximum $y$-value for all the spectrum to set the $y$-range for the subplot.  Here we will use a list comprehension to build a list of the maximum $y$-values for each spectrum.  Then, we use \lstinline{max()} to get the maximum of these values.
    \item [67--70:] We need to build the plot of the integral vs. time.  Create arrays for time and integral value using list comprehensions to pull the data out of the list of data dictionaries, and plot these arrays in the third subplot.
    \item [73--80:] Format the plot areas and show the plot.
\end{itemize}
\end{footnotesize}

\hrule
\input{Exercises/Chapter 09 NumInt}
\hrule

\end{document}
