1% BEGIN LICENSE BLOCK 2% Version: CMPL 1.1 3% 4% The contents of this file are subject to the Cisco-style Mozilla Public 5% License Version 1.1 (the "License"); you may not use this file except 6% in compliance with the License. You may obtain a copy of the License 7% at www.eclipse-clp.org/license. 8% 9% Software distributed under the License is distributed on an "AS IS" 10% basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See 11% the License for the specific language governing rights and limitations 12% under the License. 13% 14% The Original Code is The ECLiPSe Constraint Logic Programming System. 15% The Initial Developer of the Original Code is Cisco Systems, Inc. 16% Portions created by the Initial Developer are 17% Copyright (C) 2006 Cisco Systems, Inc. All Rights Reserved. 18% 19% Contributor(s): 20% 21% END LICENSE BLOCK 22 23\chapter{Working with real numbers and variables} 24\label{chapreal} 25%HEVEA\cutdef[1]{section} 26 27\setcounter{topnumber}{1} 28 29% - larger field/lake example 30 31\section{Real number basics} 32 33In general, real values cannot be represented exactly if the representation 34is explicit. As a result, they are usually approximated on computers by 35floating point numbers, which have a finite precision. This approximation 36is sufficient for most purposes; however, in some situations it can lead to 37significant error. Worse, there is usually nothing to indicate that the 38final result has significant error; this can lead to completely wrong 39answers being accepted as correct. 40 41One way to deal with this is to use \emph{interval arithmetic}. 42 \index{interval arithmetic}% 43The basic idea is that rather than using a single floating point value to 44approximate the true real value, a pair of floating point bounds are used 45which are guaranteed to enclose the true real value. Each arithmetic 46operation is performed on the interval represented by these bounds, and the 47result rounded to ensure it encloses the true result. The result is that 48any uncertainty in the final result is made explicit: while the true real 49value of the result is still not known exactly, it is guaranteed to lie 50somewhere in the computed interval. 51 52Of course, interval arithmetic is no panacea: it may be that the final 53interval is too wide to be useful. However this indicates that the problem 54was probably ill-conditioned or poorly computed: if the same computation had 55been performed with normal floating point numbers, the final floating point 56value would probably not have been near the true real value, and there would 57have been no indication that there might be a problem. 58 59\index{bounded reals|(} 60In \eclipse{}, such intervals are represented using the \emph{bounded real} 61data type. 62 63\quickref{Bounded reals}{ 64\begin{itemize} 65\item Bounded reals are written as two floating point bounds separated 66 by a double underscore (e.g.\ \texttt{1.5__2.0}, \texttt{1.0__1.0}, 67 \texttt{3.1415926535897927__3.1415926535897936}) 68\item Other numeric types can be converted to bounded reals by giving them 69 a \texttt{breal/1} wrapper, or by calling 70 \bipref{breal/2}{../bips/kernel/arithmetic/breal-2.html} directly 71\item Bounded reals are not usually entered directly by the user; normally 72 they just occur as the results of computations 73\item A bounded real represents a single real number whose value is 74 known to lie somewhere between the bounds and is uncertain only 75 because of the limited precision with which is has been calculated 76\item An arithmetic operation is only performed using bounded reals if at 77 least one of its arguments is a bounded real 78\end{itemize} 79} 80 81An example of using bounded reals to safely compute the square root of 2: 82 83\begin{quote}\begin{verbatim} 84?- X is sqrt(breal(2)). 85X = 1.4142135623730949__1.4142135623730954 86Yes 87\end{verbatim}\end{quote} 88 89To see how using ordinary floating point numbers can lead to inaccuracy, try 90dividing 1 by 10, and then adding it together 10 times. Using floats the 91result is not 1.0; using bounded reals the computed interval contains 1.0 92and gives an indication of how much potential error there is: 93 94\begin{quote}\begin{verbatim} 95?- Y is float(1) / 10, X is Y + Y + Y + Y + Y + Y + Y + Y + Y + Y. 96X = 0.99999999999999989 97Y = 0.1 98Yes 99?- Y is breal(1) / 10, X is Y + Y + Y + Y + Y + Y + Y + Y + Y + Y. 100X = 0.99999999999999978__1.0000000000000007 101Y = 0.099999999999999992__0.1 102Yes 103\end{verbatim}\end{quote} 104 105 106\section{Issues to be aware of when using bounded reals} 107 108When working with bounded reals, some of the usual rules of arithmetic no 109longer hold. In particular, it is not always possible to determine whether 110one bounded real is larger, smaller, or the same as another. This is 111because, if the intervals overlap, it is not possible to know the 112relationship between the true values. 113 114{\enableunderscores 115An example of this can be seen in Figure~\ref{interval-compare}. If the 116true value of \texttt{X} is $\texttt{X}_1$, then depending upon whether the 117true value of \texttt{Y} is (say) $\texttt{Y}_1$, $\texttt{Y}_2$ or 118$\texttt{Y}_3$, we have \texttt{X > Y}, \texttt{X =:= Y} or \texttt{X < Y}, 119respectively. 120} 121 122\begin{figure} 123\begin{center} 124%\epsfbox{interval-compare.eps} 125\resizebox{0.13\textwidth}{!}{\includegraphics{interval-compare.eps}} 126\end{center} 127\caption{Comparing two bounded reals} 128\label{interval-compare} 129\end{figure} 130 131Different classes of predicate deal with the undecidable cases in different 132ways: 133 134\index{bounded reals!comparison} 135 136\begin{description} 137\item 138[Arithmetic comparison] (</2, =:=/2, etc.) 139If the comparison cannot be determined definitively, the comparison succeeds 140but a delayed goal is left behind, indicating that the result of the 141computation is contingent on the relationship actually being true. 142Examples: 143 144\begin{quote}\begin{verbatim} 145?- X = 0.2__0.3, Y = 0.0__0.1, X > Y. 146X = 0.2__0.3 147Y = 0.0__0.1 148Yes 149 150?- X = 0.2__0.3, Y = 0.0__0.1, X < Y. 151No 152 153?- X = 0.0__0.1, Y = 0.0__0.1, X < Y. 154X = 0.0__0.1 155Y = 0.0__0.1 156Delayed goals: 157 0.0__0.1 < 0.0__0.1 158Yes 159 160?- X = Y, X = 0.0__0.1, X < Y. 161No 162\end{verbatim}\end{quote} 163 164\item 165[Term equality or comparison] (=/2, ==/2, compare/3, @</2, etc.) 166These predicates consider bounded reals from a purely syntactic point of 167view: they determine how the bounded reals compare syntactically, without 168taking into account their meaning. Two bounded reals are considered equal 169if and only if their bounds are syntactically the same (note that the 170floating point numbers 0.0 and -0.0 are considered to be syntactically 171different). A unique ordering is also defined between bounded reals which 172do not have identical bounds; see the documentation for compare/3 for 173details. This is important as it means predicates such as sort/2 behave in 174a sensible fashion when they encounter bounded reals (in particular, they do 175not throw exceptions or leave behind large numbers of meaningless delayed 176goals) --- though one does need to be careful when comparing or sorting 177things of different types. 178Examples: 179 180\begin{quote}\begin{verbatim} 181?- X = 0.2__0.3, Y = 0.0__0.1, X == Y. 182No 183 184?- X = 0.0__0.1, Y = 0.0__0.1, X == Y. 185X = 0.0__0.1 186Y = 0.0__0.1 187Yes 188 189?- X = 0.2__0.3, Y = 0.0__0.1, compare(R, X, Y). 190R = > 191X = 0.2__0.3 192Y = 0.0__0.1 193Yes 194 195?- X = 0.1__3.0, Y = 0.2__0.3, compare(R, X, Y). 196R = < 197X = 0.1__3.0 198Y = 0.2__0.3 199Yes 200 201?- X = 0.0__0.1, Y = 0.0__0.1, compare(R, X, Y). 202R = = 203X = 0.0__0.1 204Y = 0.0__0.1 205Yes 206 207?- sort([-5.0, 1.0__1.0], Sorted). 208Sorted = [1.0__1.0, -5.0] % 1.0__1.0 > -5.0, but 1.0__1.0 @< -5.0 209Yes 210\end{verbatim}\end{quote} 211 212\end{description} 213 214Note that the potential undecidability of arithmetic comparisons has 215implications when writing general code. For example, a common thing to do 216is test the value of a number, with different code being executed depending 217on whether or not it is above a certain threshold; e.g.\ 218 219\begin{code} 220( X >= 0 -> 221 % Code A 222; 223 % Code B 224) 225\end{code} 226 227When writing code such as the above, if \texttt{X} could be a bounded real, 228one ought to decide what should happen if \texttt{X}'s bounds span the 229threshold value. In the above example, if \texttt{X = -0.1__0.1} then a 230delayed goal \texttt{-0.1__0.1 >= 0} will be left behind and Code A 231executed. If one does not want the delayed goal, one can instead write: 232 233\begin{code} 234( not X >= 0 -> 235 % Code B 236; 237 % Code A 238) 239\end{code} 240 241The use of \texttt{not} ensures that any actions performed during the test 242(in particular the set up of any delayed goals) are backtracked, regardless 243of the outcome of the test. 244 245Finally, if one wishes Code B to be executed instead of Code A in the case 246of an overlap, one can reverse the sense of the test: 247 248\begin{code} 249( not X < 0 -> 250 % Code A 251; 252 % Code B 253) 254\end{code} 255 256\index{bounded reals|)} 257 258 259\section{IC as a solver for real variables} 260 261\index{library!ic|(} 262 263The IC solver is a hybrid solver which supports both real and integer 264variables. 265 266\See{See Chapter~\ref{chapicintro} for an introduction to IC and how to use 267it with integer variables.} 268 269\See{See the IC chapter in the Constraint Library Manual for a full list of 270the arithmetic operators which are available for use in IC constraint 271expressions.} 272 273 274\quickref{Real variables and constraints}{ 275\begin{itemize} 276\item Real variables may be declared using 277 \biptxtrefni{reals/1}{reals/1!ic}{../bips/lib/ic/reals-1.html}, 278 \biptxtrefni{\$::/2}{\$::/2!ic}{../bips/lib/ic/SNN-2.html}, 279 \biptxtrefni{::/2}{::/2!ic}{../bips/lib/ic/NN-2.html} 280 \index{reals/1@\texttt{reals/1}!ic}\index{::/2@\texttt{::/2}!ic} 281 \index{\$::/2@\texttt{\$::/2}!ic}(specifying 282 non-integer bounds) or just by using them in an IC constraint 283\item Basic constraints available for real variables are 284 \biptxtrefni{\$=/2}{\$=/2!ic}{../bips/lib/ic/SE-2.html}, 285 \biptxtrefni{\$>=/2}{\$>=/2!ic}{../bips/lib/ic/SGE-2.html}, 286 \biptxtrefni{\$=</2}{\$=</2!ic}{../bips/lib/ic/SEL-2.html}, 287 \biptxtrefni{\$>/2}{\$>/2!ic}{../bips/lib/ic/SG-2.html}, 288 \biptxtrefni{\$</2}{\$</2!ic}{../bips/lib/ic/SL-2.html} 289 \index{\$=/2@\texttt{\$>=/2}!ic}\index{\$=/2@\texttt{\$>=/2}!ic} 290 \index{\$=</2@\texttt{\$=</2}!ic}\index{\$>/2@\texttt{\$>/2}!ic} 291 \index{\$</2@\texttt{\$</2}!ic}and 292 %\biptxtrefni{\$\bsl=/2}{\$\=/2@\$$\backslash$=/2!ic}{../bips/lib/ic/ERE-2.html}, 293 %\biptxtrefni{\$\bsl=/2}{\$\=/2@\$\verb+\+=/2!ic}{../bips/lib/ic/ERE-2.html}, 294 {\bf \$\bsl=/2} 295 %\html{\htmladdnormallink{\$\bsl=/2}{../bips/lib/ic/SRE-2.html}}, 296 %\protect\index{\$=@\$\verb+\+=/2!ic}% 297 %\protect\index{\$\=@=\bsl=/2!ic}% 298 as well as their reified versions and the reified connectives 299\item Real constraints also work with integer variables and a mix of integer 300 and real variables 301\item Solutions to real constraints can be found using 302 \bipref{locate/2}{../bips/lib/ic/locate-2.html}, 303 \bipref{locate/3}{../bips/lib/ic/locate-3.html}, 304 \bipref{locate/4}{../bips/lib/ic/locate-4.html} or 305 \bipref{squash/3}{../bips/lib/ic/squash-3.html} 306\end{itemize} 307} 308% This belongs inside the quickref, but the \bsl gets garbled in the index 309% if it's there. 310\index{\$\=@\$\bsl=/2!ic}% 311 312IC's real constraints perform bounds propagation in the same way as the 313integer versions; indeed, most of the basic integer constraints are 314transformed into their real counterparts, plus a declaration of the 315integrality of the variables appearing in the constraint. 316 317Note that the interval reasoning performed to propagate real bounds is the 318same as that used for bounded reals; that is, the inferences made are safe, 319taking into account potential floating point errors. 320 321 322\section{Finding solutions of real constraints} 323 324In very simple cases, just imposing the constraints may be sufficient to 325directly compute the (unique) solution. For example: 326 327\begin{quote}\begin{verbatim} 328?- 3 * X $= 4. 329X = 1.3333333333333333__1.3333333333333335 330Yes 331\end{verbatim}\end{quote} 332 333Other times, propagation will reduce the domains of the variables to 334suitably small intervals: 335 336\begin{quote}\begin{verbatim} 337?- 3 * X + 2 * Y $= 4, X - 5 * Y $= 2, X $>= -100. 338Y = Y{-0.11764705946382902 .. -0.1176470540212896} 339X = X{1.4117647026808551 .. 1.4117647063092196} 340There are 2 delayed goals. 341Yes 342\end{verbatim}\end{quote} 343 344In general though, some extra work will be needed to find the solutions of a 345problem. The IC library provides two methods for assisting with this. 346Which method is appropriate depends on the nature of the solutions to be 347found. If it is expected that there a finite number of discrete solutions, 348\bipref{locate/2}{../bips/lib/ic/locate-2.html} and 349\bipref{locate/3}{../bips/lib/ic/locate-3.html} would be good choices. If 350solutions are expected to lie in a continuous region, 351\bipref{squash/3}{../bips/lib/ic/squash-3.html} may be more appropriate. 352 353Locate works by nondeterministically splitting the domains of the variables 354until they are narrower than a specified precision (in either absolute or 355relative terms). Consider the problem of finding the points where two 356circles intersect (see Figure~\ref{locatefig}). Normal propagation does not 357deduce more than the obvious bounds on the variables: 358 359\begin{figure} 360\begin{center} 361%\epsfbox{locate2.eps} 362\resizebox{0.13\textwidth}{!}{\includegraphics{interval-compare.eps}} 363\end{center} 364\caption{Example of using locate/2} 365\label{locatefig} 366\end{figure} 367 368\begin{quote}\begin{verbatim} 369?- 4 $= X^2 + Y^2, 4 $= (X - 1)^2 + (Y - 1)^2. 370X = X{-1.0000000000000004 .. 2.0000000000000004} 371Y = Y{-1.0000000000000004 .. 2.0000000000000004} 372There are 12 delayed goals. 373Yes 374\end{verbatim}\end{quote} 375 376Calling \bipref{locate/2}{../bips/lib/ic/locate-2.html} quickly determines 377that there are two solutions and finds them to the desired accuracy: 378 379\begin{quote}\begin{verbatim} 380?- 4 $= X^2 + Y^2, 4 $= (X-1)^2 + (Y-1)^2, locate([X, Y], 1e-5). 381X = X{-0.8228756603552696 .. -0.82287564484820042} 382Y = Y{1.8228756448482002 .. 1.8228756603552694} 383There are 12 delayed goals. 384More 385 386X = X{1.8228756448482004 .. 1.8228756603552696} 387Y = Y{-0.82287566035526938 .. -0.82287564484820019} 388There are 12 delayed goals. 389Yes 390\end{verbatim}\end{quote} 391 392Squash works by deterministically cutting off parts of the domains of 393variables which it determines cannot contain any solutions. In effect, it 394is like a stronger version of bounds propagation. Consider the problem of 395finding the intersection of two circular discs and a hyperplane (see 396Figure~\ref{squashfig}). Again, normal propagation does not deduce more 397than the obvious bounds on the variables: 398 399\begin{figure} 400\begin{center} 401%\epsfbox{squash2.eps} 402\resizebox{0.47\textwidth}{!}{\includegraphics{squash2.eps}} 403\end{center} 404\caption{Example of propagation using the squash algorithm} 405\label{squashfig} 406\end{figure} 407 408\begin{quote}\begin{verbatim} 409?- 4 $>= X^2 + Y^2, 4 $>= (X-1)^2 + (Y-1)^2, Y $>= X. 410Y = Y{-1.0000000000000004 .. 2.0000000000000004} 411X = X{-1.0000000000000004 .. 2.0000000000000004} 412There are 13 delayed goals. 413Yes 414\end{verbatim}\end{quote} 415 416Calling \bipref{squash/3}{../bips/lib/ic/squash-3.html} results in the 417bounds being tightened (in this case the bounds are tight for the feasible 418region, though this is not true in general): 419 420\begin{quote}\begin{verbatim} 421?- 4 $>= X^2 + Y^2, 4 $>= (X-1)^2 + (Y-1)^2, Y $>= X, 422 squash([X, Y], 1e-5, lin). 423X = X{-1.0000000000000004 .. 1.4142135999632601} 424Y = Y{-0.41421359996326 .. 2.0000000000000004} 425There are 13 delayed goals. 426Yes 427\end{verbatim}\end{quote} 428 429\See{For more details, see the IC chapter of the Library Manual or the 430documentation for the individual predicates.} 431 432 433\section{A larger example} 434\label{farm-example} 435 436Consider the following problem: 437 438\begin{quote} 439George is contemplating buying a farm which is a very strange shape, 440comprising a large triangular lake with a square field on each side. The 441area of the lake is exactly seven acres, and the area of each field is an 442exact whole number of acres. Given that information, what is the smallest 443possible total area of the three fields? 444\end{quote} 445 446A diagram of the farm is shown in Figure~\ref{lake-fields}. 447 448\begin{figure} 449\begin{center} 450%\epsfbox{lake-fields.eps} 451\resizebox{0.37\textwidth}{!}{\includegraphics{lake-fields.eps}} 452\end{center} 453\caption{Triangular lake with adjoining square fields} 454\label{lake-fields} 455\end{figure} 456 457This is a problem which mixes both integer and real quantities, and as such 458is ideal for solving with the IC library. A model for the problem appears 459below. The \texttt{farm/4} predicate sets up the constraints between the 460total area of the farm \texttt{F} and the lengths of the three sides of the 461lake, \texttt{A}, \texttt{B} and \texttt{C}. 462 463\begin{code} 464:- lib(ic). 465 466farm(F, A, B, C) :- 467 [A, B, C] :: 0.0 .. 1.0Inf, % The 3 sides of the lake 468 triangle_area(A, B, C, 7), % The lake area is 7 469 470 [F, FA, FB, FC] :: 1 .. 1.0Inf, % The square areas are integral 471 square_area(A, FA), 472 square_area(B, FB), 473 square_area(C, FC), 474 F #= FA+FB+FC, 475 476 FA $>= FB, FB $>= FC. % Avoid symmetric solutions 477 478triangle_area(A, B, C, Area) :- 479 S $>= 0, 480 S $= (A+B+C)/2, 481 Area $= sqrt(S*(S-A)*(S-B)*(S-C)). 482 483square_area(A, Area) :- 484 Area $= sqr(A). 485\end{code} 486 487A solution to the problem can then be found by first instantiating the area 488of the farm, and then using \bipref{locate/2}{../bips/lib/ic/locate-2.html} 489to find the lengths of the sides of the lakes. Instantiating the area of 490the farm first ensures that the first solution returned will be the minimal 491one, since \bipref{indomain/1}{../bips/lib/ic/indomain-1.html} always 492chooses the smallest possible value first: 493 494\begin{code} 495solve(F) :- 496 farm(F, A, B, C), % the model 497 indomain(F), % ensure that solution is minimal 498 locate([A, B, C], 0.01). 499\end{code} 500 501\section{Exercise} 502 503\begin{enumerate} 504 505\item 506 507Consider the ``farm'' problem in section~\ref{farm-example}. (Source code 508may be found in \texttt{farm.ecl}, if you have access to it.) Try running 509this program to find the answer. Note that other, larger solutions are 510available by selecting \emph{more}. 511 512This implementation sums three integer variables (\texttt{FA}, \texttt{FB} 513and \texttt{FC}), and then constrains their order to remove symmetries. 514Would this be a good candidate for the global constraint 515\texttt{ordered_sum/2}? Modify the program so that it does use 516\texttt{ordered_sum/2}. How does the run time compare with the original? 517 518\end{enumerate} 519 520\index{library!ic|)} 521 522%HEVEA\cutend 523