©Copyright JASSS

JASSS logo ----

J. Gary Polhill, Luis R. Izquierdo and Nicholas M. Gotts (2005)

The Ghost in the Model (and Other Effects of Floating Point Arithmetic)

Journal of Artificial Societies and Social Simulation vol. 8, no. 1

To cite articles published in the Journal of Artificial Societies and Social Simulation, reference the above information and include paragraph numbers if necessary

Received: 09-May-2004    Accepted: 01-Aug-2004    Published: 31-Jan-2005

* Abstract

This paper will explore the effects of errors in floating point arithmetic in two published agent-based models: the first a model of land use change (Polhill et al. 2001; Gotts et al. 2003), the second a model of the stock market (LeBaron et al. 1999). The first example demonstrates how branching statements with floating point operands of comparison operators create a high degree of nonlinearity, leading in this case to the creation of 'ghost' agents -- visible to some parts of the program but not to others. A potential solution to this problem is proposed. The second example shows how mathematical descriptions of models in the literature are insufficient to enable exact replication of work since mathematically equivalent implementations in terms of real number arithmetic are not equivalent in terms of floating point arithmetic.

Agent Based Modelling, Floating Point Arithmetic, Interval Arithmetic, Replication

* Introduction

Bistromathics … is … a revolutionary new way of understanding the behaviour of numbers. … Numbers written on restaurant bill pads within the confines of restaurants do not follow the same mathematical laws as numbers written on any other pieces of paper in any other part of the Universe.

Douglas Adams, Life, the Universe, and Everything (1982)

For some time now, Agent Based Modelling has been used to simulate and explore complex systems, which have proved intractable to other modelling approaches such as mathematical modelling. More generally, computer modelling offers a greater flexibility and scope to represent phenomena that do not naturally translate into an analytical framework. Agent Based Models however, by their very nature, require more rigorous programming standards than other computer simulations. This is because researchers are cued to expect the unexpected in the output of their simulations: they are looking for the 'surprise' that shows an interesting emergent effect in the complex system. It is important, then, to be absolutely clear that the model running in the computer is behaving exactly as specified in the design. It is very easy, in the several thousand lines of code that are involved in programming an Agent Based Model, for bugs to creep in. Unlike mathematical models, where the derivations are open to scrutiny in the publication of the work, the code used for an Agent Based Model is not checked as part of the peer-review process, and there may even be Intellectual Property Rights issues with providing the source code in an accompanying web page.

Even when the code is correct, certain assumptions can be made about the computer platform on which the software will be run that do not necessarily apply, and which can be responsible for effects that are 'surprising' for all of the wrong reasons. The assumption causing the issues raised in this paper is that floating point operations in mathematical processors on computer chips can be treated as though they are operations with real numbers. (Henceforth 'mathematics' and 'arithmetic' will refer to the use of real numbers by default.) That this assumption is false is well-known (Higham 2002), but nothing seems to be done about it in the computer models that use floating point numbers. Indeed, many programming languages, including Java and the Smalltalk platform on which SDML (Strictly Declarative Modelling Language; Moss et al. 1998) is built, do not even supply the library functions stipulated by the IEEE standard for floating point arithmetic that allow the programmer to detect and deal with errors in floating point arithmetic as and when they occur. Issues with floating point arithmetic are not restricted to Agent-Based Social Simulation. McCullough and Vinod (1999) review problems with floating point arithmetic in econometric software, for example, and floating point arithmetic errors have even been responsible for more serious incidents such as the failure of Ariane 5's flight 501 (Lions 1996) and the Patriot missile defence system (GAO 1992; Skeel 1992).

The errors involved in a single floating point arithmetic operation are tiny. Double precision floating point arithmetic guarantees a minimum of fifteen significant decimal figures of accuracy in the result of any one operation (Sun Microsystems 2001). This impressive statistic is perhaps the reason why most programmers make the mistake of assuming that floating point errors can be ignored. A social system on which an Agent-Based Model might be founded is unlikely to involve measurements that would meaningfully require more than three significant figures of accuracy. Indeed, Lawson and Hanson (1995) are cited by Winkler (2003) as claiming that measurement errors are more significant than floating point rounding errors. Beyond that, the abstractions and simplifications of the real system involved in creating the model must surely have a far more significant impact. Appealing though such reasoning may be, this paper shows two examples that challenge it. The first, a simulation of land use change developed by the authors (Polhill et al. 2001; Gotts et al. 2003), is shown to enter invalid states because of floating point errors. The second example shows how some very simple modifications to the code of the well-known Artificial Stock Market (LeBaron et al. 1999) that should be mathematically equivalent, lead to different behaviour in the output variables. A solution is proposed and evaluated using interval arithmetic.

* A model of land use change

The following describes the set up of a model of land use change that has emergent effects arising from errors in floating point arithmetic. The model is based on a particular parameterisation, P, of FEARLUS (Framework for Evaluation and Assessment of Regional Land Use Scenarios), described in Polhill et al. (2001) and Gotts et al. (2003). In what follows, terms referring to objects in the model are capitalised, and italicised on first use. FEARLUS consists of a set of Land Parcels arranged on a 2D grid of squares, each of which is assigned a Land Use by Land Managers (the agents in the model) using a Land Use Selection Algorithm each Year (a cycle in the model). After Land Uses have been assigned, the Yield from each Land Parcel is calculated according to the Land Use selected, and used to derive an amount of Wealth accrued by the Land Manager owning the Land Parcel. This is followed by a process of exchange of Land Parcels between Land Managers. Those Land Managers with negative Wealth put up enough Land Parcels for sale to bring their Wealth to zero or more. Once all Land Managers who need to have put Land Parcels up for sale, a new owner for each such Land Parcel is selected at random from a set consisting of the neighbouring Land Managers with sufficient wealth plus one potential new Land Manager. The cost of the Land Parcel is determined by a model parameter, the Land Parcel Price which is constant over space and time. If a Land Manager has sold all of its Land Parcels, it is deemed to be bankrupt, and plays no further part in the simulation.

Land Managers are divided into Subpopulations according to their Land Use Selection Algorithm. In P, there are two Subpopulations, RS, whose members choose a Land Use at random for each Land Parcel they own, and HR, whose members select a Land Use for each Land Parcel they own in turn as follows: Choose a Land Use at random if the previously selected Land Use lost money in the last Year, but otherwise retain the same non-loss-making Land Use for the current Year. Newly created Land Managers have an equal probability of belonging to HR or RS.

There are two Land Uses to choose from in P, lu1 and lu2. The Wealth accrued to a Land Manager who puts lu1 on a Land Parcel when choosing a Land Use for it is -0.6, whilst for lu2, it is +0.4. The Land Parcel Price is set to 10000 in P, which is a large number intended to prevent Land Managers from ever accruing sufficient Wealth to become the new owner of another Land Parcel, should they be selected to do so if their neighbour goes bankrupt during the 200 Years for which the model is run. The process is summarised in Figure 1.

All Land Managers begin with Wealth 0 after having been assigned their first Land Parcel. At the start of the run, Land Parcels are initially assigned Land Uses at random, with an equal probability of lu1 or lu2 being applied, and Land Managers begin with a single Land Parcel.

Figure 1. UML-style diagrams based on Booch et al. (1999) showing the yearly cycle. At the simplest level of detail, the coloured boxes represent a statechart diagram showing the flow from one state during the yearly cycle to the other, with the dashed line showing the transition from one Year to the next. Inside each coloured box is an activity diagram showing more detail on the processing involved in each state. Arrows with solid lines represent flow of program control from the start (black circle) to the end (black circle with concentric white ring). Diamonds represent decision points, with one out-flowing arrow labelled with text in square brackets indicating the condition under which that branch is used, and the other out-flowing arrow indicating the 'else' branch. The dotted lines between the arrows to/from the concurrent process (thick) lines indicate many objects engaged in the same activity. Arrows with dashed lines are object flows, indicating the involvement of an object in a particular action. Objects are represented in grey boxes divided into possibly three sections. The topmost shows the object and class, underlined and separated by a colon, with the state optionally in square brackets underneath. The next section, which appears optionally, shows certain instance variables of the object, and the bottom-most section, also optional, shows methods. Comments are indicated in a box with a folded down corner, and connected by a dashed line without an arrowhead to the item with which they are associated. In the 'Purchase' activity diagram, the random() method in the Stack class returns a random member of the stack with equal probability of any member being selected.

The expected outcome, given enough time, is for all Land Parcels to be owned by members of HR on a one-to-one basis, each using the profit-making Land Use (lu2); though there is a small non-zero probability that a member of RS will survive for any finite number of Years. However, whilst some runs do yield the expected result in the 200 Year time period, more runs than expected end up with at least one member of RS apparently persisting for many years, and in some cases even acquiring an extra Land Parcel. What is particularly strange is that these members of RS continue using the loss-making Land Use (lu1), but never go bankrupt.

To understand what is happening requires an in-depth analysis of the output of the model, focusing on Land Managers acquiring an extra Land Parcel. In one run (environment size 25 by 25), one such Land Manager, lm879, receives a Land Parcel, lp254 at the end of Year 2 and has initial wealth 0. Over the next 25 Years, lm879 chooses lu2 fifteen times, and lu1 ten times. The accumulated Wealth should be (15 × 0.4) - (10 × 0.6) = 0 at the end of Year 27 and lm879 should not have to sell lp254 -- the code stipulates that this happens only if the Wealth is less than 0 (see Figure 1, in the 'Sale' box). However, accumulated floating point errors in adding 0.4 fifteen times and subtracting 0.6 ten times from the Wealth of lm879 means that instead of 0, lm879 has a Wealth of -2.22045 × 10-16. As a consequence, lm879 puts lp254 up for sale in the 'Sale' step, and receives 10000 units of Wealth. Since 2.22045 × 10-16 is too small a quantity to subtract accurately from 10000 in floating point arithmetic, lm879 now has a Wealth of 10000.

At the same time in the 'Sale' step of Year 27, the owner of a neighbouring Parcel to lp254, lp195, has Wealth -0.6 after choosing lu1 in the 'Choose Land Use' step. Neighbouring Parcel lp195 is then also put up for sale. Once all Land Managers that need to have put their Land Parcels up for sale, new owners are sought. The code for the model assumes that no Land Manager who had put a Parcel up for sale could possibly be eligible to buy one since: (a) A Manager only sells a Parcel if its Wealth is less than zero, (b) on putting a Land Parcel up for sale, a Manager receives the Land Parcel Price, and must therefore have Wealth less than the Land Parcel Price, (c) the only Land Managers eligible to buy Land Parcels are those with Wealth greater than or equal to the Land Parcel Price. Mathematically, this is correct, but as the analysis has shown so far, (a) and (b) do not quite hold when errors in floating point arithmetic are involved. Since lp254 has not yet been fully transferred, lm879 is still the owner of a neighbouring Parcel to lp195,[1] and since lm879 has Wealth 10000, it can also afford to buy it. (See the 'Purchase' box in Figure 1.) With probability 0.5, lm879 is selected over a potential new Land Manager as the new owner of lp195. A further assumption, that only new Land Managers have no Land Parcels when becoming the owner of a Land Parcel (also mathematically correct for real numbers but not for floating point), means that lm879 loses no Wealth on acquiring lp195. Land Managers that have had to sell all their Land Parcels in step 3 are removed from the 'Choose Land Use' schedule step and from the 'Accrue Wealth' schedule step. The software objects that store these Land Managers are not destroyed by the program since they may be used for reporting purposes. The situation at the end of Year 27 is therefore that lm879 owns lp195, has Wealth 10000, is still a valid object in the system, but does not appear in the list of Land Managers to be called in any of the steps of the simulation schedule where 'each Land Manager' appears in the UML diagram in Figure 1.

In Year 28, the new owner of lp254, lm2042, chooses lu1 in the 'Choose Land Use' step, and loses 0.6 units of Wealth in the 'Accrue Wealth' step. Meanwhile, lm879 has not made any selection of Land Use for lp195, and the latter's Land Use from Year 27, lu1, remains unchanged. In the 'Accrue Wealth' step, lm879 is not asked to harvest, since it is not scheduled to do so, and its Wealth remains unchanged despite lu1 being applied to its Land Parcel. In the 'Exchange Land Parcels' step, lm2042 puts lp254 up for sale, and lm879, with Wealth 10000, is an eligible new owner. Once again, the random selection of a new owner is in lm879's favour, and it is chosen as the new owner of lp254 rather than a potential new Land Manager. This time, lm879 does pay 10000 units of Wealth to own lp254, and the situation at the end of Year 28 is that lm879 owns Land Parcels lp195 and lp254, has zero Wealth, and is not on the schedule.

From then on, lm879 remains in this 'ghostly' state until the simulation is terminated, never making any Land Use decision on its accumulated estate of lp254 and lp195, never gaining or losing any Wealth, and never putting either Parcel up for sale. Both Parcels retain the same loss-making Land Use lu1 chosen by the owners who lost them to lm879.

The creation of 'ghost' Land Managers thus requires a number of coincidences to occur:
  1. The real-valued Wealth of a Land Manager, A, must reach exactly zero. For this to happen A must choose lu2 one and a half times as often as lu1 despite both having equal probability of selection by a member of RS. Such sequences of Land Use selection are of length 5i, with 3i selections of lu2 and 2i selections of lu1, where i is a positive integer. A member of HR will not get into this situation because when it chooses lu2 for the first time it will apply it from then on. Since a member of HR choosing lu1 the first time it makes a Land Use decision will go bankrupt, the surviving members of HR are all those that chose lu2 in their first decision.
  2. The floating point error in the Wealth calculation of A must result in a negative number small enough in magnitude that when subtracted from 10000 the result is 10000.
  3. A neighbouring Land Manager, B, must go bankrupt in the same Year.
  4. Land Manager A must then be selected (with probability 0.5) over a potential new Land Manager as the new owner of B's Land Parcel.

To estimate the likelihood of (b), checks were made on all possible sequences of length 25 involving 15 selections of lu2 and 10 selections of lu1 for which the real-valued Wealth is greater than zero until the 25th selection. Of the 112227 possible such sequences, 7496 had the floating point Wealth exactly equal to zero, 7195 had floating point errors resulting in a positive Wealth, and 97536 (about 87%) had floating point errors resulting in a negative Wealth. All of the cases with negative Wealth had the potential to create ghost Land Managers. For shorter sequences, the proportion that have the potential to create ghost Land Managers is smaller, with 75% for sequences of length 20, 54% for length 15, 21% for 10, and 0% for 5.

It is then reasonable to ask, given that some runs have no such problems, whether the problem is sufficiently rare that it can be ignored. To explore this, a model was run consisting solely of members of RS. Given that these Land Managers choose at random, with an equal probability of selecting lu1 and lu2, a time-series graph showing the number of times each Land Use is applied in each Year should show roughly equal applications of each Land Use. Instead the outcome, for any seed, is much like that shown in Figure 2, with the application of lu2 gradually declining as more and more ghost Land Managers occupy the space. This demonstrates that floating point errors have the potential to create a systematic bias in the behaviour of a model.

JASSS-GhostModel-3-fig2a JASSS-GhostModel-3-fig2b
(a) (b)
Figure 2. Systematic bias in a run of the FEARLUS model outlined in Figure 1 using Subpopulation RS only. The space should contain a roughly equal proportion of lu1 (red) and lu2 (green) Land Uses, but by Year 750, the situation is as shown on the left (a). Black lines separate Land Parcels owned by different Land Managers, and the image shows a number of Managers owning two rather than one Land Parcels. The time-series graph in (b) shows the gradual decline in the use of lu2 as more and more ghost Land Managers are created.

This effect has arisen purely because of errors in floating-point arithmetic accumulating through repeated additions of numbers (0.4 and -0.6) that cannot be exactly represented in binary, which allowed the model to enter states that were not anticipated because they are not possible mathematically. Fortunately, none of the results published using this model (Polhill et al. 2001; Gotts et al. 2003) have used such settings, and their validity is unaffected by problems with floating-point arithmetic documented here. This example does raise the issue, however, of how sensitive an agent-based model can be to floating point errors. Any branching statement in the agent's decision mechanism (here in the decision to sell land) containing a comparison expression involving a floating-point variable creates a high degree of nonlinearity with respect to that variable.

It could be argued that, rather than being an issue with the implementation of the model in a digital computer, this is instead an issue with the design. Perhaps a different abstraction of the subset of reality that inspired the model shown here would not suffer from problems with floating point arithmetic, and if so, perhaps the alternative design should be seen as being better than the one used here. Whatever faults the design used here might have, the only issue that relates to the problems it has with floating point arithmetic concerns the fact that it relies on certain mathematical laws applying in the implementation. Surely we would not wish to assert that all models that rely on the truth of mathematical laws are badly designed? Neither does it seem sensible to assert that if two models A and B are abstractions of the same subset of reality and A relies on the truth of mathematical laws whilst B does not, then B is necessarily superior to A. Given that certain mathematical laws do not always hold in floating point arithmetic, perhaps a better approach would be to check in the implementation of a model that does rely on mathematical laws whether those laws have been broken, or to make some measurement of the extent to which they might have been broken. The IEEE 754 (1985) standard for floating point arithmetic, discussed later, stipulates facilities that permit either of these approaches to be adopted.

As this example has shown, the 15 decimal places of accuracy provided by double precision variables was not enough. It is also a salutary lesson in the dangers of assuming that the laws of arithmetic with real numbers apply in floating point arithmetic. For now, we can avoid these problems through using parameters that can be exactly represented in binary, but remedial steps will need to be taken to manage errors in floating point arithmetic should a scenario require parameters to be used that do not meet this demanding requirement. One approach could be to explicitly represent the state of a Land Manager in the code (e.g. 'Newly created', 'Bankrupt', 'Scheduled') and perform sanity checks each time the state changes. Whilst such an approach would prevent the creation of ghost Land Managers, it would not prevent the Wealth from dropping below zero prematurely due to accumulated floating point errors. Indeed, had the code ensured that only Land Managers whose state is 'Scheduled' or 'Newly created' could become the new owners of Land Parcels rather than relying only on a comparison with the Wealth, it is unlikely that the issues with floating point arithmetic would ever have been detected, even though some Land Managers would be bankrupted prematurely.

Another approach to addressing these issues would be to apply small constants to the comparisons, effectively in this case meaning that a Land Manager only becomes bankrupt if their Wealth is less than -k, where k is the small constant in question. Knuth (1998, p. 233) compares various approaches to this, and makes his own recommendations, which have been implemented by Theodore Belding in the fcmp package available at http://fcmp.sourceforge.net. Whilst this would probably work in the case above, it is not a generally safe approach in that it does not guarantee any model will behave correctly. This is explored further in Polhill et al. (in press). For this reason we attempted a conversion to interval arithmetic, discussed later.

* The Artificial Stock Market

The Artificial Stock Market (ASM; LeBaron et al. 1999) is a well-known agent-based model that has been used by other authors as a basis for their work (e.g. Chen & Yeh 2002). In this model, agents have a choice between investing in a safe bond with a fixed interest rate, and a risky stock in limited supply with an autoregressive stochastic dividend dt calculated as per equation (1) (LeBaron et al. 1999):

dt = D + r(dt - 1 - D) + mt (1)

where D is the average dividend (a model parameter), mt is a stochastic component sampled from a Normal distribution with zero mean and constant variance, and r is the autocorrelation parameter: if 0 then dt and dt - 1 are uncorrelated, and if 1, dt and dt - 1 are maximally correlated.

Agents adapt their strategies for choosing how much to invest in a series of repeated cycles using a Genetic Algorithm. The ASM has been shown to replicate features in time-series from real stock markets, such as weak forecastability and volatility persistence (LeBaron et al. 1999, p. 1512). These phenomena cannot be explained by the well-established Efficient Market Hypothesis, under which prices only change when there is new relevant information, meaning that they are not forecastable at all and are less volatile than prices observed in real markets; from this Shleifer (2000) argues that, "More than news seems to move stock prices." (p. 20).

The source code for the artificial stock market, in its original Objective-C format, is available in a package called 'sfsm' at http://artstkmkt.sourceforge.net/. A version of that model, dubbed ASM, has been created using the Swarm simulation system (as discussed in Johnson 2002). An earlier version of this paper (Polhill et al. 2003) used version 2.2.1 of ASM, but since then version 2.4 has been released, with the assistance of Badegruber (2003). This paper uses version 2.4, which Johnson claims on the above website to be the first truly accurate Swarm version of the original artificial stock market that LeBaron et al. published.

With version 2.2.1, we had to create a baseline version that had a separate seed for the stochastic component mt of the dividend to prevent any other changes we made from having a disproportionate effect on the output. With version 2.4 this alteration is not necessary, as Badegruber and Johnson have made provision for it. However, in order to enable compilation of ASM-2.4 on a PC running Windows XP with Swarm 2.2 pretest 10, we did have to comment out line 12 of Dividend.m, which declares an external function random() to have a type that conflicts with the declaration in the stdlib.h supplied with Cygwin. This, then, is taken as the baseline version of the ASM.

Two new versions were created from the baseline to test the sensitivity of the ASM to errors in floating point arithmetic. Version 1 modifies the implementation of equation (4) from LeBaron et al.'s paper (p. 1491) (Figure 3), and Version 2 modifies the implementation of equation (14) (p. 1496) (Figure 4). Neither modification should make any difference to the behaviour of the software according to the laws of mathematics. Version 1 computes (A/B) - (C/B) rather than (A - C)/B, and Version 2 computes Z + Y + … + A, rather than A + B + … + Z.

Figure 3. Modifications required to create Version 1 of the ASM, involving a change to line 733 of the Objective-C source file BFagent.m. Since division is distributive over subtraction, this change should have no effect according to the laws of arithmetic.

Figure 4. Modifications made to the ASM for Version 2 highlighted in bold. (Some irrelevant lines, which are the same in both versions, have been removed.) Essentially, the sum computed in equation (14) in LeBaron et al. (1999) to clear the market for the risky stock is totalled in reverse order of the agents. Since addition is associative, the order the agents are called to make the sum should not make any difference to the behaviour of the model according to the laws of arithmetic.

The results are summarised in Figure 5, which apply to running the ASM with seed 123654789, and the original parameter file asm.scm supplied with the source code. Using a graph plotted at a resolution of 25 time steps, the baseline, Version 1 and Version 2 show a difference in the behaviour of the volume at about the 2000 cycle point, and despite having the same dividend, differences continue to be observed thereafter. To get a precise picture of where the differences occurred, we ran the ASM in batch mode, creating an output data file recording the price, risk free price, and volume each time step to six decimal places. The first difference in output between the baseline and Version 1 occurs at step 1926, whilst that between the baseline and Version 2 occurs at step 2007. Using the recorded output files also enabled us to confirm that the behaviour is initially the same (to six decimal places), and hence that the modifications made in Versions 1 and 2 have not substantially altered the mathematics of the model. The differences in output observed are thus attributable purely to the cumulative effects of the changed floating point rounding errors caused by rearranging the terms in the expressions in alternative but mathematically equivalent ways.

Figure 5. Divergence of behaviour in various versions of the ASM. The same output is seen for the first 2000 cycles or so (though the graphs begin at cycle 1000), but thereafter there are differences in the number, height and width of peaks. The area highlighted in red shows the first differences observed between the versions. Version 1 shows the greater difference from the baseline version in this region, with two thin peaks where the baseline has a single wide peak, whereas Version 2 also features a wide peak but with a flatter top than the baseline.

The statistical properties of the output that the authors use to claim a match with properties of real markets seem by eye not to be affected, as qualitatively the graphs have a similar structure. Formal verification that the statistical properties are unaffected is the potential subject of future work. Though unlikely, it is not inconceivable for such a verification to fail: rounding errors in floating point numbers are not random (Higham 2002), leaving open the possibility of introducing a bias through implementing an equation one way rather than another. This may be related to the systematic bias seen in the RS-only run of the FEARLUS model discussed above.

The fact that implementing the mathematical formulae in LeBaron et al.'s paper in different but mathematically equivalent ways changes the output of the model raises a question mark over any exact replication of their work. At the platform level, compiler optimisation settings can affect the way that expressions are evaluated (Goldberg 1991), meaning that potentially, different outcomes could occur on different machines even when using the same source code. However, the modifications introduced were intended to illustrate the possible effect of implementing the ASM from its description in the literature rather than by using the authors' source code. If nothing else, the differences in behaviour highlight how critical making source code available to other researchers can be, since otherwise exact replication of results would be difficult if not impossible.[2]

* Floating point arithmetic and the IEEE 754-1985 standard

A real number n can be represented in any base b without loss of precision using two variables, a fraction, f, also a real number, known as the mantissa or significand, and an integer exponent, e, such that n = fbe, with f being smaller in magnitude than b.

In a computer, there are a discrete set of values that the binary significand can take, meaning that most fractions are not exactly representable. For example, 0.1 (in decimal) is not an exactly representable number (Wallis 1990a). Indeed, of all the fractions 0.1i, i = 1, …, 9, only 0.5 is exactly representable. An arithmetic operation (plus, minus, multiply, divide) on two floating point numbers is also not necessarily going to result in a number that is exactly representable (e.g. 1 / 3). There are thus three potential sources of error in a single arithmetic calculation: conversion from a decimal to a binary number for each of the two operands, and a rounding error from the exact result to a representable result.

Historically, computer manufacturers implemented floating point arithmetic in various ways, creating difficulties with porting code. Programs that worked in one machine could cause exceptions to be raised on another, even with simple operations such as setting y = -x or y = x / x (Wallis 1990b). The IEEE 754 standard (IEEE 1985) changed this, and now it is unusual to find a non-conforming floating-point unit in a computer.

The IEEE standard provides strict definitions for 32-bit (single) and 64-bit (double) precision floating point formats, and a more flexible definition for extended formats. It stipulates that the accuracy of arithmetic operations be such that the nearest representable number to the infinitely precise result be delivered to the destination.[3] It requires that the user be able to set the rounding mode to one of four possibilities: round-to-nearest even (the default), round towards positive infinity, round towards negative infinity, and round-to-zero (truncate). It further requires that conforming architectures generate trappable exceptions in various cases where the result of a calculation has not been exactly representable in the destination format. This last requirement enables users to check when mathematical laws might not have been strictly adhered to.

The standard is not designed to enable precise reproducibility on all platforms, however. The format of the destination is not stipulated to be the format suggested by the type of the variable. For example, Intel floating-point units work using 80-bit double extended format by default, even if a variable is defined by a programming language to be a 64-bit variable. Although the floating-point unit can be configured to do calculations in 64-bit arithmetic (Intel 2003), compilers do not necessarily issue instructions to do this (e.g. Borland C compiler, and the GNU C compiler in Cygwin). For example, the calculation 253 / (253 - 1) is, in binary, the repeating fraction 1.0…52…010…52…01 …, where "0…x…0" is replaced with x zeros. The double precision 64-bit format has 53 bits for the significand, and therefore the nearest representable number is 1.0…51…01, rounding up the least significant bit. However, in Intel's 80-bit arithmetic, with 64 bits for the significand, the result delivered to the destination is 1.0…52…010…10…0, which is then rounded to the 64-bit variable as 1.0…51…00. Nevertheless, Intel's floating point architecture is fully IEEE 754 compliant.

One platform that doesn't fully comply with IEEE 754 is Java (Kahan & Darcy 2001). Java provides no facilities for checking exception flags for any errors in floating point computations, and no facilities for changing the rounding mode. Since many popular agent-based modelling platforms are written in Java (Ascape (Parker 1998), RePast (Collier 1996)) or provide Java interfaces (Swarm (Minar et al. 1996)), this is bound to be contentious. Without these facilities, Java programmers are severely impaired: they cannot easily check whether or not a series of calculations has caused an inaccurate result, and they cannot directly implement interval arithmetic (see below). It is possible to get equivalent functionality in Java, but not without considerable effort on the part of the programmer (Darcy pers. comm.). For interval arithmetic, this involves tweaking the least significant bits of the significand to generate the closest floating point numbers on either side of the floating point result. This approach does not achieve bounds as tight as those obtained using hardware rounding, as mentioned by Abrams et al. (1998) in a paper comparing various techniques for implementing interval arithmetic in C. Detecting rounding errors can be achieved using Dekker's (1971) algorithms to make computations that effectively double the available precision using a low order and a high order floating point number of the available precision to represent the greater precision number. If, after performing a computation, the low order number has value zero, then, according to Darcy (pers. comm.), the computation in the available precision would have been computed exactly.

* Interval arithmetic

A comparison of various techniques for handling errors in floating point errors concluded that interval arithmetic offered a safe approach (Polhill et al. in press) in that it could be used to provide a warning when a comparison operator might deliver the wrong result due to accumulated floating point errors. The idea behind using interval arithmetic is always to ensure that the true value of a calculation is within known bounds. Rather than just checking whether mathematical laws have been adhered to, it can be used to provide an indication of the extent to which they have not. Of the available forms that interval arithmetic can take, Ullrich & von Gudenberg (1990) suggest that those representing the bounds explicitly are most successfully implemented in a computer. In particular, it is possible to use the IEEE 754 stipulated rounding mode functions to compute a minimum and maximum result for each calculation: one computed using round towards negative infinity, and the other using round towards positive infinity.

Let an interval be represented by [r, s], where r <= s. Let <(x) represent the largest floating point number that is not more than x, and >(x) represent the smallest floating point number that is not less than x, then the arithmetic operators are defined (Alefeld & Herzberger, 1983):

[r1, s1] + [r2, s2] = [<(r1 + r2), >(s1 + s2)] (2)

[r1, s1] - [r2, s2] = [<(r1 - s2), >(s1 - r2)] (3)

[r1, s1] × [r2, s2] = [min{<(r1 × r2), <(r1 × s2), <(s1 × r2), <(s1 × s2)},
max{>(r1 × r2), >(r1 × s2), >(s1 × r2), >(s1 × s2)}]

[r1, s1] / [r2, s2] = [min{<(r1 / r2), <(r1 / s2), <(s1 / r2), <(s1 / s2)},
max{>(r1 / r2), >(r1 / s2), >(s1 / r2), >(s1 / s2)}]

There are then two forms for the comparison operators, echoing modal logic (Kripke 1963): a 'necessarily' form where the comparison is true iff it is true for all pairs of members of each operand, and a 'possibly' form, where the comparison is true iff it is true for one or more pairs of members of each operand (Alefeld & Herzberger 1983).

Let >L represent 'necessarily greater-than' and <M represent 'possibly less-than', and similarly for the other comparison operators. Then:

[r1s1] >L [r2s2] iff r1 > s2; and [r1s1] >M [r2s2] iff s1 > r2 (6)

[r1s1] <L [r2s2] iff s1 < r2; and [r1s1] <M [r2s2] iff r1 < s2 (7)

[r1s1] =L [r2s2] iff (r1 = s2) & (r2 = s1); and [r1s1] =M [r2s2] iff (s1 >= r2) & (r1 <= s2) (8)

In a section of code where x > y leads to an action, it is possible to ensure that, regardless of floating point errors, the action only takes place when it should by replacing x and y with intervals containing the true values, and performing the action if x >L y. To check that the action is only not executed when it should not be, then after x >L y returns false, there should be a further check whether x >M y. If the latter is true, a warning should be issued to the effect that there is no longer certainty that the program has executed the correct instructions in accordance with the model design. A similar exercise for other comparison operators ensures that no action takes place unless it necessarily should and no action does not take place that possibly should without a warning being given. A model run can then be known to have executed without deviating from the design due to floating point errors given the absence of a warning.

Implementing intervals in FEARLUS should be a relatively trivial exercise, since FEARLUS only makes use of the basic arithmetic operators in calculations. However, this has not turned out to be the case. There are a number of places in the FEARLUS code where a weighted selection is made between various alternatives, in which the chance of selecting an alternative is given by its weight divided by the sum of the weights of all the alternatives.

A trivial example is in selecting which Subpopulation a newly created Land Manager will belong to. This is not an issue in the example above because the probability of belonging to HR or RS is 0.5, which is exactly representable in binary floating point. If there were three subpopulations to chose from with equal probability, then an issue would arise because 1/3 would have to be represented using the interval [<(1/3), >(1/3)]. More generally, however, consider a choice between two alternatives U and V, where the weight for U is u = [1, 2], and that for V is v = [2, 3]. Looking at the numbers, the probability of choosing U, P(U), has then a minimum of 0.25 when u = 1 and v = 3, and a maximum of 0.5 when u = v = 2. Similarly P(V), is in the interval [0.5, 0.75]. Since there is uncertainty in their relative probabilities, it is not clear how to choose between U and V in an unbiased manner.

Interestingly, an issue with weighted selection also occurs when adding small constants to comparison operations as discussed at the end of Section 2, at least using the algorithm to implement weighted selection in FEARLUS. This algorithm works as follows: If there are n options to choose from, each option i = 1, …, n having probability P(i) of being selected such that the sum of all P(i)s is equal to 1, then construct a vector c with elements ci set to the sum of all P(j)s for j = 1, …, i. Next, generate a random number r in the range [0, 1] inclusive, and select the option s with minimum value such that cs > r. However, when a small constant is added to the comparison operation, this becomes effectively cs > r + k, which would effectively bias against earlier options. Though this could be addressed by not applying the constant in this case, it serves as an illustration that applying techniques such as Belding's fcmp to floating comparisons without careful consideration of the implications would be unwise.

Since the weighted selection issue does not apply to the use of FEARLUS in this paper, an interval arithmetic version (0-7) of FEARLUS was created. An Objective-C class was used to implement interval arithmetic operations, and all floating point variables were replaced with instances of this class. Expressions containing floating point variables used method calls instead of mathematical operations. When a comparison is made, warnings are issued as described above. Using this version of FEARLUS no ghost agents are ever created before a warning is issued, and to this extent, interval arithmetic has proved successful.

In the Artificial Stock Market, 250,000 cycles of the model are used to enable agents to learn before generating the time-series on which the analysis was based (LeBaron et al. 1999, p. 1499). Given that simply rearranging terms in a single mathematical expression changed the behaviour of the volume as early as 2000 or so cycles into the model, it is conceivable that after 250,000 cycles, there could be a considerable number of comparisons in which the outcome is uncertain were interval arithmetic to be used.

Approaches to managing floating point errors have various properties that may need to be considered in evaluating them. Approaches such as interval arithmetic are safe in that they ensure the model behaves as though it is using real numbers until a warning is given. This might impose an unacceptable constraint on the model in that it will not be possible to run it for the required time period without a warning being given. A less stringent requirement could be robustness, in which the floating point error management approach could provide some guaranteed protection against the model entering unanticipated or invalid states through floating point errors. In FEARLUS, for example, ghost Land Managers could be prevented by explicitly representing the state of a Land Manager as 'bankrupt' and then preventing bankrupt Land Managers from buying Land Parcels. However, this would not protect against Land Managers with zero real-valued wealth from going bankrupt when they shouldn't.

Another potential property an approach might have is platform independence, which would facilitate exact repetition of results. At the very least, however, floating point error management techniques should not allow a systematic bias in the behaviour of the model.

* Conclusion

In an ideal world, floating point arithmetic should not be something that agent-based modellers need to trouble themselves with. This paper has demonstrated that the world is far from ideal, through two examples from published agent-based models. The first example illustrates how branching in agent decision mechanisms based on comparisons involving floating point variables can lead to emergent effects with an entirely unwelcome element of 'surprise' because of the high degree of nonlinearity this introduces, and the seeming consequent demand for infinite precision. The second example shows the necessity for authors to release their source code in public domain if there is to be any chance of other researchers repeating their results with any degree of precision. Mathematical equations in models can be implemented in a number of different ways in a computer program, and although these may be mathematically equivalent, they are not equivalent in terms of floating point arithmetic.

Traditionally, a number of heuristics have been used for dealing with errors in floating point arithmetic, such as using small tolerances when making comparisons, or writing the result of an operation to a string buffer with a degree of accuracy not greater than that provided by the floating point numbers and then copying the result out of the buffer into the floating point variable. Whilst such measures provide some ability to cope with the problems of floating point errors, they do not provide any certainty that a given run has strictly followed the model design (Polhill et al. in press). Such certainty is provided using interval arithmetic: it can be known at the point of comparison that, even allowing for floating point errors, the correct course of action is being taken according to the design of the model. However, this is no panacea: use of intervals creates difficulties for certain algorithms such as weighted selection between alternatives and sorting; it also limits the number of cycles a model can be run before accumulated floating point errors mean there is uncertainty about whether the model is behaving according to the design.

David Hales, in a presentation to the 2003 Model-2-Model workshop (Edmonds & Hales 2003), noted that reimplementation of agent-based models should not use the original source code, as they would then introduce the same artefacts. A distinction needs to be made between reimplementation and repetition in this context. A supposed advantage of using computer simulation is that experiments can be exactly repeated, allowing the confirmation of the results obtained by other authors, including errors and artefacts. Such functionality is not typically available in the natural sciences -- it is like being able to reuse the apparatus and samples to take exactly the same measurements as in the original experiment. Reimplementation is, of course, a more rigorous test of a reported effect, but this does not necessarily mean that repetition is without value. Repetition is an exercise researchers can undertake themselves with relatively little effort. Compiling the model on different platforms (e.g. a dual-boot Linux/Windows PC), at least provides a check that results are not dependent on such things as operating system, compiler and software library versions.

* Acknowledgements

This work is funded by the Scottish Executive Environment and Rural Affairs Department. The authors would also like to express their gratitude to Paul Johnson for some very helpful comments on an earlier draft of this paper, and Joe Darcy, Java's floating point guru, for his advice. This work was first presented at ESSA-1/SIMSOC-VI Conference, 19-21 September 2003, Groningen, Netherlands.

* Notes

1Technically, lm879 has no members on its Land Parcels owned stack (they have all been removed during the 'Sale' step), but the owner property of lp254 is still set to lm879 (and is not changed until the Land Parcel is found a new owned in the 'Purchase' step). See Figure 1.

2In an earlier version of this paper (Polhill et al. 2003) we commented that Epstein & Axtell's (1996) SugarScape features a number of mathematical expressions with potential for floating point issues and yet the source code is not publicly available. The latter may be incorrect. Although Leigh Tesfatsion states on the Agent-Based Computational Economics webpages (http://www.econ.iastate.edu/tesfatsi/acecode.htm) that “the Object Pascal source code for SugarScape is not publicly available”, in a working paper (Axtell, Axelrod, Epstein & Cohen 1995), the authors do state that the code is available from Robert Axtell, and Tobias & Hofmann (2004) classify SugarScape as having an open source licence in the appendix to their paper. In the final draft we therefore decided to give them the benefit of the doubt and have removed the comment here.

3However, this is the precise result of operating on the floating point operands, not on the original numbers the operands may have been rounded from. Thus, for example, the result of 0.4 × 0.1 in IEEE 754 floating point arithmetic is not the nearest representable number to 0.04, but the nearest representable number to the product of the nearest representable number to 0.4 and the nearest representable number to 0.1.

* Software

The software used to generate the results in this paper is written in Objective-C for the Swarm libraries. It is released under the GNU General Public Licence, and is available for download at http://www.macaulay.ac.uk/fearlus/floating-point/ghost/.

To reproduce the results in section 2, download FEARLUS model 0-6-2 source code and the associated parameters from the table in the above webpage. This is known to work on a Windows 2000 PC with Swarm 2.1.1, and a Sun Sparc running Solaris 8 with Swarm test version 2001-12-18. Create a directory in which you put the parameter and source code compressed tar files, and decompress them using the following commands from a Cygwin or X terminal:

	gunzip -c FEARLUS-model0-6-2.tar.gz | tar xf -
	gunzip -c FEARLUS-fpdemo.tar.gz | tar xf -

Compile the source code with the following commands, ensuring that your SWARMHOME environment variable points to the appropriate location:

	cd model0-6-2

To run the FEARLUS model as described in Figure 1, issue the following commands:

	cd ../demo/fp
	../../model0-6-2/model0-6-2 -S 1552078961 -p RH-fp.model -o demo.obs

Not all seeds generate ghost land managers. To try other seeds, replace the -S 1552078961 in the command above with -s. If you are using Swarm 2.1.1, note that -S should be replaced with -X.

To get a result similar to that shown in Figure 2, do:

	../../model0-6-2/model0-6-2 -s -p R-fp.model -o demo.obs

The interval arithmetic version of FEARLUS (0-7) is also available from the above website. The default version to download is FEARLUS-model0-7.tar.gz, but on PCs, the Cygwin environment supplied with Swarm 2.1.1 does not implement the C API hardware rounding function fpsetround(). For PCs with a Pentium CPU, you can download FEARLUS-model0-7ci.tar.gz which contains some extra code to implement this function.

The code used to estimate the likelihood of condition (b) in Section 2 is supplied in lu1lu2sequences.c. It can be compiled and run thus:

	gcc -o lu1lu2sequences lu1lu2sequences.c

The results given here apply to running this program on a Sun Sparc platform running Solaris 2.8 and quoting the line of output beginning 'i = 5' (the program loops rather optimistically until i = 40, so will need to be interrupted after this line has been displayed). If you run the program on an Intel PC, it is likely you will get different results for the reasons discussed in Section 4. For example, on one Pentium 4 PC we got 7496 having Wealth exactly equal to zero and 7195 with floating point errors resulting in a positive Wealth (as per the Sun), but only 27144 (about 19%) with a negative Wealth small enough that when subtracted from 10000 the result is 10000, as opposed to the 97536 on the Sun. Here then, Intel's 80-bit arithmetic has worked out favourably.

The ASM version 2.4 requires a later version of Swarm. We used Swarm 2.2 pretest 10 on a Windows XP PC to generate the results in Figure 5. To generate the results shown here, either download the collection of premodified versions of the source code from the website given at the beginning of this section, or download the original ASM-2.4 from http://artstkmkt.sourceforge.net and follow the instructions in Section 3 for modifying the code. Once you have the three versions of the software (baseline, Version 1 and Version 2) in separate directories, compile them with:


and run each one from its own directory with the command:
	./asm -S 123654789

In the ObserverSwarm pop-up window, you can change the time interval by setting the 'displayFrequency' parameter. The graphs in Figure 5 were generated using an interval of 25.

To generate output data files with 1 step time intervals, modify the parameter file asm.scm in each version's directory so that the BatchSwarm section reads:

	(cons 'asmBatchSwarm
	        (make-instance 'ASMBatchSwarm
	  #:loggingFrequency 1 
	  #:experimentDuration 10000

Then run each version model using:

	./asm -b -S 123654789

The output will appear in a file beginning with the name output.data. This can be compared with the output from other versions using the diff command.

* References

AXTELL, R, AXELROD, R, EPSTEIN, J M, COHEN, M D (1995) Aligning simulation models: A case study and results. Santa Fe Institute Working Paper 95-07-065. Available on-line at http://www.santafe.edu/sfi/publications/wpabstract/199507065.

ABRAMS, S L, CHO, W, HU, C-Y, MAEKAWA, T, PATRIKALAKIS, N M, SHERBROOKE, E C, & YE, X (1998) Efficient and reliable methods for rounded-interval arithmetic. Computer-Aided Design 30 (8) pp. 657-665.

ALEFELD, G, & HERZBERGER, J (1983) Introduction to Interval Computation. New York: Academic Press.

BADEGRUBER, T (2003) Agent-Based Computational Economics: New Aspects in Learning Speed and Convergence in the Santa Fe Artificial Stock Market. Inaugural dissertation, Universitat Graz, Graz, Austria.

BOOCH, G, RUMBAUGH, J, & JACOBSON, I (1999) The Unified Modeling Language User Guide. Boston, MA: Addison-Wesley.

CHEN, S-H, & YEN, C-H (2002) On the emergent properties of artificial stock markets: the efficient market hypothesis and the rational expectations hypothesis. Journal of Economic Behavior & Organization 49 pp. 217-239.

COLLIER, N (1996-) RePast: An extensible framework for agent simulation. http://repast.sourceforge.net/.

DEKKER, T J (1971) A floating point technique for extending the available precision. Numerische Mathematik 18 pp. 224-242.

EDMONDS, B, & HALES, D (2003) Replication, replication and replication: Some hard lessons from model alignment. Model to Model International Workshop, GREQAM, Vielle Charité, 2 rue de la Charité, Marseille, France. 31 March - 1 April 2003. pp. 133-150.

EPSTEIN, J M, & AXTELL, R (1996) Growing Artificial Societies: Social Sciences from the Bottom Up. Cambridge, MA: MIT Press.

GAO. 1992. Patriot Missile defense: Software problem led to system failure at Dhahran, Saudi Arabia. GAO/IMTEC-92-26, Washington, DC: U.S. Government Accountability Office. Available on-line at

GOLDBERG, D (1991) What every computer scientist should know about floating point arithmetic. ACM Computing Surveys 23 (1) pp. 5-48. Reproduced in Sun's Numerical Computing Guide (Appendix D) and available on-line at http://portal.acm.org/citation.cfm?doid=103162.103163

GOTTS, N M, POLHILL, J G, LAW, A N R (2003) Aspiration levels in a land use simulation. Cybernetics & Systems 34 (8) pp. 663-683

HIGHAM, N J (2002) Accuracy and Stability of Numerical Algorithms, 2nd Edition. Philadelphia, USA: SIAM.

IEEE (1985) IEEE Standard for Binary Floating-Point Arithmetic, IEEE 754-1985, New York, NY: Institute of Electrical and Electronics Engineers.

INTEL (2003) IA-32 Intel Architecture Software Developer's Manual Volume 1: Developer's Manual. http://www.intel.com/design/Pentium4/manuals/.

JOHNSON, P E (2002) Agent-Based Modeling: What I learned from the Artificial Stock Market. Social Science Computer Review 20 pp. 174-186.

KAHAN, W, & DARCY, J D (2001) How Java's floating point hurts everyone. Originally presented at the ACM 1998 Workshop on Java for High-Performance Network Computing, Stanford University. Revised and updated version (5 Nov 2001, 08:26 version used here) available online at http://www.cs.berkeley.edu/~wkahan/JAVAhurt.pdf.

KNUTH, D E (1998) The Art of Computer Programming Vol. 2. Seminumerical Algorithms. Third Edition. Boston, MA: Addison-Wesley.

KRIPKE, S (1963) Semantical considerations on modal logic. Acta Philosophica Fennica 16, pp. 83-94.

LAWSON, C L, & HANSON, R J (1995) Solving Least Squares Problems. Philadelphia, PA: SIAM.

LEBARON, B, ARTHUR, W B, & PALMER, R (1999) Time series properties of an artificial stock market. Journal of Economic Dynamics & Control 23 pp. 1487-1516.

LIONS, J-L (1996) Ariane 5 Flight 501 Failure: Report by the Inquiry Board. Paris: European Space Agency. Available online at http://ravel.esrin.esa.it/docs/esa-x-1819eng.pdf.

MCCULLOUGH, B D, & VINOD (1999) The numerical reliability of econometric software. Journal of Economic Literature 37 (2) pp. 633-665.

MINAR, N, BURKHART, R, LANGTON, C, ASKENAZI, M (1996-) The Swarm simulation system: A toolkit for building multi-agent systems. http://wiki.swarm.org/

MOSS, S, GAYLARD, H, WALLIS, S, & EDMONDS, B (1998) SDML: A multi-agent language for organizational modelling. Computational & Mathematical Organization Theory 4(1) pp. 43-69.

PARKER, M (1998-) Ascape. http://www.brook.edu/dybdocroot/es/dynamics/models/ascape/. The Brookings Institution.

POLHILL, J G, GOTTS, N M, & LAW, A N R (2001) Imitative and nonimitative strategies in a land use simulation. Cybernetics & Systems 32 (1-2) pp. 285-307.

POLHILL, J G, IZQUIERDO, L R, GOTTS, N. M. (2003) The ghost in the model (and other effects of floating point arithmetic). Proceedings of the first conference of the European Social Simulation Association, held in Groningen, Netherlands, 18-21 September 2003, published online at http://www.uni-koblenz.de/~kgt/ESSA/ESSA1/.

POLHILL, J G, IZQUIERDO, L R, GOTTS, N. M. (in press) What every agent based modeller should know about floating point arithmetic. Environmental Modelling and Software.

SHLEIFER, A (2000) Inefficient Markets: An Introduction to Behavioral Finance. Oxford UP.

SKEEL, R. (1992) Roundoff error cripplies Patriot missile. SIAM News 25 (4) p. 11. Available online at http://www.siam.org/siamnews/general/patriot.htm.

SUN MICROSYSTEMS (2001) Numerical Computing Guide. http://docs.sun.com/source/806-3568/.

TOBIAS, R & HOFMANN, C (2004) Evaluation of free Java-libraries for social-scientific agent based simulation. Journal of Artificial Societies and Social Simulation 7 (1). http://jasss.soc.surrey.ac.uk/7/1/6.html

ULLRICH, C, & VON GUDENBERG, J W (1990) Different approaches to interval arithmetic. In Wallis, P. J. L. (Ed.) Improving Floating Point Programming, ch. 5. Chichester, UK: John Wiley & Sons.

WALLIS, P J L (1990a) Basic concepts. In Wallis, P J L (Ed.) Improving Floating Point Programming, ch. 1. Chichester, UK: John Wiley & Sons.

WALLIS, P J L (1990b) Machine arithmetic. In Wallis, P J L (Ed.) Improving Floating Point Programming, ch. 2. Chichester, UK: John Wiley & Sons.

WINKLER, J R (2003) A statistical analysis of the numerical condition of multiple roots of polynomials. Computers & Mathematics with Applications 45 (1-3) pp. 9-24.


ButtonReturn to Contents of this issue

© Copyright Journal of Artificial Societies and Social Simulation, [2005]