text:
In order to achieve the established goal, we make use of indirect measurements.
Specifically, instead of \g, we measure an auxiliary parameter \h, which we
refer to as the quantity of measurement. The observations of \h are then
processed via Bayesian inference (see \xref{bayesian-statistics}) in order to
derive the on-wafer distribution of \g. The quantity \h is chosen such that
\one~\h is convenient and cheap to be measured; \two~\h depends on \g, which is
signified by $\h = \f(\g)$; and \three~there is a way to compute \h given \g.
The last requirement means that \f should be known. However, \f does not have to
be specified analytically, since our framework treats \f as a ``black box.'' For
example, \f can be a piece of code.

Without loss of generality, we adhere to the following convention. Each die is a
potential measurement site, and $\hnd < \nd$ denotes the number of those sites
that are actually measured. Each site comprises \np measurement points, and
there are \ns data instances per point. For instance, in the example given in
\sref{bayes-motivation}, an observation at a site is a temperature profile \mq,
which is a matrix capturing the temperatures of \np processing elements at \ns
moments in time as defined in \eref{temperature-profile}. Denote the collected
data by $H$ and assume that the locations of the measurement points are recorded
along with $H$.

It is worth noting that, if \f is the identity function (that is, $\h = \f(\g) =
\g$), the proposed technique focuses primarily on the reconstruction of any
missing observations in $H$, that is, on the unobserved sites on the wafer. From
this standpoint, our approach is a generalization of those developed in
\cite{reda2009, zhang2010}.

\inputfigure{bayes-overview}
In the rest of this section, we present our framework for characterizing process
variation. The technique revolves around Bayes's theorem \cite{gelman2013}
\begin{equation} \elab{bayes-theorem}
  p(\u | H) \propto p(H | \u) p(\u).
\end{equation}
In our context, the theorem relates the posterior density function $p(\u |H)$ of
\u given $H$ with the likelihood function $p(H | \u)$ of $H$ given \u and the
prior density function $p(\u)$ of \u; these concepts are introduced in
\xref{bayesian-statistics}, and they will be discussed further in the following.
The framework is divided into four major stages depicted in
\fref{bayes-overview}. Stage~1 is the data-harvesting stage where the designer
collects a set of observations of the quantity of measurement \h, thereby
forming the data set $H$. At Stage~2, we perform an optimization procedure that
is designed to assist the sampling procedure at Stage~3. The latter produces a
collection of samples of the quantity of interest \g, such as the effective
channel length, denoted by $G$. This data set $G$ is then processed at Stage~4
in order to estimate the desired statistics about \g, such as the probability of
\g being smaller than a certain threshold; recall the example in
\sref{bayes-motivation}.

It can be seen in \fref{bayes-overview} that Stage~2 and Stage~3 actively
communicate with the two models shown on the right-hand side, which are called
the data model and statistical model. We begin by elaborating on these models.

\subsection{Data Model}
\slab{bayes-data-model}

The data model relates the quantity of interest \g with the quantity of
measurement \h, which is denoted as follows:
\[
  \h = \f(\g).
\]
The function \f depends on the choice of \h and is specified by the designer.
The data model is utilized in order to predict the values of \h at the same
measurement sites, at the same measurement points, and with the same number and
meaning of data instances as the ones in $H$ obtained at Stage~1. The resulting
data are stacked into a single vector denoted by $\vh \in \real^{\hnd \np \ns}$.
Let also $\hvh \in \real^{\hnd \np \ns}$ be a stacked version of the data in $H$
such that the respective elements of \vh and \hvh correspond to the same
locations on the wafer.

In order to acquire a better understanding of the data model, let us return to
the example given in \sref{bayes-motivation}. In this case, \g stands for the
effective channel length, and \h for the temperature profile \mq corresponding
to a fixed workload. The data model can be roughly divided into two transitions:
\one~from the effective channel length \g to the static power consumed by the
platform at hand and \two~from this static power to the corresponding
temperature profile \h. At this point, it is worth recalling the power model
presented in \sref{power-model}. The first transition is due to the dependence
of the leakage current on the effective channel length, which is implicitly
present in \eref{static-power}. The transition can then be made by means of the
model shown in \eref{static-power} or any of its variations; see, for instance,
\cite{chandrakasan2000, srivastava2010, yang2010, juan2012}. In particular, an
adequate model of static power can be constructed via fitting to \up{SPICE}
simulations of reference electrical circuits. The only requirement for such a
model is that it should be parameterized by \g. In addition, it can be
parameterized by temperature in order to take account of the power-temperature
interdependence described in \sref{power-model}. The second transition is made
by combining the static power component with the dynamic power that corresponds
to the workload being considered. The resulting total power and the relevant
temperature-related information---such as the floorplan and thermal parameters
of the platform---are fed to a temperature simulator in order to acquire the
corresponding temperature \h, which is discussed in detail in
\cref{certainty-development}.

\subsection{Statistical Model}
\slab{bayes-statistical-model}

\inputfigure{bayes-statistical-model}
Once the wafer has been fabricated, the values of \g across the wafer are fixed;
however, they remain unknown to the designer. In order to infer them, we employ
the model developed in the current subsection, which can also be seen in
\fref{bayes-statistical-model}. The development consists of the five steps
described below.

Step~1 is to assign an adequate model to the unknown \g. We model \g as a
Gaussian process \cite{rasmussen2006}, since \one~it is flexible in capturing
the variation patterns induced by the fabrication process, \two~it is
computationally efficient, and \three~Gaussian distributions are often adequate
models of uncertainty due to process variation \cite{reda2009, srivastava2010,
juan2012}. The model is denoted by
\begin{equation} \elab{bayes-prior-quantity}
  \g | \u_\g \sim \mathrm{Gaussian\ Process}(\mu, v)
\end{equation}
where $\mu: \real^2 \to \real$ and $v: \real^2 \times \real^2 \to \real$ are the
mean function and covariance function of \g, respectively, which take locations
on the wafer as arguments. The notation also indicates that \g depends on a set
of parameters $\u_\g$, which we shall identify later on. Prior to taking any
measurements, \g is assumed to be spatially unbiased. Therefore, we let $\mu$ be
a single location-independent parameter $\mu_\g$, which means that $\mu(\v{r}) =
\mu_\g$ for any $\v{r} \in \real^2$. The covariance function $v$ is chosen to be
the following:
\begin{equation} \elab{bayes-covariance}
  v(\v{r}_1, \v{r}_2) = \sigma_\g^2 \, k(\v{r}_1, \v{r}_2)
\end{equation}
for $\v{r}_1 \in \real^2$ and $\v{r}_2 \in \real^2$ where
\begin{equation} \elab{bayes-correlation}
  k(\v{r}_1, \v{r}_2) = w \, k_\SE(\v{r}_1, \v{r}_2) + (1 - w) k_\OU(\v{r}_1, \v{r}_2)
\end{equation}
is the correlation function, and
\begin{align*}
  k_\SE(\v{r}_1, \v{r}_2) & = \exp\left(-\frac{\norm[2]{\v{r}_1 - \v{r}_2}^2}{\ell_\SE^2}\right) \text{ and} \\
  k_\OU(\v{r}_1, \v{r}_2) & = \exp\left(-\frac{\absolute{\norm[2]{\v{r}_1} - \norm[2]{\v{r}_2}}}{\ell_\OU}\right)
\end{align*}
are the squared-exponential and Ornstein--Uhlenbeck kernels
\cite{rasmussen2006}, respectively. In these formulae, $\sigma_\g^2$ is the
variance of \g, $w \in [0, 1]$ is a weight coefficient, $\ell_\SE > 0$ and
$\ell_\OU > 0$ are so-called length-scale parameters, and $\norm[2]{\cdot}$
stands for the Euclidean distance. The choice of $v$ is based on the
observations of the correlation structures induced by the fabrication process
\cite{chandrakasan2000, cheng2011}. Specifically, $k_\SE$ imposes similarities
on locations that are close to each other on the wafer, and $k_\OU$ imposes
similarities on locations that are at the same distance from the center of the
wafer. The parameters $\ell_\SE$ and $\ell_\OU$ control the extent of these
similarities, that is, the range where the correlation between two locations is
significant. Although all the above parameters of \g can be inferred from data,
for simplicity, we focus only on $\mu_\g$ and $\sigma_\g^2$. The rest of the
parameters---namely $w$, $\ell_\SE$, and $\ell_\OU$---are assumed to be
determined prior to our analysis based on knowledge of the variation patterns
that are typical for the fabrication process being considered; see
\cite{marzouk2009} and the references therein.

Step~2 is to make the above model of \g computationally tractable. The model is
an infinite-dimensional object, as it characterizes a continuum of locations.
For practical computations, it should be reduced to a finite-dimensional one.
First, \g is discretized with respect to the union of two sets of locations. The
first one is composed of the $\hnd \np$ points where the observations in $H$ are
made (\hnd measurement sites with \np measurement points each), and the other
one of the locations where the designer wishes to characterize \g. For
simplicity, assume that the designer is interested in all the sites, which
results in $\nd \np$ locations in total. Let $\vg \in \real^{\nd \np}$ store the
values of \g at these locations. Second, the dimensionality of the problem is
reduced further via the \ac{KL} decomposition, which is introduced in
\xref{probability-transformation}. Concretely, we perform the transformation
shown in \eref{model-order-reduction} with respect to the correlation matrix of
\vg computed via \eref{bayes-correlation}. The result is
\begin{equation} \elab{bayes-reduction}
  \vg = \mu_\g \v{1} + \sigma_\g \m{U} \tm{\Lambda}^\frac{1}{2} \vz
\end{equation}
where $\v{1} = (\range{1}{1}) \in \real^{\nd \np}$, and $\vz = (\z_i) \in
\real^\nz$ is a vector of independent random variables that obey the standard
Gaussian distribution. The number \nz is the final dimensionality of the model
of \g; typically, $\nz \ll \nd \np$. In addition, the parameters $\u_\g$ in
\eref{bayes-prior-quantity} are now known: $\u_\g = \set{\vz, \mu_\g,
\sigma^2_\g}$; see \fref{bayes-statistical-model}. The model is now ready for
practical computations.

Step~3 is to specify the likelihood function in \eref{bayes-theorem}, which is
where the observed information is taken into account; see
\xref{bayesian-statistics}. In our case, the observed information is the
measurements $H$ stacked into \hvh as described in \sref{bayes-data-model}.
Since the measurement process is not perfect, we also have to take measurement
noise into consideration. To this end, the observed \hvh is assumed to deviate
from the data model's prediction \vh as follows:
\[
  \hvh = \vh + \v{\epsilon}
\]
where $\v{\epsilon}$ is an $\hnd \np \ns$-dimensional vector of noise, which is
typically assumed to be white Gaussian noise \cite{rasmussen2006, marzouk2009}.
Without loss of generality, the noise is assumed to be independent of \g and to
have the same magnitude for all measurements. Hence, the model of the noise is
as follows:
\[
  \v{\epsilon} | \sigma^2_\epsilon \sim \mathrm{Gaussian}(\v{0}, \sigma^2_\epsilon \m{I})
\]
where $\sigma^2_\epsilon$ is the variance of the noise. At this point, all the
parameters of the inference are identified, and they are $\u = \u_\g \cup
\set{\sigma_\epsilon^2} = \set{\vz, \mu_\g, \sigma_\g^2, \sigma_\epsilon^2}$;
see \fref{bayes-statistical-model}. Taking the above into account, we obtain
\begin{equation} \elab{bayes-likelihood}
  \hvh | \u \sim \mathrm{Gaussian}(\vh, \sigma_\epsilon^2 \m{I}).
\end{equation}
The density function of this distribution is the likelihood $p(H | \u)$ of our
statistical model, which is the first element of the posterior in
\eref{bayes-theorem}.

Step~4 is to decide on the second element of the posterior in
\eref{bayes-theorem}, that is, on the prior $p(\u)$. We put the following priors
on the parameters in \u:
\begin{align} \elab{bayes-prior}
  \begin{split}
    \vz               & \sim \mathrm{Gaussian}(\v{0}, \m{I}), \\
    \mu_\g            & \sim \mathrm{Gaussian}(\mu_0, \sigma^2_0), \\
    \sigma^2_\g       & \sim \mathrm{Scaled\ Inverse\ } \chi^2(\nu_\g, \tau^2_\g), \text{ and} \\
    \sigma^2_\epsilon & \sim \mathrm{Scaled\ Inverse\ } \chi^2(\nu_\epsilon, \tau^2_\epsilon).
  \end{split}
\end{align}
The prior of \vz is due to the decomposition in \eref{bayes-reduction}. The
other three priors, that is, a Gaussian and two scaled inverse chi-squared
distributions, are a common choice for a Gaussian model whose mean and variance
are unknown. The parameters $\mu_0$, $\tau^2_\g$, and $\tau^2_\epsilon$
represent the presumed values of $\mu_u$, $\sigma^2_\g$, and
$\sigma^2_\epsilon$, respectively, and are set by the designer based on prior
knowledge of the technological process and measurement equipment. The parameters
$\sigma_0$, $\nu_\g$, and $\nu_\epsilon$ reflect the precision of this
information. When the prior knowledge is weak, less specific priors can be
considered \cite{gelman2013}. Finally, $p(\u)$ in \eref{bayes-theorem} is
obtained by multiplying the densities of the priors in \eref{bayes-prior}.

Step~5 is to calculate the posterior $p(\u | H)$ in \eref{bayes-theorem}. To
this end, the likelihood function $p(H | \u)$, which is the density of the
distribution shown in \eref{bayes-likelihood}, and the prior $p(\u)$, which is
the product of the densities of the distributions shown in \eref{bayes-prior},
are put together. The density of the resulting posterior distribution is as
follows:
\begin{equation} \elab{bayes-posterior}
  p(\u | H) \propto p(\hvh | \vz, \mu_\g, \sigma^2_\g, \sigma^2_\epsilon) p(\vz) p(\mu_\g) p(\sigma^2_\g) p(\sigma^2_\epsilon).
\end{equation}
Provided that there is a way to draw samples from \eref{bayes-posterior}, \g can
be readily analyzed, as we shall see in \sref{bayes-processing}. The problem,
however, is that direct sampling of the posterior is difficult due to the data
model involved in the likelihood function via \vh; see \eref{bayes-likelihood}.
In order to circumvent this problem, we utilize Markov chain Monte Carlo
sampling and, more specifically, the Metropolis--Hastings algorithm
\cite{gelman2013}, which is outlined in \xref{bayesian-statistics}. The
algorithm operates on an auxiliary distribution called the proposal
distribution. The construction of an adequate proposal is discussed next.

\subsection{Optimization Procedure}
\slab{bayes-optimization}

In this subsection, we describe the objective of Stage~2 in
\fref{bayes-overview}. Although the requirements for the proposal distribution
are mild, it is often difficult to pick an efficient proposal, that is, a
proposal that would yield a good approximation with as few evaluations of the
posterior---and thus of the data model in \sref{bayes-data-model}---as possible.
This choice is especially difficult in high-dimensional problems, and our
problem---involving around 30 parameters, as we shall see in
\sref{bayes-results}---is one them. Therefore, a careful construction of the
proposal is an essential component of the proposed framework.

A common technique for constructing a high-quality proposal is to optimize the
posterior given in \eref{bayes-posterior}. Specifically, we seek a value $\u^*$
of \u that maximizes \eref{bayes-posterior} and hence has the maximum posterior
probability. In addition, we calculate the negative of the Hessian matrix at
$\u^*$, which is called the observed information matrix and denoted by $\m{J}$;
see the output of Stage~2 in \fref{bayes-overview}. Using $\u^*$ and $\m{J}$, we
can construct a proposal that allows the Metropolis--Hastings algorithm \one~to
start producing samples directly from the desired regions of high probability
and \two~to explore those regions more rapidly. The usage of $\u^*$ and $\m{J}$
is explained next.

\subsection{Sampling Procedure}

Let us turn to Stage~3 in \fref{bayes-overview}. In order to construct an
adequate proposal and utilize it for sampling, we have at our disposal $\u^*$
and $\m{J}$ from Stage~2. A commonly used proposal is a multivariate Gaussian
distribution where the mean is the current location of the chain of samples
started at $\u^*$, and the covariance matrix is the inverse of $\m{J}$
\cite{gelman2013}. In order to speed up the sampling process, we would like to
make use of parallel computing. The aforementioned proposal, however, is purely
sequential, since the mean for the next sample draw is dependent on the previous
sample. Therefore, we appeal to a variant of the Metropolis--Hastings algorithm
known as the independence sampler \cite{gelman2013}. In this case, a typical
choice of the proposal is a multivariate t-distribution independent of the
current position of the chain as follows:
\begin{equation} \elab{bayes-proposal}
  \u \sim t_\nu(\u^*, \alpha^2 \m{J}^{-1})
\end{equation}
where $\u^*$ and $\m{J}$ are as in \sref{bayes-optimization}, \nu is the number
of degrees of freedom, and $\alpha$ is a tuning constant controlling the
standard deviation of the proposal. Now sampling the proposal in
\eref{bayes-proposal} and evaluating the posterior in \eref{bayes-posterior} can
be done in parallel. The obtained samples are then accepted or rejected as in
the usual Metropolis--Hastings algorithm.

Having completed the sampling procedure, we obtain a collection of samples of
the parameterization $\u = \set{\vz, \mu_\g, \sigma^2_\g}$. The first portion of
the samples is typically discarded as being unrepresentative; this portion is
known as the burn-in period. The preserved samples of \u are then passed through
\eref{bayes-reduction} in order to compute samples of \g, which are $\nd
\np$-dimensional. Denote the corresponding set of samples by $G$ and let its
cardinality be \no.

\subsection{Post-Processing}
\slab{bayes-processing}

At Stage~4 in \fref{bayes-overview}, using $G$, the designer computes the
desired statistics about the quantity of interest \g, such as the most probable
value of \g at some location on the wafer and the probability of a certain area
on the wafer being defective. These computations are performed in the same way
as it is typically done when \ac{MC} sampling is utilized. Specifically, they
reduce to the estimation of expected values with respect to the posterior
distribution of \u given in \eref{bayes-posterior}: in order to calculate a
certain quantity dependent on \g, one evaluates it for each sample in $G$ and
then takes the average value.

The strength of the Bayesian approach to inference starts to shine when one is
interested in assessing the trustworthiness of the measured data and thus the
credibility of estimates and decisions based on these data. Such an assessment
can be readily undertaken using our framework, since the delivered posterior
contains all the necessary information about the quantity of interest \g. As
discussed in \sref{bayes-motivation}, this is especially helpful in
decision-making.

ext: tex

max_stars_repo_name: IvanUkhov/thesis

max_issues_repo_name: IvanUkhov/thesis

max_forks_repo_name: IvanUkhov/thesis


text:
\lect{3}{30}
\marginnote{Guest lecture}
% \marginnote{Lecture X: Wednesday March 30, 2016. Guest lecture.}

\section*{Density version of VdW}
Consider $[0,N]\subset \Z$. We will color these numbers with $m$ colors. At least one color-class will contain more than $\frac{N}{m}$ elements by pigeonhole. Now, is the following density version of VdW theorem true?

\begin{theorem}[Density VdW]
\label{thm:density_VdW}
For every $k\in \N$ and for every $\delta>0$, there exists $N_0$ such that if $A\subset\{1,\dotsc,N\}$ and $N>N_0$, and $|A|\geq \delta N$, then $A$ contains a $k$-AP (arithmetic progression of length $k$).
\end{theorem}

In 1963, Roth proved this theorem for $k=3$ using the ``circle method'' of number theory. Later, Szemeredi proved the statement for arbitrary $k$. His proof, however, gives terrible bounds for $N_0$ and was extremely technical. Next, F\"urstenberg used ergodic theory (developing a whole theory) to prove the result in a clear way, but also has very soft bounds. In 1998, Tim Gowers developed a different proof using higher-order Fourier analysis (giving more reasonable bounds); this won him a fields medal.

Let $A\subset \{1,\dotsc,N\}$ be chosen randomly: every element is chosen independently with probability $\delta$\marginnote{For each $x\in \{1,\dotsc,N\}$ we include $x$ in $A$ with probability $\delta$ and don't include it with probability $1- \delta$ (say by flipping a biased coin).}.
Then $|A|\approx \delta N$. How many $3$-AP's are in $A$? $\delta^3 \times$ \# of $3$-AP's in $\{1,\dotsc,N\}$. This is then $c \delta^3 N^2$ for some $c>0$.

So random sets contain many $3$-APs. If $A$ does not have any 3-AP's, then ``it is far from looking random.'' Then we will try to find some structure in $A$ and exploit it. We will come back to this later.

\paragraph{Triangles in graphs.} How many triangles are in $G(n, \frac{1}{2})$? This is the \erdos-Ramsey (?) random graph, where each edge is decided independently by a fair coin flip. The probability of a triangle existing is $\left( \frac{1}{2} \right)^3 = \frac{1}{8}$, and there are ${n\choose 3}$ possible triangles, so we expect $\frac{1}{8}{n\choose 3}$ triangles.

The edge density is $\frac{1}{2}$ meaning the number of edges is $\frac{1}{2}{n\choose 2}$.

If we consider the adjacency matrix of $G(n,\frac{1}{2})$, besides symmetry, it just looks ``gray,'' meaning half filled, half not.

If we consider the complete bipartite graph on $n$ vertices, we have the same edge density, and the adjacency matrix looks like
\[
\begin{pmatrix}
0 & 1 \\
1 & 0
\end{pmatrix}
\]
as a block matrix. We can think of each block as uniform, either with density $0$ or $1$, but the whole matrix is certainly not uniform.

 For $G(n,\frac{1}{2})$, the whole matrix is uniform, with density $\frac{1}{2}$.

We could also consider a graph with adjacency matrix
\[
\begin{pmatrix}
\text{choose edge with prob. }\frac{2}{3} &\text{choose edge with prob. } \frac{1}{3}\\
\phantom{\text{choose edge with prob. }}\frac{1}{3} & \phantom{\text{choose edge with prob. }}\frac{2}{3}
\end{pmatrix}
\]
in some division into blocks. This graph would yield a different number of triangles.
\begin{definition}
Bipartite graph on $(A,B)$ is called $\epsilon$-regular if for all $A'\subset A$, for all $B'\subset B$, $|A'|\geq \epsilon |A|$, $|B'|\geq \epsilon B$ implies
\[
|d(A',B') - d(A,B)| \leq \epsilon,
\]
where $d(X,Y) = \frac{\edges (X,Y)}{|X|\cdot |Y|} $ is the edge density between $X$ and $Y$.\marginnote{$\edges(X,Y)$ denotes the number of edges between vertex sets $X$ and $Y$.}
\end{definition}
\begin{remark}
This means if we restrict to any large enough subsets, the edge density of those subsets is close to the edge density of the whole bipartite graph.
\end{remark}

\begin{lemma} \label{lem:eregular_yields_triangles}
If $A,B,C$ are vertex sets so that $d(A,B) = a$, $d(A,C)=b$, and $d(B,C)=c$, and $a,b,c\geq 2 \epsilon$, and $(A,B)$, $(A,C)$, $(B,C)$ are $\epsilon$-regular, then $(A,B,C)$ contains at least \[
(1-2 \epsilon)(a- \epsilon)(b- \epsilon)(c- \epsilon)\cdot |A|\cdot |B|\cdot|C|
\]
triangles.
\end{lemma}
\begin{proof}
First, $A$ has at most $\epsilon |A|$ vertices whose number of neighbors in $B$ is less than $(a- \epsilon)|B|$. Otherwise, take $B'=B$ and $A' = \{v\in A: |N_B(v)| < (a - \epsilon)|B|\}$. If $|A'|\geq \epsilon |A|$, then $d(A',B') \geq a - \epsilon$, contradicting the definition of $A'$.

In the same way, $A$ has at most $\epsilon|A|$ vertices whose neighbors in $C$ is less than $(b - \epsilon)|C|$.

Likewise, $B$ has at most $\epsilon |B|$ vertices whose neighbors in $C$ is less than $(c- \epsilon)|C|$.

\missingfigure{Oval for each of $A,B,C$ far away from each other. Cross out two patches of $A$ corresponding to vertices with few neighbors in $B$ and few in $C$. Connect $A$ to paches in $B$ and $C$. Patch in $B$ has size $(a - \epsilon) |B|$ and pach in $C$ has size at least $(b- \epsilon)|C|$.  }
\end{proof}

\begin{theorem}[Regularity lemma] \label{thm:regularity_lemma}
For every $\epsilon>0$ there exists $T = T(\epsilon)$ such that the vertices of every graph $G$ can be partitioned into almost equal size parts $V_1,\dotsc, V_T$ such that at least $(1-\epsilon)$ fraction of pairs $(V_i,V_j)$ are $\epsilon$-regular.
\end{theorem}

\begin{remark}
$T$ does not depend on the number of vertices, only on $\epsilon$. So instead of considering an arbitrarily large graph, we may simply think of a $T\times T$ matrix of pairs $(V_i,V_j)$ such that most pairs are regular. The bad news is $T(\epsilon)\geq 2^{2^{2^{\iddots^2}}}$ where there are approximately $\frac{1}{\epsilon^2}$ $2$'s. This bound can't be substationally improved, as shown by examples.
\end{remark}

\begin{definition}
A graph is \emph{$\epsilon$-far} from being triangle free if one needs to remove $\geq \epsilon n^2$ edges to make $G$ triangle free.
\end{definition}
\begin{theorem}[Triangle removal lemma] \label{thm:triangle_removal_lemma}
If $G$ is $\epsilon$-far from being triangle-free, then it contains at least $\delta( \epsilon)n^3$ triangles.
\end{theorem}
\begin{proof}
Let $G$ be $\epsilon$-far from being triangle free. Consider an $\frac{\epsilon}{4}$-regularization of $G$ via \cref{thm:regularity_lemma}. We get a partition of $V(G)$ into $T(\epsilon/4)$ parts such that most pairs are $\epsilon/4$-regular.
% \begin{tabular}{c|c|c|c}
% &&&\\
%  \hline\\
%  \hline\\
%  \hline\\
%  \hline
% \end{tabular}

Let us clean up $G$. Remove all the edges from irregular $(V_i,V_j)$'s. Remove all the edges from sparse\sidenote{where here, sparse means the density is less than $\epsilon/2$.} pairs $(V_i,V_j)$, and remove the edges from the diagonal pairs $(V_i,V_i)$. This is less than $\frac{3}{4}\epsilon n^2$  edges in total.

Since $G$ was $\epsilon$-far from being triangle free, the cleaned-up version still has triangles. Since we have at least one triangle, there is $V_i,V_j,V_k$ with a triangle\sidenote{Recall, we removed all edges from $(V_i,V_i)$, so each vertex of the triangle is in a different $V_i$.}. But then each pair has density at least $\epsilon/2$, is $\epsilon/4$-regular and \cref{lem:eregular_yields_triangles} yields $ \epsilon^3 \left( \frac{n}{T} \right)^3$ triangles.
\end{proof}

Let $A\subset \{1,\dotsc,N\}$ with density at least $\epsilon$. Let $U,V,W$ be...

For every $s\in A$ we add a triangle $x\in U$, $x+s \in V$, $x+2s\in W$ for each $x$. Label the edges by $s$.
\missingfigure{3 ovals, $U=\{1,\dotsc,N\},V,W$. Label $x\in U$, connect to $x+s\in V$, $x+2s \in W$ to make a triangle. Label the edges with $s$.}

These are $N\cdot |A|$ edge-disjoint triangles, since each $(x,s)$ yields a new triangle: any one edge uniquely determines the triangle because it determines $x$ and $s$. We need to remove at least $|A|N \geq \epsilon N^2$ edges to make this graph triangle free. So by \cref{thm:triangle_removal_lemma}, there are $\sim N^3$ triangles, so we can find a new triangle different from any determined by $(x,s)$.

Then the labels $(s_1,s_2,s_3)$ of this new triangle form a 3-AP in $A$. We can check by the way we constructed these triangles that $s_1+s_3=2s_2$.

Check: we have an edge $(x,x+s_1)$, $(x+s_1,y+2s')$, and $(y+2s',x)$. Then $s_2 = y+2s'- x - s_1$, and $s_3 = \frac{1}{2}(y+2s'-x)$. So $s_1 + s_2 = y+ 2s' - x = 2s_3$, as we wanted.

max_stars_repo_name: ericphanson/CombinatoricsNotes

max_issues_repo_name: ericphanson/CombinatoricsNotes

max_forks_repo_name: ericphanson/CombinatoricsNotes


text:
\documentclass[a4paper]{scrartcl}

\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}

\title{The TTC 2019 TT2BDD Case (v1.1)}
\author{
  Antonio García-Domínguez\\
  Aston University\\
  B4 7ET, Birmingham, United Kingdom\\
  a.garcia-dominguez@aston.ac.uk
}

\usepackage{tikz}
\usepackage{graphicx}
\graphicspath{{}{figures/}}

\usepackage[pdftex,colorlinks=true]{hyperref}

\usepackage{listings}
\lstset{columns=flexible}

\newcommand*{\class}[1]{\textsc{#1}}
\newcommand*{\feature}[1]{\emph{#1}}
\newcommand*{\file}[1]{\texttt{#1}}

\begin{document}

\maketitle

\begin{abstract}
  Model transformation tools have reached a considerable level of maturity in
  the core features, and are currently developing in many directions. Some tools
  are focusing on providing higher performance for large models or complex
  transformations. Others focus on bidirectionality, visualisation,
  traceability, or verifiability, among other research directions. Whereas past
  cases in TTC have focused on specific research directions, this case study
  presents a well-known simple transformation and welcomes researchers to apply
  their research to it. The aim of this case is to serve as a showcase of the
  various directions that model transformation research is going towards at the
  moment.
\end{abstract}

\section{Introduction}

Past editions of the Transformation Tool Contest have focused on a variety of
topics:
\begin{itemize}
\item In 2018, the Quality-based Software Selection and Hardware-Mapping
  case~\cite{gotz_quality-based_2018} discussed optimisation-oriented model
  transformations (with a combination of performance and solution quality). The
  Social Media Live Case~\cite{hinkel_ttc_2018} considered performance in
  updating model views as models changed (with a strong preference for
  approaches supporting incrementality).

\item In 2017, the Smart Grid case focused on
  incrementality~\cite{hinkel_ttc_2017}, the Families to Persons case discussed
  bidirectional transformations~\cite{anjorin_families_2017}, State Elimination
  focused on performance and the live case on Transformation
  Reuse~\cite{live2017} discussed mechanisms to share complex logic across
  multiple versions of a transformation.

\item In 2016, optimisation-oriented model transformations were discussed in
  considerable breadth through the Class Responsibility Assignment
  case~\cite{fleck_class_2016}, and an alternative dataflow-based notation for
  model transformation was evaluated in the live case study~\cite{live2016}.
\end{itemize}

While these were notable examples of realistic transformations, they were
narrowly focused on a specific topic,and their inherent complexity discouraged
some attendees from trying their hand with their own research agenda on the
transformation.

This year, we propose a broader contest that welcomes all active lines of work
on model transformation. It is based on a simpler, well-known transformation
from the ATL Zoo~\cite{atlzoo}: TT2BDD (Truth Tables to Binary Decision
Diagrams). Striving for raw performance is an option, but the case welcomes
approaches that focus on other attributes of interest of a high-quality model
transformation: for example, verifiability, traceability, bidirectionality, or
understandability. Solution providers are welcome to propose their own
attributes of interest. In general, this case is proposed as a showcase of the
current variety of model transformation tools.

All resources for this case are available on
Github\footnote{\url{https://github.com/TransformationToolContest/ttc2019-tt2bdd}}.
Please follow the link in the footnote and create a pull request with your own
solution after you have submitted your description to EasyChair.

The rest of the document is structured as follows:
Section~\ref{sec:transf-descr} describes the TT2BDD transformation.
Section~\ref{sec:task-suggestions} several tasks of interest that
should be tackled in a solution (authors are free to propose their own tasks of
interest). Section~\ref{sec:benchmark-framework} mentions the benchmark
framework for those solutions that focus on raw performance. Finally,
Section~\ref{sec:evaluation} mentions an outline of the initial audience-based
evaluation across all solutions, and the approach that will be followed to
derive additional prizes depending on the attributes targeted by the solutions.

\section{Transformation Description}
\label{sec:transf-descr}

This section introduces the Truth Tables to Binary Decision Diagram
transformation, with a description of the input and output metamodels, and an
outline of an implementation.

\subsection{Input Metamodel: Truth Tables}
\label{sec:input-metam-truth}

\begin{figure}
  \centering
  \includegraphics[width=.9\textwidth]{tt}
  \caption{Class diagram for the input Truth Tables metamodel}
  \label{fig:tt-metamodel}
\end{figure}

\begin{table}
  \centering
  \begin{tabular}{llll|l}
    $A$ & $B$ & $C$ & $D$ & $S$ \\
    \hline
    0 & 0 & - & - & 0 \\
    0 & 1 & 0 & 0 & 1 \\
    0 & 1 & 0 & 0 & 1 \\
    0 & 1 & 0 & 1 & 0 \\
    0 & 1 & 1 & - & 0 \\
    1 & 0 & 0 & 0 & 0 \\
    1 & 0 & 1 & 0 & 1 \\
    1 & - & - & 1 & 0 \\
    1 & 1 & 0 & 0 & 1 \\
    1 & 1 & 1 & 0 & 0 \\
  \end{tabular}
  \caption{Example truth table: $A$ to $D$ are input ports and $S$ is an output port. ``-'' means ``ignored''.}
  \label{tab:tt-example}
\end{table}

The input metamodel is shown on Figure~\ref{fig:tt-metamodel}. The
\class{Truth\-Table} class is as the root of the model, and contains a
collection of \class{Port}s and \class{Row}s. \class{Port}s come in two types:
\class{Input\-Port}s and \class{Output\-Port}s. \class{Row}s contain sequences
of \class{Cell}s, which assign values to the \class{Input\-Port}s and
\class{Output\-Port}s of the table.

Automated EMF opposite references are used liberally in the metamodel to
simplify the specification of the transformation. \feature{owner} is used to
access the container from several child objects (e.g. the \class{Truth\-Table}
of a \class{Row}), and it is possible to see which \feature{cells} referenced a
specific \class{Port}.

A sample model is shown on Table~\ref{tab:tt-example}. The truth table has four
\class{Input\-Port}s named $A$ to $D$, and one \class{Output\-Port} named $S$.
The first \class{Row} contains only 3 \class{Cell}s, specifying that if $A$ and
$B$ are 0, then $S$ should be 0.

\subsection{Output Metamodel: Binary Decision Trees}
\label{sec:outp-metam-binary}

\begin{figure}
  \centering
  \includegraphics[width=.9\textwidth]{bdd}
  \caption{Class diagram for the output Binary Decision Diagram metamodel, in the original tree-based version.}
  \label{fig:bdd-metamodel}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=.9\textwidth]{bddGraph}
  \caption{Changed classes in the graph-based variant of the Binary Decision Diagram metamodel (2019-05-26). Omitted clases remain the same.}
  \label{fig:bddg-metamodel}
\end{figure}

\begin{figure}
  \centering
  \begin{tikzpicture}[
      level 1/.style={level distance=12em},
      level 2/.style={level distance=4em},
      level 6/.style={level distance=3em},
      level 7/.style={level distance=2em},
      port/.style={circle, draw},
      value/.style={draw},
      assignment/.style={fill=blue!20},
  ]
    \node[port] {$A$}
    child[grow=left] {
      node[value] {0}
      [grow=down]
      child {
        node[port] {$B$}
        child[grow=left] {
          node[value] {0}
          child[grow=down] { node[assignment] {$S = 0$} }
        }
        child[grow=right] {
          node[value] {1}
          child[grow=down] {
            node[port] {$C$}
            child[grow=left] {
              node[value] {0}
              child[grow=down] {
                node[port] {$D$}
                child[grow=left] {
                  node[value] {0}
                  child[grow=down] { node[assignment] {$S=1$} }
                }
                child[grow=right] {
                  node[value] {1}
                  child[grow=down] { node[assignment] {$S=0$} }
                }
              }
            }
            child[grow=right] {
              node[value] {1}
              child[grow=down] {node[assignment] {$S=0$}}
            }
          }
        }
      }
    }
    child[grow=right] {
      node[value] {1}
      [grow=down]
      child {
        node[port] {$D$}
        child[grow=left] {
          node[value] {0}
          child[grow=down] {
            node[port] {$B$}
            child[grow=left] {
              node[value] {0}
              child[grow=down] {
                node[port] {$C$}
                child[grow=left] {
                  node[value] {0} child[grow=down] { node[assignment] {$S=0$} }
                }
                child[grow=right] {
                  node[value] {1} child[grow=down] { node[assignment] {$S=1$} }
                }
              }
            }
            child[grow=right] {
              node[value] {1}
              child[grow=down] {
                node[port] {$C$}
                child[grow=left] {
                  node[value] {0} child[grow=down] { node[assignment] {$S=1$} }
                }
                child[grow=right] {
                  node[value] {1} child[grow=down] { node[assignment] {$S=0$} }
                }
              }
            }
          }
        }
        child[grow=right] {
          node[value] {1} child[grow=down] { node[assignment] {$S=0$} }
        }
      }
    }
    ;
  \end{tikzpicture}
  \caption{Equivalent BDD for the truth table on Table~\ref{tab:tt-example}}
  \label{fig:bdd-equivalent}
\end{figure}

The original output metamodel from the ATL Zoo version is shown on
Figure~\ref{fig:bdd-metamodel}. The \class{BDD} class is the root of the model,
and contains the root of the \class{Tree} and a collection of \class{Port}s.
Similarly to the Truth Tables metamodel, there are \class{Input\-Port}s and
\class{Output\-Port}s.

Figure~\ref{fig:bddg-metamodel} shows the changes introduced into a revised
version of the original metamodel, which allow the \class{BDD} to be a rooted
acyclic graph and not just a tree. This would allow the \class{BDD} to reuse
common subtrees, which would produce more compact circuits in some situations.
The case provides both versions as separate metamodels: solution implementers
are free to choose either version.

\class{Tree} is the common superclass for any node in the tree. Inner nodes
check the value of an \class{Input\-Port}: if it is a false value, evaluation
will proceed through the \feature{tree\-For\-Zero} \class{Tree}; otherwise,
evaluation will go through the \feature{tree\-For\-One} \class{Tree}. Leaf nodes
are \class{Leaf} objects, which provide an \class{Assignment} of a boolean value
to each of the available \class{Output\-Port}s.

The equivalent BDD for the truth table on Table~\ref{tab:tt-example} is shown on
Figure~\ref{fig:bdd-equivalent}. \class{Subtree}s are represented by the circle
referencing an \class{Input\-Port} and their subtrees for when the port takes a
0 or 1 value. \class{Assignment}s are represented by the highlighted nodes that
provide values to the \class{Output\-Port}s.

\subsection{Process Outline}
\label{sec:process-outline}

The transformation essentially needs to construct a (preferably minimal) binary
decision diagram that produces the same values for the output ports, given the
same values in the input ports. Given the high interest about this problem in
circuit design, many algorithms have been proposed in the literature.

Solution authors are welcome to implement a more optimal algorithm in their
transformation tool. This section will only outline a simple approach that can
be readily implemented without too much complication.

There are some basic mappings which are immediately obvious:

\begin{itemize}
\item Each \class{TruthTable} object should correspond to a \class{BDD} object,
  with the same name and equivalent \class{Port}s.

\item Each \class{Input\-Port} and \class{Output\-Port} should be mapped to an
  object of the BDD type with the same name, and should be assigned the same
  \feature{name}.

\item Each \class{Row} should become a \class{Leaf} node: the \class{Cell}s for
  the \class{Output\-Port}s will become \class{Assignment}s.
\end{itemize}

The complexity is in deriving the inner nodes: the \class{Subtree} objects. One
simple approach is to find a TT \class{Input\-Port} which is (ideally) defined
in all the \class{Row}s, and turn it into an inner node (a \class{Subtree})
which points to the equivalent BDD \class{Input\-Port} and has two
\class{Tree}s:
\begin{itemize}
\item The zero subtree, produced from the \class{Row}(s) where the port was 0.
  This will be a \class{Subtree} if there are at least two rows in that
  situation: the transformation should proceed recursively in this case with
  those rows, excluding the input ports that have already been considered. If
  there is only one such row, this would simply point to the \class{Leaf}
  created above.

\item The one subtree, produced in a parallel way to the zero subtree.
\end{itemize}

This simple approach does not necessarily ensure a minimal subtree, as in some
points there may be multiple ports to choose from. It may require improvements
for cases where there are no input ports which are defined in all available
rows. Authors are welcome to implement a more efficient or optimal approach in
their solutions.

\section{Main Task}
\label{sec:task-suggestions}

The obvious main task in this case is to implement the TT2BDD transformation in
your favourite model transformation language, ensuring that the resulting BDDs
are equivalent to the original truth table. To simplify the work involved, the
case includes an Eclipse Modelling Framework implementation of the metamodels in
the \file{ttc2019.metamodels} project within the \file{metamodels} folder, and a
set of sample XMI input models conforming to the metamodels in the \file{models}
folder.

The \file{models} folder includes two additional tools, which can be run with
\file{java -jar tool.jar}:
\begin{itemize}
\item \file{generator.jar} can produce truth tables of an arbitrary number of
  input and output ports, given seed. Solution implementers may want to be
  careful with the values given, as truth tables grow exponentially as more
  input ports are added.

\item \file{validator.jar} can validate if a BDD model is equivalent to a source
  truth table model, by evaluating all the rows in the TT model through the BDD
  and comparing the values of the \class{Output\-Port}s. This tool is required
  to check that your transformation is correct. The tool accepts both tree-based
  and graph-based BDD models.
\end{itemize}

Solutions are free to focus on any specific quality attribute of the
transformation beyond optimality and performance. For instance, it may be useful
to be able to prove that the transformation does indeed produce a BDD which is
equivalent to the TT (even if suboptimal). It may also be interesting to be able
to visualize the BDD being built from the TT, or to be able to trace which
subtree came from which collection of rows.

The conciseness of the transformation may also be of interest, given that its
recursive nature may be difficult to map to some existing transformation
languages.

\section{Benchmark Framework}
\label{sec:benchmark-framework}

If focusing on performance, the solution authors should integrate their solution
with the provided benchmark framework. It is based on that of the TTC 2017 Smart
Grid case~\cite{hinkel_ttc_2017}, and supports the automated build and execution
of solutions. For this specific case study, the visualisation of the results is
currently disabled.

The benchmark consists of three phases:

\begin{enumerate}
\item \textbf{Initialization}, which involves setting up the basic
  infrastructure (e.g. loading metamodels). These measurements are optional.
\item \textbf{Load}, which loads the input models.
\item \textbf{Run}, which runs the transformation, produces the BDD and saves
  the file.
\end{enumerate}

\subsection{Solution requirements}
\label{sec:solut-requ}

Solutions should be forks of the main Github
project\footnote{\url{https://github.com/TransformationToolContest/ttc2019-tt2bdd}},
and should be submitted as pull requests after the descriptions have been
submitted to EasyChair. For details about the submission format, please consult
the call for
solutions\footnote{\url{https://www.transformation-tool-contest.eu/cfs.html}}.

Each solution wishing to use the benchmarking framework should print to the
standard output a line with the following fields, separated by semicolons
(``;''):

\begin{itemize}
\item \textbf{Tool}: name of the tool.
\item \textbf{Model}: name of the model being transformed.
\item \textbf{RunIndex}: the index of the repetition of the transformation.
\item \textbf{Phase}: the name of the phase.
\item \textbf{MetricName}: the name of the metric. It may be the used
  \textbf{Memory} in bytes, or the wall clock \textbf{Time} spent in integer
  nanoseconds.
\end{itemize}

\lstinputlisting[
  float,frame=tb,numbers=left,
  caption={\file{solution.ini} file for the reference ATL solution},
  label=lst:ini-atl
]{../../solutions/EMFSolutionATL/solution.ini}

To enable automatic execution by the benchmark framework, solutions should add a
subdirectory to the \file{solutions} folder of the benchmark with a
\file{solution.ini} file stating how the solution should be built and how it
should be run. As an example, the \file{solution.ini} file for the reference ATL
solution is shown on Listing~\ref{lst:ini-atl}. In the \file{build} section, the
\file{default} option specifies the command to build and test the solution, and
the \file{skipTests} option specifies the command to build the solution while
skipping unit tests. In the \file{run} section, the \file{cmd} option specifies
the command to run the solution.

The repetition of executions as defined in the benchmark configuration is done
by the benchmark. For 5 runs, the specified command will be called 5 times,
passing any required information (e.g. run index, or input model name) through
environment variables. Solutions must not save intermediate data between
different runs: each run should be entirely independent.

The name and absolute path of the input model, the run index and the name of the
tool are passed using environment variables \file{Model}, \file{ModelPath},
\file{RunIndex} and \file{Tool}. Solution authors are suggested to study the
example ATL solution on how to use these values to run their transformation.

\subsection{Running the benchmark}
\label{sec:running-benchmark}

The benchmark framework only requires Python 3.3 to be installed. Furthermore,
the solutions may imply additional frameworks. We would ask solution authors to
explicitly note dependencies to additional frameworks necessary to run their
solutions.

If all prerequisites are fulfilled, the benchmark can be run using Python with
the command \file{python scripts/run.py}. Additional options can be queried
using the option \file{{-}{-}help}. The benchmark framework can be configured
through the \file{config/config.json} file: this includes the input models to be
evaluated (some of which have been excluded by default due to their high cost
with the sample solution), the names of the tools to be run, the number of runs
per tool+model, and the timeout for each command in milliseconds.

\section{Evaluation}
\label{sec:evaluation}

Given the call for a broader set of research interests in this transformation,
the evaluation will operate on two dimensions:

\begin{itemize}
\item Is it a high-quality model transformation? The transformation should be
  complete, correct, easy to understand, efficient, and produce optimal results.

  The \file{validator.jar} will be used to check the produced BDDs against the
  source truth tables, and the authors will need to provide a convincing
  argument about the correctness of the solution.

  Attendees to the contest will evaluate the understandability of the solution,
  and the benchmark framework will provide independent measurements of memory
  and time usage.

  Tree sizes for the BDDs will be considered for the optimality of the
  transformations: the smaller, the better.

\item Does it highlight a promising research direction? Although the
  transformation may not be entirely complete or may be harder to understand, it
  may serve as an example of an active research area within model
  transformations that the community may wish to showcase.

  This may include aspects such as incrementality, bidirectionality,
  traceability, verifiability, or the ability to visualize interactively the
  transformation, among other areas of interest in the field.

  Depending on the submitted solutions, the awards for this dimension may
  change. If there are no common areas of interest in the solutions, an
  audience-driven ``Most Promising'' award will be awarded. If there are common
  themes between solutions, dedicated 1st/2nd/3rd place awards for those
  research areas may be awarded.
\end{itemize}

\bibliographystyle{plain}
\bibliography{bibliography}

\end{document}


text:
\subsection{Tumor Necrosis Factor Receptor Associated Syndrome}

Tumor necrosis factor receptor-associated periodic syndrome (TRAPS) is a rare (1 : 1,000,000) hereditary autosomal dominant disease , formerly known as Familial Hibernian Fever (FHF) \cite{roth-wojcicki_tumor_2019},  that can cause recurring fevers with durations from days up to months.
Symptoms during these fevers are: skin rash, swelling, inflammatory reactions across the whole body and pain in the abdomen, muscles and/or joints. A long term and lasting effect is the accumulation of amyloids within the kidneys and may result in other diseases \cite{nih_traps_2016}.
%Tumor Necrosis Factor Receptor Associated Periodic Syndrome (Juvenile), Roth-Wojcicki
%TRAPS, NIH
TRAPS is known to be caused by mutations within the gene turmor necrosis factor receptor 1 (TNFRSFA1/TNRF1) (Section \ref{subsec:CD_TNFRSF1A}), of which few mutations are known to hinder the transport of proteins to the cell surface and trap them in cell causing the activation of inflammatory response \cite{kimberley_falling_2007, nih_traps_2016}.
% Falling into TRAPS--receptor misfolding in the TNF receptor 1-associated periodic fever syndrome, Kimberly
%TRAPS, NIH
So far 158 mutations have been associated with the disease \cite{aksentijevich_tumor-necrosis-factor_2001}, but more mutations have been identified in TNFRSF1A wherein some might be pathogenic (Sections \ref{subsec:MM_GAVIN_data_table}, \ref{subsec:MM_GnomAD}).
%Tabular list ,Aksentijevich
\label{subsec:CD_TRAPS}

\subsection{Tumor Necrosis Factor Receptor Super Family Member 1A}
Tumor Necrosis Factor Receptor Super Family Member 1A (TNFRSF1A/TNFR1) is a gene located on chromosome 12 region 1 band 3 and sub-band 31. The gene produces a trans-membrane receptor consisting of 445 residues divided into 221 residue cytoplasmic section and a 171 extracellular part that consists of 4 conserved cysteine rich domains \cite{gray_cloning_1990,walter_encyclopedia_2007,banner_crystal_1993}.
% Cloning of human tumor necrosis factor (TNF) receptor cDNA and expression of recombinant soluble TNF-binding protein, Gray
%Encyclopedia of Molecular Pharmacology page 1247-1251 book, 1240-1244 pdf, Walter Rosenthal
%Crystal structure of the soluble human 55 kd TNF receptor-human TNFβ complex: Implications for TNF receptor activation, Banner
The receptor is ubiquitous across most cell surfaces , but is not found on erythrocytes \cite{segueni_innate_2016}, and can form two different types of unbound hexagonal clusters depending on the dimer formation \cite{naismith_crystallographic_1995}.
%Innate myeloid cell TNFR1 mediates first line defence against primary Mycobacterium tuberculosis infection, Segueni
%Crystallographic Evidence for Dimerization of Unliganded Tumor Necrosis Factor Receptor, Naismith
When the structures are dimers the binding sites are exposed making it possible for tumor necrosis factor (TNF) $\alpha$ and -$\beta$  (Section \ref{subsec:CD_TNF_A_B}) to bind in homotrimeric form, by binding of TNF the dimers disconnect and three TNFR1s interact with the TNF homotrimer \cite{naismith_crystallographic_1995}.
%Crystallographic Evidence for Dimerization of Unliganded Tumor Necrosis Factor Receptor, Naismith
With the interaction of TNF homotrimers with TNFR1 it can activate several pathways such as the nuclear factor kappa-light-chain-enhancer of activated B cells (NF-$\kappa$B), which enhances the transcription of various genes during inflammation, infection or other forms of external stress. Also it is able to activate the extrinsic pathway of apoptosis by releasing the silencer of death domain (SODD) proteins, on the cytoplasmic site of the cell, which will activate the formation of Tumor Necrosis Factor Receptor type 1-Associated DEATH Domain protein (TRADD) \cite{bender_adaptor_2005}
%The adaptor protein TRADD activates distinct mechanisms of apoptosis from the nucleus and the cytoplasm, Bender
 that forms a complex that attracts Fas associated death domain (FADD) \cite{muppidi_life_2004}.
%Life And Death Decisions: Secondary Complexes and Lipid Rafts in TNF Receptor Family Signal Transduction,Muppidi
On binding of FADD initiator caspase 8  starts a cascade wherein caspase 3 is activated an will cleave aspartate out of proteins and thereby disrupting the metabolism\cite{vinay_robbins_2014,chen_tnf-r1_2002,hengartner_biochemistry_2000}.
%Robbins and Cotran Pathologic Basis of Disease, Professional Edition, chapter 2 Necroptosis,  Page 59 (Book) 78 (pdf), Vinay
%TNF-R1 Signaling: A Beautiful Pathway, Chen and Goeddel
%The biochemistry of apoptosis, Hengartner
\label{subsec:CD_TNFRSF1A}

\subsection{Tumor Necrosis Factor Alpha and Beta}
The proteins TNF $\alpha$ and $\beta$ are both pro-inflammatory cytokines that are produced as response to an infection or when a cell is damaged. Both are transcribed from their genes that reside in chromosome 6 in the p-arm at region 2 band 1 and sub-band 3. TNF $\alpha$ and $\beta$ are 35\% identical and 50\% homologous to each other consisting out of 233 and 205 amino acid residues. Both are able to form a homotrimeric structures that can bind to the dimeric form TNFR1 (Section \ref{subsec:CD_TNFRSF1A}) to activate the extrinsic pathway\cite{aggarwal_characterization_1985,hamosh_omim_2014,kriegler_novel_1988,walter_encyclopedia_2007}.
%Characterization of receptors for human tumour necrosis factor and their regulation by γ-interferon, Aggarwal
%OMIM Entry - * 153440 - LYMPHOTOXIN-ALPHA; LTA, Hamosh
%A novel form of TNF/cachectin is a cell surface cytotoxic transmembrane protein: Ramifications for the complex physiology of TNF, Kriegler
%Encyclopedia of Molecular Pharmacology page 1247-1251 book, 1240-1244 pdf, Walter Rosenthal
\label{subsec:CD_TNF_A_B}

max_stars_repo_name: Sylt-CSI/variant-protein-prediction

max_issues_repo_name: Sylt-CSI/variant-protein-prediction

max_forks_repo_name: Sylt-CSI/variant-protein-prediction


text:
%\documentstyle[epsf,12pt]{article}
\documentstyle[epsf,11pt]{article}
\begin{document}

\title{
Reducing quasi-ergodicity in a double well potential by Tsallis Monte
Carlo simulation
}

\author{
Masao Iwamatsu$^{\dagger}$\thanks{Permanent address.
E-mail:iwamatsu@ce.hiroshima-cu.ac.jp} and Yutaka Okabe$^{\dagger}$
\\$^{*}$Department of Computer Engineering, Hiroshima City University\\
Hiroshima 731-3194, Japan\\ and\\
$^{\dagger}$Department of Physics, Tokyo Metropolitan University\\
Hachioji, Tokyo 192-0397, Japan }

\date
{
%\today
}

\maketitle

\begin{abstract}
A new Monte Carlo scheme based on the system of Tsallis's generalized
statistical mechanics is applied to a simple double well potential
to calculate the canonical thermal average of potential energy.
Although we observed serious quasi-ergodicity when using the standard
Metropolis Monte Carlo algorithm, this problem is largely reduced by
the use of the new Monte Carlo algorithm.  Therefore the ergodicity
is guaranteed even for short Monte Carlo steps if we use this new
canonical Monte Carlo scheme. \\

\begin{flushleft}

PACS: 02.70.Lq; 05.70.-a

Key words: Monte Carlo; Tsallis statistics; double well potential\\

\end{flushleft}
\end{abstract}

\newpage
\section{Introduction}

The ergodic hypothesis is fundamental to statistical mechanics.
This hypothesis states that the time average of an observable
event equals the phase-space average. In practical application,
however, problems can arise in various types of simulations
if the system must overcome high energy barriers to reach other
regions of phase space. In that case, the length of a simulation
needed in order to obtain enough statistical samples of all
regions of phase space may be extremely long. In the Monte Carlo
simulation, the errors arise as a consequence of the finite
length of the Monte Carlo walk. This error can be serious in
canonical Monte Carlo sampling, especially at low temperatures.
This problem, referred as "quasi-ergodicity" by Valleau and
Whittington~\cite{VW}, appears even in the simplest double well
problem, where the two wells are separated by a large
barrier~\cite{FFD}.

Recently several authors~\cite{Pe,TS,AS1} have proposed a
generalized simulated annealing method to locate the global
minimum based on the generalized statistical mechanics
proposed by Tsallis~\cite{Ts}, which can overcome the slow
convergence of the traditional simulated annealing
method~\cite{KGV} based on the standard Monte Carlo scheme.
Subsequently, Andricioaei and Straub have pointed out that
their generalized Monte Carlo (Tsallis Monte Carlo)~\cite{AS2}
method may be used to overcome the quasi-ergodicity appearing
in the canonical average at constant temperatures.

In this short note we use a new Monte Carlo scheme~\cite{AS2,AS3,SA}
based on Tsallis's generalized statistical mechanics~\cite{TS}
for the calculation of the thermal average of the potential energy
of a simple double well potential at constant temperatures.
We examine the performance of this new algorithm designed to
overcome the quasi-ergodicity.

\section{A generalized Monte Carlo scheme}

In the generalized statistical mechanics proposed by
Tsallis~\cite{Ts}, a crucial role is played by the generalized
entropy $S_{q}$ defined as
\begin{equation}
S_{q}=k\frac{1-\sum p_{i}^{q}}{q-1}
\label{eq:1}
\end{equation}
where $q$ is a real number which characterizes the statistical
mechanics, and $p_{i}$ is the probability of states $i$. This
entropy $S_{q}$ becomes the usual Gibbs-Shannon entropy
$S_{1}=-k\sum p_{i}\ln p_{i}$ when $q\rightarrow 1$. By
maximizing the generalized entropy with the constraints
\begin{equation}
\sum p_{i}=1\;\;\;\mbox{and}\;\; \sum p_{i}^{q}\epsilon_{i}=\mbox{const}
\label{eq:2}
\end{equation}
where $\epsilon_{i}$ is the energy spectrum, the statistical weight
for the generalized canonical ensemble characterized by the
parameter $q$ is given by the Tsallis weight
\begin{equation}
p_{i}^{q}=\exp\left(-\overline{\beta}\overline{\epsilon}_{i}\right)/Z_{q}^{q}
\label{eq:3}
\end{equation}
with
\begin{equation}
\overline{\epsilon}_{i}=\frac{q}{\overline{\beta}(q-1)}
\ln\left[1-(1-q)\overline{\beta}\epsilon_{i }\right]
\label{eq:4}
\end{equation}
where $\overline{\epsilon}_{i}$ is the generalized energy and
$\overline{\beta}$ is the Lagrange multiplier, which plays the
role of the generalized temperature in the Tsallis statistical
mechanics.  We should note here that the statistical weight is
{\it not} given by the probability distribution $p_{i}$ {\it but}
rather by $p_{i}^{q}$ in order to calculate thermal averages, as
can be seen in (\ref{eq:2}). The generalized canonical partition
function $Z_{q}$ in (\ref{eq:3}) is given by
\begin{equation}
Z_{q}=\sum \exp\left(-\overline{\beta}\overline{\epsilon}_{i}/q\right)
\label{eq:5}
\end{equation}
In the limit $q\rightarrow 1$, $\overline{\epsilon}_{i}\rightarrow
\epsilon_{i}$
and $\overline{\beta}\rightarrow\beta=1/kT$ where $T$ is the usual
temperature~\cite{CT}, then the Tsallis weight (\ref{eq:3}) becomes the
Boltamann weight $p_{i}=\exp(-\beta\epsilon_{i})/Z_{1}$ of the usual
canonical ensemble.

In this generalized statistical mechanics, Andricioaei and
Straub~\cite{AS1} have noted that since the average of an
observable $O$ is defined by
\begin{equation}
<O>_{q}=\sum p_{i}^{q}O_{i},
\label{eq:6}
\end{equation}
the detailed balance condition should be written as
\begin{equation}
p_{i}^{q}W_{ij}=p_{j}^{q}W_{ji}
\label{equation}
\end{equation}
where $W_{ij}$ is the element of the transition matrix. Based
on this observation, they~\cite{AS1} pointed out that the
generalized Monte Carlo algorithm~\cite{TS} originally proposed
by Penna~\cite{Pe} in his generalized simulated annealing method
cannot satisfy the detailed balance.  They~\cite{AS1} then
proposed a new generalized Monte Carlo scheme, where the
acceptance probability $p$ of the Monte Carlo move is given by
\begin{equation}
p=\min\left[1,\exp(-\overline{\beta}\Delta \overline{\epsilon})\right]
\label{eq:7}
\end{equation}
which has a similar form to the familiar Metropolis
algorithm~\cite{MRRTT}
\begin{equation}
p=\min\left[1,\exp(-\beta\Delta \epsilon)\right] \label{eq:8}
\end{equation}
where $\Delta\overline{\epsilon}$ is the increase of the
generalized energy (\ref{eq:4}) while $\Delta\epsilon$ is that
of the usual energy. Andricioaei and Straub~\cite{AS1} combined
this Monte Carlo algorithm with the usual simulated annealing
technique~\cite{KGV} to optimize the conformation of tetrapeptides.
They showed that their algorithm is more effective than the
standard simulated annealing based on the usual molecular
dynamics or Monte Carlo methods.

Recently, Andricioaei and Straub~\cite{AS2,AS3} pointed out that their
generalized Monte Carlo scheme~\cite{AS1} for simulated annealing can
also be used to calculate the standard canonical ensemble averages at
constant temperatures
\begin{equation}
<O>_{1}=\frac{\sum O_{i}\exp(-\beta \epsilon_{i})}
{\sum \exp(-\beta \epsilon_{i})}
\label{eq:9}
\end{equation}
using the general rule of the importance sampling~\cite{VW}
\begin{equation}
<f>=\sum f_{i}=\sum \left(\frac{f_{i}}{g_{i}}\right)g_{i}
\label{eq:8b}
\end{equation}
Then, equation (\ref{eq:9}) can be calculated using the generalized
Monte Carlo scheme from
\begin{eqnarray}
<O>_{1}&=&\frac{ \sum
O_{i}\exp\left(-\left(\beta\epsilon_{i}-\overline{\beta}
\overline{\epsilon}_{i}
\right)\right)\exp\left(-\overline{\beta}
\overline{\epsilon}_{i}\right) }{\sum
\exp\left(-\left(\beta\epsilon_{i}-\overline{\beta}
\overline{\epsilon}_{i}\right)
\right)\exp(-\overline{\beta}\overline{\epsilon}_{i})}
\label{eq:10}\\
&=& \frac{< O \exp\left(-\left(\beta\epsilon-
\overline{\beta}\overline{\epsilon}
\right)\right)>_{q}}{<\exp\left(-\left(\beta\epsilon-
\overline{\beta}\overline{\epsilon}\right)\right)>_{q}}
\label{eq:11}
\end{eqnarray}
where the generalized average $<>_{q}$ is defined by (\ref{eq:6}).
In practice, the equilibrium thermal average in the standard
canonical ensemble can be calculated by conducting these
"generalized" Monte Carlo steps and accumulating the samples by
\begin{equation}
<O>_{1}=\frac{\sum O_{j} \exp\left(-\left(\beta\epsilon_{j}-
\overline{\beta}\overline{\epsilon}_{j}\right)\right)}
{\sum \exp\left(-\left(\beta\epsilon_{j}-
\overline{\beta}\overline{\epsilon}_{j}\right)\right)}
\label{eq:12}
\end{equation}
where the sum $j$ runs over to the generalized Tsallis Monte
Carlo steps produced by eq.(\ref{eq:7}) instead of the usual
Metropolis algorithm (\ref{eq:8}). The equation (\ref{eq:12})
is a generalized expression, and we may choose
$\overline{\beta}=\beta$.  Then, the equation (\ref{eq:12})
reduces to
\begin{equation}
<O>_{1}=\frac{\sum O_{j} \exp\left(-\beta\left(\epsilon_{j}-
\overline{\epsilon}_{j}\right)\right)}
{\sum \exp\left(-\beta\left(\epsilon_{j}-
\overline{\epsilon}_{j}\right)\right)},
\label{eq:12x}
\end{equation}
which has been used in \cite{AS2,AS3}.  From now on, the actual
calculation will be done by choosing $\overline{\beta}=\beta$.

We would like to point out here that a similar reweighting
technique is used in the histogram method~\cite{FS} and
the multicanonical ensemble method~\cite{BN}. Andricioaei and
Straub applied this Tsallis Monte Carlo algorithm to the
two-dimensional Ising system~\cite{AS2}, a classical one-dimensional
potential and a 13-atom Lennard-Jones cluster~\cite{AS3}, and
showed that this new Monte Carlo algorithm is more effective
than the standard Metropolis algorithm.  They also proposed the
generalized molecular dynamics based on Tsallis generalized statistical
mechanics~\cite{AS3,SA}

\section{Application to a classical double well potential}

In order to check the performance of this new Monte Carlo
algorithm~\cite{AS2,AS3} in a classical system such as molecules
and liquids, we look at the problem of a classical particle
in a one-dimensional double well potential defined by
\begin{equation}
V(x)=\delta\left(3 x^{4}+4(\alpha-1)x^{3}-6\alpha x^{2}\right)+1
\label{eq:13}
\end{equation}
where
\begin{equation}
\delta=\frac{1}{2\alpha+1}
\label{eq:14}
\end{equation}
which was examined by Frantz {\it et al.}~\cite{FFD} to demonstrate
the quasi-ergodicity of the standard Metropolis algorithm. The
potential given by (\ref{eq:13}) has the double minimum at $x=1$
and $x=-\alpha$; therefore, it represents a symmetrical double well
when $\alpha=1$ and a single well when $\alpha=0$. This potential
is more conveniently characterized by the parameter $\gamma$~\cite{FFD}
\begin{equation}
\gamma=\frac{V(0)-V(-\alpha)}{V(0)-V(1)}
=\alpha^{3}\left(\frac{\alpha+2}{2\alpha+1}\right)
\label{eq:15}
\end{equation}
In figure 1, we show the asymmetric double well potential $V(x)$
when $\gamma=0.9$ as well as the equilibrium Tsallis weight
$p_{\rm T}\propto \exp(-\beta \overline{V}(x))$, which is renormalized,
and Boltzmann weight $p_{\rm B}=\exp(-\beta V(x))/Z_{1}$. It is obvious
that the statistical weight $p_{\rm T}$ of the Tsallis generalized
statistical mechanics has a greater chance of crossing the barrier
and, therefore, a greater possibility of avoiding quasi-ergodicity.

Alternatively this Tsallis's statistical mechanics can be regarded
as a search-space smoothing method which deforms the rugged potential
landscape by a smoother one.  In figure 2, we show the smoothed
potential $\overline{V}(x)$ calculated from (\ref{eq:4}), which has
a smoother landscape and a lower energy barrier than the
original landscape $V(x)$.

We use this double well potential to look at the quasi-ergodicity
of the standard as well as the generalized Monte Carlo scheme by
examining the classical average potential energy~\cite{FFD}
\begin{equation}
<V>=\frac{\int V(x)\exp\left(-\beta V(x)\right)dx}
{\int \exp\left(-\beta V(x)\right)dx }
\label{eq:16}
\end{equation}
which can be calculated rigorously by numerical quadrature as a
function of the temperature $\beta$ and the asymmetry parameter
$\gamma$. We employ the same Monte Carlo procedure of Frantz
{\it et al.}~\cite{FFD} and use the uniform sampling distribution
from previous position $x$ to new position $x^{'}$ given by
\begin{eqnarray}
T(x^{'}|x)&=&\frac{1}{\Delta}\;\;\;\mbox{for}
\;\;|x-x^{'}|<\frac{\Delta}{2}\ \nonumber \\
&=& 0\;\;\;\mbox{otherwise}
\label{eq:17}
\end{eqnarray}
The step size $\Delta$ can be chosen to be large enough for the
one-dimensional problem so that the quasi-ergodicity can
be reduced~\cite{FFD}.  However, we follow Frantz {\it et
al.}~\cite{FFD} and use the scaling $\Delta=2.5/\sqrt{\beta}$,
which maintains an approximately 50\% acceptance rate at all
temperatures $\beta$, by considering the application to
multidimensional systems.  This 50\% acceptance is commonly
employed in the classical Monte Carlo community~\cite{AT}.

Since the error caused by the quasi-ergodicity depends upon
the walk initialization, we start the Monte Carlo walk at
the global minimum $x=1$ and at the metastable minimum $x=-\alpha$.
If the walk starts at the global minimum, then the average
potential $<V>$ will be low if the walk is quasi-ergodic, as
the distribution associated with the higher energy well around
$\alpha$-minimum will be insufficiently sampled. While the walk
is initialized at $-\alpha$, the walk may be trapped at this
metastable $\alpha$-well and the average potential $<V>$
will be too high.

In figure 3, we show the average potential energy $<V>$ calculated
from (i) the standard Metropolis Monte Carlo and (ii$\sim$iv) the
generalized Tsallis Monte Carlo algorithm as a function of
temperature $\beta$. Similar curves can be found in figure 2 of
reference \cite{FFD}; however, those authors plotted the results
from the usual Metropolis algorithm initialized at the global
minimum $x=1$ and at the random positions. In order to check
the quasi-ergodicity of the walk, we instead show two curves
initialized at the global minimum $x=1$ and the metastable
$\alpha$-minimum at $x=-\alpha$ for each scheme. The average
potential $<V>$ is obtained from 100 independently initialized
walks, each consisting of 500 warm-up steps followed by $10^{4}$
steps with data accumulation.  The average energy $<V>$ calculated
from the standard Metropolis algorithm shows large errors due to
the quasi-ergodicity. In particular, the walk initialized at
the metastable $\alpha$-minimum greatly overestimates the potential
energy. The walk starting from the global minimum also shows
noticeable error and underestimates the average energy, as expected.

In contrast to the results from the standard Metropolis algorithm,
the results from the generalized Monte Carlo algorithm~\cite{AS2,AS3}
show excellent agreement with the numerically exact result. The
walks started from the global minimum $x=1$ produce the average
energy which is almost indistinguishable from the exact result up
to the lowest temperature (high $\beta$) examined. Even the walk
started from the metastable $\alpha$ minimum produces the energies
that are fairly accurate up to rather low temperature $\beta \sim 15$
for $q=1.5$, which shows that this generalized algorithm does
a much better job than the standard Monte Carlo scheme. The agreement
with the exact result is further improved if we increase the value
of parameter $q$ from $1.5$ (see Figure 3).  It should be noted
that the probability distribution associated with the parameter
$q=2$ is given by the Cauchy distribution and with $q>5/3$ by the
super-diffusive L\'evy distribution~\cite{TS}. The theoretical upper
boundary is $q=3$. The quasi-ergodicity can be largely circumvented
by the use of Tsallis's generalized statistical mechanics with
$q>1$, as shown in figure 3.

In order to overcome the quasi-ergodicity, the Monte Carlo
walk has to pass many times through the energy barrier.
Figure 4 compares the time series of the Monte Carlo walk
generated by the standard Monte Carlo scheme ($q=1$) and that
of the generalized Monte Carlo scheme.  Similar diagrams
can be found in reference ~\cite{AS3}.  Although the
trajectory of the standard Monte Carlo scheme is mostly
confined within the metastable $\alpha$ minimum, the number
of barrier crossing of the Tsallis Monte Carlo scheme is
significantly greater.  This barrier crossing occurs more
often when the parameter $q$ is larger.  The quasi-ergodicity
can be removed by Tsallis's generalized Monte Carlo scheme
by increasing the magnitude of the parameter $q$ from $1$.

In figure 5, we show the average energy $<V>$ as a function of
the asymmetry parameter $\gamma$. Again a similar diagram can
be found in reference~\cite{FFD}. Here again we can
clearly observe a large error for the initialization at the
metastable $\alpha$ minimum when the standard Monte Carlo
procedure is used. In this case, the superiority of Tsallis's
generalized Monte Carlo scheme over the traditional one is
also obvious.

In table 1 we showed the effect of quasi-ergodicity in the
convergence when $\gamma=0.9$ and $\beta=10$ for both
the standard Monte Carlo and the generalized Monte Carlo
algorithms with $q=1.5$, each originating at the metastable
$\alpha$-minimum. We show the average energy $<V>$ and the
standard deviation (STD) obtained from 100 independently
initialized walks. We found that the problem of quasi-ergodicity
obtained with the standard Metropolis algorithm is still present,
even if we increase the number of Monte Carol steps. On the other
hand the result from the generalized Monte Carlo simulation is
significantly improved by increasing the number of Monte Carlo
steps, and ergodicity can be recovered.  Furthermore, STD
decreases as $\sim 1/\sqrt{N}$ when the Monte Carlo step $N$ is
increased for the generalized Monte Carlo simulation, while such
a systemic decrease is not observed for the standard Monte
Carlo simulation.

This convergence can be clearly visualized by plotting the
visiting probability $P(x)$ of the Monte Carlo walks, which
should be given by the Boltzmann weight for the standard
Monte Carlo scheme and by the Tsallis weight (\ref{eq:3}) for
the generalized Monte Carlo scheme.  In figure 6, we compare
$P(x)$ generated by the standard and generalized Monte Carlo
walks with the Boltzmann and the Tsallis weights.  The visiting
probability $P(x)$ converges rapidly to the Tsallis weight as
the number of the Mone Carlo walks is increased for the
generalized Monte Carlo scheme.  On the other hand, $P(x)$
recovers the Boltzmann weight only when the Monte Carlo steps
are increased upto $\sim 10^{6}$ in the standard Monte
Carlo scheme.

It is also interesting to discuss the convergence in the context of
the "generalized ergodic measure" of Thirumalai, Mountain, and
Kirkpatrick~\cite{TMK,AS3}. Using the energy metric, the ergodic measure
$d_{V}(n)$ is defined by
\begin{equation}
d_{V}(n)=\left(V^{a}(n)-V^{b}(n)\right)^{2},
\label{eq:18}
\end{equation}
where $V^{a}(n)$ and $V^{b}(n)$ are the move average of the potential
energy $V$ along the trajectories $a$ and $b$, which
are defined from (\ref{eq:12x}) by
\begin{eqnarray}
V^{a}(n)&=&\frac{\sum_{j=0}^{n} V_{j}^{a} \exp\left(-\beta\left(V_{j}^{a}-
\overline{V}_{j}^{a}\right)\right)}
{\sum_{j=0}^{n} \exp\left(-\beta\left(V_{j}^{a}-
\overline{V}_{j}^{a}\right)\right)}.
\label{eq:19}
\end{eqnarray}
for the generalized Tsallis Monte Carlo scheme, and
\begin{eqnarray}
V^{a}(n)&=&\frac{\sum_{j=0}^{n} V_{j}^{a} \exp\left(-\beta V_{j}^{a}\right)}
{\sum_{j=0}^{n}1}.
\label{eq:19x}
\end{eqnarray}
for the standard Monte Carlo scheme.
For an ergodic system, Thirumalai {\it et al}.~\cite{TMK} suggested that
the ergodic measure converges as $d_{V}(n)\rightarrow \infty$ if
$n\rightarrow \infty$.  They found that
\begin{equation}
d_{V}(n)\simeq d_{V}(0)\frac{1}{D_{V}n}
\label{eq:20}
\end{equation}
where the "diffusion" constant $D_{V}$ depends on temperature.  In figure
7, we show $d_{V}(0)/d_{V}(n)$ for $\gamma=0.9$ potential obtained from 100
independent pairs ($a$, $b$) of Monte Carlo walks as a function of Monte Carlo
steps $n$.  We chose the trajectory starting from the metastable minimum
$x=-\alpha$  (=-0.9163) and that from the stable minimum $x=1$,
respectively, as the two independent trajectories $a$ and $b$.  It is
clear from figure 7 that the normalized inverse of the ergodic measure
$d_{V}(0)/d_{V}(n)$ grows linearly with the Monte Carlo steps $n$.
Therefore, the ergodic measure $d_{V}(n)$ decreases according to
(\ref{eq:20}), and the diffusion constant $D_{V}$ can be determined.

Figure 8 shows the diffusion constant $D_{V}$ as a function of temperature
$\beta$.  The diffusion constant is a decreasing function of the inverse
temperature $\beta$; it is more difficult to recover the ergodicity when
the temperature is lowered.  However, at lower temperatures $\beta>20$ with
$\beta\Delta V>1$ where $\Delta V=V(-\alpha)-V(1)=0.1$, the diffusion
constant starts to increase again for the Tsallis Monte Carlo walks ($q\neq
1$).  This is due to the fact that the equilibrium thermal distribution
around the metastable minimum at $x=-\alpha$ becomes negligibly small at
such low temperatures, and the ergodicity of the trajectory among the
stable and the metastable minima is
irrelevant once trajectories $a$ and $b$ both fall into the stable
minimum around $x=1$.

It is apparent from figure 8 that the diffusion constant $D_{V}$ obtained
from the standard Monte Carlo algorithm obeys the usual activation form
\begin{equation}
D_{V} \sim \exp\left(-\beta E_{b}\right)
\label{eq:21}
\end{equation}
where the activation energy $E_{b}$ is estimated to be $E_{b}\simeq 0.4$,
which is the same order of magnitude as the barrier height ($\sim 1$)
between two wells, as shown in figure 1.   On the other hand, the
diffusion
constant $D_{V}$ obtained from the Tsallis generalized Monte Carlo scheme is
characterized by a power-law form
\begin{equation}
D_{V} \sim \beta^{-\eta}
\label{eq:22}
\end{equation}
with the exponent $\eta$.  This result is conceivable, based on the
acceptance
probability (\ref{eq:7}) using the Tsallis weight, because
\begin{eqnarray}
D_{V} &\sim& \exp(-\beta\bar{E}_{b})
\nonumber \\
&\sim& \left[1-(1-q)\beta E_{b} \right]^{-q/(q-1)}
\label{eq:23}
\end{eqnarray}
and $D_{V}\sim \beta^{-q/(q-1)}$.  We note in figure 8 that $\eta\sim 1.6$
for $q=1.5$ and $\eta\sim 0.85$
for $q=2.5$.  These values are close to $q/(q-1)=3$ for $q=1.5$ and
$q/(q-1)=1.66$ if $q=2.5$.  Therefore, the diffusion constant decreases
more slowly as a function of the inverse temperature $\beta$ as
$\beta^{-\eta}$ for the
generalized Monte Carlo scheme, though it decreases exponentially
faster for
the standard Monte Carlo scheme.  The exponent $\eta$ becomes smaller as
$q$ is increased and the diffusion is possible even at low temperatures.
This result can be anticipated based on data presented in figures 3
and 6.  This power-law temperature dependence of the diffusion constant
$D_{V}$ for $q\neq 1$ has been pointed out by Straub and
Andricioaei~\cite{SA2}.

\section{Concluding remarks}

In conclusion, we found that the so-called "quasi-ergodicity"
in a double well potential encountered in the standard Metropolis
Monte Carlo algorithm is largely circumvented by the
use of the generalized Monte Carlo algorithm of Andricioaei
and Straub~\cite{AS2,AS3} based on Tsallis's~\cite{TS} generalized
statistical mechanics. Therefore this algorithm will be useful
in the Monte Carlo study of any system, in particular, when
the first order phase transition occurs. We believe that
this algorithm will be useful, for example, in the study of
the melting transition of clusters for which the J-walk
algorithm~\cite{FFD} and the multiple histrogram
method~\cite{WA} have been utilized. It is also interesting
to apply this new algorithm to the problems where the slow
dynamics due to randomness or frustration is serious, for
example, the spin-glass problems.  This generalized
Monte Carlo algorithm based on Tsallis's statistical mechanics
use the reweighting technique similar to the multicanonical
method~\cite{BN}.  But this algorithm is much simpler than
the multicanonical methods and appears to be more promising
because it does not require an iterative determination of the
weighting factor by preliminary Monte Carlo runs.

Recently, a new formulation of the Tsallis statistics
has been proposed~\cite{TMP}.  With a new choice of energy
constraint, unfamiliar consequences in the previous formulation
have been erased.  However, the Lagrange parameter used
in the new formulation can be related to that in the previous
formulation. Therefore, in performing Monte Carlo simulations
with reweighting technique, we may use the previous formulation
of the Tsallis statistics.

%\newpage

\begin{flushleft}
{\Large \bf Acknowlegement}
\end{flushleft}

Authors are grateful to the reviewer for several useful comments
which were incorporated into the final version of the manuscript.
M. I. acknowledges support from the Hiroshima City University
Grant for Special Academic Research. He is also grateful
to the Department of Physics, Tokyo Metropolitan University,
where this work was initiated, for its hospitality to him.

%\end{flushleft}

\begin{thebibliography}{16}

\bibitem{VW} J. P. Valleau and S. G. Whittington, in:Statistical Mechanics,
edited by B. J. Berne (Plenum Press, New York, 1977), Ch. 4.
\bibitem{FFD} D. D. Frantz, D. L. Freeman and J. D. Doll, J. Chem. Phys.
{\bf 93} (1990)
2769.
\bibitem{Pe} T. J. P. Penna, Phys. Rev. E {\bf 51} (1995) R1.
\bibitem{TS} C. Tsallis and D. A. Stariolo, Physica A {\bf 233} (1996) 395.
\bibitem{AS1} I. Andricioaei and J. E. Straub, Phys. Rev. E {\bf 53}
(1996) R3055.
\bibitem{Ts} C. Tsallis, J. Stat. Phys. {\bf 52} (1988) 479.
\bibitem{KGV} S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, Science {\bf
220} (1983) 671.
\bibitem{AS2} I. Andricioaei and J. E. Straub, Physica A {\bf 247} (1997) 553.
\bibitem{AS3} I. Andricioaei and J. E. Straub, J. Chem. Phys. {\bf 107}
(1997) 9117.
\bibitem{SA} J. E. Straub and I. Andricioaei, Braz. J. Phys. {\bf 29}
(1999) 179.
\bibitem{CT} E. M. F. Curado and C. Tsallis, J. Phys. A {\bf 24} (1991) L69.
\bibitem{MRRTT} N. Metropolis, A. W. Rosenbluth, N. N. Rosenbluth, A. H. Teller
 and E. Teller, J. Chem. Phys. {\bf 21} (1953) 1087.
\bibitem{FS} A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. {\bf
61} (1988) 2635; R. H.
Swendsen, J. S. Wang, and A. M. Ferrenberg, The Monte Carlo Methods in
Condensed Matter
Physics, edited by K. Binder (Springer Verlag, Berlin, 1992), pp. 75-91.
\bibitem{BN} B. Berg and T. Neuhaus, Phys. Lett. B {\bf 267} (1991) 249.
\bibitem{AT} M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids
(Clarendon Press, Oxford, 1987), pp.121-123.
\bibitem{TMK} D. Thirumalai, R. D. Mountain and T. R. Kirkpatrick, Phys.
Rev. A {\bf 39} (1989) 3563.
\bibitem{SA2} J. E. Straub and I. Andricioaei, "Exploiting Tsallis
statistics,'' in: Algorithms for Macromolecular Modeling, Edited by P.
Deuflhard, J. Hermans, B. Leimkuhler, A. Mark, S. Reich and R. D. Skeel,
Lecture Notes in Computational Science and Engineering (Springer Verlag
Berlin, 1998), vol. 4, pp. 189-204.
\bibitem{WA} S. Weerasinghe and F. G. Amar, J. Chem. Phys. {\bf 98} (1993)
4967.
\bibitem{TMP} C. Tsallis, R.S. Mendes, and A.R. Plastino, Physica A {\bf
261} (1998) 534.
\end{thebibliography}

\newpage

%\section*{Figure Captions}

\begin{figure}[th]
\epsfxsize=8cm \centerline{\epsfbox{tfig1.eps}}
\caption{A classical double well
potential (16) (thin solid curve),for $\gamma=0.9$, the Tsallis weight
$p_{\rm T} \propto \exp(-\beta \overline{V}(x))$ with $q=1.5$
(dotted curve), $q=2.5$ (chain curve) and the Boltzmann
weight $p_{\rm B}=\exp(-\beta V(x))/Z_{1}$ (thick solid curve)
both at temperature $\beta=10$.  The Tsallis weights are
renormalized.}
\label{fig:1}
\end{figure}

\begin{figure}[th]
\epsfxsize=8cm \centerline{\epsfbox{tfig2.eps}}
\caption{A smoothed potential landscape $\overline{V}(x)$
when $q=1.5$ (dotted curve) and $q=2.5$ (chain curve) for
$\gamma=0.9$ and $\beta=10$, compared with the original
potential landscape $V(x)$ ($q=1$, thin solid curve).
The energy barrier is lowered by the smoothing defined
by (3) of Tsallis's statistical mechanics.}
\label{fig:2}
\end{figure}

\begin{figure}[th]
\epsfxsize=8cm \centerline{\epsfbox{tfig3.eps}}
\caption{The classical average potential energy $<V>$ for
the $\gamma=0.9$ potential as a function of $\beta$. The
thin solid curves are initialized at the
metastable $\alpha$-minimum $x_{0}=-0.9613$
with the standard Monte Carlo
algorithm with the parameter $q=1$ and the generalized Monte
Carlo algorithm with $q=1.5$, $q=2.0$ and $q=2.5$.
The exact result is represented by the thick solid
curve.  The results initialized at the global minimum $x_{0}=1$
and calculated from the Metropolis algorithm with $q=1$ and that from
the generalized Monte Carlo algorithm with $q=1.5$ are shown
by the dotted curve. All results are obtained from 100 independently
initialized Monte Carlo walks, each consisting of 500 warm-up steps
followed by $10^{4}$ steps with data accumulation.  Some
representative error bars are
also shown.}
\label{fig:3}
\end{figure}

\begin{figure}[th]
\epsfxsize=8cm \centerline{\epsfbox{tfig4a.eps}}
\vspace{4mm}
\epsfxsize=8cm \centerline{\epsfbox{tfig4b.eps}}
\vspace{4mm}
\epsfxsize=8cm \centerline{\epsfbox{tfig4c.eps}}
\caption{ Time series of the Monte Carlo walks generated by
(a) the standard Monte Calro ($q=1$), (b) the Tsallis Monte
Carlo with $q=1.5$, and (c) $q=2.5$ when $\beta=10$ and
$\gamma=0.9$.  All walks start from the metastable
$\alpha$-minimum, and 500 warm up steps are also included.
The number of barrier crossing is greater
for larger $q$.}
\label{fig:4}
\end{figure}

\begin{figure}[th]
\epsfxsize=8cm \centerline{\epsfbox{tfig5.eps}}
\caption{The average potential energy $<V>$ for the temperature
$\beta=10$ as a function of the asymmetry parameter $\gamma$ of
the double well potential. The solid and dotted curves are
used the same as in figure 3.  We only show the result and some
representative error bars for the cases $q=1$ and $q=1.5$
of figure 3.}
\label{fig:5}
\end{figure}

\begin{figure}[th]
   \def\halftext{.471\textwidth}
   \parbox{\halftext}{
      \epsfxsize=\halftext
      \epsfbox{tfig6a.eps}
   }
   \hspace{4mm}
   \parbox{\halftext}{
      \epsfxsize=\halftext
      \epsfbox{tfig6b.eps}
   }
\caption{The Visiting probabilities $P(x)$ obtained from
(a) the standard Monte Carlo walk compared with the
Boltzmann weight, and (b) the generalized Monte Carlo
walks compared with the Tsallis weight (c.f. figure 1).
Probability distribution $P(x)$ generated from the
generalized Monte Carlo scheme can recovers the
equilibrium Tsallis weight within short Monte Carlo
steps (a), while it cannot recover the equilibrium
Boltzmann weight until rather long ($\sim 10^{6}$) Monte
Carlo steps are executed (b).}
\vspace*{20mm}
\label{fig:6}
\end{figure}

\begin{figure}[th]
\epsfxsize=8cm \centerline{\epsfbox{tfig7.eps}}
\caption{The normalized inverse of the generalized ergodic
measure $d_{V}(0)/d_{V}(n)$ as a function of the number of
Monte Carlo steps $n$ calculated for $\gamma=0.9$
potential at $\beta=10$ using the standard Monte
Carlo scheme ($q=1$) and the generalized Monte Carlo
scheme ($q=1.5, 2.5$).  $d_{V}(0)/d_{V}(n)$ grows linearly
with the Monte Carlo steps $n$ and the diffusion is faster
for larger $q$.}
\label{fig:7}
\end{figure}

\begin{figure}[th]
   \def\halftext{.471\textwidth}
   \parbox{\halftext}{
      \epsfxsize=\halftext
      \epsfbox{tfig8a.eps}
   }
   \hspace{4mm}
   \parbox{\halftext}{
      \epsfxsize=\halftext
      \epsfbox{tfig8b.eps}
   }
\caption{The diffusion constant $D_{V}$ as a function of
temperature $\beta$.  (a) A straight line of the semi-log plot
of $D_{V}$ vs $\beta$ for the standard Monte Carlo scheme
indicates activation-type diffusion, while (b) straight
lines of the log-log plot for the Tsallis Monte Carlo scheme
in the high temperature region indicate slow power-law
decrease of the diffusion constant as the temperature is lowered (larger
$\beta$).}
\vspace*{20mm}
\label{fig:8}
\end{figure}

\newpage
\begin{table}[th]
\caption{The convergence of the Metropolis and the generalized
Tsallis Monte Carlo algorithms at $\beta=10$. The average energy
$<V>$, the standard deviation (STD) obtained from 100
independently initialized walks, each originating at the
metastable $\alpha$-minimum.}

\vspace{4mm}
\begin{tabular}{c|c|c}
\hline
&Metropolis Monte Carlo & Tsallis Monte Carlo ($q=1.5)$ \\
steps & $<V>$ (STD) & $<V>$ (STD) \\
\hline
$10^{2}$ & $0.1500$ ($0.0177$) &$0.1065$ ($0.0518$) \\
$10^{3}$ & $0.1520$ ($0.0049$) &$0.0937$ ($0.0394$) \\
$10^{4}$ & $0.1475$ ($0.0192$) &$0.0798$ ($0.0162$) \\
$10^{5}$ & $0.1052$ ($0.0320$) &$0.0798$ ($0.0049$) \\
$10^{6}$ & $0.0805$ ($0.0124$) &$0.0800$ ($0.0017$) \\
Exact & $0.0799$	 & $0.0799$ \\
\hline
\end{tabular}
\end{table}

\end{document}

max_stars_repo_name: e-sim/pdf-text-extraction-benchmark

max_issues_repo_name: e-sim/pdf-text-extraction-benchmark

max_forks_repo_name: e-sim/pdf-text-extraction-benchmark


text:
\subsection{Make targets}

\begin{description}

\item[\emcode{make clean}]

Deletes all files created by \ccode{make}, such as \ccode{.o} files
and executables; also deletes files normally created in editing the
source, such as emacs \ccode{~} files and \ccode{TAGS}.