![](media/image1.png){width="1.375in" height="0.74375in"}**SDM Center
Report**

**Jan 2002 - June 2006**

http://sdmcenter.lbl.gov

**[Introduction]{.underline}**

Managing scientific data has been identified as one of the most
important emerging needs by the scientific community because of the
sheer volume and increasing complexity of data being collected.
Effectively generating, managing, and analyzing this information
requires a comprehensive, end-to-end approach to data management that
encompasses all of the stages from the initial data acquisition to the
final analysis of the data. Fortunately, the data management problems
encountered by most scientific domains are common enough to be addressed
through shared technology solutions. Based on the community input, we
have identified three significant requirements. First, more efficient
access to storage systems is needed. In particular, parallel file system
improvements are needed to write and read large volumes of data without
slowing a simulation, analysis, or visualization engine. These processes
are complicated by the fact that scientific data are structured
differently for specific application domains, and are stored in
specialized file formats. Second, scientists require technologies to
facilitate better understanding of their data, in particular the ability
to effectively perform complex data analysis and searches over large
data sets. Specialized feature discovery and statistical analysis
techniques are needed before the data can be understood or visualized.
To facilitate efficient access it is necessary to keep track of the
location of the datasets, effectively manage storage resources, and
efficiently select subsets of the data. Finally, generating the data,
collecting and storing the results, data post-processing, and analysis
of results is a tedious, fragmented process. Tools for automation of
this process in a robust, tractable, and recoverable fashion are
required to enhance scientific exploration.

Our approach is to employ an evolutionary development and deployment
process: from research through prototypes to deployment and
infrastructure. Accordingly, we have organized our activities in three
layers that abstract the end-to-end data flow described above. We
labeled the layers (from bottom to top):

-   Storage Efficient Access (SEA)

-   Data Mining and Analysis (DMA)

-   Scientific Process Automation (SPA)

The SEA layer is immediately on top of hardware, operating systems, file
systems, and mass storage systems, and provides parallel data access
technology, and transparent access to archival storage. The DMA layer,
which builds on the functionality of the SEA layer, consists of
indexing, feature identification, and parallel statistical analysis
technology. The SPA layer, which is on top of the DMA layer, provides
the ability to compose scientific workflows from the components in the
DMA layer as well as application specific modules.

###### 

###### This report is organized according to the three layers: SEA, DMA, and SPA

###### 

###### *[1. Storage Efficient Access]{.underline}*

Today's scientific applications demand that high performance I/O be part
of their operating environment. These applications access datasets of
many gigabytes (GB) or terabytes, checkpoint frequently, and create
large volumes of visualization data. Such applications are hamstrung by
bottlenecks anywhere in the I/O path, including the storage hardware,
file system, low-level I/O middleware, application level interface, and
in some cases the mechanism used for Grid I/O access. This work
addresses inefficiencies in all the software layers by carefully
balancing the needs of scientists with implementations that allow the
expression and exploitation of parallelism in access patterns.

![](media/image2.wmf){width="3.498425196850394in"
height="2.09251968503937in"}Just above the I/O hardware in a
high-performance machine sits software known as the parallel file
system. This software maintains the directory hierarchy and manages file
data distribution across a large number of I/O components. Our **PVFS2
parallel file system** can provide multiple GB/second parallel access
rates, is freely available, and is in use at numerous academic,
laboratory, and industry sites, including the Massachusetts Institute of
Technology, the University of Utah's Center for High Performance
Computing, the Ohio Supercomputer Center, Argonne National Laboratory,
and Acxiom Corp.

PVFS2 operates on a wide variety of Linux platforms, making it perfect
for cluster systems, and it has also been ported to the IBM Blue Gene/L
supercomputer. PVFS2 incorporates key optimizations that prepare it for
ultra-scale systems, including scalability optimizations that closely
connect PVFS2 to MPI-IO (see Figure 1) \[LRT04\].

Above the parallel file system is software designed to aid applications
in more efficiently accessing the parallel file system. Implementations
of the MPI-IO interface are arguably the best example of this type of
software. MPI-IO provides optimizations that help map complex data
movement into efficient parallel file system operations. Our **ROMIO
MPI-IO interface implementation** is freely distributed and is the most
popular MPI-IO implementation for both clusters and a wide variety of
vendor platforms \[TGL99\]. Our previous work on PVFS and ROMIO forms a
solid technology base for new capabilities that address issues in
scientific computing in general and in specific communities such as
enhanced I/O request capabilities for PVFS and ROMIO and the Parallel
netCDF (Network Common Data Form) interface.

Some applications now operate on data located across the wide area
network (WAN). Unfortunately for scientists, most of the APIs that they
are accustomed to using to access data stored locally will not currently
operate on data stored remotely. Grid-enabling ROMIO allows existing
applications to access remote data sources with no significant work on
the part of scientists, significantly improving the usability of Grid
I/O resources. We have implemented prototypes that are capable of
accessing both Storage Resource Managers and Logistical Network data
sources through the MPI-IO interface.

MPI-IO is a powerful but low-level interface that operates in terms of
basic types, such as floating point numbers, stored at offsets in a
file. Scientific applications desire more structured formats that map
more closely to the structures applications use, such as
multidimensional datasets. File formats that include attributes of the
data, such as the input parameters and date of creation, and are
portable between platforms, are also desirable. The Hierarchical Data
Format (HDF5) interface, popular in the astrophysics community among
others, is one such high level API. HDF5 uses MPI-IO for parallel I/O
access as well.

![](media/image3.wmf){width="6.688188976377953in"
height="2.1228346456692915in"}NetCDF is another widely used API and
portable file format that is popular in the climate simulation and data
fusion communities. Our **Parallel netCDF** (PnetCDF) work provides a
new interface for accessing netCDF data sets in parallel. This new
parallel API closely mimics the original API, but is designed with
scalability in mind and is implemented on top of MPI-IO. Performance
evaluations using micro-benchmarks as well as application I/O kernels
have shown major scalability improvements over previous efforts
\[LLC+03\].

Figure 2 compares the performance results of the FLASH astrophysics I/O
benchmark using HDF5 and PnetCDF APIs, in which over 50% more I/O
bandwidth is observed for PnetCDF. The current release of PnetCDF is
version 1.0.1 and is available from
<http://www-unix.mcs.anl.gov/parallel-netcdf>. It has been tested on a
number of platforms, including Linux clusters, SGI Origin, Cray X1E, NEC
SX, IBM SP, and BlueGene/L. Language supports are provided for C and
Fortran. The self-test suite from Unidata netCDF package has been ported
to validate against single-process results using PnetCDF APIs. In
addition, a new parallel test suite has also been developed to validate
the results from a various parallel access patterns.

Current PnetCDF users come from several major research centers,
including:

-   The FLASH application team from ASCI/Alliances Center at the
    University of Chicago

-   The Center for Atmospheric Chemical Transport Models (LLNL)

-   The Scientific Data Technologies group at NCSA, developing I/O
    components for the Regional Ocean Model System (ROMS)

-   The ASPECT center (ORNL)

-   The Portable, Extensible Toolkit for Scientific Computation (PETSc)
    team (ANL)

-   The PRogram for Integrated Earth System Modeling (PRISM) at C&C
    Research Laboratories, NEC Europe Ltd.

-   The Earth System Modeling Framework (ESMF) at the National Center
    for Atmospheric Research (NCAR)

Our earlier works made considerable contribution in **caching techniques
for very large files and datasets** typically encountered in the area of
high performance data intensive computing. In particular we showed that
disk caching techniques, used for staging very large files from either
tape storage or other alternate replica locations, onto local disks for
direct file accesses during computations, have different characteristics
than the caching techniques in virtual memory paging and web caching. We
defined an appropriate metric termed *the average cost per reference*,
for comparing various cache replacement algorithms in this contest. We
show through analytical derivation that the best replacement policy is
one that evicts the candidate with the *least cost beneficial* function
~i~(t), at the instant in time t. This is referred to as the *LCB-K*
replacement policy. The formula for computing ~i~(t) for a replaceable
file i in cache was derived as,

$$$\varphi_{i}(t) = \frac{k_{i}(t)}{t - t_{- k_{i}}} \ast \frac{g_{i}(t) \ast c_{i}(t)}{s_{i}}$

where k~i~(t) denotes the number of references made to file i up to a
maximum of K, $t_{- k_{i}}$denotes the time of k^th^ backward reference
to file i, g~i~(t) is the cumulative references to the file over its
active period, and c~i~(t) is the estimated cost of caching the file if
not in disk. Results of a comparative study *LCB-K* with other methods
such as LRU, LRU-K, GDS (or Greedy Dual Size), etc., are shown in
Figures 3a and 3b for synthetic and real workloads respectively. The
details may be found in \[OOS02, OtSh03\]. Further work on caching
combined the replacement policies with file access scheduling for
optimal cache loading. The result of this work, referred to as
*OptCacheLoad*, combines optimal scheduling with cache replacement
policy. This produces excellent response times for data intensive jobs,
compared with first come first serve (FCFS) scheduling policy, under
limited cache sizes. Figures 4a and 4b also show these for synthetic and
real workloads respectively. The details are presented in \[ORS03\].

![](media/image4.png){width="3.001388888888889in"
height="2.0520833333333335in"}![](media/image5.png){width="3.0069444444444446in"
height="2.0680555555555555in"}

Figure 3: Comparison of average cost per retrieval for various cache
replacement algorithms

![](media/image6.png){width="3.0083333333333333in"
height="2.136111111111111in"}![](media/image7.png){width="3.011111111111111in"
height="2.1125in"}

Figure 4: Average response times of jobs serviced using OptCacheLoad and
LRU replacement

###### 

###### *[2. Data Mining and Analysis]{.underline}*

***Using FastBit indices for Data Analyses***

**Summary**

The most important piece of searching software that we have developed
under SDM center funding is the FastBit indexing package. FastBit is a
tool for searching large read-only datasets. It organizes user data in a
column-oriented structure that is efficient for on-line analytical
processing (OLAP), and utilizes compressed bitmap indices to further
speed up query processing. Analyses have proven that the compressed
bitmap index used in FastBit is theoretically optimal for
one-dimensional queries. Compared with other optimal indexing methods,
bitmap indices are superior because they can be efficiently combined to
answer multi-dimensional queries whereas other optimal methods cannot.
In this section, we briefly describe the searching capability of
FastBit, then highlight two applications that make extensive use of
FastBit, namely Grid Collector and DEX.

**Introduction**

It is a significant challenge to search for key insight in the huge
amount of data being produced by many data-intensive science
applications. For example, a high-energy physics experiment called STAR
is producing nearly a petabyte of data a year and has accumulated many
millions of files in last five years of operation. One of the core
missions of the STAR experiment is to verify the existence of a new
state of matter called the Quark Gluon Plasma (QGP) \[BGS+1999\]. An
effective strategy for this task is to find the high-energy collisions
that contain signatures unique to QGP, such as a phenomenon called jet
quenching. Among the hundreds of millions of collision events captured,
a very small fraction of them, maybe only a few hundreds contain clear
signatures of jet quenching. Efficiently identifying these events and
transferring the relevant data files to analysis programs are a great
challenge. Many data-intensive science applications are facing similar
challenges in searching their data.

In the past few years, we have been working on a set of strategies to
address this type of searching problem. Usually, the data to be searched
are read-only[^1]. Our approach takes advantage of this fact. Since most
database systems (DBMS) are built for frequently modified data, FastBit
can perform searching operations significantly faster than those DBMS.

Conceptually, most data can be thought of as tables, where each row of
the table represents an object or a record, and each column represents
one attribute of the record. To accommodate frequent changes in records,
a typical DBMS stores each record together on disk. This allows easy
update of the records, but in many operations, the DBMS effectively
reads all attributes from disk in order to access a few that are
relevant for a particular query. FastBit stores each attribute together
on disk, which allows one to easily access the relevant columns without
involving any other columns. Even though an update may take longer to
execute, but because a typical update operation usually add a large
number of objects, the extra overhead introduced by vertical
partitioning is negligible. In the database theory, separating out the
values of a particular attribute is referred to as a projection. For
this reason, using column-wise organized data to answer user queries is
also known as the projection index \[OQ1997\].

Each column of a table is a dimension of the data. Many scientific
datasets have tens or hundreds of dimensions; they are called
high-dimensional data. User queries usually involve conditions on
several attributes; they are known as multi-dimensional queries. To
answer multi-dimensional queries on high-dimensional data, it is well
known that the projection index performs better than most indexing
methods including B-Tree. Since FastBit uses column-wise organization
for user data, without any additional indices it is using the projection
index, which is already very efficient. Our indexing technology further
speeds up the searching operations. We have analyzed our bitmap index
and showed it to be optimal for one-dimensional queries \[WOS2002,
WOS2004, and WOS2006\]. Some of the best indexing methods including
B+-tree and B\*-tree have this optimality property as well. However,
bitmap indices are superior because they can be efficiently combined to
answer multi-dimensional queries.

**FastBit Indices**![](media/image8.jpeg){width="4.178472222222222in"
height="3.1347222222222224in"}

The key to the effectiveness of the FastBit searching software is a
special compression scheme called the Word-Aligned Hybrid (WAH) code.
This compression scheme enables FastBit to compress the bitmap indices
to modest sizes \[WOS2004, WOS2006\]. At the same time, it enables
common operations on the compressed bitmaps to be extremely efficient
\[WOS2002, WOS2004\]. The performance measurements and the analyses are
published in well-known conferences and journal in the database
community \[WOS2002, WOS2004, and WOS2006\].

Through extensive analyses \[WOS2004, WOS2006\], we are able to prove
that the time required by a WAH compressed index to answer a range query
is a linear function of the number of hits in the worst case. Through
the analyses, we also identified that the worst case is achieved with
uniform random data. Figure 1 plots the query response time against the
number of hits. It is easy to see that the time to answer range queries
on uniform random attribute is indeed a linear function of the number of
hits. For attributes that do not have uniform distribution, the query
response time is short than those need to answer queries involving
uniform random attributes. In Figure 1, we also plotted the time used by
the next most efficient compression scheme for bitmap indices, the
Byte-aligned Bitmap Code (BBC). It is clearly that WAH compressed
indices always use less time.

A very significant advantage of compressed bitmap index over commonly
used B-tree variants is that bitmap indices can be easily and
efficiently combined to answer multidimensional queries. Figure 2 shows
the time required by different indexing methods to answer some random
multi-dimensional range queries on a set of data from the STAR
experiment. In this figure, the horizontal axes are the query box sizes,
which are defined to be the fraction of rows that are expected to be the
hits. It is clear from this figure that WAH compressed bitmap indices
are usually ten times faster than the next most efficient method.

  ----------------------------------------------------- ------------------------------------------------------
  ![](media/image9.jpeg){width="2.9298611111111112in"   ![](media/image10.jpeg){width="2.9298611111111112in"
  height="2.19375in"}                                   height="2.19375in"}

  Figure 2: Time required to answer some random         
  multi-dimensional range queries using different       
  indexing methods.                                     
  ----------------------------------------------------- ------------------------------------------------------

**Grid Collector**

Grid Collector is a client-server system that plugs into the STAR
analysis framework to provide a new ways for users to select events of
interest for analyses. Figure 3 shows the overall architecture of the
system. It uses the FastBit software to implement an Event Catalogue for
all events produced in the STAR experiment. It also enables analyses of
archived data and remote data. This capability of performing large-scale
analysis of data across the data Grid was recognized with a Best Paper
Award from International Supercomputer Conference 2005 \[WGL+2005\].

**DEX**

Our work on feature tracking in combustion simulation data \[WKS+2003\]
has led to the development of a set of Connected Component Labelling
algorithms for data from regular meshes as well as image data. This has
further led into fruitful collaborations with the Visualization group at
LBNL combining VTK and FastBit techlogies \[SSBW2005, SSWB2005\], and
applied to efficient network traffic analysis \[Sto+05\]. DEX is
shorthand for Dexterous Data Explorer. It is a visualization system for
scientific data produced on regular meshes. Figure 4 shows two sample
images produced by DEX. This work was recognized no only by the database
community \[SSBW2005\] by also by the visualization community
\[SSWB2005\].![](media/image11.png){width="4.29375in" height="1.0in"}

  ----------------------------------------------------- -----------------------------------------------------
  ![](media/image12.png){width="2.4541666666666666in"   ![](media/image13.png){width="2.6979166666666665in"
  height="2.598611111111111in"}                         height="2.4743055555555555in"}

  Exploding supernova                                   Methane jet flame
  ----------------------------------------------------- -----------------------------------------------------

Figure 4: Two sample images produced by DEX on simulation data.

### ***High Performance Statistical Computing***

###### Motivation

Terascale computing enables simulations of complex natural phenomena, on
a scale not possible just a few years ago. With this opportunity, comes
a new problem -- the massive quantities of data produced by these
simulations. However, answers to fundamental questions about the nature
of the universe remain largely hidden in these data. The goal of this
work is to provide a scalable high performance statistical data analysis
framework to help simulation scientists perform interactive analyses of
these raw data to extract knowledge. To achieve this goal we developed
an open source parallel statistical analysis package, called *Parallel
R*, that lets scientists employ a wide range of statistical analysis
routines on high performance shared and distributed memory architectures
without having to deal with the intricacies of parallelizing these
routines.

***"Parallel R" for Parallel Statistical Analysis with R***

Present data analysis tools such as Matlab, IDL, and R, even though
highly advanced in providing various statistical analysis capabilities,
are not apt to handle large data-sets. Most of the researchers' time is
spent on addressing data preparation and management needs of their
analyses. Our goal has been to scale up the existing data analysis
software to meet the requirements posed by application scientists. Under
*Parallel-R* project \[S05\], we have tried to address the
aforementioned issues using an open source *R* statistical package
\[VSR05\]. As with any other widely used data analysis software, *R* is
also not designed to work with huge datasets and performs poorly in
handling the same. With *Parallel-R*, we have identified data-parallel
and task-parallel strategies through which the existing data-analysis
system can be infused with certain parallel computation techniques to
meet some of the challenges in analysis of huge data-sets. Under
*Parallel-R* project we have designed, developed, benchmarked and
released two add-on software libraries for R, namely *RScaLAPACK* and
*taskPR*. These packages are currently distributed across 24 countries
through R's CRAN network. The RPMs for the same packages have been built
and distributed by RedHat, SUSE, and Quantian Linux distributions.

**Enabling Data-Parallel Statistical Analysis with RScaLAPACK.** With
*RScaLAPACK* we have replaced serial linear algebra problem solving
routines that used LAPACK library, with its corresponding parallel
algorithms provided by the ScaLAPACK library. A single *RScaLAPACK*
function call carries out a complete cycle of data-parallel computation
to obtain the solution for a requested linear algebra problem.
Furthermore the function interface remains similar to its serial
predecessor (see Accomplishments in section 3.2.3). In addition, the
performance benchmarks carried out on the *RScaLAPACK* functions
revealed that they are scalable in terms of both the problem size and
the number of processors. A consistent increase in the speedup with the
increase in the number of processors used in *RScaLAPACK* function
execution can be observed; further increase in the speed up with the
increase in the input data size was also observed (Figure 1). With three
update releases following the original release of *RScaLAPACK*, we have
tuned the supported functions to perform well in diverse environments. A
technical paper discussing the *RScaLAPACK*'s architecture and related
performance gains achieved was published, and presented at the 18^th^
International Conference on Parallel and Distributed Computing, 2005
\[SYBK05\].![](media/image14.wmf){width="3.3631944444444444in"
height="2.158333333333333in"}

**Providing a Middleware for Task-Parallel Statistical Analysis with
taskPR.** With *taskPR* package we identified a means using which the
existing serial code written in R, could be parallelized. The *taskPR*
system provides a mechanism to detect the out-of-order independent R
tasks and intelligently delegates each task to the remote worker
processes, to evaluate them in parallel.

Figure 2, depicting the snapshot of a serial R code in biology
demonstrates the ease with which a serial code could be parallelized.
Further, a significant speedup was recorded for these parallelized
protein abundance ration calculations with *taskPR* (Figure 3).

***Applying Parallel R across Various Scientific Applications***

**Novel Feature Extraction for Quantitative High-Throughput
Proteomics.** The cellular activities of biological organisms are
carried out by a complex network of thousands of proteins. To
characterize this network, biologists can perturb an organism's growth
with genetic or environmental treatments and then observe how the
organism responses with the adjustment of its different proteins'
abundance. Quantitative proteomics is an emerging high throughput
technology that enables biologists to measure the abundance change of
thousands of proteins in the treatment organisms relative to the
reference organisms. The results of quantitative proteomics would give a
global view of the protein network's response to a treatment. However,
the quantitative proteomics has posed unprecedented informatics
challenges: massive and complex data management system and development
of novel algorithms for very noisy data.

The mass spectral data are generated in heterogeneous propriety data
formats, such as Xcalibur, MassLynx and Analyst. The raw data in
different formats are then converted into an open mass spectral data
format, mzXML. The mass spectral data, combined with the genome
sequences, are used for peptide identification and protein
identification. With the identification results, the mass spectral data
are then used for peptide quantification and protein quantification.

We have collaborated with Dr. Robert Hettich's GTL proteomics team to
automate the data processing pipeline and integrated all modules for
quantification (shown in blue in Figure 4) into the ***ProRata***
program \[PHM+05, PKT+05, PKT+06, PKM+06\]. The data visualization is
especially important for quantitative proteomics, as biologists often
need to manually explore the data and verify the results. *ProRata* is
equipped with a sophisticated interactive graphical user interface to
provide the data visualization capability (Figure 5). The two core
modules in the ProRata program, peptide quantification and protein
quantification, are powered by the novel algorithms that we developed
for processing the noisy data. Principal component analysis (PCA) is
used for peptide quantification. The first principal component captures
the signal component in the data and the second principal component
captures the noise component. The PCA algorithm not only de-noises the
data, but also gives an estimation of the signal-to-noise ratio (S/N) of
the data. The estimated relative abundance for peptides is then combined
for protein quantification with maximum likelihood estimation (MLE). By
probabilistically weighing the peptides with their S/Ns, the MLE
algorithm gives a more accurate point estimation of the protein relative
abundance than the averaging method used in previous programs. More
importantly, the MLE algorithm also estimates the confidence interval
for the proteins' relative abundance. The confidence interval of a
protein's relative abundance captures the quantification precision for
the protein. The confidence interval estimation is especially important
in proteomics, in that the thousands of proteins are quantified at
vastly different precisions.

![](media/image2.png){width="2.9458333333333333in" height="2.1in"}

Managing proteomics data is a huge challenge especially in the context
of keeping pace with the high-throughput data generation rate of modern
mass spectrometers. In ORNL GTL mass spectrometry group, two gigabytes
of data are generated per day per instrument and there are four
instruments running seven days per week. A forthcoming DOE GTL
Proteomics Facility will have hundreds of mass spectrometry devices.
Some experiments run for 1 day and some other experiments run for 10
days, depending on the number of experimental conditions and the number
of replicates. For example, with two biological replicates and two
technical replicates (typical run), there will be four days of
instrument time, which would generate about 8 Gb data. There are
usually 12 raw files in a MudPIT experiment. But as the data pass
through the pipeline, the type of files and the number of files changes.
The SEQUEST generates million s of OUT files per 24 hours, which are
combined into 1 file by DTASelect. ProRata generates thousands of
xic.xml files, which are also processed to 1 file. The R implementation
of the ProRata's algorithms has been parallelized with *taskPR* has
greatly enhanced the data processing throughput (Figure 3).

**Parallel Analyses of Climate Data.** We are collaborating with the
SciDAC ORNL climate team including Dr. John Drake, Dr. Marcia
Branstetter and Dr. Auroop Ganguly to provide them with a capability to
efficiently handle terabyte climate simulation data spread across
multiple netCDF files. Currently, due to the size and complexity of the
data sets, even performing a simple statistical procedure on the
selected data involves laborious work and poses a great challenge.

In performing uni-variate and multi-variate dependence analysis in
climate extremes, they work with the climate data spanning 230 years
starting from year 1870 to 2100. The climate extremes refer to
extremities in the climate variables like, heat waves, cold waves,
rainfall etc. By studying the extremes researchers understand the global
climate pattern. Extreme value statistics that quantifies the stochastic
behavior of a process at unusually large (or small) values is used for
this purpose.  Currently, NCAR\'s Extremes Toolkit (extRemes) \[GKY05\],
an interactive program for analyzing extreme value data using the R
statistical programming language, is used by the climate researchers for
their analysis purposes.

One of the problems that the climate researchers at ORNL face was the
interactive selection and extraction of a required set of data spanning
over thousands (18,000) of netCDF files. Each file measuring  about 393
MB contain data for around 64 different climate data variables like
*precipitation rate* (PRECC), each such variable is dependent on one or
more of over 10 dimension values like, *latitude* (lat), *longitude*
(lon) etc. Further, each file contains data for 20 successive time-steps
with each having a periodic difference of 6 hours. At the initial stage
of our collaboration, we have been able to apply parallel R libraries to
parallelize the process of reading specific data of interest from
multiple netCDF files in parallel, given the details of the dataset that
needs to be extracted using the API as shown in Figure 6.

 \
The linear speed-up performance gains were observed in reading the data
in serial as compared to reading them in parallel using Parallel-R
libraries.

![](media/image3.jpeg){width="2.739763779527559in"
height="2.2133858267716535in"}**Bringing Parallel R to Geographical
Information Science (GIS).** We worked closely with Dr. Budhendra
Badhuri, Dr. George Fann and Dr. John Drake at ORNL to improve their
data analysis and visualization workflow. To address their large scale
data analysis needs we developed a software bridge to integrate
ParallelR with GRASS GIS application tool (Figure 7). GRASS (Geographic
Resources Analysis Support System) is a raster/vector GIS, image
processing system, and graphics production. It contains over 350
programs and tools to render maps and images on monitor and paper;
manipulate raster, vector, and sites data; process multi spectral image
data; create, manage, and store spatial data. However, it was lacking
advanced statistical data analysis capabilities that were required to
analyze GIS and climate data. The integration of GRASS with Parallel R
brought both advanced statistical analysis for large-data and image
rendering capabilities in a single environment.

###### *Feature selection in scientific applications*

###### Understanding changes in the global climate is challenging, since simulated and observed data include signals from many sources. To make meaningful comparisons between climate models, and to understand human effects on global climate, we need to isolate the effects of these different sources. For example, volcano eruptions and El Nino and Southern Oscillation (ENSO) variations both influence global temperatures, but recent volcanic eruptions coincided temporally with large ENSO events, making it difficult to separate the volcano effect from the ENSO effect. We collaborated with Dr. B. D. Santer from the Program for Climate Modeling and Intercomparison (PCMDI) at LLNL to investigate solutions that combine the existing techniques of principal component analysis (PCA) and independent component analysis (ICA). Neither method in isolation solved our problem adequately, but we obtained much improved results by first using PCA to reduce the dimension of the data, and then applying the more compute intensive ICA to the reduced data. Figure 1 displays three independent component basis functions obtained by applying our method to observational data from the National Center for Environmental Prediction. We determined the bases to represent the ENSO, the North Atlantic Oscillation (NAO), and the volcano components in the data. Figure 2 shows in red the monthly time series component estimates corresponding to the three bases presented in Figure 1. The blue curves display the ENSO 3.4 index, the NAO index, and an idealized volcano signal signature, correspondingly. The excellent match between the two series in the top panel of Figure 2 indicates that our approach isolated successfully the ENSO component. The agreements between the estimated NAO and volcano components, and their corresponding idealized signatures, are not quite as high, but they are scientifically significant nonetheless. 

###### ![](media/image4.png){width="4.440972222222222in" height="1.5875in"}

###### Figure 1. The three ICA bases representing the left) ENSO, middle) NAO, and right) volcano components.

###### ![](media/image5.png){width="6.504166666666666in" height="1.8041666666666667in"}

###### Figure 2. The three monthly time series corresponding to the three bases shown in Figure 1: the top) ENSO, middle) NAO, and bottom) volcano components. Red lines are our estimates. Blue lines indicate the ENSO 3.4 index, the NAO index, and an idealized volcano signal, respectively.

###### Following the success of our work with climate data, we were approached by scientists at General Atomics analyzing the large volumes of data being collected from the DIII-D Tokamak. They are interested in determining which of the numerous parameters being measured are relevant to the presence of edge harmonic oscillations (EHOs). EHOs are associated with a "Quiescent H-mode" of operation. H-mode operation is the choice for next-generation tokamaks, as it offers superior energy confinement, but it comes at a significant cost due to effects of edge localized modes (ELMs), which could cause rapid erosion of divertor components in tokamaks. During quiescent operation there are no ELMs. EHOs appear to provide an enhanced particle transport at the edge of the plasma that provides the needed density control. 

###### Currently, EHOs are identified mostly by visual means using the Fourier spectrum of the data from a magnetic probe. If an experiment seems to contain EHOs, data from other sensors (plasma velocity and the distance between the plasma edge and the tokamak wall) are consulted to verify the existence of the EHO. A summer student developed a program that implements the rules that the scientists have derived from their visual observations to identify EHOs. The program uses a sliding time window to analyze the data and assigns an \"EHOness\" value to each time window. The program seems to perform satisfactorily, but the results cannot be used to explain why or when EHOs appear.

###### Our goal was to help the scientists identify which of the numerous measurements is related to the occurrence of EHOs in their data. Each experiment lasts for about six seconds and data from numerous sensors are recorded and stored in a combination of MDSPlus and home-grown ptdata files. We collaborated with Dr. Keith H. Burrell of General Atomics to combine domain knowledge with our feature selection techniques to identify the key features. Dr. Burrell identified approximately 40 candidate variables, which we extracted from the database into text files with the format used by our tools. We faced some challenges in this task. (1) The only existing interface to the database is for IDL (Interactive Data Language), so existing IDL programs were modified to output data in the format we use. (2) The data from multiple sensors is not sampled at the same rate or may start or end at different points in time, and thus the data must be registered. (3) Some of the variables of interest have to be calculated from the raw data.

######  We extracted the features that describe approximately 500 experiments that have been analyzed (visually) by Dr. Burrell. His analysis will be useful to validate the results. We also did some additional preprocessing to remove entries with missing values, with feature values which appeared to be outliers, and to remove possibly irrelevant entries near the beginning and the end of experiments where some feature values appear very noisy. We then applied several dimension reduction techniques (partially funded by SciDAC-1) and found that we were able to identify EHOs with a low error rate using just a few of the features (See Fig. 3). Thus, a physicist can focus on a few variables instead of the original 40, in their attempts to understand EHOs. 

###### The software for dimension reduction and for anomaly detection in streaming data has since been licensed to General Atomics for further use in their work. In addition, our work on the DIII-D data caught the attention of fusion scientists at Princeton Plasma Physics Laboratory (PPPL), which led to an on-going collaboration on the analysis of Poincaré plots and the research and development of software for blob tracking. 

###### ![](media/image6.jpeg){width="3.0833333333333335in" height="2.057638888888889in"}

###### 

######  Figure 3. The performance of different dimension reduction techniques, with the features on the x-axis, and the error rate for the identification of EHOs on the y-axis .

###### Details of our analyses and publications are at [**www.llnl.gov/casc/sapphire/pubs.html**](http://www.llnl.gov/casc/sapphire/pubs.html). 

###### *[3. Scientific Process Automation (SPA)]{.underline}* 

###### [3.1 Background and Motivation]{.underline}

Science in many disciplines increasingly requires data-intensive and
compute-intensive information technology (IT) solutions for scientific
discovery. Scientific applications with these requirements range from
the understanding of biological processes at the sub-cellular level of
"molecular machines" (as is common, e.g., in genomics and proteomics),
to the level of simulating nuclear fusion reactions and supernova
explosions in astrophysics. A practical bottleneck for more effective
use of available computational and data resources is often in the IT
knowledge of the end-user; in the design of resource access and use of
processes; and the corresponding execution environments, i.e., in the
*scientific workflow environment* of end user scientists. The goal of
the Kepler/SPA thrust of the SDM Center is to provide solutions and
products for effective and efficient modeling, design and execution of
scientific workflows.

The scientific community now expects on-demand access to services which
are not bound to a fixed location (such as a specific lab) and fixed
resources (e.g., particular operating system), but which can be used
over (mobile) personal access device of a scientist (e.g., laptop or a
PDA or a cell phone) \[NSF03, DOE04\]. The need is for environments
where **application workflows are dramatically easier to develop and
use** (from *art* to *commodity*), thus expanding the feasible scope of
applications possible within budget and organizational constraints, and
**shifting the scientist's effort away from data management and
application development to concentrating it on scientific research and
discovery** (Figure 1)**[.]{.underline}**

![](media/image7.png){width="3.0in"
height="1.8444444444444446in"}***Reusability*** (elements can be re-used
in other workflows), ***substitutability*** (alternative implementations
are easy to insert, very precisely specified interfaces are available,
run-time component replacement mechanisms exist, there is ability to
verify and validate substitutions, etc), ***extensibility*** (ability to
readily extend system component pool, increase capabilities of
individual components, have an extensible architecture that can
automatically discover new functionalities and resources, etc),
***customizability*** (ability to customize generic features to the
needs of a particular scientific domain and problem), and
***composability*** (easy construction of more complex functional
solutions using basic components, reasoning about such compositions,
etc.)

###### [3.2 Accomplishments]{.underline}

As part of original SciDAC SDM Center, the Scientific Process Automation
(SPA) thrust has made several significant contributions to improve
scientific workflow technology, leading to the adoption of scientific
workflows within the DOE community:

-   **Co-founded the *Kepler* project**. Based on extensive evaluation
    of multiple, competing technologies (e.g., \[Cha02\], \[ABB+03\],
    \[LAG+03\]), we co-founded the Kepler project \[Kep06, LAB+05\]; a
    multi-site open source effort to extend the Ptolemy system \[Pto06\]
    and create an integrated scientific workflow infrastructure. We have
    also started to incorporate the SCIRun \[SCI06\] system into our
    environment, to better handle the needs of tightly-coupled and
    visualization-intensive workflows.

-   **Deployed workflows in multiple scientific domains.** We have
    worked closely with application scientists to design, implement, and
    deploy workflows that address their real-world needs (e.g.,
    \[Cha02\], \[PAL+03\], \[MVB+04\], \[Bha05\], \[KBB+05\]). In
    particular, we have active users on the SciDAC Terascale Supernova
    Initiative (TSI) team and within the LLNL Biotechnology Directorate.
    The success of these deployments has resulted in the SPA team being
    approached to participate in several SAPs.

-   **Developed specialized actors**. We designed and implemented many
    of the actors commonly needed by scientists. For example a web
    service actor, which supports invocation of distributed web serves
    and a secure shell actor supporting remote task execution.

-   **Extended existing capabilities**. We have added asynchronous
    service-based process status tracking, as well as a simple
    provenance tracking capability to the Kepler infrastructure. While
    significant work is required to fully-implement these capabilities
    for general workflows, the current implementations meet our
    customers' near-term requirements.

-   **Initial deployment infrastructure**. We have developed a web-based
    deployment infrastructure for distributing our system, including
    updates and patches, to users (available on-line at
    <https://www-casc.llnl.gov/sdm/>). This infrastructure includes
    documentation automatically extracted from the actor class
    definition.

Over the first three years of the project, the SPA team (Georgia Tech,
LLNL, NCSU, and SDSC) has focused on engineering and developing a
loosely coupled scientific workflow environment that uses data and
computational resources accessible both locally and remotely over the
network. In the last two years, we (LLNL, NCSU, SDSC, UC Davis and U of
Utah) worked on integrating loosely and tightly coupled workflows, and
understanding and implementing production level workflows for a number
of users. The ***current*** version of the SPA/Kepler workflow
construction and support system provides a range of capabilities, from
those that support ***low to medium intensity dataflows*** (e.g.
bioinformatics workflows) to those that support ***large data-intensive
workflows*** (e.g., large-scale astrophysics and fusion simulations).
Significant progress was possible over the last five years, because we
heavily leveraged existing technology instead of "re-inventing the
wheel" where that was not necessary. We initially conducted a
comprehensive survey of existing technologies and open source tools, and
developed a few small evaluation prototypes of the possible solutions.
This allowed us to mitigate the risks and make crucial decisions in an
informed way and on time. We have now a very vibrant development effort,
not least because we have adopted an excellent workflow construction and
execution engine code base (the Ptolemy II system[^2] from UC Berkeley)
that was designed to be **extensible** with **reusable**,
**substitutable, customizable** and **composable** components (actors),
and was provably sufficiently reliable and maintainable. We then built
our own extension libraries for the system and co-founded Kepler^3^, and
open source workflow support environment that has since been adopted by
a fairly broad scientific workflow community. Since then, the Kepler
community has grown into a vibrant new source of innovation, extension
and solutions that directly benefit our projects. We were very
successful in selecting the right pilot scenario (promoter
identification workflow discussed below) and a very capable and willing
collaborating domain scientist (Dr. Coleman). In addition, through a
collaborative effort within the Kepler project, we have successfully
brought in external developers in a cross-project open source
development activity.[^3]

Detailed information about this project and its results can be found in
the publications and reports produced (see publications list).

**3.2.1 SPA/Kepler Workflow Framework and Environment (Ptolemy
II/Kepler)**

The whole workflow process can be encapsulated in the (currently
Ptolemy-based) workflow framework to allow easy task
definition/description, process initiation, progress tracking,
collection of process metrics, and process repetition -- all from a
variety of end-user platforms (Windows, Unix, MacOS). The method for
this formal recording of the process is available, and when applied to
TSI or Fusion Simulation Project (FSP) or other workflows it eases
exchange of the workflow information, as well as its modification,
understanding, and verification and validation. This work required
development of a library of scientific SPA "actors" that - when designed
properly generalize and are thus re-usable for definition and management
of other scientific workflows. Figures below show the design view of the
several implementations operational workflows. The SPA prototypes and
production workflows were demonstrated at a number of occasions,
including the *Supercomputing Conferences 2002,2003, 2004 and 2005*.

One of the hallmarks of this project is a rich library of Kepler/Ptolemy
actors -- both general and domain specific system components and
extensions developed by SPA -- these include:

-   **custom actors** for a number of bioinformatics and TSI web
    services and provenance;

-   a **generic web service actor** for seamless plug-in of any standard
    web service \[AJL+04\];

-   a **web service harvester** to the system, which, after pointing the
    system to a web page or repository of web services, ingests all
    available service operations at once and automatically **creates a
    library of web services** to be used by the scientist.

-   **data transformation and data movement actors** were added to the
    system to handle many different data-type sources (e.g., \[Oot04\],
    \[OVV03a\]);

-   a **browser/user-interaction actor** was developed that allows the
    workflow designer to include standard browser inputs (e.g., forms,
    but also dynamic Javascript features, clickable images etc.) and
    outputs (e.g., HTML tables) within the user's workflow;

-   a **command-line actor** (for interaction with local platform
    commands), and a set of **remote access actors** (e.g., ssh actor);

-   Additional enhancements to the system were made in collaboration
    with external contributors to the Kepler initiative: A **generic
    design actor** (a SEEK contribution) allows a scientist to develop
    and document a design of a scientific workflow without necessarily
    having an implementation for the individual tasks/actors. SEEK
    colleagues also provided a 10-times faster build environment for
    Kepler as compared to the standard Ptolemy II build. A **database
    source actor** (a GEON contribution) allows connecting and querying
    any JDBC data source as part of the scientific workflow; and so on.

**3.2.2 Workflows**

Over the last two to three years we leveraged our experience by
partnering with several additional and very successful SciDAC projects
that helped us extend our system to data-intensive tasks and to
tightly-coupled workflows. Specifically, additions to the project were
interactions with Terascale Supernova Initiative (TSI), Fusion
Simulation Projects (FSP), and SCIRun. In addition to end-users from TSI
(Drs. Blondin and Swesty) and FSP (Dr. Klasky), the project welcomed a
member from the SCIRun[^4] team. This added to our experience with
large-scale tightly-coupled scientific data-flows as well as integration
with applications available through CCA. Active cooperation with the
TSI[^5] team was especially fruitful since their workflows served not
only as a template for high-intensity interactions but also to recognize
and explore next-generation issues and improvements such as detachable
workflows, fault-tolerance, process tracking and provenance, workflow
security, and service-based asynchronous solutions.

Figure 2 shows the "science driver" that was used in the first phase of
the project and its translation into Kepler framework: A starting point
for discovery is to link genomic biology approaches such as microarrays
with bioinformatics to identify and characterize eukaryotic promoters.
We call this the promoter identification workflow or PIW. To clearly
identify co-regulated groups of genes, high throughput computational
molecular biology tools are first needed that are scalable for carrying
out a variety of tasks such as identifying DNA sequences of interest,
comparison of DNA sequences, and identification of transcription factor
binding sites, etc. Some of these steps can be executed by querying
web-accessible databases and computation resources. However, using web
sources \"as-is\" to enact scientific workflows required many manual and
thus time-consuming and error-prone steps. It was thus desirable to
automate scientific workflows such as the PIW as much as possible. A
number of information technology and database challenges have to be
overcome. In the initial stages of the project, we had intensive
interactions with our domain science drivers (Dr. Matt Coleman) to
understand the requirements of the specific PIW workflow. Only after
having built an initial "custom prototype" did we have sufficient
information for the detailed design of the actual workflow. We also
investigated the available workflow technologies and user environments
for design and executing scientific applications, and eventually settled
for the Ptolemy II system (see SPA publications for details).

![](media/image8.png){width="4.833333333333333in"
height="2.7493055555555554in"}[^6] and FSP, and built first custom, then
more general workflow support solutions. Results have been excellent and
very well received. Figures 3 and 4 illustrate respectively translation
of TSI (Dr. Swesty) and FSP (Dr. Klasky) workflows into Kepler. Of
course, behind each of the steps shown in the figures are a number of
more sophisticated and details processes that perform actual
computations, move, aggregate, and transform the data, and deliver the
final end-user product. In the context of TSI, basic scientific
challenges center around the development of a standard model of
supernovae core collapse and related processes. In the case of FSP the
science is about understanding efficiency and stability of fusion
solutions and processes.

![](media/image9.png){width="3.125in"
height="2.4916666666666667in"}![](media/image10.png){width="3.5in"
height="2.401388888888889in"}

Figures 3 and 4 -- TSI (Swesty) workflow (left) and pilot for an FSP
workflow (Klasky).

Underlying challenges related to simulations, data analysis and data
manipulation include scalable parallel numerical algorithms for solution
of large, often sparse linear systems, flow equations, and large
eigen-value problems, running of simulation on supercomputers, movement
of large amounts of data over large distances, collaborative
visualization and computational steering, and collection of appropriate
process and simulation related status and provenance information. This
requires interdisciplinary teams of astrophysicists, nuclear physicists,
fusion physicist, applied mathematicians, computer scientists, etc. The
general underlying "template" (and the potential role-model for future
workflow construction and management "wizards") is amazingly similar:
large-scale parallel computations and steering (hundreds of processors,
gigabytes of memory, hours to weeks of CPU time), data-movement and
reduction (terabytes of data), visualization and analytics (interactive,
retrospective, and auditable). An abstraction of this and its Kepler
translation are illustrated in Figures 5 and 6 (Blondin workflows). In
creating successful abstractions of the TSI and FSP workflows we have
identified a number of next generation challenges. Some of these
challenges were collected in recent publications. (\[ABB+05, AKA05\].

###### [3.3 Next Set of Challenges]{.underline}

Significant work still remains to be done before scientific workflow
technology reaches 'utility' state of ease of use expected by the
majority of application scientists. The challenge is to generalize the
lessons learned so far, and to extend the framework support for whole
scientific discovery domains. The work we are proposing will make
another set of significant advances in this direction by developing
additional automation that will yield flexible, dependable, extensible,
customizable, scalable, and (re-)usable workflows for application
scientists. While the focus is on workflows that require
high-performance computing (HPC) support, the technology is applicable
across the spectrum, from hand-held PDAs to HPC.

For example, end-users need to trust the workflow solution before they
would depend on it. This requires a high confidence level that the
"mechanics" of their workflows will operate as they expect, it requires
some run-time indicators that this is the case, it requires some
information about the history of the computations and results, (for
example Figures 7 to 10) and it requires an environment where workflow
modifications are easy, and operation of workflows can happen in both
synchronous and asynchronous mode. It is obvious that the next
generation solution will span many scientific workflow domains and
layers. Hence, workflow framework and solutions must be functional,
dependable (reliable) \[e.g., VR03, Col04. Moul05, Vou05\], scalable,
domain-sensitive, verifiable (including semantics), efficient and user
friendly (at many levels), and must provide for both workflow status
assessments at run-time and data provenance, as well as debugging
options. The scope ranges from infrastructure development and deployment
of a highly dependable framework and components for solutions, to data
provenance and "dashboarding", to high-end semantic-level end-user
assistance and workflow verification and validation.

![](media/image11.png){width="3.3958333333333335in"
height="2.292361111111111in"}![](media/image12.png){width="3.234722222222222in"
height="2.2819444444444446in"}

Figures 5 and 6 -- TSI (Blondin) workflow abstraction with asynchronous
service-based control and tracking (left) and synchronous pilot for the
same workflow (Blondin).

![](media/image13.png){width="3.375in"
height="2.129861111111111in"}![](media/image14.png){width="3.25in"
height="2.127083333333333in"}

Figures 7 and 8 -- Layering approach to implementation of detachable TSI
(Blondin) workflows with run-time monitoring and provenance (left), and
example of a run-time information being collected on the current pilot.

As already mentioned, despite our significant successes, much work
remains before scientific workflow technology reaches state where it is
a commodity tool, easily used by application scientists. There are three
major, customer-driven, areas:

![](media/image15.png){width="3.125in"
height="2.7069444444444444in"}![](media/image16.png){width="3.4972222222222222in"
height="2.671527777777778in"}

Figures 9 and 10 -- Run-time logs and tracking information taken from
TSI (Swesty) workflow implementation.

-   **User interactions**. We will extend the current user interfaces to
    support web-based portals, domain specific "dashboards", template
    workflows for common tasks, and generalized actors that perform
    common, high-level functions while removing the specific
    implementation details from the user environment. For example, as
    illustrated in our pilot implementations related to TSI workflows
    (Blondin related in Figures 7 and 8 and Swesty related in Figures 9
    and 10), end-users do want to track progress the workflows are
    making. In order to do that effectively, we propose to extend our
    framework to enable both. Workflow construction and operation
    "wizards" are in the plan too.

-   **Data collection and management.** We will extend the current
    system to support data provenance, workflow provenance, and fault
    tolerances. In terms of data that we plan to collect, part of the
    needs are illustrated in Figures 5 through 8. We conducted several
    reliability and fault-tolerance studies and plan to incorporate the
    most promising mechanisms into the next generation SPA/Kepler
    framework.

-   ![](media/image17.png){width="3.6244094488188976in"
    height="2.538976377952756in"}**Workflow infrastructure.** It is
    interesting to note that one of the "next generation" issues is
    seamless integration of loosely coupled sub-workflows (e.g., Kepler
    based) with more tightly coupled ones (e.g., within problem solving
    environments such s SCIRun). Our pilots (Figure 11) indicate that
    such coupling is feasible. By applying SCIRun, we plan to make the
    visualizing process much more of an interactive plug-and-play
    environment. This will allow end-users to stay in the domain of the
    tools they are used to (e.g., Ensight) while making the process much
    more robust and efficient. SCIRun provides tools for modeling,
    computation, analysis and visualization of scientific datasets.
    Another desirable feature is interoperability with other potential
    workflow framework. Hence, will plan to extend our framework to
    support easy bridging between tightly-coupled and loosely coupled
    workflows; to provide distributed, asynchronous workflow execution
    and monitoring; and to improve workflow construction and deployment
    infrastructure (e.g., via "wizards").

[^1]: Another name that may characterize the data more accurately is
    Writing-Once Read-Many (WORM).

[^2]: <http://ptolemy.eecs.berkeley.edu/ptolemyII/>

[^3]: [http://kepler.ecoinformatics.org](http://kepler.ecoinformatics.org/)

[^4]: http://software.sci.utah.edu/scirun.html

[^5]: http://www.phy.ornl.gov/tsi/

[^6]: The Terascale Supernova Initiative (TSI,
    http://www.phy.ornl.gov/tsi/) is a DoE/OS/SciDAC sponsored
    "multidisciplinary collaboration of one national lab and eight
    universities to develop models for core collapse supernovae and
    enabling technologies in radiation transport, radiation
    hydrodynamics, nuclear structure, linear systems and eigenvalue
    solution, and collaborative visualization."