text:
\clearpage
\subsection{Pascal Boolean Data} % (fold)
\label{sub:pas_boolean_data}

Pascal includes a \texttt{Boolean} type, as well as \texttt{True} and \texttt{False} values.

\begin{table}[h]
\begin{minipage}{\textwidth}
\centering
\begin{tabular}{|l|c|c|}
\hline
\multicolumn{3}{|c|}{\textbf{Boolean Type}} \\
\hline
\emph{Name} & \emph{Size} & \emph{Values} \\
\hline
\texttt{Boolean} & 1 byte/8 bits & \texttt{True} or \texttt{False} \\
\hline
\end{tabular}
\caption{Pascal Boolean Type}\label{tbl:pas-boolean}
\end{minipage}
\end{table}

\begin{table}[h]
  \centering
  \begin{tabular}{|c|c|c|}
    \hline
     & \textbf{Description} & \textbf{C} \\
    \hline
    \textbf{Equal} & Are the values the same? & \texttt{a = b} \\
    \hline
    \textbf{Not Equal} & Are the values different? & \texttt{a <> b} \\
    \hline
    \textbf{Larger Than} & Is the left value larger than the right? & \texttt{a > b}  \\
    \hline
    \textbf{Less Than} & Is the left value smaller than the right? & \texttt{a < b}  \\
    \hline
    \textbf{Larger Or Equal} & Is the left value equal or larger than the right? & \texttt{a >= b}  \\
    \hline
    \textbf{Less Or Equal} & Is the left value smaller or equal to the right? & \texttt{a <= b}  \\
    \hline
  \end{tabular}
  \caption{Pascal Comparison Operators}
  \label{tbl:pas_comparison_op}
\end{table}

\begin{table}[h]
  \centering
  \begin{tabular}{|c|c|c|}
    \hline
     & \textbf{Description} & \textbf{Pascal} \\
    \hline
    \textbf{And} & Are both values True? & \texttt{a and b} \\
    \hline
    \textbf{Or} & Is at least one value True? & \texttt{a or b} \\
    \hline
    \textbf{Xor} & Is one value True, and the other False? & \texttt{a xor b} \\
    \hline
    \textbf{Not} & Is the value False? & \texttt{not a} \\
    \hline
  \end{tabular}
  \caption{Logical Operators}
  \label{tbl:pas-logical-operators}
\end{table}

\mynote{
The Pascal logical operators have precedence over the comparison operators, this means you will need to use parenthesis to group comparisons when combining them using logical operators. For example, \texttt{(a $>$ b) or (c $<$ d)}.
}

\clearpage

\passection{\pascode{plst:bool-test}{Pascal Boolean Test Code}{code/pascal/control-flow/TestBools.pas}}

% subsection c_boolean_data (end)

ext: tex

max_stars_repo_name: thoth-tech/programming-arcana

max_issues_repo_name: thoth-tech/programming-arcana

max_forks_repo_name: macite/programming-arcana


text:
%%%%%%%%%%%%%%%%%%%% author.tex %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% sample root file for your "contribution" to a proceedings volume
%
% Use this file as a template for your own input.
%
%%%%%%%%%%%%%%%% Springer %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\documentclass{svproc}
%
% RECOMMENDED %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

% to typeset URLs, URIs, and DOIs
\usepackage{url}
\def\UrlFont{\rmfamily}

\usepackage{listings}
\usepackage{color}

\definecolor{dkgreen}{rgb}{0,0.6,0}
\definecolor{gray}{rgb}{0.5,0.5,0.5}
\definecolor{mauve}{rgb}{0.58,0,0.82}

\lstset{frame=tb,
  language=Java,
  aboveskip=3mm,
  belowskip=3mm,
  showstringspaces=false,
  columns=flexible,
  basicstyle={\small\ttfamily},
  numbers=none,
  numberstyle=\tiny\color{gray},
  keywordstyle=\color{blue},
  commentstyle=\color{dkgreen},
  stringstyle=\color{mauve},
  breaklines=true,
  breakatwhitespace=true,
  tabsize=3
}

\usepackage{graphicx}

\usepackage{booktabs}

\usepackage[table,xcdraw]{xcolor}

\begin{document}
\mainmatter              % start of a contribution
%
\title{Deep Learning with Keras}
%
\titlerunning{Deep Learning with Keras}  % abbreviated title (for running head)
%                                     also used for the TOC unless
%                                     \toctitle is used
%
\author{Mohd Zamri Murah}
%
%\authorrunning{Ivar Ekeland et al.} % abbreviated author list (for running head)
%
%%%% list of authors for the TOC (use if author list has to be modified)
%\tocauthor{Ivar Ekeland, Roger Temam, Jeffrey Dean, David Grove,
%Craig Chambers, Kim B. Bruce, and Elisa Bertino}
%
\institute{Center of Artificial Intelligence Technology\\
Fakulti Teknologi Sains Maklumat\\
Universiti Kebangsaan Malaysia\\
\email{zamri@ukm.edu.my}}

\maketitle              % typeset the title of the contribution

\begin{abstract}
This technical report is a preliminary study of deep learning using Keras. We experiments deep six learning models on MNIST datasets. From the experiments, we observe the following conclusions. First, deep learning is very flexible and we can create many type of models and layers. Second, the experiments require GPU-based machine to shorten the time for each experiments. Third, the models are domain specified, not a general purpose model.

\keywords{deep learning, Keras, Python}
\end{abstract}
%
\section{Introduction}
%

Deep learning refers to artificial neural networks that are composed of many layers.
It's a growing trend in ML due to some favorable results in applications where the target function is very complex and the datasets are large\cite{hinton2012improving}. An algorithm is deep if the input is passed through several non-linearities before being output. Most modern learning algorithms (including decision trees and SVMs and naive bayes) are "shallow". Deep learning is motivated by intuition, theoretical arguments from circuit theory, empirical results, and current knowledge of neuroscience.

Deep Learning is about learning multiple levels of representation and abstraction that help to make sense of data such as images, sound, and text.

\section{Literature Review}

In machine learning and cognitive science, artificial neural networks (ANNs) are a family of models inspired by biological neural networks (the central nervous systems of animals, in particular the brain) and are used to estimate or approximate functions that can depend on a large number of inputs and are generally unknown. Artificial neural networks are generally presented as systems of interconnected "neurons" which exchange messages between each other. The connections have numeric weights that can be tuned based on experience, making neural nets adaptive to inputs and capable of learning.

\begin{figure}
% https://en.wikibooks.org/wiki/LaTeX/Importing_Graphics
% Use the relevant command to insert your figure file.
% For example, with the graphicx package use
\centering{\fbox{\includegraphics[scale=.35]{ann.png}}}
% figure caption is below the figure
\caption{Artificial Neural Network}
\label{fig:1}       % Give a unique label
\end{figure}

Deep learning is called "deep" because of the structure of ANNs. Earlier 40 years back, neural networks were only 2 layers deep as it was not computationally feasible to build larger networks. Now it is common to have neural networks with 10+ layers and even 100+ layer ANNs are being tried upon.

\begin{figure}
% https://en.wikibooks.org/wiki/LaTeX/Importing_Graphics
% Use the relevant command to insert your figure file.
% For example, with the graphicx package use
 \centering{\fbox{\includegraphics[scale=0.6]{dl_layers.png}}}
% figure caption is below the figure
\caption{Layers of neurons}
\label{fig:1}       % Give a unique label
\end{figure}

DL is essentially stack layers of neurons on top of each other. The lowest layer takes the raw data like images, text, sound, etc. and then each neurons stores some information about the data they encounter. Each neuron in the layer sends information up to the next layers of neurons which learn a more abstract version of the data below it. So the higher you go up, the more abstract features you learn.  DL systems are ANN that are constructed with multiple layers (sometimes called Multi-level Perceptrons).

\section{Keras}

Keras\cite{chollet2015keras} is a high-level neural networks library, written in Python\cite{van2007python} and capable of running on top of either TensorFlow\cite{abadi2015tensorflow} or Theano\cite{bergstra2011theano}.

%\begin{lstlisting}
%// Hello.java
%import javax.swing.JApplet;
%import java.awt.Graphics;
%
%public class Hello extends JApplet {
%    public void paintComponent(Graphics g) {
%        g.drawString("Hello, world!", 65, 95);
%    }
%}
%\end{lstlisting}

\lstinputlisting[language=Python,caption=Python code for multi-leyer perceptron]{mnist_mlp.py}
%\begin{lstlisting}[language=Python, caption=Python example]
%
%\end{lstlisting}

A deep learning model based on multi-layer perceptron\cite{cirecsan2012deep} on the MNIST dataset is given in figure~\ref{fig:1}.
A simple convnet\cite{sermanet2012convolutional} on the MNIST dataset is given in figure~\ref{fig:2}.

A Net2Net deep learning is given in figure~\ref{fig:3} and figure~\ref{fig:4}.  Net2Net is a group of methods to transfer knowledge from a teacher neural net to a student net,so that the student net can be trained faster than from scratch. This enable 'lifelong learning system' by gradually adjusting model complexity
    to data availability,and reusing transferable knowledge.

\begin{figure}
% https://en.wikibooks.org/wiki/LaTeX/Importing_Graphics
% Use the relevant command to insert your figure file.
% For example, with the graphicx package use
 \centering{\fbox{\includegraphics[scale=0.4,angle=90]{model_mlp.png}}}
% figure caption is below the figure
\caption{Deep learning architecture for multi-layers perceptrons}
\label{fig:1}       % Give a unique label
\end{figure}

\begin{figure}
% https://en.wikibooks.org/wiki/LaTeX/Importing_Graphics
% Use the relevant command to insert your figure file.
% For example, with the graphicx package use
 \centering{\fbox{\includegraphics[scale=0.35,angle=90]{model_cnn.png}}}
% figure caption is below the figure
\caption{Deep learning architecture for convolution network}
\label{fig:2}       % Give a unique label
\end{figure}

\begin{figure}
% https://en.wikibooks.org/wiki/LaTeX/Importing_Graphics
% Use the relevant command to insert your figure file.
% For example, with the graphicx package use
 \centering{\fbox{\includegraphics[scale=0.4,angle=90]{model_net2net_teacher.png}}}
% figure caption is below the figure
\caption{Deep learning architecture for teacher in Net2Net network}
\label{fig:3}       % Give a unique label
\end{figure}

\begin{figure}
% https://en.wikibooks.org/wiki/LaTeX/Importing_Graphics
% Use the relevant command to insert your figure file.
% For example, with the graphicx package use
 \centering{\fbox{\includegraphics[scale=0.4,angle=90]{model_net2net_student.png}}}
% figure caption is below the figure
\caption{Deep learning architecture for student in Net2Net network}
\label{fig:4}       % Give a unique label
\end{figure}

\section{Results}

The experiments is done on MNIST dataset\cite{lecun-mnisthandwrittendigit-2010}. The MNIST dataset is a dataset of handwritten digits, comprising 60 000 training examples and 10 000 test examples. The dataset can be downloaded from http://yann.lecun.com/exdb/mnist/.

All experiments are done on a Linux Ubuntu 16.04 with Intel i5 and 8GB RAM. No GPU is used. These experiments indicates that to do deep learning models, we need GPU-supported machines. Without GPU, the time requires to implement the models is long.

% Please add the following required packages to your document preamble:
% \usepackage{booktabs}
% \usepackage[table,xcdraw]{xcolor}
% If you use beamer only pass "xcolor=table" option, i.e. \documentclass[xcolor=table]{beamer}
% http://www.tablesgenerator.com/
% Please add the following required packages to your document preamble:
% \usepackage{booktabs}
% \usepackage[table,xcdraw]{xcolor}
% If you use beamer only pass "xcolor=table" option, i.e. \documentclass[xcolor=table]{beamer}
\begin{table}[]
\centering
\begin{tabular}{@{}lllll@{}}
\toprule
\textbf{method} & \textbf{epoch} & \textbf{accuracy} & \textbf{params} & \textbf{time} \\ \midrule
\textit{\begin{tabular}[c]{@{}l@{}}CNN (convolution \\ neural network)\end{tabular}} &  12 & 0.99 & 600810  &  35m\\
\rowcolor[HTML]{ECF4FF}
\textit{\begin{tabular}[c]{@{}l@{}}MLP (multilayer \\ perceptron)\end{tabular}} &  20 &  0.9828 &  669706&  5m \\
\textit{Net2Net} &  3 & 0.9958 &  238986 & 132m \\
\rowcolor[HTML]{ECF4FF}
\textit{\begin{tabular}[c]{@{}l@{}}IRNN (Initialize \\ Recurrent Networks)\end{tabular}} &  200 &  0.93 &  11210 & 30h  \\
\textit{\begin{tabular}[c]{@{}l@{}}Siamese multi-layer \\ perceptron\end{tabular}} &  20 & 0.995 &  133504 &  4m \\
\rowcolor[HTML]{ECF4FF}
\textit{Hierarchical RNN} &  5 &  0.9858 &  199434 & 150m  \\ \bottomrule
\end{tabular}
\caption{Experimental results using six types of deep learning networks on MNIST dataset.}
\label{my-label}
\end{table}
%
% ---- Bibliography ----

\pagebreak

\bibliography{keras}
%\bibliographystyle{splncs_srt}
\bibliographystyle{spmpsci}

\end{document}


text:
\section{Proposal Discussion}
\label{sec:sec007}

Through iterative development with four clinicians, we consolidated the following fundamental features to cover an assertiveness-based interaction.
In our prior work (Figure~\ref{fig:fig002}), we considered a medical image as the main input and attempted to exhaust all possible diagnoses.
However, clinicians often have a specific question for a breast image.
Sometimes, with important contextual information among the several clinical co-variables.
For instance, C11 (Senior) and C15 (Junior) mentioned that ``In breast cancer, some image patterns will result in a more severe prognosis depending if the patient has [co-variables] personal and/or family history''.
Also, C11 suggested that the assistant agent could alert radiologists about this available co-variables, which will substantiate the case severity.

By taking it into account, we can affirm ({\it i.e.}, in an assertive-based strategy) those co-variables.
Therefore, the assistant will justify the diagnostic and advise the clinician to make the right decision.
However, it is important to balance the affirmation levels depending on the group of medical experience.
In fact, during our focus groups (Section~\ref{sec:sec00502}), C11 (Senior) tend to prefer a more {\bf Non-Assertive} communication.
On the contrary, less experienced clinicians, such as C5 (Intern) and C15 (Junior), preferred a more {\bf Assertive} communication.

\subsection{Initial User Interface}
\label{sec:sec00701}

Concerning the assertiveness-based interactions, we will adapt the old (Figure~\ref{fig:fig002}) version of our assistant converging to the new (Figure~\ref{fig:fig024}) one.
For the new (Figure~\ref{fig:fig024}) version, we will have two prototypes.
One prototype to simulate the {\bf Assertive} scenario and another to simulate the {\bf Non-Assertive} scenario.

To fulfill this document design goals (Section~\ref{sec:sec00303}), raised from a later study, we will need to keep the following feature characteristics.
First of all, when advising clinicians for the correct use of the various modalities ({\bf MID}), we have to provide the {\it 5.3.~Modality~Selection} feature (Figure~\ref{fig:fig002}) on a {\it Multimodality} perspective.
Second, we must provide clinicians' free choice with both {\it 6.2.~Accept} and {\it 6.3.~Reject} options, so that the clinician can control ({\bf CRD}) the final answer.
Third and final, the assistant agent must substantiate the decision with co-variable explanations ({\bf BED}).

\subsection{Proposed Assistant}
\label{sec:sec00702}

As shown in Figure~\ref{fig:fig026}, the main interface of {\it BreastScreening-AI} (Figure~\ref{fig:fig024}) will embrace two scenarios: (1) {\bf Assertive}; or (2) {\bf Non-Assertive}; and two agent behaviours: (i) {\bf Proactive}; or (ii) {\bf Reactive}.
For Figure~\ref{fig:fig026}.a, the agent presents proactively the set of available co-variables outputted by our AI models and the respective accuracy.
The communication is {\bf Assertive}, meaning that the assistant substantiate the information with affirmative sentences.
Immediately, the clinician can {\bf Accept} or {\bf Reject} the assistant suggestion.
For Figure~\ref{fig:fig026}.b, the agent first present a summary of the suggested result.
After that, for this example, if the clinician press the {\bf Accept} button, it will present ({\it e.g.}, {\bf Reactive}) a more detailed information in a less affirmative ({\it i.e.}, {\bf Non-Assertive}) way.
From there, the clinician can confirm the suggestion or go {\bf Back} to {\bf Reject} the result.
With this strategies, we will try to convince the clinician to provide the best answer.
Therefore, we expect to mitigate the medical error by reducing the rates of FPs and FNs.

\subsection{Contextualizing Impressions}
\label{sec:sec00703}

All four clinicians reacted positively to our initial idea of adapting the agent according to the respective medical experience.
More precisely, C11 (Senior) mentioned during our focus group (Section~\ref{sec:sec00502}) that, the more quick and effective the communication is, the better.
Giving us an idea that Senior clinicians typically prefer ({\bf RQ2.}) a {\bf Non-Assertive} and {\bf Proactive} scenario.
On the contrary, C15 (Junior) argued that giving the last chance ({\it e.g.}, {\bf Reactive}) when validating the assistant classification with affirmative ({\it i.e.}, {\bf Assertive}) substantiation, would be important.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure*}[htbp]
\centering
\includegraphics[width=1.00\textwidth]{fig026}
\caption{Assistant agents with several scenarios and behaviours. We have the (a) first scenario, {\it i.e.}, an {\bf Assertive} agent with a {\bf Proactive} behaviour. Next, we have the (b) second scenario, {\it i.e.}, a {\bf Non-Assertive} agent with a {\bf Reactive} behaviour. In the {\bf Non-Assertive} scenario and {\bf Reactive} behaviour, the agent first show the (b1) summarized suggestion and then shows the (b2) detailed information.}
\label{fig:fig026}
\end{figure*}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Comparison Across Clinicians}
\label{sec:sec00704}

In a future study, we will introduce the 4x4 experimental settings that will be used to test the performance ({\bf RQ1.}), perception ({\bf RQ2.}) and UX ({\bf RQ3.}) of the proposed agent.
Thereafter, each scenario ({\it i.e.}, {\bf Assertive} or {\bf Non-Assertive}) will have two ({\it i.e.}, {\bf Proactive} or {\bf Reactive}) behaviour conditions.
To answer the research questions (Section~\ref{sec:sec00305}), each category of medical experience will test four conditions.
Hence, we can understand the impact of an assertive-based interaction on this medical workflow comparing it across clinicians.

\subsection{RQ1 - Improving Medical Performance}
\label{sec:sec00705}

To understand the impact ({\bf RQ1.}) of an assertiveness-based interaction, we will measure the clinician's provided BI-RADS.
Each clinician will have the respective category of medical experience associated.
We will compare the results with a later study~\cite{https://doi.org/10.13140/rg.2.2.16566.14403/1} and understand in what manner ({\it i.e.}, {\bf H1.1.} and {\bf H1.2.}) will our proposed solution improve the rates of FPs and FNs.

\subsection{RQ2 - Assertiveness Perception}
\label{sec:sec00706}

To measure perception ({\bf RQ2.}), we will ask clinicians to fill a post-task questionnaire.
The post-task questionnaire will have questions to measure the assertiveness levels of each condition.
Mainly, the questionnaire used was composed of demographic and assertiveness assessment questions.
From that, we can measure clinician's perception ({\it i.e.}, {\bf H2.1.} and {\bf H2.2.}) depending on the professional category.

\subsection{RQ3 - Improving Diagnostic Experience}
\label{sec:sec00707}

Finally, we will measure diagnostic experience ({\it i.e}, UX) of clinicians ({\bf RQ3.}).
For that, we will apply NASA-TLX, SUS and DOTS.
The scales will be answered after each subject condition.
That said, we can understand if an assertiveness-based interaction will improve ({\it i.e.}, {\bf H3.1.} and {\bf H3.2.}) the final diagnostic experience in comparison to our early studies.

max_stars_repo_name: mida-project/ates-reports

max_issues_repo_name: mida-project/ates-reports

max_forks_repo_name: mida-project/ates-reports


text:
\section{Introduction\label{sec:intro}}
%The ability to simulate complex biological flows from first principles
%has the potential to provide  insight into complicated physiological processes.
%Simulation of blood flow, in particular, is of paramount biological and clinical importance.
%Blood vessel constriction and dilation affects blood pressure, forces between \rbcs can cause clotting, various cells migrate differently through microfluidic devices.
%
%
%However, direct simulation of blood flow is an extremely challenging task.
%Even simulating the blood flow in smaller vessels requires modeling millions of cells (one microliter of blood contains around five million \rbcs) along with a complex blood vessel.
%\rbcs are highly deformable and cannot be well-approximated by rigid particles.
%The volume fraction of cells in human blood flow reaches 45\%, which means that a very large fraction of cells are in close contact with other cells or vessel walls at any given time.
%These constraints preclude a large number of discretization points per cell and make an evolving mesh of the fluid domain impractical and costly at large scale.
%
%Simulations capable of capturing these various types of flows faithfully must be
%\begin{itemize}
%\item \emph{numerically accurate}, to solve the model equations without
%  concern for numerical error;
%\item \emph{robust}, to handle high-volume-fraction flows, close contact between cells and vessel walls, complex geometries, and long simulation times;
%\item \emph{efficient and scalable}, to support a realistic number of cells in
%  flows through complex blood vessels.
%\end{itemize}

Achieving a fast, numerically accurate, and robust blood flow simulation requires that the system meets a number of stringent requirements.
While previous work has made significant progress \cite{Malhotra2017,lu2018parallel,rahimian2010petascale}, we focus on several new infrastructure components essential for handling confined flows and arbitrarily long-time, high volume fractions \rbc flows; in particular, our work is able to realize each of these goals.

We formulate the viscous flow in blood vessels as an integro-differential equation and make use of fast scalable summation algorithms for efficient implementation, as in prior \rbc simulations \cite{Veerapaneni2011}.
This is the only approach to date that maintains high accuracy at the microscopic level while avoiding expensive discretization of fluid volume: all degrees of freedom reside on the surfaces of \rbcs and blood vessels.

%To achieve high accuracy with few degrees of freedom per cell, and a compact boundary representation, we use spherical harmonic representations for cell boundaries and high-order representations (polynomials patches) for the blood vessels. As in previous work, a semi-implicit time stepping scheme is used.
%To achieve high accuracy with minimal degrees of freedom per cell, we required a smooth yet compact boundary representation; we use spherical harmonic representations for cell boundaries and high-order polynomials for the blood vessels.
%We update \rbc positions with a semi-implicit time stepping scheme.

The most important novel aspects of our system include:
(a) handling the \rbc-blood vessel interaction with a fully parallel, high-order
boundary integral equation solver;
(b) explicit handling of collisions with a parallel constraint-based resolution and detection algorithm.
The former is essential for modeling confined flows, while the latter is essential for handling high-volume fraction flows at long time scales without excessively small time steps or fine spatial discretizations.

\para{Our contributions}
\begin{enumerate}
  \item We present a parallel platform for long-time simulations of \rbcs through complex blood vessels.
    The extension to suspensions of various particulates (fibers, rigid bodies etc.) is straightforward from the boundary integral formulation.
    Flows through several complicated geometries are demonstrated.
  \item We have parallelized a boundary solver for elliptic PDEs on smooth complex geometries in \threed.
    By leveraging the parallel fast-multipole method of \cite{malhotra2015} and the parallel forest of quadtrees of \cite{BursteddeWilcoxGhattas11}, we are able to achieve good parallel performance and load balancing.
  \item We have extended the parallel collision handling of \cite{lu2018parallel} to include rigid \threed \linebreak boundaries composed of patches.
  \item We present weak and strong scalability results of our
    simulation on the Skylake cluster and weak scaling results on the
    Knights Landing cluster on Stampede2 at the Texas Advanced
    Computing Center along with several visualizations of long-time,
    large-scale blood cell flows through vessels.
    We observe %a 15$\times$ speed-up in strong scaling at
    49\% strong scaling efficiency for a 32-fold increase of compute
    cores.
    In our largest test on 12288 cores, we simulate 1,048,576 \rbcs in
    a blood vessel composed of 2,097,152
    patches with weak scaling efficiency of 71\% compared to 192 cores
    (\cref{fig:wscale-large-grain}). In each time step, this test uses
    over three billion degrees of freedom
    and over four billion surface elements (triangles) for collision.
  %\item We achieve to simulate the realistic blood flows of high-volume fraction($47\%$)\cref{fig:high-vol-snap}.
  \item We are able to simulate realistic human blood flows with \rbc volume
    fractions over $47\%$ (\cref{fig:high-vol-snap}).

\end{enumerate}
\textbf{Limitations. }
Despite the advantages and contributions of the computational framework presented
here, our work has some limitations.
We have made several simplifications in our model for \rbcs.
We are restricted to the low Reynolds number regime, i.e., small arteries and capillaries.
We use a simplified model for \rbcs, assuming the cell membranes to be inextensible
and with no in-plane shear rigidity.
It has been shown that flows in arterioles and capillaries with
diameter of <50 $\mu$m and \rbcs with 5 $\mu$m diameter have a
Reynolds number of $<5 \times 10^{-3}$ \cite{wang2013simulation}\cite[Section 5.4]{caro2012mechanics} with roughly $2$\% error in approximating confined flows \cite{al2008motion}.
This is sufficient for our interest in the qualitative behavior of particulate flows, with the possibility of investigating rheological dynamics in larger channels.

Regarding algorithms, each \rbc is discretized with an equal number of points, despite the
varied behavior of the velocity through the vessel.
Adaptive refinement is required in order to resolve the velocity accurately.
%We only use first-order time stepping in this work, but high-order time
%stepping, like spectral deferred correction, are easily integrated as shown in \cite{lu2017}.
Finally, the blood vessel is constructed to satisfy certain geometric
constraints that allow for the solution of \cref{eq:double_layer_int_eq} via
singular integration.
This can be overcome through uniform refinement, but a parallel adaptive
algorithm is required to maintain good performance.

\textbf{Related work: blood flow. }
Large-scale simulation of \rbc flows typically fall into four
categories: (a) \textit{Immersed boundary (\ib)} and \textit{immersed
  interface methods}; (b) particle-based methods such as
%\textit{lattice Boltzmann} (\lb),
\textit{dissipative particle dynamics (\dpd)} and \textit{smoothed
  particle hydrodynamics (\sph)} (c) multiscale network-based
approaches and (d) \textit{boundary integral equation (\bie)}
approaches.
%\note[MJM]{a reviewer said LB isn't a particle-based method?}
For a comprehensive review of general blood flow simulation methods, see \cite{freund2014numerical}.
\ib methods can produce high-quality simulations of heterogeneous particulate flows in complex blood vessels \cite{balogh2017direct,balogh2017computational, xu2013large}.
These methods typically require a finite element solve for each \rbc to compute membrane tensions and use \ib to couple the stresses with the fluid.
This approach quickly becomes costly, especially for high-order elements, and although reasonably large simulation have been achieved \cite{saadat2018immersed,saadat2019simulation}, large-scale parallelization has remained a challenge.
A different approach to simulating blood flow is with multiscale reduced-order models.
By making simplifying assumptions about the fluid behavior throughout the domain
and transforming the complex fluid system into a simpler flow problem, the
macroscopic behaviors of enormous capillary systems can be characterized
\cite{peyrounette2018multiscale,perdikaris2016multiscale} and
scaled up to thousands of cores \cite{perdikaris2015effective}.
This comes at a cost of local accuracy; by simulating the flows directly, we are able to accurately resolve local \rbc dynamics that are not captured by such schemes.

Particle-based methods have had the greatest degree of success at large-scale blood flow simulations \cite{gounley2017computational,grinberg2011new,randles2015massively, rossinelli2015silico}.
These types of approaches are extremely flexible in modeling the fluid and immersed particles, but are computationally demanding and usually suffer from numerical stiffness that requires very small time steps for a given target accuracy. For a comprehensive review, see \cite{ye2016particle}.
There have also been recent advances in coupling a particle-based
\dpd-like scheme with \ib in parallel \cite{ye2017hybrid,ye2018three},
but the number of \rbcs simulated and the complexity of the boundary seems to be limited.

\bie methods have successfully realized large-scale simulations of millions of \rbcs \cite{rahimian2010petascale} in free space.
Recently, new methods for robust handling of collisions between \rbcs in high-volume fraction simulations have been introduced \cite{lu2018parallel,Malhotra2017}.
This approach is versatile and efficient due to only requiring discretization of
\rbcs and blood vessel surfaces, while achieving high-order convergence and optimal complexity implementation due to fast summation methods \cite{Veerapaneni2009b,Veerapaneni2011,rahimian2015,sorgentone2018highly,sorgentone20193d,af2016fast,Zhao2010}.
To solve elliptic partial differential equations, \bie approaches have been successful in several application domains \cite{YBZ,wala20183d,wala2018optimization,bruno2013high}.
However, to our knowledge, there has been no work combining a Stokes
boundary solver on arbitrary complex geometries in \threed with a
collision detection and resolution scheme to simulate \rbc flows at large scale.
This work aims to fill this gap, illustrating that this can  be
achieved in a scalable manner.

\textbf{Related work: collisions. }
Parallel collision detection methods are a well-studied area in computer graphics for both shared memory and GPU parallelism \cite{Liu2010,Mazhar2011,Kim2009}.
\cite{Iglberger2009,Du2017} detect collisions between rigid bodies in a distributed memory architecture via domain decomposition.
\cite{Pabst2010} constructs a spatial hash to cull collision candidates and explicitly check candidates that hash to the same value.
The parallel geometry and physics-based collision resolution scheme detailed in \cite{yan2019computing} is most similar to the scheme used in this work.
However, such discrete collision detection schemes require small time steps to guarantee detections which can become costly for high-volume fraction simulations.


text:
\chapter{Sequence motif analysis using Bio.motifs}
\label{chapter:motifs}

This chapter gives an overview of the functionality of the
\verb|Bio.motifs| package included in Biopython. It is intended
for people who are involved in the analysis of sequence motifs, so I'll
assume that you are familiar with basic notions of motif analysis. In
case something is unclear, please look at Section~\ref{sec:links}
for some relevant links.

Most of this chapter describes the new \verb|Bio.motifs| package included
in Biopython 1.61 onwards, which is replacing the older \verb|Bio.Motif| package
introduced with Biopython 1.50, which was in turn based on two older former
Biopython modules, \verb|Bio.AlignAce| and \verb|Bio.MEME|. It provides
most of their functionality with a unified motif object implementation.

Speaking of other libraries, if you are reading this you might be
interested in \href{http://fraenkel-nsf.csbi.mit.edu/TAMO/}{TAMO}, another python library
designed to deal with sequence motifs. It supports more \emph{de-novo}
motif finders, but it is not a part of Biopython and has some restrictions
on commercial use.

\section{Motif objects}
\label{sec:object}

Since we are interested in motif analysis, we need to take a look at
\verb|Motif| objects in the first place. For that we need to import
the Bio.motifs library:

%doctest ../Tests/motifs
\begin{minted}{pycon}
>>> from Bio import motifs
\end{minted}

and we can start creating our first motif objects. We can either create
a \verb+Motif+ object from a list of instances of the motif, or we can
obtain a \verb+Motif+ object by parsing a file from a motif database
or motif finding software.

\subsection{Creating a motif from instances}

Suppose we have these instances of a DNA motif:

%cont-doctest
\begin{minted}{pycon}
>>> from Bio.Seq import Seq
>>> instances = [Seq("TACAA"),
...              Seq("TACGC"),
...              Seq("TACAC"),
...              Seq("TACCC"),
...              Seq("AACCC"),
...              Seq("AATGC"),
...              Seq("AATGC"),
...             ]
\end{minted}
then we can create a Motif object as follows:

%cont-doctest
\begin{minted}{pycon}
>>> m = motifs.create(instances)
\end{minted}
The instances are saved in an attribute \verb+m.instances+, which is essentially a Python list with some added functionality, as described below.
Printing out the Motif object shows the instances from which it was constructed:

%cont-doctest
\begin{minted}{pycon}
>>> print(m)
TACAA
TACGC
TACAC
TACCC
AACCC
AATGC
AATGC
<BLANKLINE>
\end{minted}
The length of the motif is defined as the sequence length, which should be the same for all instances:

%cont-doctest
\begin{minted}{pycon}
>>> len(m)
5
\end{minted}
The Motif object has an attribute \verb+.counts+ containing the counts of each
nucleotide at each position. Printing this counts matrix shows it in an easily readable format:

%cont-doctest
\begin{minted}{pycon}
>>> print(m.counts)
        0      1      2      3      4
A:   3.00   7.00   0.00   2.00   1.00
C:   0.00   0.00   5.00   2.00   6.00
G:   0.00   0.00   0.00   3.00   0.00
T:   4.00   0.00   2.00   0.00   0.00
<BLANKLINE>
\end{minted}

You can access these counts as a dictionary:

%cont-doctest
\begin{minted}{pycon}
>>> m.counts["A"]
[3, 7, 0, 2, 1]
\end{minted}
but you can also think of it as a 2D array with the nucleotide as the first
dimension and the position as the second dimension:

%cont-doctest
\begin{minted}{pycon}
>>> m.counts["T", 0]
4
>>> m.counts["T", 2]
2
>>> m.counts["T", 3]
0
\end{minted}
You can also directly access columns of the counts matrix
%Don't doctest this as dictionary order is platform dependent:
\begin{minted}{pycon}
>>> m.counts[:, 3]
{'A': 2, 'C': 2, 'T': 0, 'G': 3}
\end{minted}
Instead of the nucleotide itself, you can also use the index of the nucleotide
in the alphabet of the motif:

%cont-doctest
\begin{minted}{pycon}
>>> m.alphabet
'ACGT'
>>> m.counts["A",:]
(3, 7, 0, 2, 1)
>>> m.counts[0,:]
(3, 7, 0, 2, 1)
\end{minted}
The motif has an associated consensus sequence, defined as the sequence of
letters along the positions of the motif for which the largest value in the
corresponding columns of the \verb+.counts+ matrix is obtained:

%cont-doctest
\begin{minted}{pycon}
>>> m.consensus
Seq('TACGC')
\end{minted}
as well as an anticonsensus sequence, corresponding to the smallest values in
the columns of the \verb+.counts+ matrix:

%cont-doctest
\begin{minted}{pycon}
>>> m.anticonsensus
Seq('CCATG')
\end{minted}
Note that there is some ambiguity in the definition of the consensus and anticonsensus sequence if in some columns multiple nucleotides have the maximum or minimum count.

You can also ask for a degenerate consensus sequence, in which ambiguous
nucleotides are used for positions where there are multiple nucleotides with
high counts:

%cont-doctest
\begin{minted}{pycon}
>>> m.degenerate_consensus
Seq('WACVC')
\end{minted}
Here, W and R follow the IUPAC nucleotide ambiguity codes: W is either A or T,
and V is A, C, or G \cite{cornish1985}. The degenerate consensus sequence is
constructed following the rules specified by Cavener \cite{cavener1987}.

We can also get the reverse complement of a motif:

%cont-doctest
\begin{minted}{pycon}
>>> r = m.reverse_complement()
>>> r.consensus
Seq('GCGTA')
>>> r.degenerate_consensus
Seq('GBGTW')
>>> print(r)
TTGTA
GCGTA
GTGTA
GGGTA
GGGTT
GCATT
GCATT
<BLANKLINE>
\end{minted}

The reverse complement and the degenerate consensus sequence are
only defined for DNA motifs.

\subsection{Creating a sequence logo}
If we have internet access, we can create a \href{https://weblogo.berkeley.edu}{weblogo}:
\begin{minted}{pycon}
>>> m.weblogo("mymotif.png")
\end{minted}
We should get our logo saved as a PNG in the specified file.

\section{Reading motifs}
\label{sec:io}

Creating motifs from instances by hand is a bit boring, so it's
useful to have some I/O functions for reading and writing
motifs. There are not any really well established standards for storing
motifs, but there are a couple of formats that are more used than
others.

\subsection{JASPAR}
One of the most popular motif databases is \href{http://jaspar.genereg.net}{JASPAR}. In addition to the motif sequence information, the JASPAR database stores a lot of meta-information for each motif. The module \verb+Bio.motifs+ contains a specialized class \verb+jaspar.Motif+ in which this meta-information is represented as attributes:
\begin{itemize}
    \item \verb+matrix_id+ - the unique JASPAR motif ID, e.g. 'MA0004.1'
    \item \verb+name+ - the name of the TF, e.g. 'Arnt'
    \item \verb+collection+ - the JASPAR collection to which the motif belongs, e.g. 'CORE'
    \item \verb+tf_class+ - the structural class of this TF, e.g. 'Zipper-Type'
    \item \verb+tf_family+ - the family to which this TF belongs, e.g. 'Helix-Loop-Helix'
    \item \verb+species+ - the species to which this TF belongs, may have multiple values, these are specified as taxonomy IDs, e.g. 10090
    \item \verb+tax_group+ - the taxonomic supergroup to which this motif belongs, e.g. 'vertebrates'
    \item \verb+acc+ - the accession number of the TF protein, e.g. 'P53762'
    \item \verb+data_type+ - the type of data used to construct this motif, e.g. 'SELEX'
    \item \verb+medline+ - the Pubmed ID of literature supporting this motif, may be multiple values, e.g. 7592839
    \item \verb+pazar_id+ - external reference to the TF in the PAZAR database, e.g. 'TF0000003'
    \item \verb+comment+ - free form text containing notes about the construction of the motif
\end{itemize}

The \verb+jaspar.Motif+ class inherits from the generic \verb+Motif+ class and therefore provides all the facilities of any of the motif formats --- reading motifs, writing motifs, scanning sequences for motif instances etc.

JASPAR stores motifs in several different ways including three different flat file formats and as an SQL database. All of these formats facilitate the construction of a counts matrix. However, the amount of meta information described above that is available varies with the format.

\subsubsection*{The JASPAR \texttt{sites} format}

The first of the three flat file formats contains a list of instances. As an example, these are the beginning and ending lines of the JASPAR \verb+Arnt.sites+ file showing known binding sites of the mouse helix-loop-helix transcription factor Arnt.
\begin{minted}{text}
>MA0004 ARNT 1
CACGTGatgtcctc
>MA0004 ARNT 2
CACGTGggaggtac
>MA0004 ARNT 3
CACGTGccgcgcgc
...
>MA0004 ARNT 18
AACGTGacagccctcc
>MA0004 ARNT 19
AACGTGcacatcgtcc
>MA0004 ARNT 20
aggaatCGCGTGc
\end{minted}
The parts of the sequence in capital letters are the motif instances that were found to align to each other.

We can create a \verb+Motif+ object from these instances as follows:

%cont-doctest
\begin{minted}{pycon}
>>> from Bio import motifs
>>> with open("Arnt.sites") as handle:
...     arnt = motifs.read(handle, "sites")
...
\end{minted}
The instances from which this motif was created is stored in the \verb+.instances+ property:

%cont-doctest
\begin{minted}{pycon}
>>> print(arnt.instances[:3])
[Seq('CACGTG'), Seq('CACGTG'), Seq('CACGTG')]
>>> for instance in arnt.instances:
...     print(instance)
...
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
CACGTG
AACGTG
AACGTG
AACGTG
AACGTG
CGCGTG
\end{minted}
The counts matrix of this motif is automatically calculated from the instances:

%cont-doctest
\begin{minted}{pycon}
>>> print(arnt.counts)
        0      1      2      3      4      5
A:   4.00  19.00   0.00   0.00   0.00   0.00
C:  16.00   0.00  20.00   0.00   0.00   0.00
G:   0.00   1.00   0.00  20.00   0.00  20.00
T:   0.00   0.00   0.00   0.00  20.00   0.00
<BLANKLINE>
\end{minted}
This format does not store any meta information.

\subsubsection*{The JASPAR \texttt{pfm} format}

JASPAR also makes motifs available directly as a count matrix,
without the instances from which it was created. This \verb+pfm+ format only
stores the counts matrix for a single motif.
For example, this is the JASPAR file \verb+SRF.pfm+ containing the counts matrix for the human SRF transcription factor:
\begin{minted}{text}
 2 9 0 1 32 3 46 1 43 15 2 2
 1 33 45 45 1 1 0 0 0 1 0 1
39 2 1 0 0 0 0 0 0 0 44 43
 4 2 0 0 13 42 0 45 3 30 0 0
\end{minted}
We can create a motif for this count matrix as follows:

%cont-doctest
\begin{minted}{pycon}
>>> with open("SRF.pfm") as handle:
...     srf = motifs.read(handle, "pfm")
...
>>> print(srf.counts)
        0      1      2      3      4      5      6      7      8      9     10     11
A:   2.00   9.00   0.00   1.00  32.00   3.00  46.00   1.00  43.00  15.00   2.00   2.00
C:   1.00  33.00  45.00  45.00   1.00   1.00   0.00   0.00   0.00   1.00   0.00   1.00
G:  39.00   2.00   1.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00  44.00  43.00
T:   4.00   2.00   0.00   0.00  13.00  42.00   0.00  45.00   3.00  30.00   0.00   0.00
<BLANKLINE>
\end{minted}
As this motif was created from the counts matrix directly, it has no instances associated with it:

%cont-doctest
\begin{minted}{pycon}
>>> print(srf.instances)
None
\end{minted}
We can now ask for the consensus sequence of these two motifs:

%cont-doctest
\begin{minted}{pycon}
>>> print(arnt.counts.consensus)
CACGTG
>>> print(srf.counts.consensus)
GCCCATATATGG
\end{minted}

As with the instances file, no meta information is stored in this format.

\subsubsection*{The JASPAR format \texttt{jaspar}}

The \verb+jaspar+ file format allows multiple motifs to be specified in a single file. In this format each of the motif records consist of a header line followed by four lines defining the counts matrix. The header line begins with a \verb+>+ character (similar to the Fasta file format) and is followed by the unique JASPAR matrix ID and the TF name. The following example shows a \verb+jaspar+ formatted file containing the three motifs Arnt, RUNX1 and MEF2A:
\begin{minted}{text}
>MA0004.1 Arnt
A  [ 4 19  0  0  0  0 ]
C  [16  0 20  0  0  0 ]
G  [ 0  1  0 20  0 20 ]
T  [ 0  0  0  0 20  0 ]
>MA0002.1 RUNX1
A  [10 12  4  1  2  2  0  0  0  8 13 ]
C  [ 2  2  7  1  0  8  0  0  1  2  2 ]
G  [ 3  1  1  0 23  0 26 26  0  0  4 ]
T  [11 11 14 24  1 16  0  0 25 16  7 ]
>MA0052.1 MEF2A
A  [ 1  0 57  2  9  6 37  2 56  6 ]
C  [50  0  1  1  0  0  0  0  0  0 ]
G  [ 0  0  0  0  0  0  0  0  2 50 ]
T  [ 7 58  0 55 49 52 21 56  0  2 ]
\end{minted}
The motifs are read as follows:
\begin{minted}{pycon}
>>> fh = open("jaspar_motifs.txt")
>>> for m in motifs.parse(fh, "jaspar"))
...     print(m)
TF name	 Arnt
Matrix ID	MA0004.1
Matrix:
        0      1      2      3      4      5
A:   4.00  19.00   0.00   0.00   0.00   0.00
C:  16.00   0.00  20.00   0.00   0.00   0.00
G:   0.00   1.00   0.00  20.00   0.00  20.00
T:   0.00   0.00   0.00   0.00  20.00   0.00

TF name	 RUNX1
Matrix ID	MA0002.1
Matrix:
        0      1      2      3      4      5      6      7      8      9     10
A:  10.00  12.00   4.00   1.00   2.00   2.00   0.00   0.00   0.00   8.00  13.00
C:   2.00   2.00   7.00   1.00   0.00   8.00   0.00   0.00   1.00   2.00   2.00
G:   3.00   1.00   1.00   0.00  23.00   0.00  26.00  26.00   0.00   0.00   4.00
T:  11.00  11.00  14.00  24.00   1.00  16.00   0.00   0.00  25.00  16.00   7.00

TF name	 MEF2A
Matrix ID	MA0052.1
Matrix:
        0      1      2      3      4      5      6      7      8      9
A:   1.00   0.00  57.00   2.00   9.00   6.00  37.00   2.00  56.00   6.00
C:  50.00   0.00   1.00   1.00   0.00   0.00   0.00   0.00   0.00   0.00
G:   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   2.00  50.00
T:   7.00  58.00   0.00  55.00  49.00  52.00  21.00  56.00   0.00   2.00
\end{minted}

Note that printing a JASPAR motif yields both the counts data and the available meta-information.

\subsubsection*{Accessing the JASPAR database}

In addition to parsing these flat file formats, we can also retrieve motifs from a JASPAR SQL database. Unlike the flat file formats, a JASPAR database allows storing of all possible meta information defined in the JASPAR \verb+Motif+ class. It is beyond the scope of this document to describe how to set up a JASPAR database (please see the main \href{http://jaspar.genereg.net}{JASPAR} website). Motifs are read from a JASPAR database using the \verb+Bio.motifs.jaspar.db+ module.  First connect to the JASPAR database using the JASPAR5 class which models the the latest JASPAR schema:
\begin{minted}{pycon}
>>> from Bio.motifs.jaspar.db import JASPAR5
>>>
>>> JASPAR_DB_HOST = <hostname>
>>> JASPAR_DB_NAME = <db_name>
>>> JASPAR_DB_USER = <user>
>>> JASPAR_DB_PASS = <passord>
>>>
>>> jdb = JASPAR5(
...     host=JASPAR_DB_HOST,
...     name=JASPAR_DB_NAME,
...     user=JASPAR_DB_USER,
...     password=JASPAR_DB_PASS
... )
\end{minted}

Now we can fetch a single motif by its unique JASPAR ID with the \verb+fetch_motif_by_id+ method. Note that a JASPAR ID conists of a base ID and a version number separated by a decimal point, e.g. 'MA0004.1'. The \verb+fetch_motif_by_id+ method allows you to use either the fully specified ID or just the base ID. If only the base ID is provided, the latest version of the motif is returned.
\begin{minted}{pycon}
>>> arnt = jdb.fetch_motif_by_id("MA0004")
\end{minted}
Printing the motif reveals that the JASPAR SQL database stores much more meta-information than the flat files:
\begin{minted}{pycon}
>>> print(arnt)
TF name	Arnt
Matrix ID	MA0004.1
Collection	CORE
TF class	Zipper-Type
TF family	Helix-Loop-Helix
Species	10090
Taxonomic group	vertebrates
Accession	['P53762']
Data type used	SELEX
Medline	7592839
PAZAR ID	TF0000003
Comments	-
Matrix:
	0      1      2      3      4      5
A:   4.00  19.00   0.00   0.00   0.00   0.00
C:  16.00   0.00  20.00   0.00   0.00   0.00
G:   0.00   1.00   0.00  20.00   0.00  20.00
T:   0.00   0.00   0.00   0.00  20.00   0.00

\end{minted}

We can also fetch motifs by name. The name must be an exact match (partial matches or database wildcards are not currently supported). Note that as the name is not guaranteed to be unique, the \verb+fetch_motifs_by_name+ method actually returns a list.
\begin{minted}{pycon}
>>> motifs = jdb.fetch_motifs_by_name("Arnt")
>>> print(motifs[0])
TF name	Arnt
Matrix ID	MA0004.1
Collection	CORE
TF class	Zipper-Type
TF family	Helix-Loop-Helix
Species	10090
Taxonomic group	vertebrates
Accession	['P53762']
Data type used	SELEX
Medline	7592839
PAZAR ID	TF0000003
Comments	-
Matrix:
	0      1      2      3      4      5
A:   4.00  19.00   0.00   0.00   0.00   0.00
C:  16.00   0.00  20.00   0.00   0.00   0.00
G:   0.00   1.00   0.00  20.00   0.00  20.00
T:   0.00   0.00   0.00   0.00  20.00   0.00

\end{minted}

The \verb+fetch_motifs+ method allows you to fetch motifs which match a specified set of criteria. These criteria include any of the above described meta information as well as certain matrix properties such as the minimum information content (\verb+min_ic+ in the example below), the minimum length of the matrix or the minimum number of sites used to construct the matrix. Only motifs which pass ALL the specified criteria are returned. Note that selection criteria which correspond to meta information which allow for multiple values may be specified as either a single value or a list of values, e.g. \verb+tax_group+ and \verb+tf_family+ in the example below.
\begin{minted}{pycon}
>>> motifs = jdb.fetch_motifs(
...     collection = "CORE",
...     tax_group = ["vertebrates", "insects"],
...     tf_class = "Winged Helix-Turn-Helix",
...     tf_family = ["Forkhead", "Ets"],
...     min_ic = 12
... )
>>> for motif in motifs:
...     pass # do something with the motif
\end{minted}

\subsubsection*{Compatibility with Perl TFBS modules}

An important thing to note is that the JASPAR \verb+Motif+ class was designed to be compatible with the popular \href{http://tfbs.genereg.net/}{Perl TFBS modules}. Therefore some specifics about the choice of defaults for background and pseudocounts as well as how information content is computed and sequences searched for instances is based on this compatibility criteria. These choices are noted in the specific subsections below.

\begin{itemize}
\item\textbf{Choice of background:} \\
The Perl \verb+TFBS+ modules appear to allow a choice of custom background probabilities (although the documentation states that uniform background is assumed). However the default is to use a uniform background. Therefore it is recommended that you use a uniform background for computing the position-specific scoring matrix (PSSM). This is the default when using the Biopython \verb+motifs+ module.
\item\textbf{Choice of pseudocounts:} \\
By default, the Perl \verb+TFBS+ modules use a pseudocount equal to $\sqrt{N} * \textrm{bg}[\textrm{nucleotide}]$, where $N$ represents the total number of sequences used to construct the matrix. To apply this same pseudocount formula, set the motif \verb+pseudocounts+ attribute using the \verb+jaspar.calculate\_pseudcounts()+ function:
\begin{minted}{pycon}
>>> motif.pseudocounts = motifs.jaspar.calculate_pseudocounts(motif)
\end{minted}
Note that it is possible for the counts matrix to have an unequal number of sequences making up the columns. The pseudocount computation uses the average number of sequences making up the matrix. However, when \verb+normalize+ is called on the counts matrix, each count value in a column is divided by the total number of sequences making up that specific column, not by the average number of sequences. This differs from the Perl \verb+TFBS+ modules because the normalization is not done as a separate step and so the average number of sequences is used throughout the computation of the pssm. Therefore, for matrices with unequal column counts, the PSSM computed by the \verb+motifs+ module will differ somewhat from the pssm computed by the Perl \verb+TFBS+ modules.
\item\textbf{Computation of matrix information content:} \\
The information content (IC) or specificity of a matrix is computed using the \verb+mean+ method of the \verb+PositionSpecificScoringMatrix+ class. However of note, in the Perl \verb+TFBS+ modules the default behaviour is to compute the IC without first applying pseudocounts, even though by default the PSSMs are computed using pseudocounts as described above.
\item\textbf{Searching for instances:} \\
Searching for instances with the Perl \verb+TFBS+ motifs was usually performed using a relative score threshold, i.e. a score in the range 0 to 1. In order to compute the absolute PSSM score corresponding to a relative score one can use the equation:
\begin{minted}{pycon}
>>> abs_score =  (pssm.max - pssm.min) * rel_score + pssm.min
\end{minted}
To convert the absolute score of an instance back to a relative score, one can use the equation:
\begin{minted}{pycon}
>>> rel_score = (abs_score - pssm.min) / (pssm.max - pssm.min)
\end{minted}
For example, using the Arnt motif before, let's search a sequence with a relative score threshold of 0.8.
%TODO - Check missing ... lines, make into a doctest?
\begin{minted}{pycon}
>>> test_seq=Seq("TAAGCGTGCACGCGCAACACGTGCATTA")
>>> arnt.pseudocounts = motifs.jaspar.calculate_pseudocounts(arnt)
>>> pssm = arnt.pssm
>>> max_score = pssm.max
>>> min_score = pssm.min
>>> abs_score_threshold = (max_score - min_score) * 0.8 + min_score
>>> for position, score in pssm.search(test_seq,
                                       threshold=abs_score_threshold):
...     rel_score = (score - min_score) / (max_score - min_score)
...     print("Position %d: score = %5.3f, rel. score = %5.3f" % (
            position, score, rel_score))
...
Position 2: score = 5.362, rel. score = 0.801
Position 8: score = 6.112, rel. score = 0.831
Position -20: score = 7.103, rel. score = 0.870
Position 17: score = 10.351, rel. score = 1.000
Position -11: score = 10.351, rel. score = 1.000
\end{minted}
\end{itemize}

\subsection{MEME}

MEME \cite{bailey1994} is a tool for discovering motifs in a group of related
DNA or protein sequences. It takes as input a group of DNA or protein sequences
and outputs as many motifs as requested. Therefore, in contrast to JASPAR
files, MEME output files typically contain multiple motifs. This is an example.

At the top of an output file generated by MEME shows some background information
about the MEME and the version of MEME used:
\begin{minted}{text}
********************************************************************************
MEME - Motif discovery tool
********************************************************************************
MEME version 3.0 (Release date: 2004/08/18 09:07:01)
...
\end{minted}
Further down, the input set of training sequences is recapitulated:
\begin{minted}{text}
********************************************************************************
TRAINING SET
********************************************************************************
DATAFILE= INO_up800.s
ALPHABET= ACGT
Sequence name            Weight Length  Sequence name            Weight Length
-------------            ------ ------  -------------            ------ ------
CHO1                     1.0000    800  CHO2                     1.0000    800
FAS1                     1.0000    800  FAS2                     1.0000    800
ACC1                     1.0000    800  INO1                     1.0000    800
OPI3                     1.0000    800
********************************************************************************
\end{minted}
and the exact command line that was used:
\begin{minted}{text}
********************************************************************************
COMMAND LINE SUMMARY
********************************************************************************
This information can also be useful in the event you wish to report a
problem with the MEME software.

command: meme -mod oops -dna -revcomp -nmotifs 2 -bfile yeast.nc.6.freq INO_up800.s
...
\end{minted}
Next is detailed information on each motif that was found:
\begin{minted}{text}
********************************************************************************
MOTIF  1        width =   12   sites =   7   llr = 95   E-value = 2.0e-001
********************************************************************************
--------------------------------------------------------------------------------
        Motif 1 Description
--------------------------------------------------------------------------------
Simplified        A  :::9:a::::3:
pos.-specific     C  ::a:9:11691a
probability       G  ::::1::94:4:
matrix            T  aa:1::9::11:
\end{minted}
To parse this file (stored as \verb+meme.dna.oops.txt+), use

%cont-doctest
\begin{minted}{pycon}
>>> with open("meme.INO_up800.classic.oops.xml") as handle:
...     record = motifs.parse(handle, "meme")
...
\end{minted}
The \verb+motifs.parse+ command reads the complete file directly, so you can
close the file after calling \verb+motifs.parse+.
The header information is stored in attributes:

%cont-doctest
\begin{minted}{pycon}
>>> record.version
'5.0.1'
>>> record.datafile
'common/INO_up800.s'
>>> record.command
'meme common/INO_up800.s -oc results/meme10 -mod oops -dna -revcomp -bfile common/yeast.nc.6.freq -nmotifs 2 -objfun classic -minw 8 -nostatus '
>>> record.alphabet
'ACGT'
>>> record.sequences
['sequence_0', 'sequence_1', 'sequence_2', 'sequence_3', 'sequence_4', 'sequence_5', 'sequence_6']
\end{minted}
The record is an object of the \verb+Bio.motifs.meme.Record+ class.
The class inherits from list, and you can think of \verb+record+ as a list of Motif objects:

%cont-doctest
\begin{minted}{pycon}
>>> len(record)
2
>>> motif = record[0]
>>> print(motif.consensus)
GCGGCATGTGAAA
>>> print(motif.degenerate_consensus)
GSKGCATGTGAAA
\end{minted}
In addition to these generic motif attributes, each motif also stores its
specific information as calculated by MEME. For example,

%cont-doctest
\begin{minted}{pycon}
>>> motif.num_occurrences
7
>>> motif.length
13
>>> evalue = motif.evalue
>>> print("%3.1g" % evalue)
0.2
>>> motif.name
'GSKGCATGTGAAA'
>>> motif.id
'motif_1'
\end{minted}
In addition to using an index into the record, as we did above,
you can also find it by its name:

%cont-doctest
\begin{minted}{pycon}
>>> motif = record["GSKGCATGTGAAA"]
\end{minted}
Each motif has an attribute \verb+.instances+ with the sequence instances
in which the motif was found, providing some information on each instance:

%cont-doctest
\begin{minted}{pycon}
>>> len(motif.instances)
7
>>> motif.instances[0]
Instance('GCGGCATGTGAAA')
>>> motif.instances[0].motif_name
'GSKGCATGTGAAA'
>>> motif.instances[0].sequence_name
'INO1'
>>> motif.instances[0].sequence_id
'sequence_5'
>>> motif.instances[0].start
620
>>> motif.instances[0].strand
'+'
>>> motif.instances[0].length
13
>>> pvalue = motif.instances[0].pvalue
>>> print("%5.3g" % pvalue)
1.21e-08
\end{minted}

\subsubsection*{MAST}

\subsection{TRANSFAC}

TRANSFAC is a manually curated database of transcription factors, together
with their genomic binding sites and DNA binding profiles \cite{matys2003}.
While the file format used in the TRANSFAC database is nowadays also used
by others, we will refer to it as the TRANSFAC file format.

A minimal file in the TRANSFAC format looks as follows:
\begin{minted}{text}
ID  motif1
P0      A      C      G      T
01      1      2      2      0      S
02      2      1      2      0      R
03      3      0      1      1      A
04      0      5      0      0      C
05      5      0      0      0      A
06      0      0      4      1      G
07      0      1      4      0      G
08      0      0      0      5      T
09      0      0      5      0      G
10      0      1      2      2      K
11      0      2      0      3      Y
12      1      0      3      1      G
//
\end{minted}
This file shows the frequency matrix of motif \verb+motif1+ of 12 nucleotides.
In general, one file in the TRANSFAC format can contain multiple motifs. For
example, this is the contents of the example TRANSFAC file \verb+transfac.dat+:
\begin{minted}{text}
VV  EXAMPLE January 15, 2013
XX
//
ID  motif1
P0      A      C      G      T
01      1      2      2      0      S
02      2      1      2      0      R
03      3      0      1      1      A
...
11      0      2      0      3      Y
12      1      0      3      1      G
//
ID  motif2
P0      A      C      G      T
01      2      1      2      0      R
02      1      2      2      0      S
...
09      0      0      0      5      T
10      0      2      0      3      Y
//
\end{minted}
To parse a TRANSFAC file, use

%cont-doctest
\begin{minted}{pycon}
>>> with open("transfac.dat") as handle:
...     record = motifs.parse(handle, "TRANSFAC")
...
\end{minted}
If any discrepancies between the file contents and the TRANSFAC file format are detected, a \verb+ValueError+ is raised. Note that you may encounter files that do not follow the TRANSFAC format strictly. For example, the number of spaces between columns may be different, or a tab may be used instead of spaces. Use \verb+strict=False+ to enable parsing such files without raising a \verb+ValueError+:
\begin{minted}{pycon}
>>> record = motifs.parse(handle, "TRANSFAC", strict=False)
\end{minted}
When parsing a non-compliant file, we recommend to check the record returned by \verb+motif.parse+ to ensure that it is consistent with the file contents.

The overall version number, if available, is stored as \verb+record.version+:

%cont-doctest
\begin{minted}{pycon}
>>> record.version
'EXAMPLE January 15, 2013'
\end{minted}

Each motif in \verb+record+ is in instance of the \verb+Bio.motifs.transfac.Motif+
class, which inherits both from the \verb+Bio.motifs.Motif+ class and
from a Python dictionary. The dictionary uses the two-letter keys to
store any additional information about the motif:

%cont-doctest
\begin{minted}{pycon}
>>> motif = record[0]
>>> motif.degenerate_consensus # Using the Bio.motifs.Motif property
Seq('SRACAGGTGKYG')
>>> motif["ID"] # Using motif as a dictionary
'motif1'
\end{minted}

TRANSFAC files are typically much more elaborate than this example, containing
lots of additional information about the motif. Table \ref{table:transfaccodes}
lists the two-letter field codes that are commonly found in TRANSFAC files:
\begin{table}[h]
\label{table:transfaccodes}
\begin{center}
\caption{Fields commonly found in TRANSFAC files}
\begin{tabular}{|l|l||}
\verb+AC+ & Accession number \\
\verb+AS+ & Accession numbers, secondary \\
\verb+BA+ & Statistical basis \\
\verb+BF+ & Binding factors \\
\verb+BS+ & Factor binding sites underlying the matrix \\
\verb+CC+ & Comments \\
\verb+CO+ & Copyright notice \\
\verb+DE+ & Short factor description \\
\verb+DR+ & External databases \\
\verb+DT+ & Date created/updated \\
\verb+HC+ & Subfamilies \\
\verb+HP+ & Superfamilies \\
\verb+ID+ & Identifier \\
\verb+NA+ & Name of the binding factor \\
\verb+OC+ & Taxonomic classification \\
\verb+OS+ & Species/Taxon \\
\verb+OV+ & Older version \\
\verb+PV+ & Preferred version \\
\verb+TY+ & Type \\
\verb+XX+ & Empty line; these are not stored in the Record. \\
\end{tabular}
\end{center}
\end{table}

Each motif also has an attribute \verb+.references+ containing the
references associated with the motif, using these two-letter keys:

\begin{table}[h]
\begin{center}
\caption{Fields used to store references in TRANSFAC files}
\begin{tabular}{|l|l||}
\verb+RN+ & Reference number \\
\verb+RA+ & Reference authors \\
\verb+RL+ & Reference data \\
\verb+RT+ & Reference title \\
\verb+RX+ & PubMed ID \\
\end{tabular}
\end{center}
\end{table}

Printing the motifs writes them out in their native TRANSFAC format:

%cont-doctest
\begin{minted}{pycon}
>>> print(record)
VV  EXAMPLE January 15, 2013
XX
//
ID  motif1
XX
P0      A      C      G      T
01      1      2      2      0      S
02      2      1      2      0      R
03      3      0      1      1      A
04      0      5      0      0      C
05      5      0      0      0      A
06      0      0      4      1      G
07      0      1      4      0      G
08      0      0      0      5      T
09      0      0      5      0      G
10      0      1      2      2      K
11      0      2      0      3      Y
12      1      0      3      1      G
XX
//
ID  motif2
XX
P0      A      C      G      T
01      2      1      2      0      R
02      1      2      2      0      S
03      0      5      0      0      C
04      3      0      1      1      A
05      0      0      4      1      G
06      5      0      0      0      A
07      0      1      4      0      G
08      0      0      5      0      G
09      0      0      0      5      T
10      0      2      0      3      Y
XX
//
<BLANKLINE>
\end{minted}
You can export the motifs in the TRANSFAC format by capturing this output
in a string and saving it in a file:
\begin{minted}{pycon}
>>> text = str(record)
>>> with open("mytransfacfile.dat", "w") as out_handle:
...     out_handle.write(text)
...
\end{minted}

\section{Writing motifs}

Speaking of exporting, let's look at export functions in general.
We can use the \verb+format+ built-in function to write the motif in the simple JASPAR \verb+pfm+ format:
%the tabs in the output confuse doctest; don't test
\begin{minted}{pycon}
>>> print(format(arnt, "pfm"))
  4.00  19.00   0.00   0.00   0.00   0.00
 16.00   0.00  20.00   0.00   0.00   0.00
  0.00   1.00   0.00  20.00   0.00  20.00
  0.00   0.00   0.00   0.00  20.00   0.00
\end{minted}
Similarly, we can use \verb+format+ to write the motif in the JASPAR \verb+jaspar+ format:
\begin{minted}{pycon}
>>> print(format(arnt, "jaspar"))
>MA0004.1  Arnt
A [  4.00  19.00   0.00   0.00   0.00   0.00]
C [ 16.00   0.00  20.00   0.00   0.00   0.00]
G [  0.00   1.00   0.00  20.00   0.00  20.00]
T [  0.00   0.00   0.00   0.00  20.00   0.00]
\end{minted}

To write the motif in a TRANSFAC-like matrix format, use

%cont-doctest
\begin{minted}{pycon}
>>> print(format(m, "transfac"))
P0      A      C      G      T
01      3      0      0      4      W
02      7      0      0      0      A
03      0      5      0      2      C
04      2      2      3      0      V
05      1      6      0      0      C
XX
//
<BLANKLINE>
\end{minted}

To write out multiple motifs, you can use \verb+motifs.write+.
This function can be used regardless of whether the motifs originated from a TRANSFAC file. For example,

%cont-doctest
\begin{minted}{pycon}
>>> two_motifs = [arnt, srf]
>>> print(motifs.write(two_motifs, "transfac"))
P0      A      C      G      T
01      4     16      0      0      C
02     19      0      1      0      A
03      0     20      0      0      C
04      0      0     20      0      G
05      0      0      0     20      T
06      0      0     20      0      G
XX
//
P0      A      C      G      T
01      2      1     39      4      G
02      9     33      2      2      C
03      0     45      1      0      C
04      1     45      0      0      C
05     32      1      0     13      A
06      3      1      0     42      T
07     46      0      0      0      A
08      1      0      0     45      T
09     43      0      0      3      A
10     15      1      0     30      W
11      2      0     44      0      G
12      2      1     43      0      G
XX
//
<BLANKLINE>
\end{minted}

Or, to write multiple motifs in the \verb+jaspar+ format:
\begin{minted}{pycon}
>>> two_motifs = [arnt, mef2a]
>>> print(motifs.write(two_motifs, "jaspar"))
>MA0004.1  Arnt
A [  4.00  19.00   0.00   0.00   0.00   0.00]
C [ 16.00   0.00  20.00   0.00   0.00   0.00]
G [  0.00   1.00   0.00  20.00   0.00  20.00]
T [  0.00   0.00   0.00   0.00  20.00   0.00]
>MA0052.1  MEF2A
A [  1.00   0.00  57.00   2.00   9.00   6.00  37.00   2.00  56.00   6.00]
C [ 50.00   0.00   1.00   1.00   0.00   0.00   0.00   0.00   0.00   0.00]
G [  0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   2.00  50.00]
T [  7.00  58.00   0.00  55.00  49.00  52.00  21.00  56.00   0.00   2.00]
\end{minted}

\section{Position-Weight Matrices}

The \verb+.counts+ attribute of a Motif object shows how often each
nucleotide appeared at each position along the alignment. We can
normalize this matrix by dividing by the number of instances in the
alignment, resulting in the probability of each nucleotide at each
position along the alignment. We refer to these probabilities as
the position-weight matrix. However, beware that in the literature
this term may also be used to refer to the position-specific scoring
matrix, which we discuss below.

Usually, pseudocounts are added to each position before normalizing.
This avoids overfitting of the position-weight matrix to the limited
number of motif instances in the alignment, and can also prevent
probabilities from becoming zero. To add a fixed pseudocount to all
nucleotides at all positions, specify a number for the
\verb+pseudocounts+ argument:

%cont-doctest
\begin{minted}{pycon}
>>> pwm = m.counts.normalize(pseudocounts=0.5)
>>> print(pwm)
        0      1      2      3      4
A:   0.39   0.83   0.06   0.28   0.17
C:   0.06   0.06   0.61   0.28   0.72
G:   0.06   0.06   0.06   0.39   0.06
T:   0.50   0.06   0.28   0.06   0.06
<BLANKLINE>
\end{minted}
Alternatively, \verb+pseudocounts+ can be a dictionary specifying the
pseudocounts for each nucleotide. For example, as the GC content of
the human genome is about 40\%, you may want to choose the
pseudocounts accordingly:

%cont-doctest
\begin{minted}{pycon}
>>> pwm = m.counts.normalize(pseudocounts={"A":0.6, "C": 0.4, "G": 0.4, "T": 0.6})
>>> print(pwm)
        0      1      2      3      4
A:   0.40   0.84   0.07   0.29   0.18
C:   0.04   0.04   0.60   0.27   0.71
G:   0.04   0.04   0.04   0.38   0.04
T:   0.51   0.07   0.29   0.07   0.07
<BLANKLINE>
\end{minted}
The position-weight matrix has its own methods to calculate the
consensus, anticonsensus, and degenerate consensus sequences:

%cont-doctest
\begin{minted}{pycon}
>>> pwm.consensus
Seq('TACGC')
>>> pwm.anticonsensus
Seq('CCGTG')
>>> pwm.degenerate_consensus
Seq('WACNC')
\end{minted}
Note that due to the pseudocounts, the degenerate consensus sequence
calculated from the position-weight matrix is slightly different
from the degenerate consensus sequence calculated from the instances
in the motif:

%cont-doctest
\begin{minted}{pycon}
>>> m.degenerate_consensus
Seq('WACVC')
\end{minted}
The reverse complement of the position-weight matrix can be calculated directly from the \verb+pwm+:

%cont-doctest
\begin{minted}{pycon}
>>> rpwm = pwm.reverse_complement()
>>> print(rpwm)
        0      1      2      3      4
A:   0.07   0.07   0.29   0.07   0.51
C:   0.04   0.38   0.04   0.04   0.04
G:   0.71   0.27   0.60   0.04   0.04
T:   0.18   0.29   0.07   0.84   0.40
<BLANKLINE>
\end{minted}

\section{Position-Specific Scoring Matrices}

Using the background distribution and PWM with pseudo-counts added,
it's easy to compute the log-odds ratios, telling us what are the log
odds of a particular symbol to be coming from a motif against the
background. We can use the \verb|.log_odds()| method on the position-weight
matrix:

%cont-doctest
\begin{minted}{pycon}
>>> pssm = pwm.log_odds()
>>> print(pssm)
        0      1      2      3      4
A:   0.68   1.76  -1.91   0.21  -0.49
C:  -2.49  -2.49   1.26   0.09   1.51
G:  -2.49  -2.49  -2.49   0.60  -2.49
T:   1.03  -1.91   0.21  -1.91  -1.91
<BLANKLINE>
\end{minted}
Here we can see positive values for symbols more frequent in the motif
than in the background and negative for symbols more frequent in the
background. $0.0$ means that it's equally likely to see a symbol in the
background and in the motif.

This assumes that A, C, G, and T are equally likely in the background. To
calculate the position-specific scoring matrix against a background with
unequal probabilities for A, C, G, T, use the \verb+background+ argument.
For example, against a background with a 40\% GC content, use

%cont-doctest
\begin{minted}{pycon}
>>> background = {"A":0.3,"C":0.2,"G":0.2,"T":0.3}
>>> pssm = pwm.log_odds(background)
>>> print(pssm)
        0      1      2      3      4
A:   0.42   1.49  -2.17  -0.05  -0.75
C:  -2.17  -2.17   1.58   0.42   1.83
G:  -2.17  -2.17  -2.17   0.92  -2.17
T:   0.77  -2.17  -0.05  -2.17  -2.17
<BLANKLINE>
\end{minted}

The maximum and minimum score obtainable from the PSSM are stored in the
\verb+.max+ and \verb+.min+ properties:

%cont-doctest
\begin{minted}{pycon}
>>> print("%4.2f" % pssm.max)
6.59
>>> print("%4.2f" % pssm.min)
-10.85
\end{minted}

The mean and standard deviation of the PSSM scores with respect to a specific
background are calculated by the \verb+.mean+ and \verb+.std+ methods.

%cont-doctest
\begin{minted}{pycon}
>>> mean = pssm.mean(background)
>>> std = pssm.std(background)
>>> print("mean = %0.2f, standard deviation = %0.2f" % (mean, std))
mean = 3.21, standard deviation = 2.59
\end{minted}
A uniform background is used if \verb+background+ is not specified.
The mean is particularly important, as its value is equal to the
Kullback-Leibler divergence or relative entropy, and is a measure for the
information content of the motif compared to the background. As in Biopython
the base-2 logarithm is used in the calculation of the log-odds scores, the
information content has units of bits.

The \verb+.reverse_complement+, \verb+.consensus+, \verb+.anticonsensus+, and
\verb+.degenerate_consensus+ methods can be applied directly to PSSM objects.

\section{Searching for instances}
\label{sec:search}

The most frequent use for a motif is to find its instances in some
sequence. For the sake of this section, we will use an artificial sequence like this:

%cont-doctest
\begin{minted}{pycon}
>>> test_seq=Seq("TACACTGCATTACAACCCAAGCATTA")
>>> len(test_seq)
26
\end{minted}

\subsection{Searching for exact matches}

The simplest way to find instances, is to look for exact matches of
the true instances of the motif:

%cont-doctest
\begin{minted}{pycon}
>>> for pos, seq in m.instances.search(test_seq):
...     print("%i %s" % (pos, seq))
...
0 TACAC
10 TACAA
13 AACCC
\end{minted}
We can do the same with the reverse complement (to find instances on the complementary strand):

%cont-doctest
\begin{minted}{pycon}
>>> for pos, seq in r.instances.search(test_seq):
...     print("%i %s" % (pos, seq))
...
6 GCATT
20 GCATT
\end{minted}

\subsection{Searching for matches using the PSSM score}

It's just as easy to look for positions, giving rise to high log-odds scores against our motif:

%cont-doctest
\begin{minted}{pycon}
>>> for position, score in pssm.search(test_seq, threshold=3.0):
...     print("Position %d: score = %5.3f" % (position, score))
...
Position 0: score = 5.622
Position -20: score = 4.601
Position 10: score = 3.037
Position 13: score = 5.738
Position -6: score = 4.601
\end{minted}
The negative positions refer to instances of the motif found on the
reverse strand of the test sequence, and follow the Python convention
on negative indices. Therefore, the instance of the motif at \verb|pos|
is located at \verb|test_seq[pos:pos+len(m)]| both for positive and for
negative values of \verb|pos|.

You may notice the threshold parameter, here set arbitrarily to
$3.0$. This is in $log_2$, so we are now looking only for words, which
are eight times more likely to occur under the motif model than in the
background. The default threshold is $0.0$, which selects everything
that looks more like the motif than the background.

You can also calculate the scores at all positions along the sequence:
%Don't use a doc test for this as the spacing can differ
\begin{minted}{pycon}
>>> pssm.calculate(test_seq)
array([  5.62230396,  -5.6796999 ,  -3.43177247,   0.93827754,
        -6.84962511,  -2.04066086, -10.84962463,  -3.65614533,
        -0.03370807,  -3.91102552,   3.03734159,  -2.14918518,
        -0.6016975 ,   5.7381525 ,  -0.50977498,  -3.56422281,
        -8.73414803,  -0.09919716,  -0.6016975 ,  -2.39429784,
       -10.84962463,  -3.65614533], dtype=float32)
\end{minted}
In general, this is the fastest way to calculate PSSM scores.
The scores returned by \verb+pssm.calculate+ are for the forward strand
only. To obtain the scores on the reverse strand, you can take the reverse
complement of the PSSM:
\begin{minted}{pycon}
>>> rpssm = pssm.reverse_complement()
>>> rpssm.calculate(test_seq)
array([ -9.43458748,  -3.06172252,  -7.18665981,  -7.76216221,
        -2.04066086,  -4.26466274,   4.60124254,  -4.2480607 ,
        -8.73414803,  -2.26503372,  -6.49598789,  -5.64668512,
        -8.73414803, -10.84962463,  -4.82356262,  -4.82356262,
        -5.64668512,  -8.73414803,  -4.15613794,  -5.6796999 ,
         4.60124254,  -4.2480607 ], dtype=float32)
\end{minted}

\subsection{Selecting a score threshold}

If you want to use a less arbitrary way of selecting thresholds, you
can explore the distribution of PSSM scores. Since the space for a score
distribution grows exponentially with motif length, we are using an
approximation with a given precision to keep computation cost manageable:

%cont-doctest
\begin{minted}{pycon}
>>> distribution = pssm.distribution(background=background, precision=10**4)
\end{minted}
The \verb+distribution+ object can be used to determine a number of different thresholds.
We can specify the requested false-positive rate (probability of ``finding'' a motif instance in background generated sequence):

%cont-doctest
\begin{minted}{pycon}
>>> threshold = distribution.threshold_fpr(0.01)
>>> print("%5.3f" % threshold)
4.009
\end{minted}
or the false-negative rate (probability of ``not finding'' an instance generated from the motif):

%cont-doctest
\begin{minted}{pycon}
>>> threshold = distribution.threshold_fnr(0.1)
>>> print("%5.3f" % threshold)
-0.510
\end{minted}
or a threshold (approximately) satisfying some relation between the false-positive rate and the false-negative rate ($\frac{\textrm{fnr}}{\textrm{fpr}}\simeq t$):

%cont-doctest
\begin{minted}{pycon}
>>> threshold = distribution.threshold_balanced(1000)
>>> print("%5.3f" % threshold)
6.241
\end{minted}
or a threshold satisfying (roughly) the equality between the $-log$ of the
false-positive rate and the information content (as used in patser software by
Hertz and Stormo):

%cont-doctest
\begin{minted}{pycon}
>>> threshold = distribution.threshold_patser()
>>> print("%5.3f" % threshold)
0.346
\end{minted}

For example, in case of our motif, you can get the threshold giving
you exactly the same results (for this sequence) as searching for
instances with balanced threshold with rate of $1000$.

%cont-doctest
\begin{minted}{pycon}
>>> threshold = distribution.threshold_fpr(0.01)
>>> print("%5.3f" % threshold)
4.009
>>> for position, score in pssm.search(test_seq, threshold=threshold):
...     print("Position %d: score = %5.3f" % (position, score))
...
Position 0: score = 5.622
Position -20: score = 4.601
Position 13: score = 5.738
Position -6: score = 4.601
\end{minted}

\section{Each motif object has an associated Position-Specific Scoring Matrix}

To facilitate searching for potential TFBSs using PSSMs, both the position-weight matrix and the position-specific scoring matrix are associated with each motif. Using the Arnt motif as an example:
%TODO - Start a new doctest here?

%cont-doctest
\begin{minted}{pycon}
>>> from Bio import motifs
>>> with open("Arnt.sites") as handle:
...     motif = motifs.read(handle, "sites")
...
>>> print(motif.counts)
        0      1      2      3      4      5
A:   4.00  19.00   0.00   0.00   0.00   0.00
C:  16.00   0.00  20.00   0.00   0.00   0.00
G:   0.00   1.00   0.00  20.00   0.00  20.00
T:   0.00   0.00   0.00   0.00  20.00   0.00
<BLANKLINE>
>>> print(motif.pwm)
        0      1      2      3      4      5
A:   0.20   0.95   0.00   0.00   0.00   0.00
C:   0.80   0.00   1.00   0.00   0.00   0.00
G:   0.00   0.05   0.00   1.00   0.00   1.00
T:   0.00   0.00   0.00   0.00   1.00   0.00
<BLANKLINE>
>>> print(motif.pssm)
        0      1      2      3      4      5
A:  -0.32   1.93   -inf   -inf   -inf   -inf
C:   1.68   -inf   2.00   -inf   -inf   -inf
G:   -inf  -2.32   -inf   2.00   -inf   2.00
T:   -inf   -inf   -inf   -inf   2.00   -inf
<BLANKLINE>
\end{minted}
The negative infinities appear here because the corresponding entry in the frequency matrix is 0, and we are using zero pseudocounts by default:

%cont-doctest
\begin{minted}{pycon}
>>> for letter in "ACGT":
...     print("%s: %4.2f" % (letter, motif.pseudocounts[letter]))
...
A: 0.00
C: 0.00
G: 0.00
T: 0.00
\end{minted}
If you change the \verb+.pseudocounts+ attribute, the position-frequency matrix and the position-specific scoring matrix are recalculated automatically:

%cont-doctest
\begin{minted}{pycon}
>>> motif.pseudocounts = 3.0
>>> for letter in "ACGT":
...     print("%s: %4.2f" % (letter, motif.pseudocounts[letter]))
...
A: 3.00
C: 3.00
G: 3.00
T: 3.00
>>> print(motif.pwm)
        0      1      2      3      4      5
A:   0.22   0.69   0.09   0.09   0.09   0.09
C:   0.59   0.09   0.72   0.09   0.09   0.09
G:   0.09   0.12   0.09   0.72   0.09   0.72
T:   0.09   0.09   0.09   0.09   0.72   0.09
<BLANKLINE>
\end{minted}

%cont-doctest
\begin{minted}{pycon}
>>> print(motif.pssm)
        0      1      2      3      4      5
A:  -0.19   1.46  -1.42  -1.42  -1.42  -1.42
C:   1.25  -1.42   1.52  -1.42  -1.42  -1.42
G:  -1.42  -1.00  -1.42   1.52  -1.42   1.52
T:  -1.42  -1.42  -1.42  -1.42   1.52  -1.42
<BLANKLINE>
\end{minted}
You can also set the \verb+.pseudocounts+ to a dictionary over the four nucleotides if you want to use different pseudocounts for them. Setting \verb+motif.pseudocounts+ to \verb+None+ resets it to its default value of zero.

The position-specific scoring matrix depends on the background distribution, which is uniform by default:

%cont-doctest
\begin{minted}{pycon}
>>> for letter in "ACGT":
...     print("%s: %4.2f" % (letter, motif.background[letter]))
...
A: 0.25
C: 0.25
G: 0.25
T: 0.25
\end{minted}
Again, if you modify the background distribution, the position-specific scoring matrix is recalculated:

%cont-doctest
\begin{minted}{pycon}
>>> motif.background = {"A": 0.2, "C": 0.3, "G": 0.3, "T": 0.2}
>>> print(motif.pssm)
        0      1      2      3      4      5
A:   0.13   1.78  -1.09  -1.09  -1.09  -1.09
C:   0.98  -1.68   1.26  -1.68  -1.68  -1.68
G:  -1.68  -1.26  -1.68   1.26  -1.68   1.26
T:  -1.09  -1.09  -1.09  -1.09   1.85  -1.09
<BLANKLINE>
\end{minted}
Setting \verb+motif.background+ to \verb+None+ resets it to a uniform distribution:

%cont-doctest
\begin{minted}{pycon}
>>> motif.background = None
>>> for letter in "ACGT":
...     print("%s: %4.2f" % (letter, motif.background[letter]))
...
A: 0.25
C: 0.25
G: 0.25
T: 0.25
\end{minted}
If you set \verb+motif.background+ equal to a single value, it will be interpreted as the GC content:

%cont-doctest
\begin{minted}{pycon}
>>> motif.background = 0.8
>>> for letter in "ACGT":
...     print("%s: %4.2f" % (letter, motif.background[letter]))
...
A: 0.10
C: 0.40
G: 0.40
T: 0.10
\end{minted}
Note that you can now calculate the mean of the PSSM scores over the background against which it was computed:

%cont-doctest
\begin{minted}{pycon}
>>> print("%f" % motif.pssm.mean(motif.background))
4.703928
\end{minted}
as well as its standard deviation:

%cont-doctest
\begin{minted}{pycon}
>>> print("%f" % motif.pssm.std(motif.background))
3.290900
\end{minted}
and its distribution:

%cont-doctest
\begin{minted}{pycon}
>>> distribution = motif.pssm.distribution(background=motif.background)
>>> threshold = distribution.threshold_fpr(0.01)
>>> print("%f" % threshold)
3.854375
\end{minted}

Note that the position-weight matrix and the position-specific scoring matrix are recalculated each time you call \verb+motif.pwm+ or \verb+motif.pssm+, respectively. If speed is an issue and you want to use the PWM or PSSM repeatedly, you can save them as a variable, as in
\begin{minted}{pycon}
>>> pssm = motif.pssm
\end{minted}

\section{Comparing motifs}
\label{sec:comp}
Once we have more than one motif, we might want to compare them.

Before we start comparing motifs, I should point out that motif
boundaries are usually quite arbitrary. This means we often need
to compare motifs of different lengths, so comparison needs to involve
some kind of alignment.  This means we have to take into account two things:
\begin{itemize}
\item alignment of motifs
\item some function to compare aligned motifs
\end{itemize}
To align the motifs, we use ungapped alignment of PSSMs and substitute zeros
for any missing columns at the beginning and end of the matrices. This means
that effectively we are using the background distribution for columns missing
from the PSSM.
The distance function then returns the minimal distance between motifs, as
well as the corresponding offset in their alignment.

To give an example, let us first load another motif,
which is similar to our test motif \verb|m|:
%TODO - Start a new doctest here?

%cont-doctest
\begin{minted}{pycon}
>>> with open("REB1.pfm") as handle:
...    m_reb1 = motifs.read(handle, "pfm")
...
>>> m_reb1.consensus
Seq('GTTACCCGG')
>>> print(m_reb1.counts)
        0      1      2      3      4      5      6      7      8
A:  30.00   0.00   0.00 100.00   0.00   0.00   0.00   0.00  15.00
C:  10.00   0.00   0.00   0.00 100.00 100.00 100.00   0.00  15.00
G:  50.00   0.00   0.00   0.00   0.00   0.00   0.00  60.00  55.00
T:  10.00 100.00 100.00   0.00   0.00   0.00   0.00  40.00  15.00
<BLANKLINE>
\end{minted}

To make the motifs comparable, we choose the same values for the pseudocounts and the background distribution as our motif \verb|m|:

%cont-doctest
\begin{minted}{pycon}
>>> m_reb1.pseudocounts = {"A":0.6, "C": 0.4, "G": 0.4, "T": 0.6}
>>> m_reb1.background = {"A":0.3,"C":0.2,"G":0.2,"T":0.3}
>>> pssm_reb1 = m_reb1.pssm
>>> print(pssm_reb1)
        0      1      2      3      4      5      6      7      8
A:   0.00  -5.67  -5.67   1.72  -5.67  -5.67  -5.67  -5.67  -0.97
C:  -0.97  -5.67  -5.67  -5.67   2.30   2.30   2.30  -5.67  -0.41
G:   1.30  -5.67  -5.67  -5.67  -5.67  -5.67  -5.67   1.57   1.44
T:  -1.53   1.72   1.72  -5.67  -5.67  -5.67  -5.67   0.41  -0.97
<BLANKLINE>
\end{minted}
We'll compare these motifs using the Pearson correlation.
Since we want it to resemble a distance measure, we actually take
$1-r$, where $r$ is the Pearson correlation coefficient (PCC):

%cont-doctest
\begin{minted}{pycon}
>>> distance, offset = pssm.dist_pearson(pssm_reb1)
>>> print("distance = %5.3g" % distance)
distance = 0.239
>>> print(offset)
-2
\end{minted}
This means that the best PCC between motif \verb|m| and  \verb|m_reb1| is obtained with the following alignment:
\begin{minted}{text}
m:      bbTACGCbb
m_reb1: GTTACCCGG
\end{minted}
where \verb|b| stands for background distribution. The PCC itself is
roughly $1-0.239=0.761$.

\section{\emph{De novo} motif finding}
\label{sec:find}

Currently, Biopython has only limited support for \emph{de novo} motif
finding. Namely, we support running \verb|xxmotif| and also parsing of
MEME. Since the number of motif finding tools is growing rapidly,
contributions of new parsers are welcome.

\subsection{MEME}
\label{sec:meme}

Let's assume, you have run MEME on sequences of your choice with your
favorite parameters and saved the output in the file
\verb|meme.out|. You can retrieve the motifs reported by MEME by
running the following piece of code:

%doctest ../Tests/motifs
\begin{minted}{pycon}
>>> from Bio import motifs
>>> with open("meme.psp_test.classic.zoops.xml") as handle:
...     motifsM = motifs.parse(handle, "meme")
...
\end{minted}
\begin{minted}{pycon}
>>> motifsM
[<Bio.motifs.meme.Motif object at 0xc356b0>]
\end{minted}

Besides the most wanted list of motifs, the result object contains more useful information, accessible through properties with self-explanatory names:
\begin{itemize}
\item \verb|.alphabet|
\item \verb|.datafile|
\item \verb|.sequences|
\item \verb|.version|
\item \verb|.command|
\end{itemize}

The motifs returned by the MEME Parser can be treated exactly like regular
Motif objects (with instances), they also provide some extra
functionality, by adding additional information about the instances.

%cont-doctest
\begin{minted}{pycon}
>>> motifsM[0].consensus
Seq('GCTTATGTAA')
>>> motifsM[0].instances[0].sequence_name
'iYFL005W'
>>> motifsM[0].instances[0].sequence_id
'sequence_15'
>>> motifsM[0].instances[0].start
480
>>> motifsM[0].instances[0].strand
'+'
\end{minted}
\begin{minted}{pycon}
>>> motifsM[0].instances[0].pvalue
1.97e-06
\end{minted}

\section{Useful links}
\label{sec:links}

\begin{itemize}
\item \href{https://en.wikipedia.org/wiki/Sequence_motif}{Sequence motif} in wikipedia
\item \href{https://en.wikipedia.org/wiki/Position_weight_matrix}{PWM} in wikipedia
\item \href{https://en.wikipedia.org/wiki/Consensus_sequence}{Consensus sequence} in wikipedia
\item \href{http://bio.cs.washington.edu/assessment/}{Comparison of different motif finding programs}
\end{itemize}

max_stars_repo_name: eric-erki/biopython

max_issues_repo_name: SONNYROBOT/biopython

max_forks_repo_name: SONNYROBOT/biopython


text:
%% Copyright (C) 2010, 2011, Gostai S.A.S.
%%
%% This software is provided "as is" without warranty of any kind,
%% either expressed or implied, including but not limited to the
%% implied warranties of fitness for a particular purpose.
%%
%% See the LICENSE file for more information.

\section{Loadable}

Loadable objects can be switched on and off --- typically physical
devices.

\subsection{Example}

The intended use is rather as follows:

\begin{urbiscript}
import gsrapi.*;
class Motor: Loadable
{
  var val = 0;
  function go(var d)
  {
    if (load)
      val += d
    else
      echo("cannot advance, the motor is off")|;
  };
};
[00000002] Motor

var m = Motor.new();
[00000003] Motor_0x42364388

m.load;
[00000004] false

m.go(1);
[00000006] *** cannot advance, the motor is off

m.on();
[00000007] Motor_0x42364388

m.go(123);
m.val;
[00000009] 123
\end{urbiscript}

\subsection{Prototypes}

\begin{refObjects}
\item[Object]
\end{refObjects}

\subsection{Construction}

\lstinline|Loadable| can be constructed, but it hardly makes sense.
This object should serve as a prototype.

\subsection{Slots}

\begin{urbiscriptapi}
\item[load] The current status.

\item[off](<val>)%
  Set \refSlot{load} to \lstinline|false| and return \this.
\begin{urbiassert}
do (Loadable.new())
{
  assert
  {
    !load;
    off() === this;
    !load;
    on() === this;
    load;
    off() === this;
    !load;
  };
};
\end{urbiassert}

\item[on](<val>)%
  Set \refSlot{load} to \lstinline|true| and return \this.
\begin{urbiassert}
do (Loadable.new())
{
  assert
  {
    !load;
    on() === this;
    load;
    on() === this;
    load;
  };
};
\end{urbiassert}

\item[toggle]%
  Set \refSlot{load} from \lstinline|true| to \lstinline|false|, and
  vice-versa.  Return \var{val}.
\begin{urbiassert}
do (Loadable.new())
{
  assert
  {
    !load;
    toggle() === this;
    load;
    toggle() === this;
    !load;
  };
};
\end{urbiassert}
\end{urbiscriptapi}

%%% Local Variables:
%%% coding: utf-8
%%% mode: latex
%%% TeX-master: "../urbi-sdk"
%%% ispell-dictionary: "american"
%%% ispell-personal-dictionary: "../urbi.dict"
%%% fill-column: 76
%%% End:

max_stars_repo_name: jcbaillie/urbi

max_issues_repo_name: jcbaillie/urbi

max_forks_repo_name: jcbaillie/urbi


text:
\documentclass{ximera}
\input{../preamble}
\title{Substitution and Tables}
%%%%%\author{Philip T. Gressman}

\begin{document}
\begin{abstract}
We review substitution and the use of integral tables.
\end{abstract}
\maketitle

\section*{(Videos) Calculus: Single Variable}
\textbf{Note: The section on the Gompertz Model (7:08--11:31) relates to ideas that we will study later but have not yet seen. It is recommended that you do your best to understand it now and then come back to it again when we study ODEs.}
\youtube{q_CmUbctfzo}
\youtube{T-WgpMM813E}

\section*{Online Texts}
\begin{itemize}
\item \link[OpenStax II 1.5: Substitution]{https://openstax.org/books/calculus-volume-2/pages/1-5-substitution} and  \link[OpenStax II 3.5: Tables and Computer Systems]{https://openstax.org/books/calculus-volume-2/pages/3-5-other-strategies-for-integration}
\item \link[Ximera OSU: Substitution]{https://ximera.osu.edu/mooculus/calculus1/substitution/titlePage}
\item \link[Community Calculus 8.1: Substitution]{https://www.whitman.edu/mathematics/calculus_online/section08.01.html}
\end{itemize}

\section*{Examples}

\begin{example}%[washeralgsub2001]
The region in the plane given by \[\left|{\frac{x}{2} - \frac{1}{2 \sqrt{\frac{x^{2}}{3} + 1}}}\right| \leq y \leq \frac{x}{2} + \frac{1}{2 \sqrt{\frac{x^{2}}{3} + 1}}\] and \(0 \leq x \leq 3\) is revolved around the \(x\)-axis. Compute the volume of the resulting solid.

If the variable \(x\) is used for slicing, then slices are perpendicular to the axis of rotation, which indicates the washer method should be used.
The inequalities for \(y\) give the outer and inner radii, and \[ \left(\frac{x}{2} + \frac{1}{2 \sqrt{\frac{x^{2}}{3} + 1}}\right)^2 - \left( \left|{\frac{x}{2} - \frac{1}{2 \sqrt{\frac{x^{2}}{3} + 1}}}\right|\right)^2 = \frac{\sqrt{3} x}{\sqrt{x^{2} + 3}}. \] (Note that the absolute values go away when the radius is squared.)
This leads to the conclusion
\[ V=  \int_{0}^{3} \frac{\sqrt{3} \pi x}{\sqrt{x^{2} + 3}} dx. \]
Among the options below, the best choice for a potential substitution is
\begin{multipleChoice}
\choice{$u = \sqrt{3x}$}
\choice{$u = 1/\sqrt{x^2+3}$}
\choice{$u = \sqrt{x^2 + 3}$}
\choice[correct]{$u = x^2+3$}
\end{multipleChoice}
Using the substitution above, the differentials satisfy \(du = \answer{2x}dx\). This gives the equality
\[  \int \frac{\sqrt{3} \pi x}{\sqrt{x^{2} + 3}}\, dx  = \int \answer{\frac{\sqrt{3} \pi}{2 \sqrt{u}}} du = \answer{\sqrt{3} \pi \sqrt{u}}.  \]
Reversing the substitution gives
\[ \begin{aligned} \int_{0}^{3} \frac{\sqrt{3} \pi x}{\sqrt{x^{2} + 3}}\, dx & = \left. \left[\answer{\sqrt{3} \pi \sqrt{x^{2} + 3}} \right] \right|_{0}^{3}\\ & =  \answer{3 \pi}. \end{aligned} \]
\textbf{Note:} We would not have to reverse the substitution if we also determined the new bounds. In this case, if $u = x^2 + 3$ and $x = 0$, then $u = \answer{3}$. Likewise if $x = 3$, then $u = \answer{12}$. Thus we could also have carried out the calculation by changing bounds:
\[ \int_{0}^{3} \frac{\sqrt{3} \pi x}{\sqrt{x^{2} + 3}}\, dx = \int_{\answer{3}}^{\answer{12}} \answer{\frac{\sqrt{3} \pi}{2 \sqrt{u}}} du. \]
\end{example}

\begin{example}%%From the corresponding worksheet
Using the table of integrals below, compute the indefinite integral
\[ \int \frac{d \theta}{(1-\theta) \sqrt{\theta^2 - 2 \theta - 3}}.\]
\begin{itemize}
\item When dealing with quadratic expressions such as the ones appearing in this integrand, it is often necessary to complete the square before appealing to a table. In this case,
\[ \theta^2 - 2 \theta - 3 = \left( \theta - \answer{1} \right)^2 - \answer{4}. \]
\item Using the most appropriate entry of the table and plugging in the correct value of $a$ gives
\[ \int \frac{dx}{x \sqrt{x^2 - \answer{4}}} = \frac{1}{\answer{2}} \mathrm{arcsec} \frac{x}{\answer{2}} + C. \]
\item Based on the results of completing the square, we make a substitution $x = \theta - \answer{1}$ and conclude
\[ \int \frac{d \theta}{(1-\theta) \sqrt{\theta^2 - 2 \theta - 3}} =  \int \frac{dx}{\answer{-x} \sqrt{x^2 - \answer{4}}} = \answer{- \frac{1}{2}} \mathrm{arcsec} \answer{\frac{\theta-1}{2}} + C. \]
\end{itemize}
\end{example}

\begin{example}%[2019MT1.integraltable]
Use a table of integrals to compute the antiderivative below.
\[ \int \frac{x^2 ~ dx}{\sqrt{16 + x^6}} \]
\begin{itemize}
\item The key is to make the substitution $u = x^{\answer{3}}$ so that the expression $16 + x^6$ can be understood as a quadratic function of $u$.
\item Specifically $du = \answer{3x^2} dx$, so
\[ \int \frac{x^2 ~ dx}{\sqrt{16 + x^6}} = \frac{1}{3} \int \frac{du}{\sqrt{16+u^2}}. \]
Aside from the factor of $1/3$, the $u$ integral belongs to the table:
\[ \int \frac{du}{\sqrt{16+u^2}} = \answer{\ln | u + \sqrt{16 + u^2}| } + C \]
\item We conclude that
\[  \int \frac{x^2 ~ dx}{\sqrt{16 + x^6}} = \frac{1}{3} \int \frac{du}{\sqrt{16+u^2}} = \frac{\ln ( x^3 + \sqrt{16 + x^6})}{3} + C. \]
Absolute values are not needed in the logarithm because $x^3 + \sqrt{16 + x^6}$ can never be negative.
\end{itemize}
\end{example}

\begin{center}
\textbf{A Basic Table of Integrals} \\ \label{table_of_integrals}
(Note: To use these tables, $k$ and $a$ must represent \textit{constants} in the integral you wish to compute and cannot depend on the variable of integration.)
\begin{tabular}{|lr|}
\hline
\vspace{-10pt}&   \\
$\displaystyle \int x^k dx = \frac{x^{k+1}}{k+1} + C$ & $(k \neq -1)$ \\ $\displaystyle \int \frac{1}{x} dx = \ln |x| + C$ & \\
$\displaystyle \int a^x dx = \frac{a^x}{\ln a} + C$  &  $(a > 0, a \neq 1)$  \\   $\displaystyle \int e^x dx = e^x + C$ & \\
$\displaystyle \int x^k \ln |x| dx = \frac{x^{k+1}}{k+1} \ln |x| - \frac{x^{k+1}}{(k+1)^2} + C$  \hspace{-8pt} & $(k \neq -1)$  \\  $\displaystyle \int \frac{dx}{x^2 + a^2} = \frac{1}{a} \arctan \frac{x}{a} + C$ & $(a \neq 0)$ \\
$\displaystyle \int \frac{dx}{\sqrt{a^2 + x^2}} = \ln \left|x + \sqrt{a^2 + x^2}\right| + C$ & \\ $\displaystyle \int \frac{dx}{\sqrt{a^2 - x^2}} = \arcsin \frac{x}{a} + C$ & $(a > 0, |x| < a)$  \\
$\displaystyle \int \frac{dx}{\sqrt{x^2 - a^2}} = \ln \left| x + \sqrt{x^2-a^2} \right| + C$ & $(|x| > |a|)$ \\ $\displaystyle \int \frac{dx}{x \sqrt{x^2 - a^2}} = \frac{1}{a} \mathrm{arcsec} \frac{x}{a} + C$ & $(a > 0, x > a)$ \\
$\displaystyle \int \sin x ~dx = - \cos x + C$ & \\ $\displaystyle \int \cos x ~dx = \sin x + C$ & \\
$\displaystyle \int \tan x ~dx = \ln |\sec x| + C$ & \\ $\displaystyle \int \csc x ~dx = - \ln |\csc x + \cot x| + C$ \hspace{-18pt} & \\
$\displaystyle \int \sec x ~dx = \ln |\sec x + \tan x| + C$ & \\ $\displaystyle \int \cot x ~dx = \ln |\sin x| + C$ & \\
$\displaystyle \int \sec^2 x ~ dx = \tan x + C$ & \\ $\displaystyle \int \csc^2 x ~dx = - \cot x + C$ & \\
$\displaystyle \int \sec x \tan x dx = \sec x + C$  & \\ $\displaystyle \int \csc x \cot x dx = - \csc x + C$ &  \vspace{-8pt} \\
&   \\
\hline
\end{tabular}
\end{center}

\end{document}