text:
\chapter{Computer program analysis of data science code}
\label{ch:program-analysis}

Computer programs, though precise enough to be unambiguously executed, are
primarily written by humans in order to be intelligible to humans. Whether
expressed in their native textual format, or in the parsed format of abstract
syntax trees, computer programs are for many purposes nearly as unintelligible
to machines as is natural language text. In particular, machines do not readily
create, introspect, or manipulate code in intelligent ways.

The difficulty of computer program understanding has several sources, some more
fundamental than others. As a practical matter, the programming languages in
common use tend to be complex. Popular languages offer many features and
conveniences, such as classes, special operators, and flavors of ``syntactic
sugar,'' that are appreciated by programmers but go considerably beyond the
simple models of computation studied by theoretical computer scientists. This is
certainly true of Python and, to a lesser extent, R. Moreover, no matter how
simple or complex the language, important properties of a program are rarely
immediate from its syntactic presentation but must be inferred through
\emph{computer program analysis}. Such properties include control flow and data
flow, opportunities for parallelism, and guarantees about correctness and
termination. Foundational results, like the famous unsolvability of halting
problem, ensure that many program properties are generally impossible to decide
algorithmically, without executing the program. Finally, when viewed as human
artifacts, computer programs embody abstract concepts and domain knowledge and
so possess a meaning going beyond what is normally understood as program
semantics in programming language theory. Intelligent inspection and
manipulation of a program requires some understanding of its meaning in this
more nebulous sense.

This chapter concerns the first two problems, of managing programming language
complexity and inferring program properties, particularly data flow. It
describes the design and implementation of a software system for transforming
programs in the Python and R languages into data-flow diagrams, introduced as
\emph{raw flow graphs} in \cref{sec:introduction-semantic-enrichment}. While
assuming nothing specifically about the subject matter of the code, the system
comprises the first half of a method for the semantic analysis of data science
code. The second half, the topic of the next chapter, introduces concepts
specific to data science.

In a departure from the rest of the dissertation, the program analysis
methodology is not presented in a formal or mathematically rigorous style. This
appears to be an inevitable feature of the domain. Like many practical
programming languages, the Python language is too complex and ambiguously
specified to admit a useful formal semantics.\footnote{In his master's thesis,
  Guth gives a large, but incomplete, operational semantics of Python 3.3
  \cite{guth2013}.} The R language is simpler than Python, but still quite
complex compared to idealized mathematical models of
computation.\footnote{Operational semantics of ``core R,'' a proper subset of
  the R language, is given in \cite{morandat2012}.} No attempt is made to
formalize the process of abstraction and simplification by which Python and R
are reduced to an idealized programming language. Indeed, the main purpose of
this process to create a computational representation that is easy reason
precisely about and manipulate.

\section{Recording flow graphs using program analysis}
\label{sec:program-analysis-general}

Raw flow graphs, or in this chapter simply \emph{flow graphs}, model the flow of
data during the execution of a computer program. Several small examples have
been given in
\cref{lst:toy-kmeans-scipy,lst:toy-kmeans-sklearn,lst:toy-kmeans-r} and
\cref{fig:toy-kmeans-scipy,fig:toy-kmeans-sklearn,fig:toy-kmeans-r} in the
Introduction. Another such example, of fitting a linear regression model using
Scikit-learn, is shown in \cref{lst:toy-regression-sklearn} and
\cref{fig:toy-regression-sklearn-raw}. The correspondence between the source
listing and the flow graph should be mostly self-evident.\footnote{The names
  \verb|read_csv| in the code and \verb|_make_parser_function.<locals>.parser_f|
  in the flow graph do not match because the \verb|read_csv| and
  \verb|read_table| functions in Pandas are dynamically generated from an
  internal higher-order function called \verb|_make_parser_function|. The
  obscure name in the flow graph is the function's ``true name'' according to
  the Python interpreter.}

\begin{listing}
  \begin{tcolorbox}
    \inputminted{python}{code/regression-sklearn-metrics.py}
  \end{tcolorbox}
  \caption{Linear regression in Python, using Pandas and Scikit-learn}
  \label{lst:toy-regression-sklearn}
\end{listing}

\begin{figure}
  \centering
  \input{wiring-diagrams/program-analysis/regression-sklearn-metrics.raw}
  \caption{Raw flow graph for \cref{lst:toy-regression-sklearn}}
  \label{fig:toy-regression-sklearn-raw}
\end{figure}

Flow graphs are \emph{wiring diagrams}, also known as \emph{string diagrams}.
Formally, they represent morphisms in a cartesian closed category, a concept
partially introduced in \cref{sec:cartesian-categories} and elaborated on in the
next chapter. From a programmer's perspective, flow graphs are an idealized
model of typed functional programming, in which functions can be defined and
composed in sequence and in parallel. Data is permitted to be duplicated or
discarded.\footnote{In wiring diagrams, the duplication and discarding
  operations are represented implicitly, as ports having multiple or zero
  incident wires, rather than explicitly, as the filled circles used in
  \cref{ch:category-theory,ch:algebra-statistics,ch:zoo-statistics}.} The type
system has basic or primitive types, product types, a unit or singleton type,
and function types.

The wires in a flow graph represent typed objects and the boxes represent
function calls in the target programming language. The types can be primitive,
such as \verb|float| or \verb|str| in Python and \verb|numeric| or
\verb|character| in R, or externally defined classes, such as NumPy's core
datatype \verb|ndarray| or R's linear model type \verb|lm|. As for the boxes,
the phrase ``function call'' is used liberally to encompass essentially any
user-invoked computational action. So, in a typical object-oriented programming
language like Python, the ``functions'' include not just standalone functions,
but also static methods, instance methods, object attribute accessors, container
indexing, and special unary and binary operators. This understanding of
functions is more closely aligned with the R programming language, with its
slogan that ``everything that happens in R is a function call''
\cite{chambers2016}.

The flow graph is constructed by two stages of program analysis, the first
static and the second dynamic. Program analysis can in general be static or
dynamic or both. A \emph{static} program analysis inspects or transforms the
program without executing it, while a \emph{dynamic} program analysis involves
executing the program. In our system, the static phase transforms the original
program to emit events whenever user-level code performs certain actions, and
then the dynamic phase executes the modified program, assembling the flow graph
based on the observed events. Because both phases are aided by tight integration
with the language interpreter, the Python program analysis system is implemented
in Python and the R system in R. Nevertheless, the two systems are
architecturally similar.

Let us now consider the static and dynamic phases in greater detail.

\paragraph{Static program analysis}

The static phase of program analysis augments the source code with additional
instructions to record whenever
\begin{enumerate}[(i),nosep]
\item a function call is about to begin,
\item a function call is returned from,
\item a variable is accessed,
\item a variable in the local scope is assigned, or
\item a variable in the local scope is deleted.\footnote{Deletion is performed
    by the \verb|del| keyword in Python and the \verb|rm| function in R. Since
    these languages have garbage collectors, explicit deletion is rarely
    encountered in practice.}
\end{enumerate}
The code is transformed by compiling it to an \emph{abstract syntax tree (AST)}
using the language's built-in parser, walking the syntax tree to find
occurrences of function calls and variable uses, and modifying those nodes to
make special callbacks. In Python, the source code is transformed in a single
pass. The resulting AST can then be pretty-printed as new source code or, more
usefully, simply executed directly. In R, the code transformations are made
on-the-fly as the programs executes, which is possible due to R's unique
combination of lazy evaluation and dynamic metaprogramming.

Consider, for concreteness, the rewriting of a function call. In Python, the
function call \verb|f(x,y)| is transformed into the expression
\begin{minted}[autogobble,xleftmargin=0.5in]{python}
  __tracer.trace_return(__tracer.trace_function(f)(
      __tracer.trace_argument(x),
      __tracer.trace_argument(y)
  ))
\end{minted}
where \verb|__tracer| refers to a hidden object that will be injected into the
code's global namespace at runtime. Each of the tracer's methods returns its
argument unmodified after suitable processing. Thus, due to Python's evaluation
order for function calls, the transformed code will
\begin{enumerate}[nosep]
\item evaluate the name \verb|f|, yielding a function or other callable object,
  then call \verb|trace_function| with this value,
\item evaluate the name \verb|x|, then call \verb|trace_argument| with this
  value,
\item evaluate the name \verb|y|, then call \verb|trace_argument| with this
  value, and finally
\item call the object \verb|f| with the objects \verb|x| and \verb|y| as
  arguments, then call \verb|trace_return| with the return value.
\end{enumerate}
The effect is to evaluate the expression \verb|f(x,y)| while capturing enough
information to reconstruct what function was called, what were its arguments,
and what value it returned.

One might wonder why a simpler transformation is not performed instead. For
example, the function call \verb|f(x,y)| might be rewritten as
\verb|trace_call(f,[x,y])| or \verb|trace_call(f,x,y)|, where the function
\verb|trace_call| would evaluate \verb|f| with arguments \verb|x| and \verb|y|,
record all relevant information about the function call, and return the computed
value. However, this can have unwanted effects when the function \verb|f| being
evaluated is not \emph{referentially transparent}, that is, when it depends not
just on its inputs but on some property of its calling environment. A correct
code transformation preserves the call stack of the program by evaluating the
function \verb|f| in its original context, whereas the simpler transformations
suggested here do not. Obscure though it may seem, this issue does arise in
practice. For instance, the Python package Patsy \cite{patsy2018}, which
provides R-style model formulas for the popular package Statsmodels
\cite{statsmodels2010}, inspects the call stack to circumvent Python's lack of
metaprogramming facilities. It will break if the call stack is not preserved.

The prime directive of any program transformation, whether for optimization or
introspection, is to not alter the observed behavior of the program. This can be
surprisingly subtle, as the case above illustrates. Another difficult case,
non-standard evaluation in R, will be considered later.

\paragraph{Dynamic program analysis}

The dynamic phase of program analysis executes the transformed program,
gradually constructing the flow graph as callbacks occur. The sequence of
function calls, defining the boxes of the diagram, is immediate from the
callbacks, but the data flow between function calls, defining the wires of the
diagram, must also be recorded. Data is generally passed between functions
either directly, from one function call to another by function composition, or
indirectly, though variable assignments and accesses. Extra bookkeeping is
needed to track this data flow.

The program analysis system tracks the flow of data both across the call stack
and within a single level of the call stack. An empty wiring diagram is created
when the program begins. Whenever a function is about to be called, if the
function is defined within the program being analyzed, then a new wiring diagram
is created, pushed onto a stack maintained by the system, and the recording
process is restarted recursively. Otherwise, the function must be defined by an
external library or by the language itself and is treated as an \emph{atomic},
or indecomposable, computation. In either case, a new box is added to the outer
wiring diagram,\footnote{As a technical note, in Python, the new box can be
  created and connected \emph{before} the function is called, whereas in R this
  must happen \emph{after} the function returns, because the lazy evaluation of
  function arguments defers their availability for introspection.} the box being
a nested wiring diagram in the first case and an atomic box in the second. The
hierarchy of the wiring diagram thus mirrors the call stack of the program
throughout its execution. Once a new box has been added, the box's inputs are
wired either to outputs of previously created boxes or to inputs of the outer
box, according to the provenance of the arguments to the function call.

At each level of the call stack, the provenance of the objects in scope is
maintained by associating program events and objects with their sources in the
flow graph. For local variables, a lookup table maps each variable name to an
output port of the box representing the function call that created or most
recently mutated the variable's value. The five callbacks made by the
transformed program are then handled by:
\begin{enumerate}[(i),nosep]
\item when a function call is about to begin, create a new box for the function
  call, as described above;
\item when a function call is returned from, pass the output ports of its box as
  sources to the enclosing expression;
\item when a variable is accessed, look up variable's source in the lookup table
  and pass the source to the enclosing expression;
\item when a variable is assigned, add or replace that entry in the lookup table
  with a new source;
\item when a variable is deleted, remove that entry from the
  lookup table.
\end{enumerate}
The passing of extra information in cases (ii) and (iii) can be performed
statically, by boxing the values that pass between two statically transformed
expressions, or dynamically, by pushing them onto a queue to be emptied later.

This concludes a technical overview of the program analysis systems for Python
and R. Language-specific features and challenges are discussed in subsequent
sections, but let us first consider the general tradeoffs made by our
methodology.

\paragraph{Static versus dynamic analysis}

Although it has a static phase, the program analysis system is essentially
dynamic, for several reasons. Static analysis, especially about type inference,
is challenging for dynamically typed languages like Python and R. Moreover, in
applications to data science, it is of interest to capture values computed
during the course of the program’s execution, such as parameter estimates or the
selected hyperparameters. For this dynamic analysis is indispensable. In
general, dynamic analysis is easier to implement that static analysis, as it
skirts the computationally intractable or even undecidable problems that static
analysis easily produces \cite{landi1992}.

Of course, a disadvantage of dynamic analysis is the necessity of running the
program. Crucially, our system needs not just the code itself, but its input
data and runtime environment. These are all requirements of scientific
reproducibility, so in principle they ought to be satisfied. In practice they
are often neglected. Even assuming that all the resources are available, a
significant investment in infrastructure and curation would be needed to process
data analyses at a large scale. Thus, for metascientific purposes, it is
worthwhile to pursue purely static program analyses in addition to dynamic ones,
acknowledging the trade-offs involved in both cases.

\paragraph{User versus library code}

The program analysis system transforms and records only the code that is
directly presented to the system as input. For now, the input program is assumed
to be a single script or Jupyter notebook, as is common in data
analysis.\footnote{That said, complex data analyses may easily span multiple
  scripts or notebooks connected by intermediate files and databases. A more
  complete system would track the provenance of data both within a single
  program, as our system does, and between programs, as done by conventional
  data provenance tools \cite{simmhan2005a}.} Any code that is executed
indirectly through library or system functions is thus not recorded. This
restriction accords with the basic assumption, first stated in
\cref{ch:introduction}, that library functions are semantically meaningful
computational units.

Limiting the scope of the program analysis also reduces its runtime overhead.
For a typical data analysis script, the vast majority of its running time is
spent inside imported Python or R functions or inside the underlying C and
Fortran routine. These computations are not affected by the program analysis
system.

\paragraph{Programs versus program executions}

The presentation so far has not carefully distinguished between a computer
program and a particular \emph{execution} of that program. In general, a program
having nontrivial control flow, through conditional branching (\verb|if|
statements), looping (\verb|for| and \verb|while| statements), or recursion,
will execute a different sequence of instructions depending on the input to the
program. At present, our system ignores control flow, capturing only the
function calls made during the execution of the program, not those that might
have been made on different input data. Thus, faced with a conditional branch,
the system records only the executed branch, and faced with a loop, the system
records the unrolled sequence of iterations.

This limitation, while important, is not as severe as might initially be
supposed. Most programs, such as a graphical application or a web server, are
expected to run on many different user inputs and behave differently each time.
In contrast, data analyses are usually created for and attached to specific
datasets. As a requirement of scientific reproducibility, we expect that if the
program is run repeatedly on the same data, it will produce the same result each
time. It is therefore reasonable to ask what actually did happen on given data,
rather than what might have happened on different data. Even so, unrolling a
loop with many iterations causes a blow-up in the size of the flow graph, which
is highly impractical. Addressing this problem is left to future work.

\section{Flow graphs for Python}
\label{sec:program-analysis-python}

A clean and simple syntax belies the surprising complexity of the Python
programming language. Advertised as a multi-paradigm language, it supports
imperative, functional, and object-oriented programming styles and includes
numerous advanced features, such as generators, list comprehensions, and, in
recent versions, coroutines. Its combination of features has made Python into a
popular general-purpose programming language but causes difficulties when
translating Python code into a simpler model of computation, like that of flow
graphs. Our program analysis system resolves some but not all of these
difficulties.

\begin{listing}
  \begin{tcolorbox}
    \inputminted[fontsize=\footnotesize]{python}{code/plot_svm_regression.py}
  \end{tcolorbox}
  \caption{Support vector regression (SVR) in Scikit-learn. The code is
    reproduced without changes from Scikit-learn's official suite of examples.}
  \label{lst:plot-svm-regression}
\end{listing}

\begin{figure}
  \centering
  \resizebox{\textwidth}{!}{
    \input{wiring-diagrams/program-analysis/plot_svm_regression.raw}}
  \caption{Raw flow graph for \cref{lst:plot-svm-regression}}
  \label{fig:plot-svm-regression-raw}
\end{figure}

\paragraph{Homogenizing the syntax}

Many instructions in Python, including access and assignment of object
attributes, indexed access and assignment of containers, and special unary and
binary operators, can be interpreted as function calls yet are not officially
regarded as such by the Python grammar. During the static phase of program
analysis, such instructions are reduced to function calls so that they may be
recorded by the same methods. Conveniently for this purpose, the standard
library module \verb|operator| provides functional aliases for nearly every
special operator in the Python language.

For example, the script in \cref{lst:plot-svm-regression}, reproduced verbatim
from the Scikit-learn documentation, uses the multiplication (\verb|*|),
subtraction (\verb|-|), array indexing and slicing (\verb|[::]|), and in-place
addition (\verb|+=|) operators. In the corresponding flow graph of
\cref{fig:plot-svm-regression-raw}, these instructions appear as calls to the
\verb|mul|, \verb|sub|, \verb|getitem|, and \verb|iadd| functions from the
\verb|operator| module, as well as the built-in function \verb|slice|.

\paragraph{Object-oriented programming}

Although the foundational NumPy and SciPy packages have mainly procedural
interfaces, many Python packages for data science, such as Pandas and
Scikit-learn, are written in an object-oriented style. Some packages, such as
Matplotlib and Statsmodels, provide both procedural and object-oriented
interfaces. In comparison, flow graphs are typed, even allowing a form of
subtyping through implicit conversion (\cref{sec:concepts-as-category}), but
they are not object-oriented in the sense of class-based programming, where
objects inherit attributes and methods through classes.

To accommodate class-based objects, the program analysis system reinterprets
calls of instance methods as calls of functions having an extra first argument
for the object instance. This protocol actually agrees with Python's own syntax
for defining an instance method, where the instance object is represented by a
first argument conventionally called \verb|self|. For example, the \verb|fit|
method of a supervised model in Scikit-learn, seen in
\cref{lst:plot-svm-regression}, has the signature \verb|fit(self, X, y)|, which
appears as a box with three input ports in \cref{fig:plot-svm-regression-raw}.
Object attribute accesses and assignments are translated into \verb|getattr| and
\verb|setattr| calls, using the method described previously, but for readability
are displayed as boxes labeled by the name of the attribute.

\paragraph{Side effects}

The most severe difficulties in analyzing Python code arise from the mismatch
between the programming models of Python, which permits mutation, global state,
and other side effects, and of flow graphs, which is purely functional. To these
problems only partial solutions are available. When it is known that a function
or method mutates its inputs, it may be reinterpreted as a non-mutating function
with an additional output for every mutated input. For example, the \verb|fit|
methods in Scikit-learn return a fitted model, which is in fact the original,
unfitted model having been mutated. Regardless of whether the return value is
consumed (compare \cref{lst:toy-regression-sklearn,lst:plot-svm-regression}),
function calls to \verb|fit| become boxes with an output for the fitted model.

That is simple enough; the trouble is knowing \emph{when} a function is
mutating. Functions corresponding to certain special operations, such as
\verb|setattr| and \verb|setitem| for attribute and indexed assignment, can
always be assumed to be mutating. When they are available, function annotations
can also be used to manually mark a function as mutating (\cref{sec:dso}). But
in general the system cannot determine whether an arbitrary function is
mutating, since the mutation can occur in library code that is not statically
analyzed. Moreover, mutations can be \emph{implicit}. For example, if a column
in a data frame is mutated through a reference to that column, then the
containing data frame should be regarded as mutated as well. Implicit mutations
are currently not addressed.

Mutation of global state is a still worse problem for the functional paradigm.
As an example, \cref{lst:plot-svm-regression} generates a plot using
Matplotlib's imperative interface, which maintains global state for the active
plot. The calls made to \verb|xlabel|, \verb|ylabel|, \verb|title|,
\verb|legend|, and \verb|show| all then appear as isolated boxes in the flow
graph of \cref{fig:plot-svm-regression-raw}. From the functional perspective,
this sequence of functional calls should be reinterpreted as a chain explicitly
passing a plot object from one call to the next. However, adopting this
convention in the flow graph would require a more drastic rewriting of the
observed program events than any considered here.

\paragraph{Weak references}

For objects with a unique identity, weak references offer a dynamic alternative
to static program transformation for tracking the data flow of objects. A
\emph{weak reference} is a reference, or pointer, to an object that does not
increase the object's reference count and hence does not prevent the object from
being garbage collected. The program analysis system may therefore freely use
weak references without causing memory leaks. Most objects in Python can be
weakly referenced. Important exceptions are the primitive objects, such as
booleans, integers, and strings, and the built-in containers for lists, tuples,
and dictionaries.

Wherever possible, the Python program analysis system maintains weak references
to the objects returned by recorded function calls, allowing the objects to be
reidentified if they appear as arguments to future function calls. This method
is complementary to the default, static approach.\footnote{The original
  implementation of the Python program analysis, described in
  \cite{patterson2017ibmjrd}, relied exclusively on weak references; however,
  this is plainly inadequate.} It is not applicable to all objects, yet where it
is applicable, it provides a strong guarantee of object identity, immune to any
dynamic trickery that might fool a static analysis.

\section{Flow graphs for R}
\label{sec:program-analysis-r}

The R programming language is, in most respects, far simpler than Python.
Despite appearances, R is quite close to being a functional language. Nearly all
of its standard functions and data structures have \emph{copy-on-modify}
semantics, meaning that ostensibly mutating operations in fact copy the data
before modifying it. The only mutable data structure in base R is the
environment type, which does not figure explicitly in most data analysis
scripts. The default and most commonly used object system, known as S3, is based
on generic functions, not classes, making it a natural fit for flow graphs. The
R language also has a very simple abstract syntax, with expressions, or abstract
syntax trees, composed of only four types of nodes:\footnote{In R, all control
  flow constructs, variable assignment operators, and unary and binary operators
  are special types of functions.} scalar constants, names, function calls, and
``pairlists.'' The first three node types are accommodated by the general
methodology of \cref{sec:program-analysis-general}, while the last is a special
type for function arguments that does not appear in ordinary code.

Before discussing a few distinctive aspects of the R language, let us examine a
data analysis written in R that is more realistic than the preceding examples. A
recent DREAM Challenge \cite{saez-rodriguez2016} asked how well clinical and
genetic covariates predict patient response to anti-TNF treatment for rheumatoid
arthritis \cite{sieberts2016}. Each team of analysts was instructed to submit
two models, one using only genetic covariates and the other using any
combination of clinical and genetic covariates.

An analysis submitted by a top-ranking team \cite{kramer2014} is displayed in
\cref{lst:dream-ra}, having been lightly modified for portability. Its flow
graph is shown in \cref{fig:dream-ra-raw}. The analysts fit two predictive
models, the first using both genetic and clinical covariates and the second
using only clinical covariates. Both models use the Cubist regression algorithm
\cite[\S 8.7]{kuhn2013}, a variant of random forests based on M5 regression
model trees \cite{wang1997}. Because the genetic data is high-dimensional, the
first model is fit using the subset of the genetic covariates selected by VIF
regression, a variable selection algorithm \cite{lin2011}. The linear regression
model created by VIF regression is used only for variable selection, not for
prediction. If this code is found to be opaque, the semantic flow graph
displayed in \cref{fig:dream-ra-semantic} of the following chapter can also be
consulted.

\begin{listing}
  \begin{tcolorbox}
    \inputminted[fontsize=\footnotesize]{r}{code/dream-ra.R}
  \end{tcolorbox}
  \caption{R code for two models from the Rheumatoid Arthritis DREAM Challenge
    \cite{sieberts2016}. The code is reproduced without essential changes from
    \cite{kramer2014}.}
  \label{lst:dream-ra}
\end{listing}

\begin{figure}
  \centering
  \resizebox{\textwidth}{!}{
    \input{wiring-diagrams/program-analysis/dream-ra.raw}}
  \caption{Raw flow graph for \cref{lst:dream-ra}}
  \label{fig:dream-ra-raw}
\end{figure}

Although R conforms fairly neatly to the formalism of flow graphs, the language
has two distinctive features, lazy evaluation and non-standard evaluation, that
must be accounted for when transforming programs.

\paragraph{Laziness} The R language has \emph{call-by-need} semantics, under
which arguments to function calls are evaluated lazily. That is, when a function
is called, it does not receive previously computed values as arguments but
rather \emph{promises} to compute these values. A promise object is evaluated
only when it is accessed for the first time within the body of the function. In
particular, if a promise object is never accessed, then it is never evaluated.
Also, the order in which a function's arguments are evaluated depends on the
body of the function.

Since R functions can have side effects, a valid program transformation
generally cannot force promises to be evaluated earlier than they would
otherwise be. The R program analysis system meets this requirement through
non-standard evaluation, rewriting the code in the promise before the function
is called, without evaluating the promise. The evaluation order is thus exactly
the same as in the original program, which is not necessarily that given for
Python programs in \cref{sec:program-analysis-general}.

\paragraph{Non-standard evaluation}

Although R may be the world's most widely used lazy functional programming
language, most end-users of R do not explicitly rely on laziness in their code
and may not even be aware that R is a lazy language \cite{goel2019}. One
consequence of laziness in R is to make possible a dynamic form of
metaprogramming. In \emph{non-standard evaluation}, the interpretation of a
function argument having a fixed syntactical form can depend on the internal
logic of the function. Thus, in general, the semantics of R code cannot be known
statically. In practice, non-standard evaluation is often only a convenience,
used to abbreviate code and capture variable names for plot labels. Such uses
have no significant impact on the program semantics. However, the tidyverse
suite of packages makes much more extensive use of non-standard evaluation,
effectively implementing its own domain-specific language for data processing
within R \cite{wickham2019}.

Non-standard evaluation is the most serious obstacle to the program analysis of
R code because the meaning of an expression can depend arbitrarily on the
context in which it is evaluated. There appears to be no simple and generic
solution to this problem. At present, our program analysis system handles, on a
case-by-case basis, certain recurring patterns of non-standard evaluation in
base R. However, no attempt has been made to systematically accommodate all the
major uses of non-standard evaluation in the R ecosystem.

\section{Notes and references}

An early version of the Python program analysis system is described in previous
work by the author and collaborators \cite{patterson2017ibmjrd}. Significant
design and implementation changes have since been made to improve the robustness
of the system. Most importantly, the use of Python's system trace function
(\verb|sys.settrace|), which cannot trace built-in functions or C extensions,
has been eliminated in favor of more difficult to implement, but more reliable,
static program transformations
(\cref{sec:program-analysis-general,sec:program-analysis-python}). Preliminary
support for the R programming language has also been added
(\cref{sec:program-analysis-r}). The R program analysis system was announced,
but not described in any detail, in a more recent publication
\cite{patterson2018}.

\paragraph{Program analysis}

A general reference on computer program analysis, with a more theoretical
orientation than this chapter, is the text by Nielson et al \cite{nielson1999}.
As the book's contents suggest, a large part of the literature on program
analysis is about static analysis, with the main application being to optimizing
compilers \cite{aho2006}. A recent survey of dynamic program analysis is
\cite{gosain2015}. Most applications of dynamic analysis are to debugging and
testing, performance profiling, and security analysis.

Data flow analyses are also classified as \emph{intraprocedural} (within a
single procedure) or \emph{interprocedural} (between procedures in a larger
program) \cites[Chapter 2]{nielson1999}{khedker2009}. The program analysis of
this chapter is therefore a dynamic, interprocedural data flow analysis. A
common intermediate representation in an interprocedural data flow analysis is a
\emph{call graph} (or \emph{call multigraph}), where the vertices correspond to
functions and there is a directed edge from vertex $u$ to $v$ for each call of
the function $v$ from within $u$. Call graphs do not contain the information
about data dependency that a flow graph does. Also, a flow graph is not,
strictly speaking, a graph but a wiring diagram.

By building a custom version of the R interpreter, the authors of
\cite{morandat2012} implemented the \emph{TrackeR} system to capture detailed
execution traces of R programs. The dynamic program analysis here is
comparatively light-weight, being based on the standard R interpreter but
running transformed R programs.

\paragraph{Data provenance} The aim of recording the steps of a data analysis is
shared by the field of data provenance. The \emph{provenance} of a data resource
includes its origin and the process of transformation by which it was derived
\cite{simmhan2005a}. Insofar as the survey \cite{simmhan2005b} is
representative, the main difference between our system and the typical data
provenance system seems to be granularity. In data provenance, the finest
granularity of data resource is often files or database records, whereas our
system operates on a single file and traces arbitrary program objects. For
example, the StarFlow system targets data analyses written in Python but
operates at the level of scripts \cite{angelino2010}. The program analysis here
would be naturally complemented by a file-level provenance system to describe
data provenance at multiple levels of granularity.


TITLE: Numerical Methods - Fixed Point Iteration
QUESTION [1 upvotes]: We were set the following question before we had learned about systems of non-linear equations and I'm confused as to how we're meant to answer it using the fixed point method theorems that we have been provided with.
Consider the fixed point iteration $\underline{x}^{n+1} = \underline{\Psi}(\underline{x}^{n})$ where
$\Psi(\underline{x}) = \Psi((x_1,x_2,x_3)^T))$ =

Analyse the convergence of the method to the fixed point $\alpha$ = $(0,\frac{1}{3},0)^T$.
Are we meant to show that $|\Psi(\underline{x})|<1$ and continously differentiable? And if so, do we do this using a Jacobian matrix?

REPLY [1 votes]: When $x$ is near $\alpha$, the map $\Psi$ can be approximated by the first-order Taylor series
$$ \Psi(x) \approx \Psi(\alpha) + J(\alpha) (x - \alpha) = \alpha + J(\alpha) (x - \alpha) $$
where $J(\alpha)$ is the Jacobian matrix evaluated at $x=\alpha$.  Note that we have used the fact that $\alpha$ is a fixed point, i.e., $\Psi(\alpha)=\alpha$.
From this you can show that that the recursion
$$ x^{n+1} = \Psi(x^n) $$
is equivalent to
$$ x^{n+1} - \alpha = J(\alpha) (x^n - \alpha) $$
So if you define $y^n = x^n - \alpha$ then
$$ y^{n+1} = J(\alpha) y^n $$
so
$$ y^n = J(\alpha)^n y^0 $$
where the superscript on $J(\alpha)$ is an exponent rather than a time index.  Now the critical question is whether and how quickly $J(\alpha)^n$ shrinks with $n$.  This depends on the eigenvalues of $J(\alpha)$.
Basically, the fixed point is stable if all of the eigenvalues of $J(\alpha)$ are less than 1.  The convergence error is roughly $\lambda_1^n$, where $\lambda_1$ is the largest eigenvalue.


TITLE: If $0 < \theta < \frac{\pi}{2}$, then $\gamma$ is a logarithmic spiral
QUESTION [1 upvotes]: Let $\gamma$ be a plane curve parametrized by the arc length, having the property that its tangent vector $T(t)$ forms a fixed angle $\theta$ with $\gamma(t)$.
Before explaining where I am, let define the notations I use:

$T(t) := \dot{\gamma}(t)$;
$N(T)$ is the vector that is normal to $T(t)$;
$\kappa(t)$ is the curvature of $\gamma(t)$;
Frenet differential equation is $\ddot{\gamma}(t) = \kappa(t) N(t)$ and $N(t) = i \dot{\gamma}(t)$.

We have that $\gamma$ is parallel to $\dot{\gamma} \exp(i \theta)$, so we can write $\gamma(t) = \alpha(t) \dot{\gamma}(t) \exp(i \theta)$ where $\alpha \in C^\infty$. Taking the derivative in both sides we obtain:
$$\begin{align*}
\dot{\gamma}(t) &= \dot{\alpha}(t) \dot{\gamma}(t) \exp(i \theta) + \alpha(t) \ddot{\gamma}(t) \exp(i \theta)\\
&= (\dot{\alpha}(t) \dot{\gamma}(t) + \alpha(t) \ddot{\gamma}(t)) \exp(i \theta)\\
&= (\dot{\alpha}(t) \dot{\gamma}(t) + \alpha(t) i \kappa(t) \dot{\gamma}(t)) \exp(i \theta)\\
\Rightarrow 1 &= \exp(i \theta)(\dot{\alpha}(t) + \alpha(t) i \kappa(t))\\
\Rightarrow \exp(-i \theta) &= \dot{\alpha}(t) + \alpha(t) i \kappa(t).
\end{align*}$$
By identifying real and imaginary parts, we obtain:
$$\begin{cases}
\cos(\theta) &= \dot{\alpha}(t)\\
\sin(\theta) &= -\alpha(t) \kappa(t).
\end{cases}$$
Integrating on both side of the first equation we obtain
$$\begin{cases}
\cos(\theta)t &= \alpha(t)\\
\sin(\theta) &= -\alpha(t) \kappa(t).
\end{cases}$$
Therefore
$$\sin(\theta) = -\cos(\theta)t \kappa(t) \quad \Rightarrow \quad \kappa(t) = -\frac{1}{t} \tan(\theta).$$
Now my hint says that I should obtain $\kappa(t) = -\frac{1}{t} \cot(\theta)$, which I obviously don't although I'm not far. Then I'm a bit lost on what to do next. Could someone help on these two points?

REPLY [0 votes]: Your result is correct. Note that $\alpha(t)$ denotes the distance of the moving point from the origin. When $0<\theta\ll1$ then the tangent vector $T(t)$ is almost parallel to $\gamma(t)$, which means that almost all of the arc length goes into increasing $|\gamma(t)|=\alpha(t)$. Therefore $\dot\alpha(t)=\cos\theta\doteq1$ makes sense, etcetera.


TITLE: Question about multidimensional iid random variable
QUESTION [3 upvotes]: Let $X=(x_1,\ldots,x_d)^\top\in[0,1]^d$ be the row-wise representation of an $n\times n$ image ($d=n\times n$). Each element of $X$ is the value of a pixel, which we assume it belongs to $[0,1]$.
If we assume that $X$ is distributed normally with covariance matrix $\Sigma=\sigma^2I_d$, i.e., $X\sim\mathcal{N}(\mathbf{0},\Sigma)$, what can we say about the various sub-regions of $X$?
Starting from a single variable, say $x_i$ for some $i\in\{1,\ldots,d\}$, is it true that $x_i\sim\mathcal{N}(0,\sigma^2)$? If so, is it also true that a random selection of $x_i$'s, for instance the vector $Y=(x_3, x_7, x_{11})^\top$, is also distributed normally with the same variance? That is, in the case of the $3$-dimensional $Y$ for instance, is it true that $Y\sim\mathcal{N}(\mathbf{0},\sigma^2I_3)$?.
Finally, if we assume that the original $X$ is independent and identically distributed (iid), would that be a sufficient condition?

REPLY [1 votes]: I consider the two variable case and  take into account that the means of $X_1$ and $X_2$ are 0.
The pdf of the univariate untruncated distribution of $X_1$ is
$f_{X_1}(x_1)=\frac{1}{\sqrt{2\cdot \pi \sigma_1^2}}\cdot e^{-\frac12 \left[\frac{x_1^2}{\sigma_1 ^2}\right]} $
The cdf of the univariate untruncated distribution of $X_1$ is
$F_{X_1}(x_1)=\int_{-\infty}^{x_1}\frac{1}{\sqrt{2\cdot \pi \sigma_1^2}}\cdot e^{-\frac12 \left[\frac{s^2}{\sigma_1 ^2}\right]} \, ds$
The pdf of the univariate truncated distribution of $X_1$  on $[0,1 ] $is
$f(x_1;0,1)=\frac{\frac1{\sigma_1} \cdot f_{X_1}(x_1)}{F_{X_1}(1)-F_{X_1}(0)}
$
The cdf of the univariate truncated distribution of $X_1$  on $[0,1 ] $is
$F(x_1;0,1)=\int_{0 }^{x_1}\frac{\frac1{\sigma_1} \cdot f_{X_1}(s)}{F_{X_1}(1)-F_{X_1}(0)} \, ds$
Since $X_1$ and $X_2$ are independent the joint truncated cumulative distribution function on $[0,1 ]^2 $ is
$F(x_1,x_2;[0,1]^2)=\int_{0 }^{x_1}\frac{\frac1{\sigma_1} \cdot f_{X_1}(s)}{F_{X_1}(1)-F_{X_1}(0)} \, ds\cdot \int_{0 }^{x_2}\frac{\frac1{\sigma_2} \cdot f_{X_2}(t)}{F_{X_2}(1)-F_{X_2}(0)} \, dt$