Systematic analysis of cluster computing log data: the case of IBM BlueGene/Q

The complexity and cost of managing large computing infrastructures is on the rise. Automating management actions to minimize human intervention is an attempt at containing these costs. This in turn requires building models and making predictions through analysis of log data. In current studies, this analysis is often limited to restricted datasets, and as such, is not applicable at a system level. Building models that are accurate enough to be useful in real systems requires analysis of log data from disparate sources. In this paper we provide a characterization study with four datasets reporting on power consumption, temperature, workload and hardware/software events for an IBM BlueGene/Q cluster. We show that the system handles a very rich parallel workload, with low correlation among its components in terms of temperature and power, but higher correlation in terms of events. Power and temperature correlate strongly, while events display negative correlations with load and power. Power and workload show moderate correlations, and only at the level of components. The aim of the study is a systematic, integrated characterization of the computing infrastructure and discovery of correlation sources and levels to be used in future predictive modeling efforts.

I. INTRODUCTION
Recent technological developments have resulted in large computing infrastructures becoming critical for many applications, not only for specialized tasks, but also for governing many aspects of our daily lives. This in turn renders their continued and correct functioning of utmost importance. At the same time, increases in the size and complexity of such systems has significantly hampered their manageability. Although progress has been made in modeling and predicting the behavior of large computing infrastructures, results have not reached a level of accuracy necessary for wide-spread industry adoption. This leaves the management of large computing systems a difficult activity that continues to rely heavily on human interventions [1]. One possible cause for this disconnect is the narrow scope of most analyses limiting their focus to particular subsystems rather than considering the system as a whole, which is a necessary condition for achieving high performance [2].
Large computing systems produce large amounts of data in the form of logs tracing resource consumption, errors and various events. These data can be put to use for building predictive models of system behavior in an effort to tackle the management challenges. With recent progress in Data Science and Big Data tools, this approach to system management is becoming increasingly promising [3]. It is also becoming increasingly feasible to integrate data from multiple sources describing different components of the whole system to obtain more realistic models towards more reliable predictions.
In this paper we present the results of study aimed at integrating data from a multitude of sources in an effort to understand the behavior and relations among the components for a large computing system. We have analyzed four different datasets obtained from different subsystems of a 10-rack IBM BlueGene/Q supercomputer [4]. These include power and temperature logs, workload traces and hardware/software events during an overlapping period of about 3 months. We fist present studies considering each of the datasets independently seeking significant features that describe system behavior and correlations among system components. Next we consider multiple datasets taken together looking for correlations among various metrics so as to identify further possible relations to be used in future modeling and prediction studies.
The contributions of this paper are twofold. First, we provide a characterization of a BlueGene/Q machine from thermal, power consumption, workload and event log perspectives, stressing the relations between different components of the system. We show that power consumption and thermal behavior are independent among components, yet events exhibit significant correlations, indicating possible propagation of errors. Secondly, an integrated analysis looking at the four datasets together reveals significant positive correlation between power consumption and temperature, yet a weaker negative correlation between hardware/software events and power consumption or workload.
The rest of the paper is organized as follows. In the next section we describe of the data and the system which it comes from. We then proceed to the data analysis, both for individual datasets and correlated features, in Section III. Section IV discusses related work and Section V concludes the paper with some final remarks and ideas for future work.

II. DATA DESCRIPTION
The source of our data is Fermi [5], an IBM BlueGene/Q system run by CINECA, the non-profit consortium of 69 Italian universities and 3 Public Administration Institutions, operating the largest data center in Italy. The system has a total of 163,840 computing cores and handles HPC workload from the consortium's academic partners. Fermi consists of 10 racks, each containing 2 mid-planes. Each mid-plane has 16 node-boards with 32 16-core nodes per node-board. Each mid-plane is powered by 18 so-called bulk power modules. System logging is based on standard BlueGene/Q tools [4]. One such tool is the Mid-plane Manager Control System which is in charge of environmental monitoring. This is where we extract power and temperature information for the system. A different tool is the Machine Controller, which controls access to the different hardware components and provides socalled Reliability, Availability and Serviceability (RAS) logs. Information about the workload is obtained through a custom tool developed in-house at CINECA. Table I summarizes the four datasets. These span similar time periods, with an overlap of about three months, and include information at different levels of the system obtained from sensors on different components.
Power logs were extracted from the bulk power monitor of Fermi, which reports input and output voltages and current levels for each bulk power module, with a 5-minute resolution. We used input values to compute power levels. By summing the power levels over the different components, we obtained time series of power consumption for individual mid-planes, racks and for the entire system. Power at the node-board level cannot be reliably computed from these data since 18 bulk power modules are used to power 16 node-boards (redundant system).
Temperature logs were extracted from the node-board monitor, which reports two temperatures for each node-board, with a 15-minute resolution. Using these values we were able to compute time series of average temperatures at node-board, mid-plane, rack and system levels.
Workload data was reported by CINECA in a more aggregated form: a list of jobs with the date of completion, running time, number of cores and queue time. This allowed us to compute the CPU time used by every job, as well as time series of total daily CPU times, number of cores and queue times. Since only the date of job completions (and not the exact times) are available in the data, these totals are approximate, yet they give a very good indication of the daily load at system level. No load information at other levels (node-board, midplane, rack) was available at this time.
RAS logs were the fourth dataset available. They consist of hardware and software events from all system components. Events are labeled based on severity: our dataset includes three event categories with 163,134 FATAL, 473,982 WARN and 137,438 INFO events. Also, the exact time and exact location of each event is included. Using these data, we computed the distribution of inter-event times at system level, for each severity category. We also obtained time series monitoring the number of events in each category at 24-hour and 5-minute resolutions. These were obtained at node-board, mid-plane, rack and system levels.
III. DATA ANALYSIS While each of the above datasets may provide very useful insight into the functioning of the Fermi system, combining them for an integrated analysis has even greater potential. Consequently, in this section we first perform an analysis of each dataset individually, identifying their features and comparing them. This is followed by an integrated analysis to study how different metrics from different subsystems correlate.

A. Power logs
The system specifications for an IBM Blue/GeneQ machine declare the typical power consumption to be around 65kW, with a maximum of 100kW per rack [6]. However, real consumption varies depending on system load and state of components, e.g. how many nodes are up. Figure 1 displays the distribution of power consumption sampled at 5-minute intervals. Three different distributions are shown: total power at system level followed by power consumption by each rack and mid-plane. All distributions are centered around the official declared average values: 650kW at the system level consisting of 10 racks, 65kW at the rack level and 32.5kW at the midplane level, confirming the specifications. As we move from higher levels to lower ones, the distribution becomes broader. While total consumption is mostly between 50kW and 70kW per rack, with a bell-shaped distribution, when looking at individual racks additional peaks emerge, with some racks showing power consumptions up to 90kW for some time windows, but also frequent values under 50kW. The additional peaks become even more evident at the mid-plane level. This difference in distributions shows that power consumption is very heterogeneous. This needs to be taken into account when modeling power consumption, and indications are that while predicting overall system power consumption might be easier due to a large stability in time, finer grained predictions at rack or mid-plane levels might produce more accurate results.
Since power consumption at rack and mid-plane levels appears to be more heterogeneous, it is interesting to see if power values are similar in time when comparing different components (racks or mid-planes). Considering that power consumption depends on system load, load balancing algorithms should make it so that the different parts of the system work equally. So power consumption should be correlated among components of the system (racks or mid-planes), unless load is too low. Figure 2 shows Pearson correlations of power consumption between rack pairs and between mid-plane pairs. It is clear that at both levels, correlations are in general very low. Only a few mid-plane pairs seem to have correlations above 0.5 which are only moderate. Given that in section III-C we will observe a high general load of the system, this may indicate that load is not very well balanced.

B. Temperature logs
Distribution of average temperatures has also been analyzed. Figures 3 and 4 show histograms of average tem- peratures sampled every 15 minutes. At the overall system level, again the distribution is bell-shaped and narrow with few exceptions, with one mode around 50 • C. As we zoom in at different component levels, the distribution becomes again wider with additional peaks showing up again at very high and very low temperatures. Individual node-boards can reach up to 75 • C, significantly larger than the system average. This shows again how at different levels the system appears to behave differently, with greater heterogeneity in time at the finer grain logs.
Looking again at correlations among different components of the same type ( Figure 5), a pattern similar to power consumption is observed. With very few exceptions, temperatures among components are not correlated. In terms of thermal isolation, this is good news, since having one hot node-board does not mean all surrounding node-boards are hot as well. However since power consumption showed the same pattern, this can be considered an additional evidence that workload is not very well balanced in the system.

C. Workload logs
A first important question in terms of workload is what type of jobs are submitted to the system. For the 78,128 jobs that ran on the system during the logged period, Figure 6 displays the distribution of different job attributes: total CPU time, total running time and number of cores used. In terms of time requirements, jobs are very heterogeneous as evidenced by a long tailed CPU time distribution, with a few very heavy jobs and many short jobs present. Effective running times are bimodal, with many short jobs and many long jobs (all running times under 24 hours), and slightly fewer medium-length jobs. The number of cores per job, on the other hand, is not as heterogeneous. Jobs require only eight different values for the core count, most using over 100 cores and up to 32768. So, in general, the jobs that are run are highly parallel. The structure of the workload dataset allows for analysis of patterns in time only at system level and at 24-hour resolution. Figure 7 shows the total CPU time for jobs completed each day. Note that this does not represent the exact system load for that day. The data only contains the date of completion for each job, not the exact time, so it is not possible to check exactly how many hours each job ran in a day. Jobs completed on one day could have been started the previous day. However the total time used by jobs completed in the same day still gives an indication of system load. A roughly bell-shaped distribution is obtained again, with a mean around 13 billion seconds (around 400 years), indicating very high levels of system load for Fermi. For RAS events, we first look at inter-event times at the system level. Figure 8 shows the distribution of inter-event times, in seconds for the three event categories. The plot is a rank plot: the time intervals between events were sorted in  descending order, and then their values plotted on the vertical axis. The inter-event times do not appear to follow a known distribution. FATAL events show a few very large intervals and then many very small intervals. This indicates a pattern with spikes of many events in short periods of time but with large breaks between them. INFO and WARN events do not show very large inter-event times, and a smaller fraction of very short intervals, indicating they are more evenly spread across time.

D. RAS logs
Further analysis of the time distribution of events is performed by counting the number of events in each category at different time windows. We use two window sizes, to explore different time resolutions: 24-hour and 5-minute windows. Figure 9 displays the 24-hour resolution time series and their relative correlations. It is clear how WARN and INFO events are much more common, with most days showing some instances. While FATAL events come in spikes in a few of the monitored days. At the same time, it appears that daily INFO and WARN events are highly correlated, and so are WARN and FATAL events. However INFO and FATAL events seem to appear together less frequently. This could mean that INFO and WARN events could be used for predicting FATAL events at this time resolution. Similar data is shown for the finer grained time resolution (5-minute windows) in Figure 10. Here we display does not affect 24-hour correlations but decreases them at 5minute resolution.
A different question is whether events correlate not only across categories, but also within the same category and across different components. Figure 11 shows the distribution of correlations between racks, mid-planes and node-boards for FATAL events. Similar results were obtained for the other two categories. Unlike power and temperature, FATAL events have higher correlation across components, with a significant number of component pair correlations over 0.5. This indicates that failures may propagate across components, while temperature and power consumption are more independent.

E. The big picture
The analysis of individual datasets has shed some light into the functioning of the Fermi system, and where relations between components might exist. Here we look at the four datasets together to uncover further relations between the different components and logs.
A first analysis is at the system level, looking at the different measures for the overall system rather than at individual racks, mid-planes or nodes-boards. Figure 12 shows correlations between the different datasets for two time windows: 24 hours and 5 minutes. Workload data can only be analyzed on a daily basis, so it is not included at the 5-minute resolution.
As we have seen earlier, RAS events correlate very well at 24-hour resolution and very little at 5-minute resolution. There also seems to be a strong relation between temperature and power, at both resolutions. This is a well known correlation for computer systems, and has been observed in other systems as well [7].
Workload measures do not appear to correlate among themselves, nor with power consumption, at 24-hour resolution. Better correlations could be obtained by zooming in at rack, mid-plane or node-board levels, however the structure of our workload dataset does not allow for this analysis. This prompts for changes in workload log structure in Fermi. Low correlation is present between system temperature and CPU time. Again, this should increase at other component levels. Some negative relation appears to exist between RAS events and CPU time, which is somewhat counterintuitive: one would expect more events to appear when the system works harder. However, this could be due to the coarse-grained time resolution. When looking at daily patterns, it is quite possible for large numbers of RAS events to have resulted in system failure, which in turn resulted in fewer completed jobs. This could explain the negative relation.
A negative correlation also appears between power/temperature and RAS events, again rather counterintuitively. This means that when there are more events, the system is working and consuming less, which can be due to complete shutdown of components after frequent RAS events. So, when trying to predict power consumption or FATAL events, one needs to take into consideration the negative dependence. Similar to RAS events, when increasing resolution (5-minute windows) the correlations dampen out. This indicates possible cross-correlations, i.e. correlations obtained between the two time series after shifting one of them by a specific time, due to a time shift in the events themselves.
The system level analysis might enhance some correlations, when behavior of components propagate or dampen them out by averaging over components. For this reason, we have also analyzed correlations at rack, mid-plane and node-board levels. Figure 13 shows correlations between available datasets at each level. Again, a strong relation between temperature and power consumption is visible. This suggests that in making predictions, using only one of the two features might suffice. In the case of BlueGene/Q, this is good news since power logs at node-board level are not available.
At rack level, RAS events are well correlated, but then correlations decrease as the level of detail increases (moving from rack to node-board level). This could be an artifact of rack-level events including events related to components other than node-boards. So, as components become more specific, the number of total events decreases, resulting in a very sparse dataset for which the Pearson correlation might be unsuitable.

IV. RELATED WORK
Log analysis for large computing infrastructures has been the focus of numerous studies. Recently, Google have released two workload traces for one of their clusters, triggering a flurry of analyses trying to characterize their system. General statistics have been provided [8], [9], [10], and more heterogeneity compared to grids has been observed [11]. Some modeling work has also appeared based on these data. In terms of resource usage, classes of jobs and their prevalence were used to characterize workloads and generate new ones [12], [13], or real usage patterns were replaced by the average utilization [14]. A simulator has also been built to reproduce the Google cluster load [15]. Important insight has been obtained from these data, however they only look at one part of the system, i.e. the workload. In modeling and prediction studies, it is important to integrate data from different components and sources. Other traces have also been analyzed in the past [16], [17], [18], and tools for analysis developed [19], but again concentrating on one type of data only, such as workload or failure traces. Here we integrate various datasets for a systematic analysis.
For IBM BlueGene, some previous work on failure prediction using RAS logs exists. In [20] failure prediction in a BlueGene/Q cluster is attempted. Their aim is to combine low-level hardware failure with high level kernel logs, however they stop at classification using RAS logs only. They analyze several classification tools (SVM, customized KNN, feature selection, rule based models) and obtain best performance with the KNN classifier. Another comparison of different classification tools for failure prediction in an IBM BlueGene/L cluster is [21]. They analyze SVMs, neural networks and rule based classifiers and also compute feature significance. They try to predict whether different event categories will appear in a time window based on previous windows. These predictive studies look only at RAS events, while integrating further data related to power consumption, temperature and workload could improve prediction accuracy significantly. Here we provide the first step towards such an analysis.
Most previous methods concentrate on workload or failure events, however recent work also includes thermal analysis of systems. In [22] multi-component chiller systems used in data centers are modeled. They build a Dynamic Bayesian Network model of the system using time series data from different sensors. Their aim it to control utilization levels of such equipment in order to optimize power consumption and lifetime. Again, this study does not integrate other types of data in the analysis. Some integration is performed by a very recent study from Google [3], who model Power Usage Effectiveness of their system using thermal information (temperatures, humidity, etc) together with overall system load, using an Artificial Neural Network model. Recently, a novel monitoring system for several criteria has been built [7]. This monitors Eurora, another HPC machine from CINECA, the same organization that provided our data. Their aim is thermal characterization of the system, however they do record several types of data besides power and chiller characteristics, including workload and machine status. These data could be, in principle, used for predictive and modeling studies in the future.

V. CONCLUSIONS
Given the need for a systematic integrative analysis of large computing infrastructures, this paper has presented an analysis of four data sets describing different subsystems of an IBM BlueGene/Q machine. Temperature, power consumption, workload and RAS logs were studied independently to characterize the system systematically and then together to identify features of the data that might help in future integrated predictive or modeling studies.
The power consumption data showed that subcomponents of the system do not have correlated power consumptions, suggesting that load balancing might need attention for this system. Workload analysis showed high levels of load of the system, on average, with highly parallel jobs. RAS events of different severity showed correlation when looking at hourly patterns, but not when time resolution was increased to 5minute windows. This suggests that there might be crosscorrelations in the data that should be investigated further for modeling and prediction purposes. Across components, correlations were higher than in the case of power consumption, indicating that failures tend to propagate between components.
Not surprisingly, integrated data analysis showed that power consumption and system temperatures correlate very well, regardless of the space or time resolution. This indicates that, for the future, it might be useful to select only one of the two features for data analysis, in order to decrease dimensionality, which is good news also because power logs do not reach the same spatial resolution as temperature logs. Unexpectedly, when looking at daily levels of workload and power consumption versus RAS events, negative correlations were found, which can be explained by components turning off when errors appear and thus resulting in lower power consumption and less work completed. This has to be taken into account in predictive studies, where one would expect more events to be logged when the system is busier.
This study will continue with data external to the computing infrastructure, such as the water cooler and indoor air conditioning, together with data external to the data center, such as weather and seismic activity. RAS events will be further investigated to identify which of them indicate full shutdown of components, in order to verify our hypothesis for the negative correlations. Also, a predictive study will be undertaken, to use the correlations we already obtained to predict component failure. For this objective, cross-correlations will be also investigated.


extfieldsofstudy:
Mathematics
Computer Science

text:
Linear Ranking for Linear Lasso Programs

The general setting of this work is the constraint-based synthesis of termination arguments. We consider a restricted class of programs called lasso programs. The termination argument for a lasso program is a pair of a ranking function and an invariant. We present the---to the best of our knowledge---first method to synthesize termination arguments for lasso programs that uses linear arithmetic. We prove a completeness theorem. The completeness theorem establishes that, even though we use only linear (as opposed to non-linear) constraint solving, we are able to compute termination arguments in several interesting cases. The key to our method lies in a constraint transformation that replaces a disjunction by a sum.

Introduction
Termination is arguably the single most interesting correctness property of a program. Research on proving termination can be divided according to three (interrelated) topics, namely: practical tools [1,9,13,17,18,19,21,22], decidability questions [4,8,25], and constraint-based synthesis of termination arguments [2,3,5,6,7,10,12,14,20,23]. The work in this paper falls under the research on the third topic. The general goal of this research is to investigate how one can derive a constraint from the program text and compute a termination argument (of a restricted form) through the solution of the constraint, i.e., via constraint solving.
In this paper, we present a method for the synthesis of termination arguments for a specific class of programs that we call lasso programs. As the name indicates, the control flow graph of a lasso program is of a restricted shape: a stem followed by a loop.
Lasso programs do not appear as stand-alone programs. Lasso programs appear in practice whenever one needs a finite representation of an infinite path in a control flow graph, for example in (potentially spurious) counterexamples in a termination analysis [13,17,18,19], non-termination analysis [16], stability analysis [11,22], or cost analysis [1,15]. Importantly, the termination argument for a lasso program is a pair of a ranking function and an invariant (the rank must decrease only for states that satisfy the invariant). Figure 1 shows an example of a lasso program.
The class of lasso programs lies between two classes of programs for which constraint-based methods have been studied extensively. For the first, more specialized class, methods can be based on linear arithmetic constraint solving [2,3,10,12,20]. For the second, more general class, all known methods are based on non-linear arithmetic constraint solving [5,7]. The contribution of our method can be phrased, alternatively, as the generalization of the applicability of the 'linear methods', or as the optimization of the 'non-linear method' to a 'linear method' for a subproblem. The step from 'non-linear' to 'linear' is interesting for principled reasons (non-linear arithmetic constraint solving is undecidable in the case of integers). As we will show the step is also practically interesting.
The reader may wonder how practical tools presently handle the situation where one needs to compute termination arguments for lasso programs. One possibility is to resort to heuristics. For example, instead of computing a termination argument for the lasso program in Figure 1, one would compute the ranking function f (x) = x for the program while(x>=0){x:=x-23;}.
The key to our method is a constraint transformation that replaces a disjunction by a sum. We apply the 'or-to-plus' transformation in the context of Farkas' Lemma. Following [2,5,10,12,20], we apply Farkas' Lemma in order to eliminate the universal quantifiers in the arithmetic constraint whose solution is the termination argument. If we apply Farkas' Lemma to the constraint after the 'or-to-plus' transformation, we obtain a linear arithmetic constraint.
The effect of the 'or-to-plus' transformation to the constraint is a restriction of its solution space. The restriction seems strong; i.e., in some cases, the solution space becomes empty. We can characterize those cases. In other words, we can characterize when the 'or-to-plus' transformation leads to the loss of an termination argument, and when it does not. The characterization is formulated as a completeness theorem for which we will present the proof. This characterization allows us to establish that, even though we use only linear (as opposed to non-linear) constraint solving, we are able to compute termination arguments in several interesting cases. A possible explanation for this (perhaps initially surprising) fact is that, for synthesis, we are interested in the mere existence of a solution, and the loss of many solutions does not necessarily mean the loss of all solutions of the constraint.
We have implemented our method and we have used our implementation to illustrate the applicability and the efficiency of our method. Our implementation is available through a web interface, together with a number of example programs (including the ones used in this paper). 1 1: y := 23; 2: while( x >= 0 ) { 3: x := x -y; 4: y := y + 1; 5: } 1 2 τstem : y ′ = 23 τ loop : Example of a lasso program and its formal representation P yPositive = (τstem, τ loop ). The ranking function defined by f (x, y) = x decreases in transitions from states that satisfy the invariant y ≥ 1 (the ranking function does not decrease when y ≤ 0).

Preliminaries: Linear Arithmetic
We use x to denote the vector with entries x 1 , . . . , x n , and x ⊺ to denote the transposed vector of x. As usual, the expression A·x ≤ b denotes the conjunction of linear constraints We call a relation τ (x, x ′ ) a linear relation if τ is defined by a conjunction of linear constraints over the variables x and x ′ , i.e., if there is a matrix A with m rows and 2n columns and a vector b of size m such that the following equation holds. τ is defined by an affine linear term, i.e., there is a vector r ⊺ and a number r 0 such that the following equation holds.
We call a predicate I(x) a linear predicate, if I(x) is defined by a linear inequality, i.e., there is a vector s ⊺ and a number s 0 such that following equivalence holds.
Farkas' Lemma. We use the affine version of Farkas' Lemma [24] which is also used in [2,5,12,23,20] and states the following. Given a satisfiable conjunction of linear constraints A · x ≤ b and a linear constraint c ⊺ · x ≤ δ, the following equivalence holds.

Lasso Program
To abstract away from program syntax, we define a lasso program directly by the two relations that generate its execution sequences.
Definition 1 (Lasso Program). Given a set of states Σ, a lasso program is given by the two relations τ stem ⊆ Σ × Σ and τ loop ⊆ Σ × Σ. We call τ stem the stem of P and τ loop the loop of P . An execution of the lasso program P is a possibly infinite sequence of states σ 0 , σ 1 , . . . such that the pair of the first two states is an element of the stem, i.e., (σ 0 , σ 1 ) ∈ τ stem and each other consecutive pair of states is an element of the loop, i.e., We call the lasso program P terminating if P has no infinite execution.
We use constraints over primed and unprimed variables to denote a transition relation (see Figure 1). In order to avoid cumbersome technicalities, we consider only lasso programs that have an execution that contains at least three states. This means we consider only programs where the relational composition of τ stem and τ loop is non-empty, i.e., τ stem • τ loop = ∅.
Since Turing, a termination argument is based on an ordering which does not allow infinite decreasing chains (such as ordering on the natural numbers). Here, we use the ordering over the set of positive reals which is defined by some value δ > 0, namely Ranking Function. We call a function f from the states of the lasso program P into the reals R a ranking function for P if there is a positive number δ > 0 such that for each consecutive pair of states (x i , x i+1 ) of a loop transition (i ≥ 1) in every execution of P the value of f is decreasing by at least δ, i.e., and the value of f is non-negative, i.e., If there is a ranking function for the lasso program P , then P is terminating.
Inductive Invariant. We call a state predicate I(x) an inductive invariant of the lasso program P if the predicate holds after executing the stem, i.e., and if the predicate holds before executing the loop, then the predicate holds afterwards, i.e., Ranking Function with Supporting Invariant. We call a pair of a ranking function f (x) and an inductive invariant I(x) of the lasso program P a ranking function with supporting invariant if the following holds.
-There exists a positive real number δ > 0 such that, if the inductive invariant holds then an execution of the loop decreases the value of the ranking function by at least δ, i.e., -In states in which the inductive invariant holds and the loop can be executed, the value of the ranking function is non-negative, i.e., For example, the lasso program depicted in Figure 1 has the ranking function f (x, y) = x with supporting invariant y ≥ 1.
Linear lasso programs. Linear lasso programs. For the remainder of this paper we consider only linear lasso programs, linear ranking functions, and linear inductive invariants which we will define next. The variables of the programs will range over the reals until we come to Section 9 where we turn to programs over integers.

Definition 2 (Linear Lasso Program). A linear lasso program
is a lasso program whose states are vectors over the reals, i.e. Σ = R n , and whose relations τ stem and τ loop are linear relations.

We use the expression
Linear Ranking Function. If a ranking function f : R n → R is an (affine) linear function, we call f a linear ranking function. We use r 1 , . . . , r n as coefficients of a linear ranking function, r as their vector, Linear Invariant. If an inductive invariant I(x) is a linear predicate, we call I a linear inductive invariant. We use s 1 , . . . , s n as coefficients of the term that defines the linear predicate, s as their vector,

The Or-to-Plus Method
Our constraint-based method for the synthesis of linear ranking functions for linear lasso programs consists of three main steps: Step 1. Set up four (universally quantified) constraints whose free variables are the coefficients of a linear ranking function with linear supporting invariant.
Step 2. Apply Farkas' Lemma to the four constraints to obtain equivalent constraints without universal quantification.
Step 3. Obtain solutions for the free variables by linear constraint solving.
The particularity of our four constraints in Step 1 is that the application of Farkas' Lemma in Step 2 yields constraints that are linear.
Instead of presenting our constraints immediately, we derive them in three successive transformations of constraints. We start with the four constraints (ϕ invStem ), (ϕ invLoop ), (ϕ rkDecr ), and (ϕ rkBound ). Below, we have rephrased the four constraints for the setting where the ranking function is linear and the supporting invariant is linear. We marked them (ϕ BMS 1 ), (ϕ BMS 2 ), (ϕ BMS 3 ), and (ϕ BMS 4 ) in reference to Bradley, Manna and Sipma [5] who were the first to use them in the corresponding step of their method.

The Bradley-Manna-Sipma constraints
for the special case of lasso programs and one linear supporting invariant 2 The free variables of ϕ BMS are r, r 0 , s, and s 0 .
Transformation 1: Move supporting invariant to right-hand side. We bring the conjunct s ⊺ · x + s 0 ≥ 0 in three of the four constraints (ϕ BMS 1 ), (ϕ BMS 2 ), (ϕ BMS 3 ), and (ϕ BMS 4 ) to the right-hand side of the implication, according to the following scheme.
We obtain the following constraints.
Transformation 2: Drop supporting invariant in fourth constraint. We strengthen the fourth constraint (ψ 4 ) by removing the disjunct −s ⊺ · x − s 0 > 0. A solution for the strengthened constraint defines a ranking function whose value is bounded from below for all states (and not just those that satisfy the supporting invariant).
Transformation 3: Replace disjunction by sum. We replace the disjunction on the right-hand side of the implication in constraints (ψ 2 ) and (ψ 3 ) by a single inequality, according to the scheme below. (It is the disjunction which prevents us from applying Farkas' Lemma to the constraints (ψ 2 ) and (ψ 3 ).) In the second constraint (ψ 2 ), we replace the disjunction In the third constraint (ψ 3 ), we replace the disjunction We obtain the following four constraints.
Since we consider linear lasso programs, the relations τ stem and τ loop are given as conjunctions of linear constraints.
We have now finished the description for the three transformation steps that lead us to the the or-to-plus constraints. We are now ready to introduce our method.
The Or-to-Plus Method Input: linear lasso program P . Output: coefficients r, r 0 , s, and s 0 of a linear ranking function with linear supporting invariant 1. Set up the or-to-Plus constraints ϕ 1 , ϕ 2 , ϕ 3 , and ϕ 4 for P .
After setting up the four or-to-plus constraints ϕ 1 , ϕ 2 , ϕ 3 , ϕ 4 in Step 1, we apply Farkas' Lemma to each of the four constraints in Step 2. We obtain four linear constraints. E.g., by applying Farkas' Lemma to the constraint (ϕ 3 ) we obtain the following linear constraint.
We apply linear constraint solving in Step 3. We obtain a satisfying assignment for the free variables in the resulting constraints. The values obtained for r, r 0 , s and s 0 are the coefficients of a linear ranking function f (x) with linear supporting invariant I(x).
The or-to-plus method inherits its soundness from method of Bradley-Manna-Sipma.
Step 1 is an equivalence transformation on the Bradley-Manna-Sipma constraints, Step 2 and Step 3 strengthen the constraints, and the application of Farkas' Lemma is an equivalence transformation. Thus, a satisfying assignment of the or-to-plus constraints obtained after the application of Farkas' Lemma is also a satisfying assignment of the Bradley-Manna-Sipma constraints.

Completeness of the Or-to-Plus Method
In the tradition of constraint-based synthesis for verification, we will formulate completeness according to the following scheme: the method X applied to a program P in the class Y computes (the coefficients of) a correctness argument of the form Z whenever one exists (i.e., whenever a correctness argument of the form Z exists for the program P ). Here, X is the or-to-plus method, Y is the class of lasso programs, and Z is a termination argument consisting of a linear ranking function and an invariant of a form that we we define next.

Definition 3 (Non-decreasing linear inductive invariant). We call a linear inductive invariant s
⊺ · x + s 0 ≥ 0 of the lasso program P non-decreasing if the loop implies that the value of the term s ⊺ · x + s 0 does not decrease when executing the loop, i.e., In Section 6 we give examples which may help to convey some intuition about the meaning of 'non-decreasing', examples of those terminating programs that do have a linear ranking function with a non-decreasing linear supporting invariant, and examples of those that don't.
Theorem 1 (Completeness). The or-to-plus method applied to the linear lasso program P succeeds and computes the coefficients of a linear ranking function with non-decreasing linear supporting invariant whenever one exists.
To prove this theorem we use the following lemma.
then there exists a constant µ ≥ 0 such that Proof (of Lemma 1).
defined by a weighted sum is a rotation of H around the intersection Z of the boundary of H and the boundary of H ′ . If a polyhedron X is contained in the union H ∪ H ′ , then there is a half-space Hµ defined by a weighted sum that contains X. By assumption, (1) does not imply (3), so A·x ≤ b ∧ h ⊺ ·x+h 0 ≤ 0 is satisfiable, and by Farkas' Lemma this formula is equivalent to and thus Because A · x ≤ b is satisfiable by assumption, Farkas' Lemma can be applied again to yield Theorem 1). Let f (x) =r ⊺ · x +r 0 be a ranking function with nondecreasing supporting invariant I(x) ≡s ⊺ · x +s 0 ≥ 0 for the lasso program P . Since executions of our lasso programs comprise at least three states, there can be no supporting invariant that contradicts the loop, i.e.
is not valid. From (ϕ rkBound ) it follows that and hence the implication is valid. By (1) and Lemma 1 there is a µ 1 ≥ 0 such that is valid. If we assign r →r − µ 1 ·s, r 0 →r 0 − µ 1 ·s 0 , then (ϕ 4 ) is satisfied.
Because I(x) ≡s · x +s 0 ≥ 0 is a non-decreasing invariant, and hence, since µ 1 ≥ 0, From (ϕ rkDecr ) we know that and hence equivalently With (2) we obtain validity of the following formula.
By (1) and Lemma 1 there exists a µ 2 ≥ 0 such that We pick the assignment r →r − µ 1 ·s, r 0 →r 0 − µ 1 ·s 0 , s → µ 2 ·s, s 0 → µ 2 ·s 0 , which hence satisfies (ϕ 3 ). We already argued that it satisfies (ϕ 4 ), and from µ 2 ≥ 0 and the fact that I(x) is a non-decreasing inductive invariant it follows that the assignment also satisfies (ϕ 1 ) and (ϕ 2 ). Hence, the ranking function (r − µ 1 ·s) ⊺ · x +r 0 − µ 1 ·s 0 with supporting invariant (µ 2 ·s) ⊺ · x + µ 2 ·s 0 ≥ 0 can be found by the or-to-plus method. Our three transformations strengthened the Bradley-Manna-Sipma constraints, hence the solution space of the or-to-plus constraints is smaller than the solution space of the Bradley-Manna-Sipma constraints. This can be seen e.g., in the example depicted in Figure 4. The program P bound has the linear ranking function f (x, y) = x with linear supporting invariant y ≥ 23, but the coefficients of this ranking function and supporting invariant are no solution of the or-to-plus constraints; the constraint ϕ 4 is violated. Does this mean that our method will not succeed? No, it does not. By Theorem 1, in fact, we do know that the method will succeed. I.e., since we know of some linear ranking function with non-decreasing supporting invariant (in this case, f (x, y) = x and y ≥ 23), even if it is not a solution, we know that there exists one which is a solution (here, for example, f (x, y) = x − y with the (trivial) supporting invariant 0 ≥ 0). The prerequisite of Theorem 1 is the existence of a non-decreasing supporting invariant. There are linear lasso programs that have a linear ranking function with linear supporting invariant, but do not have a linear ranking function with a non-decreasing linear supporting invariant. E.g., for the lasso programs depicted in Figure 5 and Figure 6 our or-to-plus method is not able to synthesize a ranking function for these programs.
The linear lasso program P zeno depicted in Figure 5 has the linear ranking function f (x, y) = x with the linear supporting invariant y ≥ 1. However this inductive invariant is not non-decreasing; while executing the loop the value of the variable y converges to 1 in the following sequence. The statement havoc y; in the lasso program P wild is a nondeterministic assignment to the variable y. The relations τ stem and τ loop of this lasso program are given by the constraints y ′ ≥ 1 and x ≥ 0 ∧ x ′ = x − y ∧ y ′ ≥ 1. P wild has the ranking function f (x, y) = x with the supporting invariant y ≥ 1, however this inductive invariant is not non-decreasing in each execution of the loop the variable y can get any value greater than or equal to one.
The next example shows that nondeterministic updates are no general obstacle for our or-to-plus method. In the linear lasso program P array the loop iterates over an array of positive integers. The index accessed in the next

Lasso Programs over the Integers
In the preceding sections we considered lasso programs over the reals. In this section we discuss the applicability of the or-to-plus method to linear lasso programs over the integers, i.e., programs where the set of states Σ is a subset of Z n . We still use real-valued ranking functions. We obtain the constraints for coeffi-cients of a linear ranking function with linear supporting invariant by restricting the range of the universal quantification in the constraints ϕ 1 , ϕ 2 , ϕ 3 , and ϕ 4 to the integers. E.g., the constraint ϕ 3 for linear lasso programs over the integers is where the domain of the coefficients r, r 0 , s, and s 0 and the quantified variable δ are the reals. Now, Farkas' lemma is not an equivalence transformation, its application results in weaker formulas. This means the or-to-plus method is still sound, but we loose the completeness result of Theorem 1. An example for this is assume 2*y >= 1; while( x >= 0 ) { x := x -2*y + 1; } the program P nonIntegral , depicted in Figure 8 that has the following transition relations.
Over integer variables, P nonIntegral1 has the linear ranking function f (x, y) = x with the linear supporting invariant y ≥ 1. Over real-valued variables, P nonIntegral1 does not terminate. If we add the additional constraint y ′ ≥ 1 to τ stem , the programs' semantics over the integers is not changed, but we are able to synthesize a linear ranking function with a linear supporting invariant. Adding this additional constraint gives the constraints a property that we formally define as follows.
Integral constraints. A conjunction of linear constraints A·x ≤ b is called integral if the set of satisfying assignments over the reals S := {r ∈ R n | A · r ≤ b} coincides with the integer hull of S (the convex hull of all integer vectors in S).
For each conjunction of m linear constraints there is an equivalent conjunction of at most 2 m linear constraints that is integral [24]. We add an additional step to the or-to-plus method in which we make the constraints in the stem transition τ stem and loop transition τ loop integral.
The Or-to-Plus Method (Int) Input: linear lasso program P with integer variables Output: coefficients r, r 0 , s, and s 0 of linear ranking function with linear supporting invariant Lemma 2 (Integral version of Farkas' Lemma). Given a conjunction of linear constraints A · x ≤ b which is satisfiable and integral, and a linear constraint c ⊺ · x ≤ δ, Proof. We write this statement as a linear programming problem.
Because the constraints A · x ≤ b are integral, there is an integral vector x ∈ Z n such that c ⊺ · x is the optimum solution to (P). Thus the optimum over integers is ≤ δ if and only if the optimum of the reals is. The statement now follows from the real version of Farkas' Lemma.
⊓ ⊔ assume 2*y >= z; while( x >= 0 && z == 1 ) { x := x -2*y + 1; } However, even if τ stem and τ loop are integral, our method is not complete over the integers. In the completeness proof for the reals we applied Farkas' Lemma to conjunctions of a polyhedron A · x ≤ b and an inequality h ⊺ · x + h 0 ≤ 0. This inequality contains free variables, namely the coefficients of the supporting invariant s ⊺ · x + s 0 ≥ 0. Even if τ stem and τ loop are integral, this conjunction might not be integral and we cannot apply the integer version of Farkas' lemma in this case.
A counterexample to completeness of our integer version of the or-to-plus method is the linear lasso program P nonIntegral2 depicted in Figure 9.

Conclusion
We have presented a constraint-based synthesis method for a class of programs that was not investigated before for the synthesis problem. The class is restricted (though less restricted than the widely studied class of simple while programs) but still requires the combined synthesis of not only a ranking function but also an invariant. We have formulated and proven a completeness theorem that gives us an indication on the extent of power of a method that does without nonlinear constraint solving.
We implemented the or-to-plus method as plugin of the Ultimate software analysis framework. A version that allows one to 'play around' with lasso programs is available via a web interface at the following URL. http://ultimate.informatik.uni-freiburg.de/LassoRanker As mentioned in the introduction, the class of lasso programs is motivated by the fact that they are a natural way (and, it seems, the only way) to represent an (infinite) counterexample path in a control flow graph. It is a topic of future research to explore the different scenarios in practical tools that use a module to find a ranking function and a supporting invariant for a lasso program (e.g., in [1,13,15,16,17,21,22]) and to compare the performance of ourtheoretically motivated-synthesis method in comparison with the existingheuristically motivated-approach used presently in the module.


extfieldsofstudy:
Computer Science
Materials Science
Mathematics
Biology

text:
Feature extraction for proteomics imaging mass spectrometry data

Imaging mass spectrometry (IMS) has transformed proteomics by providing an avenue for collecting spatially distributed molecular data. Mass spectrometry data acquired with matrix assisted laser desorption ionization (MALDI) IMS consist of tens of thousands of spectra, measured at regular grid points across the surface of a tissue section. Unlike the more standard liquid chromatography mass spectrometry, MALDI-IMS preserves the spatial information inherent in the tissue. Motivated by the need to differentiate cell populations and tissue types in MALDI-IMS data accurately and efficiently, we propose an integrated cluster and feature extraction approach for such data. We work with the derived binary data representing presence/absence of ions, as this is the essential information in the data. Our approach takes advantage of the spatial structure of the data in a noise removal and initial dimension reduction step and applies $k$-means clustering with the cosine distance to the high-dimensional binary data. The combined smoothing-clustering yields spatially localized clusters that clearly show the correspondence with cancer and various noncancerous tissue types. Feature extraction of the high-dimensional binary data is accomplished with our difference in proportions of occurrence (DIPPS) approach which ranks the variables and selects a set of variables in a data-driven manner. We summarize the best variables in a single image that has a natural interpretation. Application of our method to data from patients with ovarian cancer shows good separation of tissue types and close agreement of our results with tissue types identified by pathologists.

Introduction
Mass spectrometry (MS) has become a versatile and powerful tool in proteomics for the analysis of complex biological systems, including the identification and quantification of proteins and peptides (Ong and Mann, 2005). Many different technologies have been developed under the collective field of proteomics mass spectrometry (Aebersold and Mann, 2003). The focus in this paper is the more recent development (see Groseclose et al. (2007); Aoki et al. (2007)) of matrix assisted laser desorption ionization (MALDI) imaging mass spectrometry (IMS), also known as MALDI imaging, and in particular an analysis of MALDI-IMS data acquired from tissue samples of patients with ovarian cancer.
Unlike the more common 2D gel electrophoresis (2D-GE) and liquid chromatography mass spectrometry (LC-MS) techniques in proteomics, MALDI-IMS preserves the spatial distribution inherent in the tissue; and the tens of thousands of spatially distributed spectra acquired from a single tissue sample in a MALDI-IMS experiment provide new challenges for statisticians and bioinformaticians as well as having the potential to lead to breakthroughs in biological research (see Casadonte and Caprioli (2011)). We propose a combined cluster and feature extraction method for such data which exhibits cancer specific variables whose protein associations can be inferred by parallel LC-MS experiments such as those of Meding et al. (2012).
Standard proteomics mass spectrometry methods such as 2D-GE and LC-MS have been described and reviewed in the biological, bioinformatics and statistics literature for some time, see Wasinger et al. (1995). We briefly explain LC-MS and important differences with MALDI-IMS in Section 2. For an overview and review of recent approaches in LC-MS, see America and Cordewener (2008). Statistical challenges of proteomics mass spectrometry data are outlined in Wu et al. (2003). The statistics and bioinformatics literature on the analysis of 2D-GE and LC-MS data is growing fast and covers a range of statistical methods. Testing and classification of such data are described in Morris et al. (2005); Morris (2012); Yu et al. (2006) and references therein. Other statistical approaches that have been proposed and applied to 2D-GE and LC-MS data include peak identification, alignment, and feature selection, see Yu et al. (2006); identification of proteins, see Yu et al. (2006) and Karpievitch et al. (2010); wavelet-based methods, see Morris and Carroll (2006), America and Cordewener (2008) and references therein; and methods from survival analysis for the detection of differentially expressed proteins, see Tekwe et al. (2012).
Contrasting these developments in the analyses of 2D-GE and LC-MS data, the newer MALDI-IMS methods which have been introduced into routine research practice have not yet attracted as much attention in the statistics literature, although they have been covered in proteomics/mass spectrometry journals -see Norris et al. (2007); Jones et al. (2012); Gessel et al. (2014); Stone et al. (2012) and references therein. Recent research from Alexandrov and his collaborators focuses on clustering and dimension reduction of MALDI-IMS data; see Alexandrov et al. (2010); Alexandrov and Kobarg (2011); Alexandrov et al. (2013) and references therein.
The potential of MALDI-IMS is described in Alexandrov and Kobarg (2011): 'IMS is one of the most promising innovative measurement techniques in biochemistry which has proven its potential in discovery of new drugs and cancer biomarkers. ... IMS was used in numerous studies leading to understanding chemical composition and biological processes. ... As for many modern biochemical techniques, in particular in proteomics, the development of computational methods for IMS is lagging behind the technological progress.' Related to our research are the papers by Deininger et al. (2008); Alexandrov et al. (2010) and Bonnel et al. (2011) who cluster their MALDI-IMS data using principal component analysis and hierarchical clustering, or a Gaussian mixture model, as they find that clustering of the raw data does not yield biologically interpretable results. Our proposal differs from their research in a number of important aspects. Unlike these authors, we derive suitably binned binary data instead of working with the raw or intensity data (which we describe in Section 3), and following Koch (2013, Chapter 6), who demonstrated the success in finding biologically meaningful tissue clusters for such binary data, apply k-means clustering to the binary data. Analysis of the binary data has the added advantage of being computationally more efficient since the binary data are much smaller and parallel computations are possible for all the methods we propose. Further, our approach combines clustering with explicit feature extraction -based on our Difference in Proportions of Occurrence (DIPPS) statistic which ranks and selects the 'best' variables in a data-driven way and integrates the results in a single heat map that allows a natural interpretation. The ability to visualize the selected features as a single easily interpretable image has not been a part of the above-mentioned papers and gives a significant advantage to our approach. This paper is organized in the following way. We briefly describe relevant background on proteomics in Section 2, and discuss advantages of MALDI-IMS for the biological research fields. We describe MALDI-IMS data in Section 3, indicate how to interpret such data, and derive the binary data that we will analyze. Section 4 covers the spatial smoothing and clustering steps of our method. The DIPPS approach (including a feature extraction step) is described in Section 5. In Section 6 we consider the application of our combined cluster and DIPPS approach to multiple tissue sections from three surgically excised ovarian cancers of different patients. For each cancer and each tissue section, we find a set of features that characterizes the cancer region. The results shown in Section 6 successfully separate differences between patients from within-patient variation, and illustrate how the methods can be implemented and extended in order to address biologically relevant questions such as the identification of potential tissue-specific protein markers, and classification of patients based on their response to treatment.

Proteomics Background
Ovarian cancers are virtually asymptomatic and as a result the vast majority of cases are detected when the disease is in an advanced stage and metastasis has occurred. For patients with late stage ovarian cancer, radical surgery and chemotherapy are often insufficient to address the disease adequately and many patients relapse after treatment. The combination of late-stage diagnosis and unsuccessful treatments makes ovarian cancer the most lethal gynecological cancer with only 30% of patients surviving five years following diagnosis (Ricciardelli and Oehler, 2009;Jemal et al., 2011). Compounding these clinical statistics is the heterogeneity of ovarian cancers, which makes typical proteomic analyses difficult (Deininger et al., 2008). The key to addressing ovarian cancer will be our understanding of the mechanisms driving cancer progression, the identification of molecular markers which can predict treatment success and most importantly, the identification of new treatment targets. As proteins are the main functional components of cells and tissues, determining protein distributions in cancer tissue represents a crucial step in understanding cellular dysfunction such as cancer in a holistic way.
Proteins are functional and structural components of cells. They are synthesized within cells as linear amino acid sequences and subsequently fold into more complex 3D structures that determine their functions and intra-cellular locations. The complete set of proteins which exist in a given cell, tissue or biological fluid, under defined conditions, is termed its proteome (Wilkins et al., 1996). Proteomes vary considerably between different cellular states. Understanding these differences allows for insight into the cellular processes involved in development and progression of diseases such as cancer. Proteomics is the discipline which characterizes proteome changes. Proteomic identification and quantitation, typically achieved through MS, requires the reduction of proteome complexity prior to analysis. This is often achieved with fractionation methods such as one-or two-dimensional gel electrophoresis (Gygi et al., 2000) or liquid chromatography (LC). LC is currently the predominately used fractionation technique in proteomics (Rogowska-Wrzesinska et al., 2013). In this approach molecules are affinity bound to a column material, or stationary phase, and eluted with a solvent, or mobile phase, the composition of which is usually changed as a function of time.
LC-MS analysis of a proteome uses either a top-down or bottom-up approach. Top-down approaches analyze intact proteins, whereas bottom-up approaches use proteolytic enzymes such as trypsin to digest proteins into smaller peptides prior to analysis. We are primarily interested in the bottom-up approach. In a bottom-up LC experiment peptides dissolved in a hydrophilic mobile phase are retained using hydrophobic interactions with the stationary phase. The peptides are then eluted by a linear increase in the percentage of hydrophobic solvent in the mobile phase over time. A peptide will be retained on the hydrophobic stationary phase until the peptide's affinity for the percentage organic solvent in the mobile phase exceeds its affinity for the stationary phase. The LC eluant is often directly coupled to an MS instrument, and the technique is then called LC-MS.
MS instruments contain an ion source, mass analyzer and detector. The ion source for LC-MS usually employs electro-spray ionization (ESI), but can instead use matrix assisted laser desorption ionization (MALDI) to produce gaseous peptide ions for analysis. The mass to charge ratio (m/z) of these ions is measured by the mass analyzer and detector to produce a mass spectrum, that is, measured intensity as a function of m/z. Importantly, the eluant flow from an LC system is constant, leading to the acquisition of a number of MS measurements over time. Each mass spectrum acquired in a LC-MS experiment is identified by its retention time. The peptides produced in a bottom-up proteolytic digest can also be fragmented to produce spectra which correspond to amino acid composition. Identification of these spectra is attempted using search algorithms, such as MASCOT (for more information see Koenig et al. (2008) and references therein), which match measured spectra to expected fragmentation spectra. It is the combination of LC-MS with bottom-up strategies that currently allows proteomic studies to identify and quantify thousands of unique proteins in a given biological sample.
Despite its popularity, LC-MS suffers from two method-specific limitations: 1. The length of an LC-MS run is usually measured in hours, precluding the rapid analysis of large sample numbers (≥ 20), and 2. Tissue samples used for LC-MS based proteomics are homogenized and solubilized, which removes all spatial information inherent in the tissue.
Tissues are usually a mix of different cell populations, the organization of which is directly related to their functions and hence the spatial information is essential. Different cell populations or tissue types are typically localized in distinct spatial regions, as can be seen in Figure 1(a). The loss of spatial information that occurs in LC-MS analysis has motivated the development of direct tissue analysis using MALDI-IMS Groseclose et al., 2008;Gustafsson et al., 2011).
To prepare a sample for MALDI-IMS a tissue block is thinly sectioned (2 -10µm thick slices) and mounted onto conductive microscopy slides. For an analysis of peptides, the tissue section is coated with a homogeneous layer of proteolytic enzyme and digested. The digest is overlaid with a low molecular weight organic matrix that co-crystallizes with the tissue-derived peptides. MALDI is used to ionize peptides directly from the tissue and to acquire a mass spectrum from each (x, y) position on a regular grid across the tissue section. A spacing of 20 -250µm between acquisition locations is typical and balances the requirements for high quality MS data, practical data size, and measurement time. In the ovarian cancer data a spacing of 100µm was used.
In LC-MS spectra are identified by their retention time. In contrast, MALDI-IMS spectra are curated by their spatial (x, y) coordinates, and the availability of spatial coordinates offers a unique perspective on tissue analysis. Typically, tissue structure is visualized by pathologists using haematoxylin and eosin, H&E, staining followed by light microscopy. Since tissue sections can be H&E stained post-acquisition of MALDI-IMS data (Deutskens et al., 2011), spatial changes in protein abundance measured by MALDI-IMS can be compared to the histology visible in the H&E stain of the same tissue (see Figure 1).
MALDI-IMS provides unprecedented capacities for pathology, and as a consequence has important implications for research into human cancers (Gorzolka and Walch, 2014). Pathology provides an in-depth understanding of the way in which cellular structure and appearance relate to cancer behavior and the cancer's effect on the patient. LC-MS based proteomics is able to establish identification and quantitation for large portions of the cancer proteome. MALDI-IMS, however, can combine these capacities in the following ways: 1. MALDI-IMS acquires spatially referenced data directly from cancerous tissue and has the potential to provide novel information in the form of tissue or treatment specific markers. This is not achievable with current pathology and could potentially allow clinicians to use MALDI-IMS data for individualizing patient treatments.
2. Cancers, and in particular ovarian cancers, are heterogeneous tissues (Deininger et al., 2008). Spatial analysis by MALDI-IMS allows this heterogeneity to be addressed explicitly (Gorzolka and Walch, 2014). Furthermore, if arrays of patient tumor cores are used, hundreds of patients can be analyzed in a single day. MALDI-IMS therefore provides the opportunity to screen cancer samples rapidly, target regions of interest and employ complementary LC-MS methods to characterize the tumor proteome of these tissues.
One of the key advantages of MALDI-IMS is its ability to analyze large cohorts of patients in a timely manner, for example by using a tissue micro-array or TMA (Groseclose et al., 2008;Steurer et al., 2013). To exploit this ability in future, it will be crucial to understand differences and to be able to distinguish between patients. As proof of concept, we consider data including multiple tissue types from each of three patients in this paper.

The Data and Binary Binning
In this section we describe the form of MALDI-IMS data. We also describe the steps used in transforming the initial functional data to the final 'binary binned data'. We will refer to the 'binary binned data' as 'the data' beyond this section. As discussed in Section 2, a MALDI-IMS dataset consists of mass spectra annotated by spatial (x, y) coordinates. Each mass spectrum typically contains between 5 and 200 peaks, each corresponding to a biological analyte such as a peptide or protein. A number of preprocessing steps are involved in extracting the peaks in the functional data, including: smoothing, baseline reduction, and peak picking (typically using proprietary software, e.g. flexAnalysis, Bruker Daltonik). We will only be concerned with these extracted peaks from here on.
In this paper we consider a selection of datasets from three surgically excised ovarian cancers, each from a different patient -see Gustafsson (2012). Thin crosssections of tissue are obtained from each cancer, and MALDI-IMS data are collected. We will use the term 'dataset' to refer to the data collected from a single section of tissue. Multiple sections are taken from each cancer. We expect multiple datasets from the same cancer to be similar, and aim to separate within-patient from betweenpatient variability. Initially we consider a single typifying dataset (shown in Figure 1) and will refer to it as the 'motivating dataset' throughout. A typical dataset might have 5000−100000 spectra, the motivating dataset consists of 14059. The motivating dataset will serve to illustrate our proposed method which we describe in Sections 4 and 5. In Section 6 we consider nine datasets (including the motivating dataset), three from each cancer, and compare datasets within and between patients. Although we have access to more datasets we only show nine as these suffice for comparisons both within and between patients, while still illustrating our method concisely.
As described in Section 1, each dataset contains mass spectra collected at regularly spaced points on the surface of a thin tissue section. Within a dataset, each mass spectrum is a realization of ion intensity as a function of mass-to-charge ratio (m/z). When using MALDI-MS the charge (z) is typically one and so the mass-tocharge ratio reduces to mass. Mass peaks (local maxima) can be identified and the original mass spectra reduced to 'peak list' form. This was done using the standard flexAnalysis software package (Bruker Daltonics, Bremen, Germany). Representing the data as peak lists significantly reduces the amount of data involved, often by as much as two orders of magnitude and this can be very important due to the amount of data involved (> 10GB).
The m/z locations of the peaks are discretized by binning. We call the data after binning 'binned intensity data' as in this form the variables are measures of ion intensity for each of the m/z bins. The number of variables in the binned intensity data will vary depending on the choice of bin size which turns out not to be crucial provided the chosen bin size is within a reasonable range (0.05 − 2 m/z). We will use an intermediate bin size of 0.25 m/z, which yields 4298 variables (non-empty bins) in the motivating dataset. We repeated all analyses with a number of alternate bin sizes and consistently produced very similar results, given the bin size is chosen from the above interval. The principal effect that choice of bin size has is on the total number of variables, and the number that are removed in the dimension reducing steps of our method -smaller bin sizes will result in more variables initially, and more being removed, larger bin sizes result in fewer variables initially, and fewer being removed.
Experience with exploratory analysis including k-means clustering of the binned intensity data has shown that these data do not lead to interpretable or spatially localized clusters (Koch, 2013). This is due to the effect extraneous variables (such as matrix crystal morphology) have on the intensity measurements (Garden and Sweedler, 2000). Koch (2013, Example 6.12) proposes to transform the intensity data to binary data. There exist other approaches to addressing these effects (such as the normalization approaches of Deininger et al. (2011)), but the approach of Koch (2013) bypasses the effect of the extraneous variables on peak intensities entirely. Dealing with the binary data has the added advantage of being extremely computationally efficient. The transformation to binary data is natural, as it tells us if a peak is present (takes the value one), or is absent (takes the value zero). Representing the variables in this binary form can be considered a data cleaning step, as it removes superfluous sources of noise. Subsequent analyses will refer only to the binary data, and we will refer to these simply as 'the data' from here on.

Clustering Spatially Smoothed Data
We propose a two-step method for separating spectra that belong to different tissue types. First a smoothing step, which serves dual purposes: to act as a data cleaning method, and to reduce the dimension. Second a clustering step, to separate the spectra into groups that correspond to tissue types.

Smoothing Step
The spatial information available in MALDI-IMS data can be used to clean the data and remove variables that are spatially dispersed. We incorporate this spatial information through a spatial smooth.
Let X be a d × n data matrix; the rows of X correspond to variables (m/z bins) and are denoted x i• , the columns of X correspond to spectra and are denoted x •j , and the entries of X are denoted x ij . Let 0 ≤ τ < 1 2 be a smoothing parameter and δ ≥ 0 a distance cutoff.
We iteratively update the values of X. Let X (k) denote the updated matrix at the kth iteration. Similarly let x ij denote the rows, columns, and values of X (k) respectively. At the kth iteration, the proportion of spectra, T (k) ij , in a spatial δ-neighborhood of the jth spectrum x (k) •j whose values at the ith variable x This proportion determines if the value x (k) ij should be smoothed. In (1), d j denotes a 1 × n indicator vector with entry 1 if the corresponding indexed spectrum is in a δ-neighborhood of x j• and zero otherwise. If this proportion T (k) ij is less than τ , we update the value x (k) ij as in (2). We generate the smoothed data by iteratively calculating the entries of X (k) for k = 1, 2, . . . , starting with X (0) = X. The T (k) ij term appearing in (2) is defined in (1). We stop when convergence is reached (k = k * = min{k : X (k) = X (k−1) }). The spatially smoothed data are X (k * ) .
Remarks on the smoothing process: 1. We choose δ = √ 2 which results in a range 1 Moore neighborhood. See Gray (2003). This neighborhood is used in the cellular automata literature including Conway (1970). It is worth noting that acquiring the range 1 Moore neighborhood by using the Euclidean distance and δ = √ 2 on a regular grid is equivalent to using the Tchebychev distance, and δ = 1 (without the regular grid assumption).
2. The smoothing parameter τ defines the proportion of neighboring spectra needed to agree in order for the value for an observation to remain unchanged at any given step. Small values of τ smooth less (τ = 0 leaves the data unmodified), while larger values smooth more. The limit τ → 1 2 results in maximum smoothing, and is equivalent to the intuitive median smooth. The median smooth tends to yield over-smoothed data, and often fails to converge. We chose an intermediate smoothing parameter, τ = 1 4 , for these analyses. 3. Alternative smoothing options include kernel methods (Wand and Jones, 1995) which apply to continuous data. These methods produce continuous values when applied to binary data, for which there is no clear interpretation. Our method produces binary smoothed data -maintaining the interpretability of the binary values.
4. At each smoothing iteration k, variables are smoothed independently, and within each variable all observations are smoothed simultaneously at each step. This means that it is possible to parallelize the smoothing algorithm, making efficient use of computational resources.
Variables that exhibit occurrence of peaks in a small number of spatially delocalized spectra or in almost all spectra constitute a large proportion of all variables. These variables tend not to be relevant as they are usually internal calibrants , errors, contaminants, or features that are too small (of the order 100µm) to be of interest. These variables have zero variance after the smoothing step. Following the smoothing, these zero-variance variables are removed -reducing the dimension of the data. Any such variables are not of interest as they are internal calibrants , errors, contaminants, or features that are too small (of the order 100µm) to be of interest. The motivating dataset has 4298 variables before the smoothing step. After the smoothing step 818 of these variables have non-zero variance and the remainder are discarded.

Clustering Step
The second step in our approach concerns clustering of the spatially smoothed data. We use k-means clustering with the cosine distance. Based on the information available from the histology, there are three broad tissue types present which could be labeled as fatty, connective and cancer tissue respectively. Further, there are spectra that were acquired off-tissue. Thus we perform k-means clustering with k = 4. We choose initial cluster centers as this results in the algorithm converging to a deterministic solution and allows results to be reproduced.
Even after the dimension-reducing smoothing step, the data are typically still high-dimensional. In the motivating dataset 818 variables are retained after discarding zero-variance variables at the end of the smoothing step. A natural measure of distance between high-dimensional vectors is the cosine distance -rather than the Euclidean distance, the common default in many k-means algorithms. Because we are dealing with binary data, we also repeated our analyses with the Hamming distance, which has an intuitive interpretation for binary data. The Hamming distance produces very similar results for the smoothed data, but shows slightly less stability for the unsmoothed data and so we prefer to use the cosine distance and only show the corresponding results.
For data with associated spatial information (such as MALDI-IMS data) it is natural to display the cluster membership in the form of cluster maps or cluster images: colored pixels at the (x, y) coordinates of the spectra which show the cluster membership of each spectrum using different colors to identify clusters. The H&E stained tissue cross-section and result of 4-means cluster analysis (by cosine distance) of the motivating dataset are shown in Figure 1(a) and Figure 1(b) respectively. The cluster membership in Figure 1(b) corresponds well with tissue types as determined by the histology in Figure 1(a): cyan corresponds to fatty tissue, green to cancerassociated connective tissue and orange to cancer. The cancer cluster accurately assigned the location of four cancer nodules (indicated by arrows in Figure 1(a)). The fourth cluster in pink corresponds well with off-tissue spectra. The correspondence between cluster results and histology demonstrates that the spatial smooth and cluster analysis isolate key molecular information which allows differentiation of tissue types by their mass spectra.

The Difference in Proportions of Occurrence (DIPPS) Approach
In Section 4 good agreement was demonstrated between tissue types visible in the histology of Figure 1(a) and clusters from the 4-means clustering of Figure 1(b). From the biological perspective it is of great interest to be able to quantify the differences between these tissue types in an easily interpretable way. At a mathematical level, a characterization of the differences between the tissue types translates into an identification of variables that discriminate between them, or more generally between subsets of data. We explain our approach, the DIPPS approach, for an arbitrary subset of observations with binary variables. We then apply these ideas to the motivating dataset using spectra from the cancer cluster as the subset of interest, since distinguishing these spectra from the non-cancer spectra is particularly interesting. In this section we describe the DIPPS approach which consists of 1. a ranking of the variables for their ability to characterize a subset of data; 2. a data-driven selection of the 'best' variables; and 3. a DIPPS map which combines the best variables in a single image that has a natural interpretation.

The DIPPS Approach for a General Subset of Observations
We define the DIPPS statistic for a fixed subset of the data in (3), and show how this new statistic introduces a natural ranking of the variables based on their ability to characterize the subset of interest. Further, this statistic leads to a data-driven selection of the number of ranked variables that characterize the subset. The extraction of the most highly ranked and hence most relevant variables can be interpreted as a dimension reduction step, as it effectively restricts the variables to the extracted ones. Further, we propose a way of displaying graphically the information obtained from the statistic in an easily interpretable summary image. These maps make this technique useful in exploratory analyses as combining the selected variables into a single interpretable image allows for broad conclusions to be quickly and easily drawn. When the amount of data becomes large the approach commonly used in proteomics, namely of considering each variable individually, is of limited use due to the sheer amount of information that needs to be considered simultaneously, and the ability to produce an image that summarizes all selected variables becomes particularly useful. Let S be a subset of observations of the data, i.e. a subset of the columns of X. We let p(S) denote the mean of the observations in S, and similarly let p(S c ) denote the mean of the observations in its compliment S c . Entries of the vector p(S) are the proportions of observations in S that take the value one for each variable. The interpretation of binary data as the occurrence of an event -here existence of peaks -allows the mean p(S) to be interpreted as the proportions of occurrence (for each variable) in S. These entries can be considered as the 'sensitivity' or true positive rate of each variable for S. Similarly, entries of the vector 1 d×1 − p(S c ) are the proportions of observations in S c that take the value zero for each variable and can be considered the 'specificity' or true negative rate of each variable for S. In order to characterize S, both sensitivity and specificity should be high. We sum these measures of sensitivity and specificity, and subtract one to give a range of [−1, 1] and then define the vector of DIPPS statistics for S as (3) For convenience of notation we omit the dependence on S, and write d(S) = d. We will similarly omit the S dependence for t a , n a , c, and a * below, as we are treating S as fixed. We use the entries of the vector d to rank the variables: the entry with the greatest value corresponds to the individual variable that characterizes S best. Next we determine the set of variables that collectively best characterize S. We do this by finding a cutoff value in a data-driven way as follows. For a > 0 let t a be the d × 1 vector that takes the value one if the corresponding element of d is greater than a, and is zero otherwise, and let n a be the number of ones in the vector t a . Let c be the centroid of S, calculated using the same distance, D, as in Section 4.2 -so in our case the cosine distance. Our suggested cutoff a * is where D (c, t a ) denotes the distance between c and t a . The n a * variables that best characterize S are those whose entries are one in t a * . The sum of these variables for each spectrum can then be visualized spatially in a DIPPS map.

The DIPPS Approach for the Cancer Cluster
Now we consider the DIPPS approach applied to the motivating dataset. We are particularly interested in characterizing the cancer spectra in these data, and we thus choose S to be the set of spectra belonging to the cancer cluster shown in Figure 1(b). The cutoff of (4) results in n a * = 45 variables at a * = 0.139. These variables characterize the cancer tissue, and thus warrant closer inspection. Followup analysis can be done to identify the peptides that these variables correspond to, and to draw inference as to their parent proteins. To illustrate the biological significance of these variables, we compared their m/z values to the previously published results of Gustafsson (2012). Two of the highly ranked variables (listed in Table 1, along with inferred parent proteins) were previously identified as highly expressed in cancer by manual assessment of variable distribution in datasets from this patient. The identifications were achieved through both mass matching to LC-MS/MS as well as in situ MS/MS (Gustafsson, 2012). The identity and histological distribution of these variables were successfully validated using immuno-histochemistry (unpublished results). This illustrates that the DIPPS approach can find features of known importance. Crucially, the DIPPS approach produces a list of characterizing variables more rapidly and comprehensively than a manual inspection of individual variables. An important aspect of our proposed method is the visualization of the selected variables in a single summary image which is easy to interpret, and which we call a DIPPS map. The ability to summarize results in this way is particularly attractive as considering so many variables individually can be time consuming and fail to provide an 'big-picture' perspective. We construct the DIPPS map pointwise at each (x, y) location by counting the number of extracted variables that exhibit occurrence for the spectrum at (x, y), and by visualizing this count in heat colours: cold indicating a count of zero, hot indicating that all exhibit occurrence. Figure 2(b) shows the DIPPS map for the cancer cluster.
A DIPPS map such as that shown in Figure 2(b), and the variables that define it have intuitive interpretations that the results of cluster analysis do not. This method allows variables crucial to cluster/tissue differentiation to be isolated and visualized. The visualization of such DIPPS maps becomes increasingly important in exploratory analyses when the number of patients and datasets increases as it becomes quickly infeasible to consider each of the selected variables separately. The DIPPS map highlights gradations which reveal finer detail and provides a more representative visualization of the data that is not possible in cluster maps. We discuss this point further in Section 6. Our DIPPS approach also acts as a variable reduction technique by further reducing the data from 818 to 45 variables. More importantly it successfully separates tissue type specific variables, addressing the heterogeneity of the tissue. This facilitates the comparison of datasets, which we discuss in Section 6.
6 Use of DIPPS and the Jaccard Distance in Differential Comparison of Multiple Datasets In this section we consider nine datasets, including the motivating dataset. Of these datasets three are from each of three different patients which we will refer to as patients A, B, and C, respectively. We will refer to the three datasets from patient A as A1, A2, and A3, and similarly for the datasets from patients B and C. Dataset A2 is the motivating dataset. Each of the three datasets arising from the same patient are acquired from thin (6µm) tissue sections of a single surgically excised tissue. Because of this experimental setup, we expect the cluster and DIPPS maps of datasets from the same cancer to be similar in terms of the location of the cancer clusters and the selected variables. The DIPPS map intuitively combines the variables that characterize cancer clusters, making the interpretation of large patient cohorts much easier.
The results of our analyses of the nine datasets are displayed in Figures ??, ??, and ?? for patients A, B, and C, respectively. Each figure corresponds to one cancer and shows the three datasets in rows. Each row consists of an H&E stain on the left, a cluster map in the center, and a DIPPS map (for the cancer cluster) on the right. The second row of Figure 3 repeats the results for the motivating dataset shown in Figures 1 and 2.
In all nine datasets, the clustering results correspond well with the underlying tissue morphology. Minor exceptions to this relate to dataset specific misassignment of tissue types. These exceptions do however highlight the ability of the DIPPS maps to find and extract information in the data that is not available in the cluster maps, as the DIPPS maps overcome these difficulties despite slight misalignment between cluster membership and tissue morphology. For an example of this robustness property of the DIPPS maps, consider Figure 4 -although the cluster maps show a noticeable difference in cluster shapes across the datasets from patient B, the DIPPS maps are comparatively consistent.
The DIPPS maps are also more representative of the data as they allow for gradual changes and fine detail to be seen that cannot be represented by the 'hard' boundaries in a cluster map. For example consider Figure 5, particularly C2 and C3. The DIPPS maps for these two datasets are much more representative of the cancer regions (compared with the cluster membership), as these maps indicate that the connective tissue shares variables with the cancer tissue. This is to be expected, given that the connective tissue is influenced by the presence of the cancer. However, the cancer tissue contains a higher proportion of those key cancer-specific variables, leading to a gradual change in color in the DIPPS map corresponding to fewer of these key variables occurring further away from the cancer. These gradual changes are more intuitive and informative than the cluster membership which has 'hard' boundaries and hence does not represent the data as accurately.
In addition to producing the DIPPS maps, the DIPPS approach yields the actual variables which are representative of a chosen subset. This property can be used to compare tissue types across several datasets. We calculate the DIPPS statistic and select characterizing variables as described in Section 5 for the cancer cluster in each of the nine datasets. To determine how similar these characterizing variables are, we could compare the list of variables but as the number of datasets increases, pairwise comparisons become burdensome. We use the Jaccard distance, see Jaccard (1901), for capturing similarity between subsets of observations. Other criteria could be used to calculate similarity indices, such as the χ 2 -based distances used in Koch (2013, Example 8.11) but for comparing subset overlap the Jaccard distance is a natural choice. For a pair of sets S i , S j , the Jaccard distance is We calculate the Jaccard distance for pairs of sets of characterizing variables for the cancer cluster of all nine datasets. Figure 6 shows these distances. The 'colors' range from black for a Jaccard distance of zero (indicating equality of sets) to white for a Jaccard distance of one which indicates disjoint sets. The datasets are displayed from left to right (or top to bottom) starting with those of patient A, followed by those of patient B, and finally the datasets of patient C. The block diagonal of Figure 6 is notably dark -this illustrates that differences within patients are smaller than differences between patients. This difference allows variables that are patient-specific in their characterization of the cancer to be separated from variables that are cancer specific across multiple patients. There are a number of variables that characterize cancer in all 9 of these datasets -including the m/z bin centered at 1628.75, mentioned in Section 5. The identity of this variable (listed in Table 1) has already been validated using two independent methods (in situ MS/MS and immuno-histochemistry), indicating that the protein from which this peptide derives is highly expressed in the cancer tissue of these patients. This protein could therefore be investigated further as a marker for ovarian cancer metastasis in a larger patient cohort.
On the other hand, the variables that are patient-specific in their characterization of cancer could be further investigated for their ability to classify cancers according to clinical or diagnostic criteria such as response to treatment. There appears to be a similarity between patients A and C not shared by B, as can be seen in Figure 6. This discrepancy could be a consequence of the relatively larger number of variables that characterize the cancer of patient B.' The discrepancy in variable composition between patients A/C and B can be explained by the large amount of necrotic tissue in the sections from patient B. Mass spectra from necrotic tissue are expected to be vastly different from those of other cancer tissue, which would indicate that the patient-specific variables in patient B are likely to be markers for necrotic tissue.
Other tissue types could be considered (fatty and connective tissue), and a figure similar to Figure 6 could be produced for each of them. Such figures do not show a notably darker block-diagonal in the way that Figure 6 does. This tells us that the main differences between patients are in the patients' cancer, rather than in their other tissues, and reenforces how crucial it is to address the heterogeneity of these data by separating tissue types before conducting comparisons between patients.

Conclusion
This paper proposes an integrated approach to clustering and feature extraction for functional and high-dimensional data. The research is motivated by and applied to proteomics MALDI-IMS data obtained from tissue samples and is part of a comprehensive investigation into ovarian cancer. However, our approach is generic and applies to functional or high-dimensional data with an underlying spatial structure that belong to two or more groups.
The combined cluster and feature extraction approach is based on our difference in proportions (DIPPS) statistic and includes novel visualizations which enhance the cluster maps. For the MALDI-IMS cancer data, these maps have a natural interpretation in terms of the variables that characterize cancer tissue. In other applications, the DIPPS maps overcome the limitations of single-variable visualizations and cluster maps.
A study of datasets from multiple patients shows that the DIPPS approach is able to separate variables that are patient-specific in their characterization of the cancer from variables that are cancer specific across multiple patients. This suggests that use of MALDI-IMS could lead to a differentiation of cancer types or responses of patients to potential treatments in future research. We only observed such a separation of variables for cancer tissue and not for other tissue types, which reinforced the need to address the heterogeneity of tissue data in a cluster analysis prior to feature extraction.
In proteomics, the ability to automate feature extraction and to present these features as DIPPS maps provides an opportunity for holistic appraisal of MALDI-IMS data. This is crucial due to the size and number of such datasets and their high-dimensional nature. By isolating variables important to specific tissue types and reporting similarities across patients, it will be easier to identify variables for further validation as tissue markers and to build models that could assign tissue types in new datasets by classification. Isolating variables important to specific patients and reporting differences between patients will make it easier to identify variables for further validation as predictors for addressing clinical questions such as predicting chemotherapy response and patient survival.
Yu, W., B. Wu, T. Huang, X. Li, K. Williams, and H. Zhao (2006). Statistical methods in proteomics. In Springer Handbook of Engineering Statistics, pp. 623-638. Springer. Figure 1: Left, (a) shows the H&E stained tissue section, arrows indicate the four cancer nodules. Right, (b) shows pixels plotted in their relative spatial locations, color identifies the cluster membership produced by 4-means clustering of the binary spatially smoothed data.  , and DIPPS maps for the cancer cluster (right column) for the three datasets A1, A2, and A3 (each represented in a row). Note that n a * = 43, 45, 61 variables are visualized in the DIPPS map for the datasets A1, A2, and A3, respectively.
Figure 5: H&E stains (left column), cluster maps (centre column), and DIPPS maps for the cancer cluster (right column) for the three datasets C1, C2, and C3 (each represented in a row). Note that n a * = 75, 39, 53 variables are visualized in the DIPPS map for the datasets C1, C2, and C3, respectively. Figure 6: The 9 × 9 grid shows pairwise Jaccard distances between the sets of variables characterizing the cancer clusters. Rows (and columns) correspond to the 9 datasets. The colour of each pixel indicates the value of the Jaccard distance -white corresponding to a value of one, black to zero.