text:
# Photometric Mapping with ISOPHOT using the “P32” Astronomical Observation Template Based on observations with ISO, an ESA project with instruments funded by ESA member states (especially the P/I countries France, Germany, the Netherlands and the United Kingdom) with participation of ISAS and NASA.

## 1 Introduction

From the point of view of signal processing and photometry, diffraction-limited mapping in the FIR with cryogenic space observatories equipped with photoconductor detectors poses a particular challenge. In this wavelength regime the number of pixels in detector arrays is limited in comparison with that in mid- and near-IR detectors. This means that more repointings are needed to map structures spanning a given number of resolution elements. Due to the logistical constraints imposed by the limited operational lifetime of a cryogenic mission, this inevitably leads to the problem that the time scale for modulation of illumination on the detector pixels becomes smaller than the characteristic transient response timescale of the detectors to steps in illumination. The latter timescale can reach minutes for the low levels of illumination encountered on board a space observatory. The result is that the signals from the detectors can depend as much on the illumination history as on the instantaneous illumination. Unless corrected for, the transient response behaviour of the detectors lead to distortions in images, as well as to systematic errors in the photometry, especially for discrete sources appearing on the maps. In general, these artifacts become more long lived and more difficult to correct for at fainter levels of illumination, since the transient response timescales increase with decreasing illumination.

For mapping observations in the FIR, the Infrared Space Observatory (ISO; Kessler et al. 1996) was equipped with two Gallium-doped Germanium photoconductor detectors which comprised part of the ISOPHOT instrument (Lemke et al. 1996). The ISOPHOT C200 detector was a 2$`\times `$2 (stressed) pixel array operating in the 100$``$200$`\mu \mathrm{m}`$ range, and the ISOPHOT C100 detector was a 3$`\times `$3 pixel array operating in the 50$``$100$`\mu \mathrm{m}`$ range. Compared to the IRAS survey detectors, the ISOPHOT-C detectors had relatively small pixels designed to provide near diffraction limiting imaging. ISOPHOT thus generally encountered larger contrasts in illumination between source and background than IRAS did, making the artifacts from the transient response more prominent, particularly for fields with bright compact sources and faint backgrounds.

A further difficulty specific to mapping in the FIR with ISO was that, unlike IRAS, the satellite had no possibility to cover a target field in a controlled raster slew mode. This limited the field size that could be mapped using the spacecraft raster pointing mode alone, since the minimum time interval between the satellite fine pointings used in this mode was around 8 s. This often greatly exceeded the nominal exposure time needed to reach a required level of sensitivity (or even for many fields the confusion limit). Furthermore, the angular sampling and redundancy achievable using the fine pointing mode in the available time was often quite limited, so that compromises sometimes had to be made to adequately extend the map onto the background.

A specific operational mode for ISO - the “P32” Astronomical Observation Template (AOT) - was developed for the ISOPHOT instrument to alleviate these effects (Heinrichsen et al. 1997). This mode employed a combination of standard spacecraft repointings and rapid oversampled scans using the focal plane chopper. The technique could achieve a Nyquist sampling on map areas of sky ranging up to 45$`\times `$45 arcmin in extent (ca. 70$`\times `$70 FWHM resolution elements) on timescales of no more than a few hours. In addition to mapping large sources, the P32 AOT was extensively used to observe very faint compact sources where the improved sky sampling and redundancy alleviated the effects of confusion and glitching.

In all, over 6$`\%`$ of the observing time of ISO was devoted to P32 observations during the 1995-1998 mission, but the mode could not until now be fully exploited scientifically due to the lack of a means of correcting for the complex non-linear response behaviour of the Ge:Ga detectors. Here we describe the basic concept of an algorithm which can successfully correct for the transient response artifacts in P32 observations. This algorithm forms the kernel of the “P32Tools” package, which is now publically available as part of the ISOPHOT Interactive Analysis package PIA (Gabriel et. al 1997; Gabriel & Acosta-Pulido 1999). Information on the algorithm, as well as the first scientific applications, can also be found in Tuffs et al. (2002a,b) and Tuffs & Gabriel (2002). The user interface of P32Tools that connects the P32 algorithm to PIA is described by Lu et al. (2002). After a brief overview of relevant aspects of the P32 AOT in Sect. 2, we describe the semi-empirical model used to reproduce the transient response behaviour of the PHT-C detectors in Sect. 3. Sect. 4 describes the algorithms used to correct data. Examples demonstrating the photometric and imaging performance of ISOPHOT P32 observations are given in Sect. 5, based on maps corrected using the P32 algorithm.

## 2 The P32 AOT

The basic concept of the P32 AOT, in which imaging was achieved using a combination of standard spacecraft repointings and rapid oversampled scans using ISOPHOT’s focal plane chopper, is summarised in Fig. 1.

There were two basic requirements, outlined in the original proposal for the AOT (Tuffs & Chini 1990) which led to the final design: The first goal was to give ISOPHOT the capability of achieving a Nyquist sky sampling ($`\mathrm{\Delta }\mathrm{\Theta }`$=$`\lambda `$/2D = 17 arcsec at 100$`\mu \mathrm{m}`$) on extended sources subtending up to 50 resolution elements in a tractable observation time. As well as providing a unique representation of any arbitrary sky brightnesss distribution, this also improved the effective confusion limit in deep observations, which, in directions of low background was $`100`$ and 20 mJy rms at 200 and 100$`\mu \mathrm{m}`$, respectively, for ISOPHOT-C.

The second requirement was to achieve a redundancy in the observed data. This led to the sampling of each of a given set of sky directions by each detector pixel at different times. The redundancy was achieved in part by making several (at least 4) chopper sweeps at each spacecraft pointing direction. Further redundancy was achieved by having an overlap in sky coverage of at least half the chopper sweep amplitude between chopper sweeps made at successive spacecraft pointing directions in the Y spacecraft coordinate. Thus, a given sky direction was sampled on three different time scales by each detector pixel: on intervals of the detector non-destructive read interval, on intervals of the chopper sweep, and on intervals of the spacecraft raster pointings. Although the redundancy requirement was originally made to help follow the long term trends in detector responsivity during an observation, in practice it also proved very useful in removing effects caused by cosmic particle impacts (Gabriel & Acosta-Pulido 2000), the so-called “glitches”, from the data. This was especially important to optimise the sensitivity of deep observations using the C100 detector, which are generally limited by glitching rather than by confusion.

Both these requirements led to dwell times per chopper plateau which could be as low as 0.1 s. For a given sky brightness, this also had the consequence that the detector was generally read faster than in the other AOTs using the PHT-C detectors, which allows a higher time resolution in the deglitching procedures. A further beneficial side effect of the P32 AOT was that because the transient response timescale of the detectors is longer than the dwell time of the chopper, even for bright sources, the effective limiting source brightness that could be observed without saturating the detector was higher than for the other AOTs.

### 2.1 The P32 “natural grid”

The fact that the spacecraft pointing increment in Y was constrained to be a multiple of the chopper step interval means that a rectangular grid of directions on the sky can be defined such that all data taken during P32 measurement fine pointings will, to within the precision of the fine pointings, exactly fall onto the directions of the grid. This so-called P32 “natural grid” allows maps to be made without any gridding function, thus optimising angular resolution. The P32 natural grid is also a basic building block of the algorithm to correct for the transient response behaviour of the detector, as the algorithm solves for the intrinsic sky brightness in M$`\times `$N independent variables for a grid dimension of M$`\times `$N. Here M and N are the number of columns and rows on the P32 natural grid. An example of the almost perfect registration of data onto the grid is shown in Fig. 2.

All the example maps shown in this paper will be sampled without a gridding function on the P32 natural grid. This means that the data in each of the map pixels is independent. This is unlike the situation for most other maps calculated using PIA, where a gridding function is employed and so raw data can contribute to more than one map pixel.

## 3 Transient response of the PHT-C detectors

Gallium doped germanium photoconductor detectors exhibit a number of effects under low background conditions which create severe calibration uncertainties. The ISOPHOT C200 (Ge:Ga stressed) and especially the C100

(Ge:Ga) detectors have a complex non-linear response as a function of illumination history on timescales of $`0.1100`$ sec, which depends on the absolute illumination as well as the changes in illumination. The following behaviour is seen in response to illumination steps (see also Fouks 1992, Acosta-Pulido et al. 2000 and Laureijs et al. 2000, Sect 4.2.2.):

* a rapid jump in signal, which however undershoots the asymptotic equilibrium response (often severely for the C100 array)
* a hook response (overshooting or undershooting on intermediate timescales of $`110`$ s)
* a slow convergence to the asymptotic equilibrium response of the detector on timescales of $`10100`$ s

An example showing these effects for a staring observation is given in Fig. 3. Also seen are points discarded by the deglitching algorithm in PIA, marked by stars. Many glitches are accompanied by a longer lived “tail”. It is these tails that determine the ultimate sensitivity of the ISOPHOT C100 array.

### 3.1 Semi-empirical detector model

Unlike the ISOCAM and ISOPHOT-S Si:Ge detectors on board ISO, which operate in the mid-IR, the transient response behaviour of the Ge:Ga detectors of ISOPHOT-C are not described by the analytical Fouks Schubert model (Fouks & Schubert 1995; Schubert et al. 1995; Acosta-Pulido 1998; see Coulais & Abergel 2000 for its application to ISOCAM). Although numerical simulations of the transient behaviour of Ge:Ga detectors have been made (Hagel et al. 1996), these do not cover the wide range of operating conditions experienced by the ISOPHOT detectors. However, the effects observed in ISOPHOT-C can to some extent be qualitatively understood in terms of the build up of charge on the contacts between semiconductor and readout circuit, as analytically modelled by Sclar (1984). Sclar’s model can be formulated as the superposition of three exponentials to describe the response of the detector to a step in illumination. A superposition of three exponentials has also been proposed by Fujiwara et al. (1995) and Church et al. (1996) to represent the transient behaviour. The three timescales involved are a minimum requirement to reproduce the hook response and the long term approach to the asymptotic response. However, no theory exists to predict values for the constants in this representation, and thus they must be found empirically.

In practice it is a formidable and unsolved problem to extract the constants from in flight data for as many as three temporal components to the transient response behaviour. In this paper we describe an approximate solution involving a superposition of two exponentials, each with illumination-dependent constants. This was found to adequately represent the detector response on timescales greater than a few seconds, though it only approximately model the hook response (see also Blomme & Runacres 2002). Our model can be viewed as a generalisation of the single exponential model developed for correction of longer stares by Acosta-Pulido, Gabriel & Castañeda (2000).

In the model considered here, schematically illustrated in Fig. 4, the basic detector response $`S`$ to an instantaneous step in illumination from $`S_\mathrm{}\mathrm{p}`$ to $`S_{\mathrm{}}`$ is given by the sum of a “slow” signal response $`S_1`$ and a “fast” signal response $`S_2`$ as:

$`S`$ $`=`$ $`S_1+S_2`$ (1)
$`S_1`$ $`=`$ $`(1\beta _2)S_{\mathrm{}}(1\mathrm{exp}[t/\tau _1])+S_{01}\mathrm{exp}[t/\tau _1]`$ (2)
$`S_2`$ $`=`$ $`\beta _2S_{\mathrm{}}(1\mathrm{exp}[t/\tau _2])+S_{02}\mathrm{exp}[t/\tau _2]`$ (3)
$`S_{01}`$ and $`S_{02}`$ are signals corresponding to the slow and fast response components immediately following the illumination step. $`S_{01}`$ and $`S_{02}`$ are related to the corresponding signals immediately prior to the illumination step, $`S_{1\mathrm{p}}`$ and $`S_{2\mathrm{p}}`$, which contain the information about the illumination history. $`S_{01}`$ and $`S_{02}`$ are specified by the jump conditions:

$`S_{01}`$ $`=`$ $`\beta _1(S_{\mathrm{}}S_\mathrm{}\mathrm{p})+S_{1\mathrm{p}}`$ (4)
$`S_{02}`$ $`=`$ $`S_{2\mathrm{p}}`$ (5)
The model is thus specified by four primary parameters: the slow and fast timescales $`\tau _1`$ and $`\tau _2`$, the “jump” factor $`\beta _1`$ and a constant of proportionality for the fast response $`\beta _2`$. In order to obtain satisfactory model fits to complex illumination histories, each of these primary parameters had to be specified as monotonic functions of illumination, each with three subsiduary parameters:

$`\beta _1=\beta _{10}+\beta _{11}S_{\mathrm{}}^{\beta _{12}}`$ (6)
$`\tau _1=\tau _{10}+\tau _{11}S_{\mathrm{}}^{\tau _{12}}`$ (7)
$`\beta _2=\beta _{20}+\beta _{21}S_{\mathrm{}}^{\beta _{22}}`$ (8)
$`\tau _2=\tau _{20}+\tau _{21}S_{\mathrm{}}^{\tau _{22}}`$ (9)
Thus 12 subsiduary parameters in all are needed to characterise the transient response behaviour of the detectors. In addition, the initial starting state of the detector, given by the values of $`S_{1\mathrm{p}}`$ and $`S_{2\mathrm{p}}`$ immediately prior to the start of the mapping observation, must be specified.

### 3.2 Determination of model parameters

The parameterisation of the illumination dependence of the primary parameters $`\tau _1`$, $`\beta _1`$, $`\tau _2`$ and $`\beta _2`$ was made in engineering units (V/s), rather than in astrophysical units such as MJy/sr. Thus, it was assumed that there is no dependence of the detector response on the wavelength of the IR photons reaching the detector. Since the individual pixels of the C100 and the C200 arrays behave as independent detectors, each of the 12 parameters of the detector model had to be separately determined for each detector pixel.

The “slow” primary parameters $`\tau _1`$ and $`\beta _1`$ were found from fits to the detector response to fine calibration source (FCS) exposures. To make this determination independent of the determination of the “fast” parameters, only data points corresponding to times of 5s or longer after the switch-on of the FCS was used. An example of a fit to the signal from an FCS exposure is given in Fig. 5. The model fit to the signal response on time scales ranging from a few seconds up to the duration of the FCS exposure depends mainly on the values of $`\tau _1`$ and $`\beta _1`$.

A database of values for $`\tau _1`$ and $`\beta _1`$ was built up for a sample of FCS measurements, selected so as to correspond to the complete range of sky and source illuminations encountered by ISOPHOT during the ISO mission. These allowed the dependence of $`\tau _1`$ and $`\beta _1`$ on illumination to be determined according to Eqns. 6$`\&`$7. An example showing the illumination dependence of $`\beta _1`$ for the central pixel of the C100 array is given in Fig. 6a. The solid line is the fitted illumination dependence of $`\beta _1`$ used to find the values of the parameters $`\beta _{10}`$, $`\beta _{11}`$ and $`\beta _{12}`$ used in the model. A corresponding plot for the illumination dependence of $`\tau _1`$ is given in Fig. 6b.

At shorter timescales the response is primarily determined by the “fast” primary parameters $`\tau _2`$ and $`\beta _2`$. However, the FCS exposures could not be used to find $`\tau _2`$ and $`\beta _2`$ since no data was recorded in the first two seconds following the detector switch-on. This is also the reason why the hook response is not seen in Fig. 5. Furthermore, there is very probably a non-negligible timescale needed by the FCS to reach its final operating temperature. Therefore, the “fast” parameters (with timescales $``$0.1$`\tau `$5 s) were found using a self calibration technique (described in Sect. 4.2), operating on co-added repeated chopper sweeps crossing standard point source calibrators. During the self calibration the “slow” parameters were fixed to the values determined in the fits to the FCS exposures. Observations of standard point source calibrators (12 source-wavelength combinations for both the C100 and C200 arrays) were analysed using the self-calibration technique. Each observation of a source provided several points for the determination of the ”fast” detector constants, since individual subscans through the source could be processed independently. In all, self calibrations were performed on data spanning a range in detector illumination of about 100. This allowed the dependence of $`\tau _2`$ and $`\beta _2`$ on illumination to be determined according to the parameterisation of Eqns. 8$`\&`$9. Examples showing the illumination dependence of $`\beta _2`$ and $`\tau _2`$ for the central pixel of the C100 array are given in Figs. 6c and 6d, respectively. The values of the 12 parameters determined for each pixel of the C100 and C200 detectors are given in Tables 1 and 2, respectively.

## 4 The P32 algorithm

In a standard reduction of ISOPHOT data using the interactive analysis procedures of PIA, the analysis is made in several irreversible steps, starting with input “edited raw data” at the full time resolution, and finishing with the final calibrated map. To correct for responsivity drift effects, however, it is necessary to iterate between a sky map and the input data at full time resolution. This different concept required a completely new data reduction package. In its development phase, and for the determination of the detector model parameters, this data reduction package was run as a set of IDL scripts. Running in the PIA environment, P32Tools provides access to the main functionality of this data reduction package via a GUI interface (Lu et al. 2002).

There are three basic elements to the concept of the P32 algorithm. The first step is signal conditioning, which aims to provide a stream of signal values, each with an attached time, from which artifacts such as glitches, dark current, and non-linearity effects have been removed, but which retains the full imprint of the transient response of the detector to the illumination history at the sharpest available time resolution. We will refer to such a stream of signals as the “signal timeline”.

The second step is the transient correction itself, which is an iterative process to determine the most likely sky brightness distribution giving rise to the observed signals. This is a non-linear optimisation problem with the values of sky brightness as variables. For large maps several hundred independent variables are involved.

Lastly, there is the calculation of calibrated transient corrected maps. This entails the conversion from V/s to MJy/sr, flat fielding, and the coaddition of data from different detectors pixels. Using P32Tools, this third stage requires the exporting of the transient corrected signal timeline into PIA at the SCP level, followed by the application of standard PIA processing. Here we just describe the new procedures for signal conditioning and transient correction.

### 4.1 Signal Conditioning

The processing steps for signal conditioning are as follows:

1. Edited Raw Data is imported, after having been corrected in PIA for the integration ramp non-linearity.
2. The specification of chopper phase in the data is checked and corrected where necessary.
3. The pointing is established from the instantaneous satellite pointing information, unlike PIA which uses the commanded pointing. This also establishes the pointing directions on the slews between the spacecraft fine pointings. Internally defined “on target flags” are set for each spacecraft fine pointing direction, and the central direction and sampling intervals directions for the P32 “natural grid” (see Sect. 2.1) are established.
4. The integration ramps are differentiated and the data is corrected for dark current and dependence of photometry on reset interval, using standard calibration parameters taken from PIA.
5. Systematic variations of the signal according to position in the integration ramps are removed from the data. This process is called “signal relativisation”. In contrast to the default pipeline or PIA analysis, this allows the use of all non-destructive reads, as well as providing a more complete signal timeline for use in the correction of the transient behaviour of the detector. The “signal relativisation” is applied such that signal variations due to the transient response behaviour are preserved.
6. Random noise is determined by examining the statistics of the signals on each chopper plateau. This is necessary due to the high readout rates of the detectors in many P32 observations. Otherwise, the determination of noise from individual readout ramps, as done in a standard PIA analysis, would be insufficiently precise.
7. Spikes are detected and removed from the data in a first-stage deglitching procedure. This operates at the full time resolution of the data given by the non-destructive read interval.
8. A second stage deglitching procedure is performed, operating at the time resolution of the chopper plateau (which comprise of at least 16 non-destructive read intervals). This procedure (described by Peschke & Tuffs 2002) removes any long lived glitch “tails” following the spikes detected in the first-stage deglitch.

### 4.2 Transient Correction

A data flow diagram summarising the iterative non-linear optimisation process to solve for the sky brightnesses sampled on the P32 natural grid is summarised in Fig. 7. The calculation to correct for the transient response behaviour is made separately for each detector pixel. The kernel of the procedure is to solve for the illumination corresponding to the signal on each chopper plateau. This is done in a constrained optimisation using a binary chop scheme for set values of the 12 detector parameters. The constraints set for the solution in illumination define the search range for illumination. They are set such that the solution for illumination should be positive, and less than 10 times the value of the highest uncorrected signal. In practice the illuminations found on bright sources are never more than a factor of 5 higher than the uncorrected signal, so the determined photometry is not sensitive to the exact value of the upper limit in the search range. The values for $`S_{1\mathrm{p}}`$ and $`S_{2\mathrm{p}}`$ needed to solve for illumination on each chopper plateau are calculated from the response of the model to the preceding illumination history.

After solving for the illumination on the current chopper plateau (in the representation in Fig. 7, to give the kth value $`I_K`$ of the sequence of illuminations in the “signal timeline”), the found value $`I_K`$ is divided by the vignetting corresponding to the chopper phase to obtain a new estimate of the sky brightess for the viewed direction on the P32 natural grid. This sky brightness estimate is then combined with any previous solutions for the same P32 natural grid pixel obtained with the same detector pixel from previous chopper plateaus, and the trial sky map (sampled on the P32 natural grid) is updated for this improved knowledge. The response of the detector pixel to the illumination history is then recalculated and a solution is found for the illumination viewed on the next chopper plateau. This procedure is continued until all the signal timeline has been processed. Then a second pass through the data can be made (the “second major iteration” loop in Fig. 7), and so on, until the overall solution for the sky brightness distribution no longer changes between successive major iterations. Any sky directions for which no solution can be found are flagged in an output mask, sampled on the P32 natural grid.

The goodness of the solution is quantified through a value for $`\chi ^2`$, calculated from a comparison between the model fit and the signal timeline. The random uncertainties assigned to each signal sample in the timeline can optionally be recalculated from the self consistency of solutions for the illumination found in a given sky direction. Examples of model fits to the signal timeline are shown in Fig. 8 for the C100 detector and in Fig. 9 for the C200 detector.

The algorithm can also be used for datasets for which the readout rate was too rapid to allow all the data to be transmitted to the ground. In such cases the illumination history can still be determined at all times on fine pointings, even where there are gaps in the signal timeline. On slews, the illumination history during the gaps are either interpolated in time, or taken from the nearest previously solved direction on the P32 natural grid.

For most observations using the C200 detector, and observations of moderate intensity sources made using the C100 detector, good fits to the entire unbroken signal timeline can be obtained for each detector pixel. For observations of bright sources using the C100 detector, however, it is often necessary to break up the solutions into smaller segments, according to the distribution of the sources in the field mapped. This is best achieved by solving separately for the portions of the signal timeline passing through the bright sources. The most conservative processing strategy, appropriate for observations using the C100 detector of complex brightness distributions with high source to background contrast, is to break up the calculation into time intervals corresponding to individual spacecraft fine pointings. In this latter case, no processing of the slews is required.

Several processing options are available for the interactive processing, to optimise solutions for the particular characteristics of each dataset:

Determination of detector starting state

The default operation of the program is to assume that the detector is in equilibrium at the start of the observation, viewing the sky direction corresponding to the first detector plateau. This corresponds to the case $`S_{1\mathrm{p}}=S_{01}=S_{\mathrm{}}`$ and $`S_{1\mathrm{p}}=\mathrm{\hspace{0.17em}0}`$. This will very rarely be the case, however, since an FCS measurement which may not be perfectly matched to the sky brightness will have been made just seconds prior to the start of the spacecraft raster. Therefore, the algorithm can optionally find a solution for the detector starting state, as characterised by the values of $`S_{1\mathrm{p}}`$ and $`S_{2\mathrm{p}}`$ immediately prior to the start of the first chopper plateau. This is done by solving for $`S_{1\mathrm{p}}`$ and $`S_{2\mathrm{p}}`$ for a fixed $`S_{\mathrm{}}`$ on the first plateau, where $`S_{\mathrm{}}`$ is determined from the last chopper plateau in the first spacecraft pointing to view the same sky direction. The timeline for the first spacecraft pointing is then repeatedly processed until a stable solution for $`S_{1\mathrm{p}}`$ and $`S_{2\mathrm{p}}`$ prior to the first plateau is found. This option was found to be particularly useful for data taken with the C200 detector.

“Self calibration” of detector parameters

In order to determine the primary detector parameters $`\beta _2`$ and $`\tau _2`$ in the first instance, “a self calibration” processing technique was used by which any of the 12 detector parameters (Sect. 3.1) could be determined from observations of bright sources. The self calibration functionality has also been implemented in the publically available P32Tools. This is because, even though the standard parameters of the detector model given in Tables 1 and 2 generally give good results, some fine tuning of the parameters may still yield improved fits to the signal timeline for some astronomical targets, particularly for targets with a high contrast to the background.

Sources of any arbitrary brightness distribution and position can be used for the self calibration. However, in the initial determination of the detector parameters for the “fast” response of the detector, point source calibrators of known brightness were used, and the solutions for $`\tau _2`$ and $`\beta _2`$ were additionally constrained by the requirement that the solutions for illumination were consistent with the known flux densities of the calibrators.

The self calibration procedure works by simply repeating the optimisation process for different subsets of the 12 detector parameters, specified by the user. Each combination of detector parameters produces an individual solution for the sky brightness distribution, as well as a value for $`\chi ^2`$. A search is made in the parameter space of the detector parameters until a minimum is found in $`\chi ^2`$. This is a lengthy process, however, making it advisable to perform self calibrations on limited portions of the timeline. This is often chosen to correspond to a single spacecraft pointing where the chopper sweeps for the detector pixel being investigated pass through bright structure. The process can be further accelerated in cases that the repeated chopper sweeps through a source give the same repeated signal pattern. Then, the optimisation can be performed on an average of the coadded chopper sweeps performed on a single spacecraft pointing. This latter technique is referred to in P32Tools as a “composite self calibration”. Some examples of the useage of the self calibration procedure are given by Schulz et al. (2002a).

Processing of slews

The calculation of the illumination history for data taken during the slews is less straightforward than for data taken during the spacecraft fine pointings. This is because the instantaneous pointing directions during slews are effectively at random positions with respect to the P32 natural grid. Therefore, the use of solutions for the sky brightness distribution on the P32 natural grid found from analysis of previous fine pointings will in general lead to inaccuracies in the determination of the illumination history during the slews. There are two main approaches which can be adopted to counter this problem.

The first is to solve for the illumination not over a chopper plateau, but for individual readouts in the signal timeline during the slew. This automatically provides a solution in which the model fit exactly matches the data, as there is only one fitted data point for each determined illumination. An example is given in the fit to the signal timeline for the C200 detector shown in Fig 9. This approach is useful for observations of bright structured sources in which the variations in illumination with position on the slew on timescales of a chopper plateau exceed the signal to noise on individual data samples. It is the default processing option for the C200 detector.

In cases of low signal to noise, and for all instances where not all data are transmitted to the ground, solving for the average illumination on each chopper plateau is a preferable technique. Otherwise, the values of $`S_{1\mathrm{p}}`$ and $`S_{2\mathrm{p}}`$ calculated just prior to the first chopper plateau on the next fine pointing can be wildly inaccurate, which can corrupt the solutions for the remainder of the signal timeline of the observation. Processing on a plateau by plateau basis is the default processing option to calculate the illumination history during slews for observations with the C100 detector. Another reason for processing slews on a plateau-by-plateau basis for the C100 detector is that the procedure to solve for illumination for individual readouts is very susceptible to any residual glitches, and the C100 detector is more prone to glitching than the C200 detector.

## 5 Photometric performance

In this section we give examples of the photometric performance achieved from maps made using the P32 algorithm. In all cases the standard calibration parameters used in the signal conditioning step prior to the correction of the transient response behaviour and the subsequent

conversion from engineering units ($`V/s`$) to astronomical units ($`MJy/str`$ and $`Jy`$) were taken from version 8.1 of PIA. To illustrate the quantitative effect of the transient correction we give results for datasets processed by the P32 algorithm both with and without the transient correction. In all cases, the corrections for the transient response behaviour were made using the standard detector model parameters given in Tables 1 and 2. In Sect. 5.1 we give examples of the effects of the transient correction on the derived spatially integrated flux densities. Where possible these integrated flux densities are also compared with flux densities from IRAS or from source models for the standard calibrators. In Sect. 5.2 we give examples of the imaging performance.

### 5.1 Integrated Flux Densities

In general the corrections in integrated flux densities made by the algorithm depend on the source brightness, structure, size, the source/background ratio, and the dwell time on each chopper plateau. The largest corrections are for bright point sources on faint backgrounds observed with the C100 detector.

In Fig. 10 we plot the measured integrated and colour-corrected flux densities of standard ISOPHOT calibrators (the asteroids Ceres and Vesta, and the planet Neptune) versus the nominal flux densities of these sources (determined as described by Schulz et al. 2002b). Since these same sources were used to derived the detector model, one expects a good statistical match between the nominal and the derived flux densities. Indeed, the integrated flux densities found after correction for the transient behaviour of the detector lie reasonably close to the solid line indicating an equality between the measured and nominal flux densities, though with some scatter (of around 15 percent). It can be seen that the flux densities derived from data uncorrected for the transient effects have been raised by factors ranging from $``$2 for the fainter standards to $``$5 for the brighter ones.

As an example of the photometric performance for fainter point sources we give results achieved for the faint standard star HR 1654. This source was not used in the determination of the detector model parameters, so it constitutes a test of the photometric performance of ISOPHOT in its P32 observing mode. The derived integrated flux densities, with and without correction for the transient response behaviour of the C100 detector, were compared in Fig. 11 with predicted flux densities from a stellar model. The corrected photometry is in reasonable agreement with the theoretical predictions. As expected, there is again a trend for observations with larger detector illuminations to have larger corrections in integrated photometry.

A good linear correlation is also seen between integrated flux densities of Virgo cluster galaxies (Tuffs et al. 2002a), derived from P32 ISOPHOT observations processed using the P32 algorithm, and flux densities from the IRAS survey (Fig. 12). The interpretation of these measurements constituted the first science application (Popescu et al. 2002) of the P32 algorithm. The ISOPHOT observations of Virgo cluster galaxies were furthermore used to derive the ratios of fluxes measured by ISO to those measured by IRAS. The ISO/IRAS ratios were found to be 0.95 and 0.82 at 60 and 100$`\mu \mathrm{m}`$, respectively, after scaling the ISOPHOT measurements onto the COBE-DIRBE flux scale (Tuffs et al. 2002a).

The effect of the transient response behaviour on the integrated flux densities of very extended sources is illustrated by the photometry of a few bright extended nearby galaxies. Thus, in Table 3, Col. 5, we give the integrated flux densities of NGC 6946, M 51 and M 101 taken from Tuffs et al. (2003; in preparation). For NGC 6946 the integrated flux densities are given here for the first time; those for M 51 and M 101 supercede the values given by Hippelein et al. (1996). The standard values for the detector parameters (as given in Tables 1 and 2) were used for the correction of the transient response. Table 3 also gives a comparison with the IRAS flux densities at 60 and 100$`\mu \mathrm{m}`$. In general, the ISOPHOT photometry yields somewhat lower integrated flux densities than the corresponding IRAS photometry, similar to the ISO/IRAS ratios found for the Virgo cluster galaxies (Tuffs et al. 2002a). As for the point sources, we also reduced the maps of these nearby galaxies without correcting for the transient response of the detectors. We were thus able to derive the “correction factor” (Col. 6) indicating the ratio between the integrated flux densities with and without correcting for the transient. In comparison to the point sources of similar integrated flux densities, these corrections factors are moderate, varying between 1.07 and 1.45. This is because the timescale for the detector to traverse these extended sources becomes comparable to the time constants for the transient behaviour. The result therefore is a redistribution of flux density within the area of the galaxies along the scan direction, which goes some way to offset the loss of signal. We emphasise however that the remaining corrections are systematic, rather than random effects, and that the precise value of the correction factor for integrated flux density depends on the extent and morphology of the target. Furthermore, the correction factors for any point-like sources within extended sources, such as galactic nucleii or HII regions, will be much larger, comparable to those found for the point source calibrators.

### 5.2 Imaging Performance

Fig. 13 shows an example of a brightness profile through HR 1654 at 100$`\mu \mathrm{m}`$, for data processed with and without the responsivity correction. Some 95$`\%`$ of the flux density has been recovered with the P32 algorithm. Without the correction, some 30$`\%`$ of the integrated emission is missing and the signal only reaches 50$`\%`$ of the peak illumination. The local minimum near 60 arcsec in the Y offset is a typical hook response artifact, where the algorithm has overshot the true solution. This happens for rapid chopper sweeps passing through the beam kernel. This is a fundamental limitation of the detector model, which, as described in Sect. 3.1, does not correctly reproduce the hook response on timescales of up to a few seconds. This problem is particularly apparent for downwards illumination steps. The only effective antidote is to mask the solution immediately following a transition through a bright source peak.

Another effect of the inability to model the hook response is that the beam profile becomes somewhat distorted. This also has the consequence, that for observations of bright sources, the measured FWHM can become narrower than for the true point spread function, as predicted from the telescope optics and pixel footprint. We emphasise that, because of the beam distortion, the flux density of bright point sources cannot be found from the peak surface brightness, but instead should be determined by integrating the map. An example of an extremely bright point source reduced using the P32 algorithm is given in Fig. 14, depicting maps of Ceres in the C105 filter, made with and without the correction for the transient response of the C100 detector.

Despite the limitations due to the lack of a proper modelling of the hook response, the algorithm can effectively correct for artifacts associated with the transient response on timescales from a few seconds to a few minutes. This is illustrated in Fig. 15 by the maps of the interacting galaxy pair KPG 347 in the C200 filter (again, after and before correction for the transient response of the detector, respectively). The uncorrected map shows a spurious elongation in the direction of the spacecraft Y coordinate, which is almost completely absent in the corrected map. If uncorrected, such artifacts could lead to false conclusions about the brightness of FIR emission in the outer regions of resolved sources. Also visible are artifacts at Y$`=\pm `$310 arcsec. Similar artifacts are also seen at the same offset in other P32 observations made with the C200 detector, with position and amplitude relative to the central source independent of the filter used. These features are attributable to optical effects peculiar to the P32 mode, in which a small fraction of the light from a centrally positioned compact source is deflected away from the main beam into the vignetted part of the field of view for a certain combination of the offset in spacecraft fine pointing and chopper phase. Although these artifacts have nothing to do with the transient response behaviour of the detector, they are alleviated by the application of the P32 algorithm when solutions for detector illumination are converted into sky brightnesses, due to the vignetting correction made in this processing step.

The image of M 101 in Fig. 16, made using the C100 detector with the C100 filter, is given as a state of the art example of what can be achieved with a careful interactive processing of P32 data using the P32 algorithm. In addition to the transient response corrections and a masking of residual hook response artifacts, a time dependent flat field has been applied. The spiral structure of the galaxy, with embedded HII region complexes and a component of diffuse interarm emission can clearly be seen.

###### Acknowledgements.

This work was supported by grant 50-QI-9201 of the Deutsches Zentrum für Luft- und Raumfahrt. We would like to acknowledge all those who have helped us in many ways in the development of the algorithm described here. Richard Tuffs would like to thank his colleagues at the Max-Planck-Institut für Kernphysik, in particular Prof. Heinrich Völk, for their support and encouragement. Dr. Cristina Popescu is thanked for the careful reading of the manuscript and comments. Thanks are also due to Prof. Rolf Chini for his input to the definition of the original concept of the P32 AOT. We have also benefited from many useful discussions with Drs. R. Laureijs, S. Peschke and B. Schulz and the team in the ISO data centre at Villafranca, with Prof. D. Lemke, Drs. U. Klaas, M. Haas and M. Stickel at the ISOPHOT data centre at the Max-Planck-Institut für Astronomie, and with Drs. N. Lu and I. Khan at the Infrared Processing and Analysis Center. We thank the referee for his helpful and astute comments.


format: text

text:
# Diffractive Final States with the H1 Detector at HERA

## 1 Introduction

Understanding diffraction in hadronic interactions at high energies, where at least one of the beam hadrons remains intact, losing only a small fraction $`x_{IP}`$ of its incident longitudinal momentum, represents one of the most important challenges in Quantum Chromodynamics (QCD). The $`ep`$ collider HERA offers the unique possibility to study hard diffractive processes, such as dijet and heavy quark production, over a wide range of photon virtualities $`Q^2`$.

Recent high precision measurements of the inclusive diffractive deep-inelastic scattering (DIS) cross section by the H1 collaboration f2d97  have been used to extract diffractive parton distributions (dpdf’s) of the proton by means of a DGLAP QCD fit. If QCD factorization is valid in diffractive $`ep`$ scattering, these dpdf’s are universal and can be used to predict cross sections for exclusive hard diffractive processes such as jet and heavy quark production.

The proof of QCD factorization in diffractive $`ep`$ interactions collins  is restricted to deep-inelastic scattering at large $`Q^2`$ and not valid for the case of an almost real photon (photoproduction). Furthermore, a severe breakdown of factorization has been observed when using dpdf’s obtained at HERA to predict diffractive jet production in hadron-hadron interactions at the Tevatron cdf .

Here, QCD factorization is tested in DIS to next-to-leading order (NLO) by comparing cross section measurements for dijet and $`D^{}`$ meson production in diffractive DIS with predictions obtained from convoluting the dpdf’s from f2d97  with NLO order matrix elements for dijet and charm quark production nlopaper . In addition, diffractive dijet photoproduction data are compared with leading order (LO) Monte Carlo predictions gpjets .

### 1.1 Kinematics

Fig. 1 shows an example diagram for diffractive dijet or charm quark production at HERA.

The photon (virtuality $`Q^2`$) emitted from the beam electron (not shown) interacts with the proton, which loses only a small fraction $`x_{IP}`$ of its incident momentum and stays intact. The longitudinal momentum fraction of the parton entering the hard scattering process relative to the diffractive exchange is labelled $`z_{IP}`$. A pair of high transverse momentum ($`p_T`$) jets or heavy quarks is produced. The photon-proton centre-of-mass energy is labelled $`W`$ and the invariant masses of the diffractively produced system and of the two partons emerging from the hard sub-process are denoted $`M_X`$ and $`M_{12}`$, respectively. The inelasticity variable $`y`$ is given by $`ys=W^2+Q^2`$, where $`s`$ is the squared $`ep`$ centre-of-mass energy.

### 1.2 H1 Diffractive Parton Distributions

QCD factorization permits to express the diffractive DIS $`\gamma ^{}p`$ cross section as a convolution of diffractive parton distributions $`p_i^D`$ and partonic cross sections $`\widehat{\sigma }^{\gamma ^{}i}`$:

$$\frac{\mathrm{d}^2\sigma _{\gamma ^{}p}^D}{\mathrm{d}x_{IP}\mathrm{d}t}\underset{i}{}\widehat{\sigma }^{\gamma ^{}i}(x,Q^2)p_i^D(x,Q^2,x_{IP},t),$$
(1)
where $`t`$ is the squared 4-momentum transferred at the proton vertex. The dpdf’s are universal and obey the DGLAP evolution equations. In the QCD fit to the H1 inclusive diffractive DIS data in f2d97 , Regge factorization was assumed as well (consistent with the data within the present level of precision), i.e.

$$p_i^D(x,Q^2,x_{IP},t)=f_{IP}(x_{IP},t)p_i^{IP}(z=x/x_{IP},Q^2),$$
(2)
so that the shape of the dpdf’s is independent of $`x_{IP},t`$ and the normalization is controlled by a flux factor $`f_{IP}`$, for which a value for the pomeron intercept $`\alpha _{IP}(0)=1.173`$ is used. The parton distributions $`p_i^{IP}`$, as determined from both the LO and NLO QCD fits, are dominated by the gluon density and extend to large fractional momenta $`z`$.

## 2 Dijet Production in Diffractive DIS

To test QCD factorization for diffractive dijet production in DIS, the H1 cross sections from h1jets  were used, corresponding to the kinematic range $`4<Q^2<80\mathrm{GeV}^2`$ and $`x_{IP}<0.05`$. The cross sections were corrected to asymmetric cuts on the jet transverse momentum $`p_{T,1(2)}>5(4)\mathrm{GeV}`$, to facilitate comparisons with NLO calculations.

The NLO dpdf’s obtained in f2d97  were convoluted with NLO matrix elements for dijet production by means of the DISENT disent  program, as suggested in hautmann . The renormalization and factorization scales were set to the average $`p_T`$ of the two highest $`p_T`$ partons. $`\alpha _s`$ was set via $`\mathrm{\Lambda }_{QCD}^{n=4}=200\mathrm{MeV}`$, as in the QCD fit. To allow comparisons with the measured data, the same jet algorithm as for the data, as well as hadronization corrections, were applied.

Comparisons of the LO and NLO QCD calculations with the dijet data are shown in Figs. 2 and 3. The size of the NLO corrections is on average more than a factor 2 (increasing with decreasing $`p_T`$ and $`Q^2`$). The inner error band of the NLO calculations represents the renormalization scale uncertainty, whereas the outer band includes the uncertainty in the hadronization corrections. Within the uncertainties, the data are well described in both shape and normalization by the NLO calculations based on the dpdf’s from f2d97 , in agreement with QCD factorization.

## 3 $`𝑫^{\mathbf{}}`$ Meson Production in Diffractive DIS

In h1dstar , cross sections for $`D^{}`$ meson production in diffractive DIS were published for the kinematic range $`2<Q^2<100\mathrm{GeV}^2`$, $`x_{IP}<0.04`$ and $`p_{T,D^{}}^{}>2\mathrm{GeV}`$, where the latter variable corresponds to the transverse momentum of the $`D^{}`$ meson in the photon-proton centre-of-mass frame.

Both LO and NLO QCD calculations were performed using the diffractive version dhvq  of HVQDIS hvqdis , interfaced to the H1 dpdf’s. The renormalization and factorization scales were set to $`\mu ^2=Q^2+4m_c^2`$. Further parameter values are a charm quark mass of $`m_c=1.5\mathrm{GeV}`$, a $`cD^{}`$ hadronization fraction of $`f(cD^{})=0.233`$ and $`ϵ=0.078`$ for the used Peterson fragmentation function.

A comparison of the calculations with the $`D^{}`$ data is shown in Fig. 4. The inner error band of the NLO calculation represents the renormalization scale uncertainty, whereas the outer error band includes variations of $`m_c`$ and $`ϵ`$. Within the uncertainties, the data are well described in both shape and normalization by the NLO calculations, using the dpdf’s from f2d97 , supporting the idea of QCD factorization.

## 4 Dijets in Diffractive Photoproduction

In gpjets , H1 has presented a measurement of diffractive dijet photoproduction, i.e. for $`Q^20`$. The cross sections correspond to the kinematic range $`Q^2<0.01\mathrm{GeV}^2`$, $`0.3<y<0.65`$, $`x_{IP}<0.03`$ and $`p_{T,1(2)}>5(4)\mathrm{GeV}`$, where jets are defined using the $`k_T`$ algorithm.

The data are compared in Figs. 5 (absolute cross sections) and 6 (normalized cross sections) with LO Monte Carlo predictions using both the LO version of the dpdf’s from f2d97 , as well as an earlier version from h1f2d94 . $`x_\gamma ^{jets}`$ is an estimator for the momentum fraction of the photon entering the hard process, which is used to distinguish between direct and resolved photon interactions, where in the latter case the photon develops hadronic structure, and a parton from the photon with momentum fraction $`x_\gamma `$ enters the hard process.

For both direct and resolved interactions, the dijet data are well described by the Monte Carlo predictions if the new H1 dpdf’s from f2d97  are used. Thus, the data are not consistent with a breakdown of factorization even in resolved photoproduction, which resembles hadron-hadron scattering where the presence of the second hadron may lead to a suppression of the rate of diffractive events due to spectator interactions.

## 5 Conclusions

The new H1 diffractive parton distributions have been interfaced with NLO QCD calculations and compared with dijet and $`D^{}`$ production data in diffractive DIS. The good agreement within the uncertainties lends support for the validity of QCD factorization in diffractive DIS, now tested to next-to-leading order. In addition, also photoproduction diffractive dijet measurements are described in shape and normalization by Monte Carlo predictions using these dpdf’s, not indicating any breakdown of factorization in diffractive $`ep`$ interactions at HERA.


text:
# Local-search techniques for propositional logic extended with cardinality constraints

## 1 Introduction

We propose and study local-search satisfiability solvers for an extension of propositional logic with explicit means to represent cardinality constraints.

In recent years, propositional logic has been attracting considerable attention as a general-purpose modeling and computing tool, well suited for solving search problems. For instance, to solve a graph $`k`$-coloring problem for an undirected graph $`G`$, we construct a propositional theory $`T`$ so that its models encode $`k`$-colorings of $`G`$ and there is a polynomial-time method to reconstruct a $`k`$-colorings of $`G`$ from a model of $`T`$. Once we have such a theory $`T`$, we apply to it a satisfiability solver, find a model of $`T`$ and reconstruct from the model the corresponding $`k`$-coloring of $`G`$.

Instances of many other search problems can be represented in a similar way as propositional theories and this modeling capability of the propositional logic has been known for a long time. However, it has been only recently that we saw a dramatic improvement in the performance of programs to compute models of propositional theories . These new programs can often handle theories consisting of hundreds of thousands, sometimes millions, of clauses. They demonstrate that propositional logic is not only a tool to represent problems but also a viable computational formalism.

The approach we outlined above has its limitations. The repertoire of operators available for building formulas to represent problem constraints is restricted to boolean connectives. Moreover, since satisfiability solvers usually require CNF theories as input, for the most part the only formulas one can use to express constraints are clauses. One effect of these restrictions is often very large size of CNF theories needed to represent even quite simple constraints and, consequently, poorer effectiveness of satisfiability solvers in computing answers to search problems. Researchers recognized this limitation of propositional logic. They proposed extensions to the basic language with the equivalence operator , with cardinality atoms  and with pseudo-boolean constraints , and developed solvers capable of computing models for theories in the expanded syntax.

In this paper, we focus on an extension of propositional logic with cardinality atoms, as described in . Specifically, a cardinality atom is an expression of the form $`kXm`$, where $`k`$ and $`m`$ are non-negative integers and $`X`$ is a set of propositional atoms. Cardinality atoms offer a direct means to represent cardinality constraints on sets and help construct concise encodings of many search problems. We call this extension of the propositional logic the propositional logic with cardinality constraints and denote it by $`\mathrm{𝑃𝐿}^{cc}`$.

To make the logic $`\mathrm{𝑃𝐿}^{cc}`$ into a computational mechanism, we need programs to compute models of $`\mathrm{𝑃𝐿}^{cc}`$ theories. One possible approach is to compile cardinality atoms away, replacing them with equivalent propositional-logic representations. After converting the resulting theories to CNF, we can use any off-the-shelf satisfiability solver to compute models. Another approach is to design solvers specifically tailored to the expanded syntax of the logic $`\mathrm{𝑃𝐿}^{cc}`$. To the best of our knowledge, the first such solver was proposed in . A more recent solver, aspps111The acronym for answer-set programming with propositional schemata., was described in .

These two solvers are complete solvers. In this paper, we propose and study local-search satisfiability solvers that can handle the extended syntax of the logic $`\mathrm{𝑃𝐿}^{cc}`$. In our work we built on ideas first used in WSAT, one of the most effective local-search satisfiability solvers for propositional logic 222In the paper, we write WSAT instead of WALKSAT to shorten the notation.. In particular, as in WSAT, we proceed by executing a prespecified number of tries. Each try starts with a random truth assignment and consists of a sequence of local modification steps called flips. Each flip is determined by an atom selected from an unsatisfied clause. We base the choice of an atom on the value of its break-count (some measure of how much the corresponding flip increases the degree to which the clauses in the theory are violated). In $`\mathrm{𝑊𝑆𝐴𝑇}`$, the break-count of an atom is the number of clauses that become unsatisfied when the truth value of the atom is flipped. In the presence of cardinality atoms, this simple measure does not lead to satisfactory algorithms and modifications are necessary.

In this paper, we propose two approaches. In the first of them, we change the definition of the break-count. To this end, we exploit the fact that cardinality atoms are only high-level shorthands for some special propositional theories and, as we already indicated earlier, can be compiled away. Let $`T`$ be a $`\mathrm{𝑃𝐿}^{cc}`$ theory and let $`T^{}`$ be its propositional-logic equivalent. We define the break-count of an atom $`a`$ in $`T`$ as the number of clauses in the compiled theory $`T^{}`$ that become unsatisfied after we flip $`a`$. Important thing to note is that we do not need to compute $`T^{}`$ explicitly in order to compute the break-count of $`a`$. It can be computed directly on the basis of $`T`$ alone.

Our second approach keeps the concept of the break-count exactly as it is defined in WSAT but changes the notion of a flip. This approach applies whenever a $`\mathrm{𝑃𝐿}^{cc}`$ theory $`T`$ can be separated into two parts $`T_1`$ and $`T_2`$ so that: (1) $`T_2`$ consists of propositional clauses, (2) it is easy to construct random assignments that satisfy $`T_1`$, and (3) for every truth assignment satisfying $`T_1`$, (modified) flips executed on this assignment result in assignments that also satisfy $`T_1`$. In such cases, we can start a try by generating an initial truth assignment to satisfy all clauses in $`T_1`$, and then executing a sequence of (modified) flips, choosing atoms for flipping based on the number of clauses in $`T_2`$ (which are all standard propositional CNF clauses) that become unsatisfiable after the flip.

In the paper, we develop and implement both ideas. We study experimentally the performance of our algorithms on several search problems: the graph coloring problem, the vertex-cover problem and the open latin-square problem. We compare the performance of our algorithms to that of selected SAT solvers executed on CNF theories obtained from $`\mathrm{𝑃𝐿}^{cc}`$ theories by compiling away cardinality atoms.

## 2 Logic $`\mathrm{𝑃𝐿}^{cc}`$

The language of the logic $`\mathrm{𝑃𝐿}^{cc}`$ is determined by the set $`\mathrm{𝐴𝑡}`$ of propositional atoms and two special symbols $``$ and $``$ that we always interpret as false and true, respectively. A cardinality atom (c-atom, for short) is an expression of the form $`kXm`$, where $`X`$ is a set of propositional atoms, and $`k`$ and $`m`$ are non-negative integers. If $`X=\{a_1,\mathrm{},a_n\}`$, we will also write $`k\{a_1,\mathrm{},a_n\}m`$ to denote a c-atom $`kXm`$. One (but not both) of $`k`$ and $`m`$ may be missing. Intuitively, a c-atom $`kXm`$ means: at least $`k`$ and no more than $`m`$ of atoms in $`X`$ are true. If $`k`$ (or $`m`$) is missing, the c-atom constrains the number of its propositional atoms that must be true only from above (only from below, respectively).

A clause is an expression of the form $`\neg \alpha _1\mathrm{}\neg \alpha _r\beta _1\mathrm{}\beta _s,`$ where each $`\alpha _i`$, $`1ir`$, and each $`\beta _j`$, $`1js`$, is a propositional atom or a c-atom. A theory of the logic $`\mathrm{𝑃𝐿}^{cc}`$ is any set of clauses333It is easy to extend the language of $`\mathrm{𝑃𝐿}^{cc}`$ and introduce arbitrary formulas built of atoms and c-atoms by means of logical connectives. Since clausal theories, as in propositional logic, are most fundamental, we focus on clausal theories only..

An interpretation is an assignment of truth values $`𝐭`$ and $`𝐟`$ to atoms in $`\mathrm{𝐴𝑡}`$. An interpretation $`I`$ satisfies an atom $`a`$ if $`I(a)=𝐭`$. An interpretation $`I`$ satisfies a c-atom $`k\{a_1,\mathrm{},a_n\}m`$ if $`k|\{i:I(a_i)=𝐭\}|m.`$

This notion of satisfiability extends in a standard way to clauses and theories. We will write interchangingly “is a model of” and “satisfies”. We will also write $`IE`$, when $`I`$ is a model of an atom, c-atom, clause or theory $`E`$.

We will now illustrate the use of the logic $`\mathrm{𝑃𝐿}^{cc}`$ as a modeling tool by presenting $`\mathrm{𝑃𝐿}^{cc}`$ theories that encode (1) the graph-coloring problem, (2) the graph vertex-cover problem, and (3) the open latin-square problem. We later use these theories as benchmarks in performance tests.

In the first of these problems we are given a graph $`G`$ with the set $`V=\{1,\mathrm{},n\}`$ of vertices and a set $`E`$ of edges (unordered pairs of vertices). We are also given a set $`C=\{1,\mathrm{},k\}`$ of colors. The objective is to find an assignment of colors to vertices so that for every edge, its vertices get different colors. A $`\mathrm{𝑃𝐿}^{cc}`$ theory representing this problem is built of propositional atoms $`c_{i,j}`$, where $`1in`$, and $`1jk`$. An intended meaning of an atom $`c_{i,j}`$ is that vertex $`i`$ gets color $`j`$. We define the theory $`\mathrm{𝑐𝑜𝑙}(G,k)`$ to consist of the following clauses:

1. $`1\{c_{i,1},\mathrm{},c_{i,k}\}1`$, for every $`i`$, $`1in`$. These clauses ensure that every vertex obtains exactly one color
2. $`\neg c_{p,j}\neg c_{r,j}`$, for every edge $`\{p,r\}E`$ and for every color $`j`$. These clauses enforce the main colorability constraint.

It is easy to see that models of the theory $`\mathrm{𝑐𝑜𝑙}(G,k)`$ are indeed in one-to-one correspondence with $`k`$-colorings of $`G`$.

In a similar way, we construct a theory $`\mathrm{𝑣𝑐}(G,k)`$ that represents the vertex-cover problem. Let $`G`$ be an undirected graph with the set $`V=\{1,\mathrm{},n\}`$ of vertices and a set $`E`$ of edges. Given $`G`$ and a positive integer $`k`$, the objective is to find a set $`U`$ of no more than $`k`$ vertices, such that every edge has at least one of its vertices in $`U`$ (such sets $`U`$ are vertex covers). We build the theory $`\mathrm{𝑣𝑐}(G,k)`$ of atoms $`in_i`$, $`1in`$, (intended meaning of $`in_i`$: vertex $`i`$ is in a vertex cover) and define it to consist of the following clauses:

1. $`\{in_1,\mathrm{},in_n\}k`$. This clause guarantees that at most $`k`$ vertices are chosen to a vertex cover
2. $`in_pin_r`$, for every edge $`\{p,r\}E`$. These clauses enforce the main vertex cover constraint.

Again, it is evident that models of theory $`\mathrm{𝑣𝑐}(G,k)`$ are in one-to-one correspondence with those vertex covers of $`G`$ that have no more than $`k`$ elements.

In the open latin-square problem, we are given an integer $`n`$ and a collection $`D`$ of triples $`(i,j,k)`$, where $`i`$, $`j`$ and $`k`$ are integers from $`\{1,\mathrm{},n\}`$. The goal is to find an $`n\times n`$ array $`A`$ such that all entries in $`A`$ are integers from $`\{1,\mathrm{},n\}`$, no row and column of $`A`$ contains two identical integers, and for every $`(i,j,k)D`$, $`A(i,j)=k`$. In other words, we are looking for a latin square of order $`n`$ that extends the partial assignment specified by $`D`$. To represent this problem we construct a $`\mathrm{𝑃𝐿}^{cc}`$ theory $`\mathrm{𝑙𝑠}(n,D)`$ consisting of the following clauses:

1. $`a_{i,j,k}`$, for every $`(i,j,k)D`$ (to represent the partial assignment $`D`$ given as input)
2. $`1\{a_{i,j,1},\mathrm{},a_{i,j,n}\}1`$, for every $`i,j=1,\mathrm{},n`$ (to enforce that every entry receives exactly one value)
3. $`\{a_{i,1,k},\mathrm{},a_{i,n,k}\}1`$, for every $`i,k=1,\mathrm{},n`$ (in combination with (2) these clauses enforce that an integer $`k`$ appears exactly once in a row $`i`$)
4. $`\{a_{1,j,k},\mathrm{},a_{n,j,k}\}1`$, for every $`j,k=1,\mathrm{},n`$ (in combination with (2) these clauses enforce that an integer $`k`$ appears exactly once in a column $`j`$).

One can verify that models of the theory $`\mathrm{𝑙𝑠}(n,D)`$ correspond to solutions to the open latin-square problem with input $`D`$.

The use of c-atoms in all these three examples results in concise representations of the corresponding problems. Clearly, we could eliminate c-atoms and replace the constraints they represent by equivalent CNF theories. However, the encodings become less direct, less concise and more complex.

## 3 Using SAT Solvers to Compute Models of $`\mathrm{𝑃𝐿}^{cc}`$ Theories

We will now discuss methods to find models of $`\mathrm{𝑃𝐿}^{cc}`$ theories by means of standard SAT solvers. A key idea is to compile away c-atoms by replacing them with their propositional-logic descriptions. We will propose several ways to do so.

Let us consider a c-atom $`C=k\{a_1,\mathrm{},a_n\}m`$ and let us define a CNF theory $`C^{}`$ to consist of the following clauses:

1. $`\neg a_{i_1}\mathrm{}\neg a_{i_{m+1}}`$, for any $`m+1`$ atoms $`a_{i_1},\mathrm{},a_{i_{m+1}}`$ from $`\{a_1,\mathrm{},a_n\}`$ (there are $`\left(\genfrac{}{}{0pt}{}{n}{m+1}\right)`$ such clauses); and
2. $`a_{i_1}\mathrm{}a_{i_{nk+1}}`$, for any $`nk+1`$ atoms $`a_{i_1},\mathrm{},a_{i_{nk+1}}`$ in $`\{a_1,\mathrm{},a_n\}`$ (there are $`\left(\genfrac{}{}{0pt}{}{n}{k1}\right)`$ such clauses).

It is easy to see that the theory $`C^{}`$ has the same models as the c-atom $`C`$.

Let $`T`$ be a $`\mathrm{𝑃𝐿}^{cc}`$ theory. We denote by $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}(T)`$ the CNF theory obtained from $`T`$ by replacing every c-atom $`C`$ with the conjunction of clauses in $`C^{}`$ and by applying distributivity to transform the resulting theory into the CNF. This approach translates $`T`$ into a theory in the same language but it is practical only if $`k`$ and $`m`$ are small (do not exceed, say 2). Otherwise, the size of the theory $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}(T)`$ quickly gets too large for SAT solvers to be effective.

Our next method to compile away c-atoms depends on counting. To simplify the presentation, we will describe it in the case of a c-atom of the form $`kX`$ but it extends easily to the general case. We will assume that $`k1`$ (otherwise, $`kX`$ is true) and $`k|X|`$ (otherwise $`kX`$ is false).

Let us consider a $`\mathrm{𝑃𝐿}^{cc}`$ theory $`T`$ and let us assume that $`T`$ contains a c-atom of the form $`C=k\{a_1,\mathrm{},a_n\}`$. We introduce new propositional atoms: $`b_{i,j}`$, $`i=0,\mathrm{},n`$; $`j=0,\mathrm{},k`$. The intended role for $`b_{i,j}`$ is to represent the fact that at least $`j`$ atoms in $`\{a_1,\mathrm{},a_i\}`$ are true. Therefore, we define a theory $`C^{}`$ to consist of the following clauses:

1. $`b_{0,j}`$, $`j=1,\mathrm{},k`$,
2. $`b_{i,0}`$, $`i=0,\mathrm{},n`$,
3. $`b_{i,j}b_{i1,j}(b_{i1,j1}a_i)`$, $`i=1,\mathrm{},n`$, $`j=1,\mathrm{},k`$.

Let $`I`$ be an interpretation such that $`IC^{}`$. One can verify that $`Ib_{i,j}`$ if and only if $`Ij\{a_1,\mathrm{},a_i\}`$. In particular, $`Ib_{n,k}`$ if and only if $`IC`$. Thus, if we replace $`C`$ in $`T`$ with $`b_{n,k}`$ and add to $`T`$ the theory $`C^{}`$ the resulting theory has the same models (modulo new atoms) as $`T`$. By repeated application of this procedure, we can eliminate all c-atoms from $`T`$. Moreover, if we represent theories $`C^{}`$ in CNF, the resulting theory will itself be in CNF. We will denote this CNF theory as $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑢𝑐}(T)`$, where $`uc`$ stands for unary counting. One can show that the size of $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑢𝑐}(T)`$ is $`O(R\times \mathrm{𝑠𝑖𝑧𝑒}(T))`$, where $`R`$ is the maximum of all integers appearing in $`T`$ as lower or upper bounds in c-atoms. It follows that, in general, this translation leads to more concise theories than $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}`$. However, it does introduce new atoms.

The idea of counting can be pushed further. Namely, we can design a more concise translation than $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑢𝑐}`$ by following the idea of counting and by representing numbers in the binary system and by building theories to model binary counting and comparison. For a $`\mathrm{𝑃𝐿}^{cc}`$ theory $`T`$, we denote the result of applying this translation method to $`T`$ by $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑐}`$ ($`bc`$ stands for binary counting). Due to space limitation we omit the details of this translation. We only note that the size of $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑐}(T)`$ is $`O(\mathrm{𝑠𝑖𝑧𝑒}(T)\mathrm{log}_2(R+1))`$, where $`R`$ is the maximum of all integer bounds of c-atoms appearing in $`T`$.

## 4 Local-search Algorithms for the Logic $`\mathrm{𝑃𝐿}^{cc}`$

In this section we describe a local-search algorithm $`\mathrm{𝐺𝑒𝑛𝑒𝑟𝑖𝑐}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ designed to test satisfiability of theories in the logic $`\mathrm{𝑃𝐿}^{cc}`$. It follows a general pattern of WSAT . The algorithm executes $`\mathrm{𝑀𝑎𝑥}\text{-}\mathrm{𝑇𝑟𝑖𝑒𝑠}`$ independent tries. Each try starts in a randomly generated truth assignment and consists of a sequence of up to $`\mathrm{𝑀𝑎𝑥}\text{-}\mathrm{𝐹𝑙𝑖𝑝𝑠}`$ flips, that is, local changes to the current truth assignment. The algorithm terminates with a truth assignment that is a model of the input theory, or with no output at all (even though the input theory may in fact be satisfiable). We provide a detailed description of the algorithm $`\mathrm{𝐺𝑒𝑛𝑒𝑟𝑖𝑐}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ in Figure 1.

We note that the procedure $`\mathrm{𝐹𝑙𝑖𝑝}`$ may, in general, depend on the input theory $`T`$. It is not the case in $`\mathrm{𝑊𝑆𝐴𝑇}`$ and other similar algorithms but it is so in one of the algorithms we propose in the paper. Thus, we include $`T`$ as one of the arguments of the procedure $`\mathrm{𝐹𝑙𝑖𝑝}`$.

We also note that in the algorithm, we use several parameters that, in our implementations, we enter from the command line. They are $`\mathrm{𝑀𝑎𝑥}\text{-}\mathrm{𝑇𝑟𝑖𝑒𝑠}`$, $`\mathrm{𝑀𝑎𝑥}\text{-}\mathrm{𝐹𝑙𝑖𝑝𝑠}`$ and $`p`$. All these parameters affect the performance of the program. We come back to this matter later in Section 5.

To obtain a concrete implementation of the algorithm $`\mathrm{𝐺𝑒𝑛𝑒𝑟𝑖𝑐}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$, we need to define $`\mathrm{𝑏𝑟𝑒𝑎𝑘}\text{-}\mathrm{𝑐𝑜𝑢𝑛𝑡}(x)`$ and to specify the notion of a flip. In this paper we follow two basic directions. In the first of them, we use a simple notion of a flip, that is, we always flip just one atom. We introduce, however, a more complex concept of the break-count, which we call the virtual break-count. In the second approach, we use a simple notion of the break-count — the number of clauses that become unsatisfied — but introduce a more complex concept of a flip, which we call the double-flip.

To specify our first instantiation of the algorithm $`\mathrm{𝐺𝑒𝑛𝑒𝑟𝑖𝑐}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}(T)`$, we define the break-count of an atom $`x`$ in $`T`$ as the number of clauses in the CNF theory $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}(T)`$ that become unsatisfied after flipping $`x`$. The key idea is to observe that this number can be computed strictly on the basis of $`T`$, that is, without actually constructing the theory $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}(T)`$. It is critical since the size of the theory $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}(T)`$ is in general much larger than the size of $`T`$ (sometimes even exponentially larger). We refer to this notion of the break-count as the virtual break-count as it is defined not with respect to an input $`\mathrm{𝑃𝐿}^{cc}`$ theory $`T`$ but with respect to a “virtual” theory $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}(T)`$, which we do not explicitly construct.

Further, we define the procedure $`\mathrm{𝐹𝑙𝑖𝑝}(\sigma ,a)`$ (it does not depend on $`T`$ hence, we dropped $`T`$ from the notation) so that, given a truth assignment $`\sigma `$ and an atom $`a`$, it returns the truth assignment $`\sigma ^{}`$ obtained from $`\sigma `$ by setting $`\sigma ^{}(a)`$ to the dual value of $`\sigma (a)`$ and by keeping all other truth values in $`\sigma `$ unchanged (this is the basic notion of the flip that is used in many local-search algorithms, in particular in WSAT). We call the resulting version of the the algorithm $`\mathrm{𝐺𝑒𝑛𝑒𝑟𝑖𝑐}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}(T)`$, the virtual break-count $`\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ and denote it by $`\mathrm{𝑣𝑏}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$.

The second instantiation of the algorithm $`\mathrm{𝐺𝑒𝑛𝑒𝑟𝑖𝑐}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ that we will discuss applies only to $`\mathrm{𝑃𝐿}^{cc}`$ theories of some special syntactic form. A $`\mathrm{𝑃𝐿}^{cc}`$ theory $`T`$ is simple, if $`T=T^{cc}T^{\mathrm{𝑐𝑛𝑓}}`$, where $`T^{cc}T^{\mathrm{𝑐𝑛𝑓}}=\mathrm{}`$ and

1. $`T^{cc}`$ consists of unit clauses $`C_i=k_iX_im_i`$, $`1ip`$, such that sets $`X_i`$ are pairwise disjoint
2. $`T^{\mathrm{𝑐𝑛𝑓}}`$ consists of propositional clauses
3. for every $`i`$, $`1ip`$, $`k_i<|X_i|`$ and $`m_i>0`$.

Condition (3) is not particularly restrictive. In particular, it excludes c-atoms $`kXm`$ such that $`k>|X|`$, which are trivially false and can be simplified away from the theory, as well as those for which $`k=|X|`$, which forces all atoms in $`X`$ to be true and again implies straightforward simplifications. The effect of the restriction $`m>0`$ is similar; it eliminates c-atoms with $`m=0`$, for which it must be that all atoms in $`X`$ be false. We note that $`\mathrm{𝑃𝐿}^{cc}`$ theories we proposed as encodings of the graph-coloring and vertex-cover problems are simple; the theory encoding the latin-square problem is not.

In this section, we consider only simple $`\mathrm{𝑃𝐿}^{cc}`$ theories. Let us assume that we designed the procedure $`\mathrm{𝐹𝑙𝑖𝑝}(T,\sigma ,a)`$ so that it has the following property:

(DF) if a truth assignment $`\sigma `$ is a model of $`T^{cc}`$ then $`\sigma ^{}=\mathrm{𝐹𝑙𝑖𝑝}(T,\sigma ,a)`$ is also a model of $`T^{cc}`$.

Let us consider a try starting with a truth assignment $`\sigma `$ that satisfies all clauses in $`T^{cc}`$. If our procedure $`\mathrm{𝐹𝑙𝑖𝑝}`$ satisfies the property (DF), then all truth assignments that we will generate in this try satisfy all clauses in $`T^{cc}`$. It follows that the only clauses that can become unsatisfied during the try are the propositional clauses in $`T^{\mathrm{𝑐𝑛𝑓}}`$. Consequently, in order to compute the break-count of an atom, we only need to consider the CNF theory $`T^{\mathrm{𝑐𝑛𝑓}}`$ and count how many clauses in $`T^{\mathrm{𝑐𝑛𝑓}}`$ become unsatisfiable when we perform a flip.

Since all c-atoms in $`T^{cc}`$ are pairwise disjoint, it is easy to generate random truth assignments that satisfy all these constraints. Thus, it is easy to generate a random starting truth assignment for a try. Moreover, it is also quite straightforward to design a procedure $`\mathrm{𝐹𝑙𝑖𝑝}`$ so that it satisfies property (DF). We will outline one such procedure now and provide for it a detailed pseudo-code description.

Let us assume that $`\sigma `$ is a truth assignment that satisfies all clauses in $`T^{cc}`$ and that we selected an atom $`a`$ as the third argument for the procedure $`\mathrm{𝐹𝑙𝑖𝑝}`$. If flipping the value of $`a`$ does not violate any unit clause in $`T^{cc}`$, the procedure $`\mathrm{𝐹𝑙𝑖𝑝}(T,\sigma ,a)`$ returns the truth assignment obtained from $`\sigma `$ by flipping the value of $`a`$. Otherwise, since the c-atoms forming the clauses in $`T^{cc}`$ are pairwise disjoint, there is exactly one clause in $`T^{cc}`$, say $`kXm`$, that becomes unsatisfied when the value of $`a`$ is flipped. In this case, clearly, $`aX`$.

We proceed now as follows. We find in $`X`$ another atom, say $`b`$, whose truth value is opposite to that of $`a`$, and flip both $`a`$ and $`b`$. That is, $`\mathrm{𝐹𝑙𝑖𝑝}(T,\sigma ,a)`$ returns the truth assignment obtained from $`\sigma `$ by flipping the values assigned to $`a`$ and $`b`$ to their duals. Clearly, by performing this double flip we maintain the property that all clauses in $`T^{cc}`$ are still satisfied. Indeed all clauses in $`T^{cc}`$ other than $`kXm`$ are not affected by the flips (these clauses contain neither $`a`$ nor $`b`$) and $`kXm`$ is satisfied because flipping $`a`$ and $`b`$ simply switches their truth values and, therefore, does not change the number of atoms in $`X`$ that are true.

The only question is whether such an atom $`b`$ can be found. The answer is indeed positive. If $`\sigma (a)=𝐭`$ and flipping $`a`$ breaks clause $`kXm`$, we must have that the number of atoms that are true in $`X`$ is equal to $`k`$. Since $`|X|>k`$, there is an atom in $`X`$ that is false. The reasoning in the case when $`\sigma (a)=𝐟`$ is similar.

A pseudo-code for the procedure is given in Figure 2.

## 5 Experiments, Results and Discussion

We performed experimental studies of the effectiveness of our local-search algorithms in solving several difficult search problems. For the experiments we selected the graph-coloring problem, the vertex-cover problem and the latin-square problem. We discussed these problems in Section 2 and described $`\mathrm{𝑃𝐿}^{cc}`$ theories that encode them. To build $`\mathrm{𝑃𝐿}^{cc}`$ theories for testing, we randomly generate or otherwise select input instances to these search problems and instantiate the corresponding $`\mathrm{𝑃𝐿}^{cc}`$ encodings. For the graph-coloring and vertex-cover problems we obtain simple $`\mathrm{𝑃𝐿}^{cc}`$ theories and so all methods we discussed apply. The theories we obtain from the latin-square problem are not simple. Consequently, the algorithm $`\mathrm{𝑑𝑓}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ does not apply but all other methods do.

Our primary goal is to demonstrate that our algorithms $`\mathrm{𝑣𝑏}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ and $`\mathrm{𝑑𝑓}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ can compute models of large $`\mathrm{𝑃𝐿}^{cc}`$ theories and, consequently, are effective tools for solving search problems. To this end, we study the performance of these algorithms and compare it to the performance of methods that employ SAT solvers, specifically $`\mathrm{𝑊𝑆𝐴𝑇}`$ and $`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$ . We chose $`\mathrm{𝑊𝑆𝐴𝑇}`$ since it is a local-search algorithm, as are $`\mathrm{𝑣𝑏}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ and $`\mathrm{𝑑𝑓}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$. We chose $`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$ since it is one of the most advanced complete methods. In order to use SAT solvers to compute models of $`\mathrm{𝑃𝐿}^{cc}`$ theories, we executed them on the CNF theories produced by procedures $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑐}`$ and $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}`$ (Section 3). We selected the method $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑐}`$ as it results in most concise translations444Our experiments with the translation $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑢𝑐}`$ show that it performs worse. We believe it is due to larger size of theories it creates. We do not report these results here due to space limitations.. We selected the method $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}`$ as it is arguably the most straightforward translation and it does not require auxiliary atoms.

For all local-search algorithms, including $`\mathrm{𝑊𝑆𝐴𝑇}`$, we used the same values of $`\mathrm{𝑀𝑎𝑥}\text{-}\mathrm{𝑇𝑟𝑖𝑒𝑠}`$ and $`\mathrm{𝑀𝑎𝑥}\text{-}\mathrm{𝐹𝑙𝑖𝑝𝑠}`$: 100 and 100000, respectively. The performance of local-search algorithms depends to a large degree on the on the value of the parameter $`p`$ (noise). For each method and for each theory, we ran experiments to determine the value of $`p`$, for which the performance was the best. All results we report here come from the best runs for each local-search method.

To assess the performance of solvers on families of test theories, we use the following measures.

1. The average running time over all instances in a family
2. The success rate of a method: the ratio of the number of theories in a family, for which the method finds a solution, to the total number of instances in the family for which we were able to find a solution using any of the methods we tested (for all methods we set a limit of 2 hours of CPU time/instance).

The success rate is an important measure of the effectiveness of local-search techniques. It is not only important that they run fast but also that they are likely to find models when models exist.

We will now present and discuss the results of our experiments. We start with the coloring problem. We generated for testing five families $`C_1,\mathrm{},C_5`$, each consisting of 50 random graphs with 1000 vertices and 3850, 3860, 3870 3880 and 3890 edges, respectively. The problem was to find for these graphs a coloring with 4 colors (each of these graphs has a 4-coloring). We show the results in Table 1. Columns $`\mathrm{𝑣𝑏}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$, $`\mathrm{𝑑𝑓}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ show the performance results for our local-search algorithms run on $`\mathrm{𝑃𝐿}^{cc}`$ theories encoding the 4-colorability problem on the graphs in the families $`C_i`$, $`1i5`$. Columns $`\mathrm{𝑊𝑆𝐴𝑇}`$-bc and $`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$-bc show the performance of the algorithms $`\mathrm{𝑊𝑆𝐴𝑇}`$ and $`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$ on CNF theories obtained from the $`\mathrm{𝑃𝐿}^{cc}`$-theories by the procedure $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑐}`$. Columns $`\mathrm{𝑊𝑆𝐴𝑇}`$-basic and $`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$-basic show the performance of the algorithms $`\mathrm{𝑊𝑆𝐴𝑇}`$ and $`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$ on CNF theories produced by the procedure $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}`$ (since the bounds in c-atoms in the case of 4-coloring are equal to 1, there is no dramatic increase in the size when using the procedure $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}`$). The first number in each entry is the average running time in seconds, the second number — the percentage success rate. The results for local-search algorithms were obtained with the value of noise $`p=0.4`$ (we found this value to work well for all the methods).

In terms of the success rate, our algorithms achieve or come very close to perfect 100%, and are comparable or slightly better than the combination of $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}`$ and $`\mathrm{𝑊𝑆𝐴𝑇}`$. When comparing the running time, our algorithms are slower but only by a constant factor. The algorithm $`\mathrm{𝑣𝑏}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ is only about 0.3 times slower and the algorithm $`\mathrm{𝑑𝑓}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ is about 3.5 times slower.

Next, we note that the combination $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑐}`$ and $`\mathrm{𝑊𝑆𝐴𝑇}`$ does not perform well at all. It fails to find a 4-coloring even for a single graph. We also observe that $`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$ performs well no matter which technique is used to eliminate c-atoms. It finds a 4-coloring for every graph that we tested. In terms of the running time there is no significant difference between its performance on theories obtained by $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑐}`$ as opposed to $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}`$. However, $`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$ is, in general, slower than $`\mathrm{𝑊𝑆𝐴𝑇}`$ and our local-search algorithms $`\mathrm{𝑣𝑏}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ and $`\mathrm{𝑑𝑓}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$.

Finally, we note that our results suggest that our algorithms are less sensitive to the choice of a value for the noise parameter $`p`$. In Table 2 we show the performance results for our two algorithms and for the combination $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}`$ and $`\mathrm{𝑊𝑆𝐴𝑇}`$ on theories obtained from the graphs in the family $`C_1`$ and for $`p`$ assuming values 0.1, 0.2, 0.3 and 0.4.

We also tested our algorithms on graph-coloring instances that were used in the graph-coloring competition at the CP-2002 conference. We refer to http://mat.gsia.cmu.edu/COLORING02/ for details. We experimented with 63 instances available there. For each of these graphs, we identified the smallest number of colors that is known to suffice to color it. We then tested whether the algorithms $`\mathrm{𝑣𝑏}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$, $`\mathrm{𝑊𝑆𝐴𝑇}`$ and $`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$ (the latter two in combination with the procedure $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}`$ to produce a CNF encoding) can find a coloring using that many colors. We found that the algorithms $`\mathrm{𝑑𝑓}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$, $`\mathrm{𝑊𝑆𝐴𝑇}`$ and $`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$ (the latter two in combination with $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}`$) were very effective. Their success rate (the percentage of instances for which these methods could match the best known result) was 62%, 56% and 54%, respectively. In comparison, the best among the algorithms that participated in the competition, the algorithm MZ, has success rate of 40% only and the success rate of other algorithms does not exceed 30%.

For the vertex cover problem we randomly generated 50 graphs with 200 vertices and 400 edges. For $`i=103,\mathrm{},107`$, we constructed a family $`VC_i`$ of $`\mathrm{𝑃𝐿}^{cc}`$ theories encoding, for graphs we generated, the problem of finding a vertex cover of cardinality at most $`i`$. For this problem, the translation $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}`$ is not practical as translating just a single c-atom $`\{in_1,\mathrm{},in_{200}\}i`$ requires $`\left(\genfrac{}{}{0pt}{}{200}{i+1}\right)`$ clauses and these numbers are astronomically large for $`i=103,\mathrm{},107`$. The translation $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑐}`$ also does not perform well. Neither $`\mathrm{𝑊𝑆𝐴𝑇}`$ nor $`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$ succeed in finding a solution to even a single instance (as always, within 2 hours of CPU time/instance). Thus, for the vertex-cover problem, we developed yet another CNF encoding, which we refer to as ad-hoc. This encoding worked well with $`\mathrm{𝑊𝑆𝐴𝑇}`$ but not with $`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$. We show the results in Table 3. For this problem, the value of noise $`p=0.1`$ worked best for all local-search methods.

Our algorithms perform very well. They have the best running time (with $`\mathrm{𝑣𝑏}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ being somewhat faster than $`\mathrm{𝑑𝑓}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$) and find solutions for all instances for which we were able to find solutions using these and other techniques. In terms of the success rate $`\mathrm{𝑊𝑆𝐴𝑇}`$, when run on ad-hoc translations, performed as well as our algorithms but was several (7 to 67, depending on the method and family) times slower.

As in the case of graph coloring, our algorithms again were less sensitive to the choice of the noise value $`p`$, as shown in Table 4 (the tests were run on the family $`VC_{103}`$.

We also experimented with the vertex-cover problem for graphs of an order of magnitude larger. We randomly generated 50 graphs, each with 2000 vertices and 4000 edges. For these graphs we constructed a family $`VC_{1035}`$ consisting of 50 $`\mathrm{𝑃𝐿}^{cc}`$ theories, each encoding the problem of finding a vertex-cover of cardinality at most 1035 in the corresponding graph. With graphs of this size, all compilation methods produce large and complex CNF theories on which both $`\mathrm{𝑊𝑆𝐴𝑇}`$ and $`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$ fail to find even a single solution. Due to the use of c-atoms, the $`\mathrm{𝑃𝐿}^{cc}`$ theories are relatively small. Each consists of 2000 atoms and 4001 clauses and has a total of about 10,000 atom occurrences. Our algorithms $`\mathrm{𝑣𝑏}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ and $`\mathrm{𝑑𝑓}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ run on each of the theories in under an hour and the algorithm $`\mathrm{𝑑𝑓}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ finds a vertex cover of cardinality at most 1035 for 9 of them. The algorithm $`\mathrm{𝑣𝑏}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ is about two times faster but has worse success rate: finds solutions only in 7 instances.

The last test concerned the latin-square problem. We assumed $`n=30`$ and randomly generated 50 instances of the problem, each specifying values for some 10 entries in the array. Out of these instances we constructed a family $`LS`$ of the corresponding $`\mathrm{𝑃𝐿}^{cc}`$ theories. Since these $`\mathrm{𝑃𝐿}^{cc}`$ theories are not simple, we did not test the algorithm $`\mathrm{𝑑𝑓}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ here. The results are shown in Table 5. For the local-search methods, we used the value of noise $`p=0.1`$.

These results show that our algorithms are faster than the combination of $`\mathrm{𝑊𝑆𝐴𝑇}`$ and $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}`$ ($`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑐}`$ again does not work well with $`\mathrm{𝑊𝑆𝐴𝑇}`$) and have a better success rate. The fastest in this case is, however, the combination of $`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$ and $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑐}`$. The combination of $`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$ and $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}`$ works worse and it is also slower than our algorithms.

## 6 Conclusions

Overall, our local-search algorithms $`\mathrm{𝑣𝑏}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ and $`\mathrm{𝑑𝑓}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$, designed explicitly for $`\mathrm{𝑃𝐿}^{cc}`$ theories, perform very well.

It is especially true in the presence of cardinality constraints with large bounds where the ability to handle such constraints directly, without the need to encode them as CNF theories, is essential. It makes it possible for our algorithms to handle large instances of search problems that contain such constraints. We considered one problem in this category, the vertex-cover problem, and demonstrated superior performance of our search algorithms over other techniques. For large instances (we considered graphs with 2000 vertices and 4000 edges and searched for vertex covers of cardinality 1035) SAT solvers are rendered ineffective by the size of CNF encodings and their complexity. Even for instances of much smaller size (search for vertex covers of 103-107 elements in graphs with 200 vertices and 400 edges), our algorithms are many times faster and have a better success rate than $`\mathrm{𝑊𝑆𝐴𝑇}`$ ($`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$ is still ineffective).

Also for $`\mathrm{𝑃𝐿}^{cc}`$ theories that contain only c-atoms of the form $`1X1`$, $`X1`$ and $`1X`$, the ability to handle such constraints directly seems to be an advantage and leads to good performance, especially in terms of the success rate. In the graph-coloring and latin-square problems our algorithms consistently had comparable or higher success rate than methods employing SAT solvers. In terms of time our methods are certainly competitive. For the coloring problem, they were slower than the method based on $`\mathrm{𝑊𝑆𝐴𝑇}`$ and $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}`$ but faster than all other methods. For the latin-square problem, they were slower than the combination of $`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$ and $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑐}`$ but again faster than other methods.

Finally, we note that our methods seem to be easier to use and more robust. SAT-based method have a disadvantage that their performance strongly depends on the selection of the method to compile away c-atoms and no method we studied is consistently better than others. The problem of selecting the right way to compile c-atoms away does not appear in the context of our algorithms. Further, the performance of local-search methods, especially the success rate, highly depends on the value of the noise parameter $`p`$. Our results show that our algorithms are less sensitive to changes in $`p`$ than those that employ $`\mathrm{𝑊𝑆𝐴𝑇}`$, which makes the task of selecting the value for $`p`$ for our algorithms easier.

These results provide further support to a growing trend in satisfiability research to extend the syntax of propositional logic by constructs to model high-level constraints, and to design solvers that can handle this expanded syntax directly. In the expanded syntax, we obtain more concise representations of search problems. Moreover, these representations are more directly aligned with the inherent structure of the problem. Both factors, we believe, will lead to faster, more effective solvers.

In this paper, we focused on the logic $`\mathrm{𝑃𝐿}^{cc}`$, an extension of propositional logic with c-atoms, that is, direct means to encode cardinality constraints. The specific contribution of the paper are two local-search algorithms $`\mathrm{𝑣𝑏}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ and $`\mathrm{𝑑𝑓}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$, tailored to the syntax of the the logic $`\mathrm{𝑃𝐿}^{cc}`$. These algorithms rely on two ideas. The first of them is to regard a $`\mathrm{𝑃𝐿}^{cc}`$ theory as a compact encoding of a CNF theory modeling the same problem. One can now design local-search algorithms so that they work with a $`\mathrm{𝑃𝐿}^{cc}`$ theory but proceed as propositional SAT solvers would when run on the corresponding propositional encoding. We selected the procedure $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}`$ to establish the correspondence between $`\mathrm{𝑃𝐿}^{cc}`$ theories and CNF encodings, as it does not require any new propositional variables and makes it easy to simulate propositional local-search solvers. We selected a particular propositional local-search method, $`\mathrm{𝑊𝑆𝐴𝑇}`$, one of the best-performing local-search algorithms. Many other choices are possible. Whether they lead to more effective solvers is an open research problem.

The second idea is to change the notion of a flip. We applied it designing the algorithm $`\mathrm{𝑑𝑓}\text{-}\mathrm{𝑊𝑆𝐴𝑇}^{cc}`$ for the class of simple $`\mathrm{𝑃𝐿}^{cc}`$ theories. However, this method applies whenever a $`\mathrm{𝑃𝐿}^{cc}`$ theory $`T`$ can be partitioned into two parts $`T_1`$ and $`T_2`$ so that (1) it is easy to generate random truth assignments satisfying constraints in $`T_1`$, and (2) there is a notion of a flip that preserves satisfaction of constraints in the first part and allows one, in a sequence of such flips, to reach any point in the search space of truth assignments satisfying constraints in $`T_1`$. Identifying specific syntactic classes of $`\mathrm{𝑃𝐿}^{cc}`$ theories and the corresponding notions of a flip is also a promising research direction.

In our experiments we designed compilation techniques to allow us to use SAT solvers in searching for models of $`\mathrm{𝑃𝐿}^{cc}`$ theories. In general, approaches that rely on counting do not work well with $`\mathrm{𝑊𝑆𝐴𝑇}`$, as they introduce too much structure into the theory. The translation $`\mathrm{𝑐𝑜𝑚𝑝𝑖𝑙𝑒}\text{-}\mathrm{𝑏𝑎𝑠𝑖𝑐}`$ is the best match for $`\mathrm{𝑊𝑆𝐴𝑇}`$ (whenever it does not lead to astronomically large theories). All methods seem to work well with $`\mathrm{𝑧𝑐ℎ𝑎𝑓𝑓}`$ at least in some of the cases we studied but none worked well for the vertex-cover problem. To design better techniques to eliminate c-atoms and to make the process of selecting an effective translation systematic rather than ad hoc is another interesting research direction.

Our work is related to  and , which describe local-search solvers for theories in propositional logic extended by pseudo-boolean constraints. However, the classes of formulas accepted by these two solvers and by ours are different. We use cardinality atoms as generalized “atomic” components of clauses while pseudoboolean constraints have to form unit clauses. On the other hand, pseudoboolean constraints are more general than cardinality atoms. At present, we are comparing the performance of all the solvers on the class of theories that are accepted by all solvers (which includes all theories considered here).

## Acknowledgments

The authors are grateful to the reviewers for comments and pointers to papers on extensions of propositional logic by pseudoboolean constraints. This research was supported by the National Science Foundation under Grant No. 0097278.