%% Define fonts for use in code.
\def\prf{\bf}   %% PRimitive Font
\def\ppf{\sf}   %% Pcm  Primitive Font
\def\cmf{\it}   %% CoMment Font

\def\fig#1{Figure~\ref{fig:#1}}



\subsection{Protein Folding}

A second application that we have implemented is protein folding.
The reasons for choosing this application are similar to the reasons for
choosing ray-tracing. The problems are large enough to warrant the use of
parallelism. A common problem size takes over two hours when run
sequentially, and we would like to run a series of problems. Also, an
initial attempt to parallelize the program did not make efficient use of the
machine.  
This attempt statically broke the computation into subcomputations; but the
subcomputations were too coarse, and their run times too variable, to keep
all processors busy.  An implementation using PCM avoids this
problem. 


The work on this problem was done in conjunction with Vijay Pande of the
Center for Material Sciences and Engineering at MIT. In their work
\cite{vijay} they use the lattice model \cite{prot-lat-model} for protein
design. 
%
In this approach to the problem a protein is described as a chain of
monomers (each monomer may represent a number of amino acids).
%
%It is assumed that each monomer will sit on a lattice point, and that the
%polymer will take on a compact conformation. In other words this model
%assumes the polymer will fold such that it completely fills every point of a
%cube. 
%
It is assumed that each monomer will sit on a lattice point ({\it i.e.,} a
point on a 3D grid). Each possible folding of the polymer can then be
described as some path along the set of lattice points. 
\fig{pfoldcube} shows a polymer of length 26; each shade represents a
different type of monomer.
The model assumes the polymer will take on the most compact possible paths,
so it is only concerned with paths that completely fill some cube. 
%
In a folded polymer, two neighboring monomers will exert some attractive or
repulsive force on one another. This force depends on the types of the two
monomers. The energy of a folded polymer is the sum of the forces between
all of the neighboring monomers. Of course, this energy value depends
greatly on the way in which the polymer is folded.
For this work it is necessary to consider all possible foldings of a given
polymer and compute a histogram of the energy values.
This is what we have implemented using the PCM model on the CM-5.

\begin{figure}[htbp]
\centerline{\psfigure{pfoldcube}{4.0in}}
\caption{\label{fig:pfoldcube} A folded polymer}
\end{figure}


For the rest of this section we will be concerned mainly with the
implementation of this problem using PCM. For more details on the
algorithms used and the results obtained see \cite{pfold-cm5}.
%
At its heart, this program is a search program that finds all possible
unique paths in a cube; each path must traverse each point in the cube
exactly once. Each complete path represents one possible folding of the
polymer; so each time a path is found the energy value for that folding
is calculated and a result histogram is updated with the value.

This algorithm works by incrementally building up paths through the cube
until complete paths are reached.
The function {\it Count\_Entries} performs the core of the search. An outline of
the sequential code for this function is given in \fig{pf-original}.



\begin{figure}[htbp]
\begin{center}
\begin{code}
{\prf int} Count\_Entries(STATE st, {\prf int} point) \{ \\
\>  memcpy(state,st,sizeof(STATE)); \\
\>  add\_point\_to\_path(point,state); \hspace{.5in}   /** {\cmf add point to path} **/ \\
%%\>  /** {\cmf If we found a complete path update the result histogram} **/ \\
%%\>  /** {\cmf and return 1 (the number of paths found.)}    **/ \\
\>  {\prf if} (complete\_path(state))\{ \\
\>\>    update\_result\_histogram(state); \\
\>\>    {\prf return} 1 ; \\
\>  \} \\
\>  /** {\cmf If no complete paths are possible return 0} ***/ \\
\>  give\_up = f(state,point); \\
\>  {\prf if} (give\_up) {\prf return} 0; \\
\>  /** {\cmf Otherwise call Count\_Entries recursively on each neighbor} **/ \\
\>  {\prf else}\{ \\
\>\>    sum = 0; \\
\>\>    {\prf for}(i=0;i$<$num\_neighbors;i++)\{ \\
\>\>\>      next\_neighbor=neighbor[i]; \\
\>\>\>      {\prf if} (not\_occupied(next\_neighbor,state))\{ \\
\>\>\>\>        sum += count\_entries(state, next\_neighbor); \\
\>\>\>      \} \\
\>\>    \} \\
\> \} \\
\>  {\prf return} sum; \\
\}
\end{code} 
\caption{\label{fig:pf-original} Kernel of sequential protein folding}
\end{center}
\end{figure}


The first argument to this function is a STATE structure. This structure
defines the partial path that has been constructed so far. 
This structure contains information describing which points are occupied, the
type of monomer at each occupied point, and other data used to increase the
efficiency of the search. 
%This structure must tell whether a point is empty or occupied and give the
%type of the monomer at each occupied point. 
%In addition we also store information in this structure that will improve
%the efficiency of the algorithm. For example, we keep track of the number
%of empty neighbors around each point. This is useful in determining when we
%have reached a point where it will be impossible to find a complete path.
%This allows us to greatly prune the search.
%
The size of the STATE structure is on the order of 100 bytes.
The second argument to the function is {\it point}, the lattice point to be
added to the partial path.
This function returns an integer, namely the number of paths found.
 The function first adds {\it point} to the partial path. If this
completes the path it updates the result histogram with the energy value of
this new path and returns. Otherwise it applies some tests to see if it
impossible for this state to lead to a complete path. If a complete path is
impossible it prunes the search here. 
Otherwise it calls itself recursively for each empty neighbor of {\it
point}. At the end it returns the total number of complete paths found.
Typically we run the program on several polymers at a time. Each time we
find a complete path we just calculate several energy values. This
amortizes the time spent searching over several polymers.


We start the search by calling {\it Count\_Entries} on a set of starting
paths. The obvious way to start this search is to call {\it Count\_Entries}
once for each point in the cube, each time using a starting path consisting
only of that point. The problem with this approach is that it would count
duplicate paths. Two paths are considered duplicates if one path is a
rotation and/or reflection of the other. Instead, we have a second
(sequential) program which generates a set of starting paths that guarantees
that no generated paths will be duplicated. The initial routine simply
calls {\it Count\_Entries} on each of these starting paths.

\subsubsection{Protein Folding Parallelization}

An earlier attempt at parallelization was made without using PCM. In this
attempt the starting paths were statically divided between the nodes.
Each node then ran the sequential code for its subset of the starting paths.
At the end the result histograms on all the nodes are merged.
This was trivial to implement but gave poor results because of inadequate
load balancing. The amount of work involved for calculating with different
starting paths can differ by many orders of magnitude. 

To get a more efficient parallelization, the computation needed to be broken
into finer grains. PCM was ideal for this task. The procedure that makes
use of the PCM primitives is the {\it Count\_Entries} procedure. A skeleton of
the code for this procedure is given in \fig{pf-parallel}.
The major changes to this code to run it using PCM are similar to the
changes made to the ray-tracer. In particular, most of the recursive calls
to {\it Count\_Entries} are now made by creating and posting a closure. Also,
this version takes an additional argument {\it parent}. This is the
continuation to which the result will be sent.


%
%%In this version before making any recursive calls we first determine how
%%many recursive calls we will need to make. Since we want to sum up the
%%results of all the recursive calls we want to know how many items we will
%%have to sum over.\footnote{This is not strictly necessary, but is done for
%%performance reasons. If we could not easily determine the number of calls
%%in advance then each time we made a call we would create a new summation
%%closure which would add the result of this call to the sum of the remaining calls.}
%%This is also useful because it lets us perform a ``Fast Tail Call.'' This
%%is basically a ``fallthrough'' that 

\begin{figure}[tp]
\begin{center}
\begin{code} 
{\prf THREAD} Count\_Entries({\prf Cont} parent; {\prf int} point, STATE state,) \{  \\
\>  add\_point\_to\_path(point,state);    \\
\>  {\prf if} (complete\_path(state))\{   \\
\>\>    update\_result\_histogram(state);      \\
\>\>    {\ppf send\_argument}(1, parent);  \\
\>\>    return;     \\
\>  \}     \\
        \\
\>  /** {\cmf Determine number of neighbor nodes to be visited}  **/     \\
\>  give\_up = f(state,point);    \\
\>  num\_ntv = g(...);  \>\>\>\>  /** {\cmf num\_ntv = num of neighbors to visit} **/      \\
   \>\>\>\>\>      /** {\cmf nbrs\_to\_visit[i] = 'i'th neighbor to visit}  **/       \\
        \\
\>  /** {\cmf Case 0: If no  paths to search, then return 0 (no paths found)} **/  \\
\>  {\prf if} (give\_up $||$ (num\_ntv==0)) send\_argument(0, parent);     \\
        \\
\>  /** {\cmf Case 1:  If exactly 1 neighbor to visit --  only try that one } **/  \\
\>  {\prf else if} (num\_ntv==1) \\
\>\>      {\ppf tail\_call}(Count\_Entries,parent,nbrs\_to\_visit[0],state);       \\
\>  {\prf else}\{ \\
\>\>    /** {\cmf General case -- n neighbors to try  [n$>$1]}  **/    \\
\>\>    /** {\cmf create a closure to sum results of all sub-computations} **/      \\
\>\>    /** {\cmf post num\_ntv-1 threads and perform a tail call for the last} **/ \\
\>\>    sum\_closure = {\ppf make\_closure}(sum,parent,num\_ntv,\_=[num\_ntv]);    \\
\>\>    {\prf for}(i=0;i$<$(num\_ntv-1);i++)\{        \\
\>\>\>      next\_neighbor=nbrs\_to\_visit[i];        \\
\>\>\>      k1 = {\ppf make\_closure}(Count\_Entries,        \\
\>\>\>\>     \>\>        {\prf cont}\{sum\_closure,sum:val[i]\},  \\
\>\>\>\>     \>\>        next\_neighbor, new\_st=[sizeof (STATE)]);  \\
\>\>\>      memcpy(new\_st,st,sizeof(STATE)); \\
\>\>\>      {\ppf post}(k1); \\
\>\>\>    \}   \\
\>\>    /**{\cmf Perform a tail call for final neighbor} **/ \\
\>\>    new\_parent = {\prf Cont}\{sum\_closure,sum:val[num\_ntv-1]\};       \\
\>\>    {\ppf tail\_call}(Count\_Entries,new\_parent,next\_neighbor+1,state);        \\
\>  \}     \\
\}
\end{code}
\caption{\label{fig:pf-parallel} Kernel of parallel protein folding}
\end{center}
\end{figure}


There have been other changes to the code, however, most of
these changes were made for performance reasons rather than correctness
reasons. 
%%As before this program first adds {\it point} to the partial path and checks
%%if it has a complete path. 
The first difference is that this version determines in advance the
number of neighbors that will be visited. If there is just one neighbor
that needs to be visited, then exactly one recursive call needs to be made.
%%(and therefore no summation closure is necessary).
We do this by making use of a tail-call.
In this instance the tail call eliminates two overheads: first, the posting and
scheduling of the closure, and second, the copying of the state argument
into the new closure.


If there are more than one recursive call to be made then a summation
closure is created. This closure will add up the results of the
subcomputations and send the result to the parent of this thread. The
unusual syntax in the call to make\_closure ({\it i.e.,} `` \_ = [{\it
num\_ntv}]'') signifies that a specified number of empty slots (here {\it
num\_ntv}) should be left in the closure. These slots will be filled in
later with the results of the subcomputations. Most of the recursive calls
are made by making and posting closures. The final call makes use of a
tail call, for the same reasons given above.


\begin{figure}[t]
\centerline{\hspace{-0.6in}\psfigure{pf333}{2.8in} \hspace{0.3in} \psfigure{pf334}{2.8in}}
\caption{\label{fig:pf-perf} Run times for protein folding }
\end{figure}


\subsubsection{Protein Folding Results}

Many variations of this program have been run on a range of problem and
machine sizes. Results for two problem sizes are shown in \fig{pf-perf}.
The first figure shows a smaller problem size, namely a $3\times3\times3$
cube which has 103,346 paths. The second shows a larger size, namely a
$3\times3\times4$ cube which has over 48 million paths. In both figures the
horizontal line shows the execution time for the optimized sequential code.
The sloped line shows the execution time using PCM on various machine sizes.

The first observation is that the overhead added by the PCM model is
fairly small. The PCM program running on one processor was 14 and 22
percent slower than the sequential version in these
examples\footnote{Due to its large runtime, for the second example the
sequential version was run on a standard Sparc, not on one processor
of the CM-5. The runtime on a single CM-5 node is likely to be
slightly larger than this. This causes us to slightly overstate the
overhead.}.  This overhead is due both to posting and scheduling
closures and to the copying of the state structure into each closure.

The second observation is that the speedup curves are almost
perfectly linear for both of these problems.  (Even on the smallest
problems we ran we saw nearly linear speedups.  On 64 processors we
achieved a speedup of 60 and a run time of only 0.25 seconds.)  This
is because the granularity of the threads is small enough that the
workload can be evenly distributed.
%{\em mention effect of tail-calls ... }

%% ???? Mention 3x4x4 run ????


%
%In judging the efficiency of this program we look at two items, the overhead
%imposed when running on a single processor
%
%What is the overhead due to the PCM model, and how well does the program
%scale when the number of processors is increased.
%To judge the overhead we ran 
%
%



