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) 2000 - 2006 Cisco Systems, Inc. All Rights Reserved. 18% 19% Contributor(s): Warwick Harvey and Kish Shen, IC-Parc 20% 21% END LICENSE BLOCK 22% ---------------------------------------------------------------------- 23% System: ECLiPSe Constraint Logic Programming System 24% Version: $Id: eplex_relax.pl,v 1.1 2006/09/23 01:53:27 snovello Exp $ 25% 26% Linear relaxations of various constraints for use with eplex. 27% 28% (Currently just contains the convex hull relaxation of the piecewise 29% constraint.) 30% 31% W.Harvey, IC-Parc 32% 33% Modified by K.Shen, IC-Parc, to work with standalone eplex. Based on 34% ic_eplex_relax.ecl rather than the old eplex_relax.pl code. 35% 36% TODO: 37% 38% - Add support for infinite bounds in piecewise_linear_hull/3. Currently, 39% if any bound is infinite, no convex hull is generated. 40% 41% - When the eplex module provides support for subsuming constraints, 42% exploit it by subsuming constraints from old convex hulls when they are 43% not part of the new convex hull. 44% 45% ---------------------------------------------------------------------- 46 47% ---------------------------------------------------------------------- 48:- module(eplex_relax). 49% ---------------------------------------------------------------------- 50 51:- export piecewise_linear_hull/3. 52:- export piecewise_linear_hull/4. 53 54:- lib(eplex). 55 56% ---------------------------------------------------------------------- 57% The piecewise linear (convex hull) constraint 58% ---------------------------------------------------------------------- 59 60:- lib(ic). 61:- lib(ic_constraints). % For piecewise_linear_2/4 and struct(chunk/3). 62 63% 64% Computing the convex hull: 65% 66% The convex hull is found by computing the convex hull of all the points in 67% the piecewise linear specification which are feasible with respect to the 68% bounds on X and Y, as well as all the points where a current bound 69% intersects a segment from the piecewise linear constraint. 70% 71% The actual convex hull computation is done using an approach similar to 72% that used by the Graham scan algorithm for computing convex hulls. 73% In this implementation, we compute and maintain the convex hull in two 74% parts; the "upper" part and the "lower" part. These are computed 75% simultaneously, but are treated as separate entities. Essentially, 76% we process the points we wish to "wrap" in order from left to right, 77% maintaining the upper and lower convex hulls of the points processed so 78% far. When processing a new point, some of the points on the existing 79% convex hulls may fall inside the new ones (i.e. below the upper hull or 80% above the lower one). These are always the most recent points to have 81% been added, and are fairly easy to identify (and thus remove) by examining 82% gradients. 83% 84% When all points have been processed, one is left with a list of points 85% which are on the convex hull. It is then straightforward to generate 86% appropriate constraints defining this hull. 87% 88% 89% Updating the convex hull: 90% 91% To update the convex hull when a bound has changed, we simply recompute it 92% and send any new constraints to eplex. Note that it should be possible to 93% optimise this by only recomputing the parts of the convex hull which are 94% now infeasible (and thus changed). However, the implementation seems 95% complicated enough already. :-) 96% 97 98 99 % Structure for holding information about the convex hull, namely 100 % lists of the points on the upper and lower hulls and lists of the 101 % corresponding constraints. 102:- local struct(hull_info(upoints, lpoints, uconstraints, lconstraints)). 103 104 % Structure for holding information about a constraint. Currently 105 % this is not used except for diagnostics, but when the eplex module 106 % supports subsumption of constraints, this will store the 107 % constraint identifier required to subsume the constraint. 108:- local struct(constraint(dummy)). 109 110 % The predicate to be called by the user. It just calls other 111 % predicates to set up the propagation-only constraint, and then add 112 % the convex hull extension. 113 114piecewise_linear_hull(X, Points, Y) :- 115 piecewise_linear_hull(X, Points, Y, eplex). 116 117piecewise_linear_hull(X, Points, Y, Pool) :- 118 % Make sure X and Y have IC attributes before we try to suspend on 119 % them... 120 ic:reals([X, Y]), 121 % lower priority than propagate_hull/6 122 suspend(ic_to_eplex_bounds(X,Pool,S), 7, [X->ic:min,X->ic:max], S), 123 piecewise_linear_2(X, Points, Y, ChunkList), 124 piecewise_linear_hull_2(X, Y, ChunkList, Pool). 125 126 127 % Construct the initial hull info, set up the propagation demon, and 128 % perform the initial propagation. 129piecewise_linear_hull_2(X, Y, ChunkList, Pool) :- 130 Chunks =.. [[] | ChunkList], 131 HullInfo = hull_info with 132 [upoints:[], lpoints:[], 133 uconstraints:[], lconstraints:[]], 134 suspend(propagate_hull(X, Y, Chunks, HullInfo, Pool, Susp), 6, 135 [X->min, X->max, Y->min, Y->max], Susp), 136 propagate_hull(X, Y, Chunks, HullInfo, Pool, Susp). 137 138 139 % The demon which updates the convex hull whenever bounds change. 140:- demon(propagate_hull/6). 141propagate_hull(X, Y, Chunks, HullInfo, Pool, Susp) :- 142 ( var(X), var(Y) -> 143 get_bounds(X, XMin, XMax), 144 get_bounds(Y, YMin, YMax), 145 HullInfo = hull_info with 146 [upoints:UHullPoints, lpoints:LHullPoints, 147 uconstraints:UConstraints, lconstraints:LConstraints], 148 149 % Compute the new integer hull and update the constraints. 150 compute_hull_points(Chunks, XMin, XMax, YMin, YMax, 151 NewUHullPoints, NewLHullPoints), 152 %printf("Upper Hull Points:%n%DQw.%n", [NewUHullPoints]), 153 %printf("Lower Hull Points:%n%DQw.%n", [NewLHullPoints]), 154 update_hull_constraints(upper, Pool, X, Y, UHullPoints, NewUHullPoints, 155 UConstraints, NewUConstraints), 156 update_hull_constraints(lower, Pool, X, Y, LHullPoints, NewLHullPoints, 157 LConstraints, NewLConstraints), 158 159 % Update the description of the convex hull. 160 setarg(upoints of hull_info, HullInfo, NewUHullPoints), 161 setarg(lpoints of hull_info, HullInfo, NewLHullPoints), 162 setarg(uconstraints of hull_info, HullInfo, NewUConstraints), 163 setarg(lconstraints of hull_info, HullInfo, NewLConstraints) 164 ; 165 % If X or Y have become fixed, then we have nothing to add 166 % over what the piecewise_linear/2 constraint will infer. 167 kill_suspension(Susp) 168 ). 169 170 171 % Compute the points on the convex hull. 172 % Note: returns points in reverse order (right-to-left). 173compute_hull_points(Chunks, XMin, XMax, YMin, YMax, UHullPoints, LHullPoints) :- 174 functor(Chunks, _, NumChunks), 175 176 ( XMin > -1.0Inf, XMax < 1.0Inf -> 177 % Get the convex hull started. 178 compute_initial_hull_points(Chunks, XMin, XMax, YMin, YMax, 179 LeftChunkNum, NextPointNum, 180 UHullPoints0, LHullPoints0, Complete), 181 182 ( Complete = yes -> 183 UHullPoints = UHullPoints0, 184 LHullPoints = LHullPoints0 185 ; 186 fx(Chunks, LeftChunkNum, NumChunks, XMax, YRight, 187 right, RightChunkNum, RightSegNum), 188 189 LeftChunk is Chunks[LeftChunkNum], 190 LeftChunk = chunk with points:LeftPoints, 191 functor(LeftPoints, _, NumLeftPoints), 192 193 ( LeftChunkNum = RightChunkNum -> 194 % Left and right points are both in the same chunk. 195 % Wrap the remaining complete segments in this chunk. 196 wrap_points(NextPointNum, RightSegNum, LeftPoints, 197 YMin, YMax, between, 198 UHullPoints0, UHullPoints3, 199 LHullPoints0, LHullPoints3) 200 ; 201 % Wrap the rest of this chunk. 202 wrap_points(NextPointNum, NumLeftPoints, LeftPoints, 203 YMin, YMax, between, 204 UHullPoints0, UHullPoints1, 205 LHullPoints0, LHullPoints1), 206 207 % Wrap the remaining complete chunks. 208 NextChunkNum is LeftChunkNum + 1, 209 LastChunkNum is RightChunkNum - 1, 210 wrap_chunks(NextChunkNum, LastChunkNum, Chunks, 211 YMin, YMax, 212 UHullPoints1, UHullPoints2, 213 LHullPoints1, LHullPoints2), 214 215 % Wrap the complete segments of the final chunk. 216 RightChunk is Chunks[RightChunkNum], 217 RightChunk = chunk with points:RightPoints, 218 wrap_points(1, RightSegNum, RightPoints, 219 YMin, YMax, between, 220 UHullPoints2, UHullPoints3, 221 LHullPoints2, LHullPoints3) 222 ), 223 224 % Finish it off. 225 compute_final_hull_points(Chunks, XMax, YRight, YMin, YMax, 226 RightChunkNum, RightSegNum, 227 UHullPoints3, UHullPoints, 228 LHullPoints3, LHullPoints) 229 ) 230 ; 231 % We don't support infinite domains yet 232 UHullPoints = [], 233 LHullPoints = [] 234 ), 235 %write("Upper hull points: "), 236 %writeln(UHullPoints), 237 %write("Lower hull points: "), 238 %writeln(LHullPoints), 239 true. 240 241 242 % Compute the start of the convex hull (typically up to the first 243 % point to the right of XMin). 244compute_initial_hull_points(Chunks, XMin, XMax, YMin, YMax, 245 LeftChunkNum, NextPointNum, 246 UHullPoints, LHullPoints, Complete) :- 247 functor(Chunks, _, NumChunks), 248 fx(Chunks, 1, NumChunks, XMin, YLeft, left, LeftChunkNum, LeftSegNum), 249 Chunk is Chunks[LeftChunkNum], 250 Chunk = chunk with [points:Points, right:EndRight], 251 functor(Points, _, NumPoints), 252 253 ( NumPoints = 1 -> 254 % Discontinuity --- simply check whether it should be 255 % included or not. 256 P0 is Points[1], 257 (X0, Y0) = P0, 258 ( breal_max(Y0) >= YMin, breal_min(Y0) =< YMax -> 259 % XXX - for non-zero-width intervals, might be worth 260 % trimming Y bounds to YMin/YMax if appropriate. 261 UHullPoints = [P0], 262 LHullPoints = [P0] 263 ; 264 UHullPoints = [], 265 LHullPoints = [] 266 ), 267 NextPointNum = 2, 268 Complete = no 269 ; 270 P0 is Points[LeftSegNum], 271 P1 is Points[LeftSegNum+1], 272 (X0, Y0) = P0, 273 (X1, Y1) = P1, 274 275 ( ( XMax < breal_min(X1) ; EndRight = infinite, LeftSegNum =:= NumPoints-1 ) -> 276 % We're all in the one segment, which is pretty boring... 277 % Return some fake hull points that result in the correct 278 % constraints and say we're done. 279 % XXX - if the segment extends infinitely in either 280 % direction, the constraints added by these fake points 281 % could inadvertently tighten bounds due to roundoff. 282 UHullPoints = [P1, P0], 283 LHullPoints = [P1, P0], 284 NextPointNum is LeftSegNum + 2, 285 Complete = yes 286 ; XMin =:= breal_max(X1), LeftSegNum =:= NumPoints-1 -> 287 % We're at the very last point in the segment --- simply 288 % check whether it should be included or not. 289 ( breal_max(Y1) >= YMin, breal_min(Y1) =< YMax -> 290 % XXX - for non-zero-width intervals, might be worth 291 % trimming Y bounds to YMin/YMax if appropriate. 292 UHullPoints = [P1], 293 LHullPoints = [P1] 294 ; 295 UHullPoints = [], 296 LHullPoints = [] 297 ), 298 NextPointNum is LeftSegNum + 2, 299 Complete = no 300 ; 301 % General case 302 303 % Add XMin to the hull - being careful that the added point 304 % lies between the bounds (discontinuities can mean YLeft is 305 % actually way outside YMin .. YMax). 306 breal_bounds(YLeft, YLeftMin, YLeftMax), 307 ( YLeftMin < YMin -> 308 LHullPoints0 = [(XMin, YMin)] 309 ; YLeftMin > YMax -> 310 LHullPoints0 = [(XMin, YMax)] 311 ; 312 LHullPoints0 = [(XMin, YLeftMin)] 313 ), 314 ( YLeftMax > YMax -> 315 UHullPoints0 = [(XMin, YMax)] 316 ; YLeftMax < YMin -> 317 UHullPoints0 = [(XMin, YMin)] 318 ; 319 UHullPoints0 = [(XMin, YLeftMax)] 320 ), 321 322 % If (the upper bound of) this segment slopes upwards and 323 % intersects YMax, we need to add this intersection point to 324 % the upper hull. 325 ( YLeftMin < YMax, YMax < breal_max(Y1) -> 326 interpolate(Y0, X0, Y1, X1, YMax, Xa), 327 breal_min(Xa, XaMin), 328 329 ( XMin < XaMin, XaMin < breal_min(X1) -> 330 wrap_point_upper((XaMin, YMax), 331 UHullPoints0, UHullPoints) 332 ; 333 UHullPoints = UHullPoints0 334 ) 335 ; 336 UHullPoints = UHullPoints0 337 ), 338 339 % If (the lower bound of) this segment slopes downwards and 340 % intersects YMin, we need to add this intersection point to 341 % the lower hull. 342 ( YLeftMax > YMin, YMin > breal_min(Y1) -> 343 interpolate(Y0, X0, Y1, X1, YMin, Xb), 344 breal_min(Xb, XbMin), 345 346 ( XMin < XbMin, XbMin < breal_min(X1) -> 347 wrap_point_lower((XbMin, YMin), 348 LHullPoints0, LHullPoints) 349 ; 350 LHullPoints = LHullPoints0 351 ) 352 ; 353 LHullPoints = LHullPoints0 354 ), 355 356 NextPointNum is LeftSegNum + 1, 357 Complete = no 358 ) 359 ). 360 361 % Compute the end of the convex hull (from the point before the 362 % segment XMax is deemed to be in). 363compute_final_hull_points(Chunks, XMax, YRight, YMin, YMax, 364 RightChunkNum, RightSegNum, 365 UHullPoints0, UHullPoints, LHullPoints0, LHullPoints) :- 366 Chunk is Chunks[RightChunkNum], 367 Chunk = chunk with points:Points, 368 functor(Points, _, NumPoints), 369 breal_bounds(YRight, YRightMin0, YRightMax0), 370 YRightMin is min(YRightMin0, YMax), 371 YRightMax is max(YRightMax0, YMin), 372 373 ( NumPoints = 1 -> 374 % Discontinuity --- simply check whether it should be 375 % included or not. 376 P0 is Points[1], 377 (X0, Y0) = P0, 378 ( breal_max(Y0) >= YMin, breal_min(Y0) =< YMax -> 379 % XXX - for non-zero-width intervals, might be worth 380 % trimming Y bounds to YMin/YMax if appropriate. 381 wrap_point_upper(P0, UHullPoints0, UHullPoints), 382 wrap_point_lower(P0, LHullPoints0, LHullPoints) 383 ; 384 UHullPoints = UHullPoints0, 385 LHullPoints = LHullPoints0 386 ) 387 ; 388 P0 is Points[RightSegNum], 389 P1 is Points[RightSegNum+1], 390 (X0, Y0) = P0, 391 (X1, Y1) = P1, 392 393 % General case 394 395 % If (the lower bound of) this segment slopes upwards and 396 % intersects YMin, we need to add this intersection point to 397 % the lower hull. 398 ( breal_min(Y0) < YMin, YMin < breal_min(YRightMax) -> 399 interpolate(Y0, X0, Y1, X1, YMin, Xa), 400 breal_max(Xa, XaMax), 401 402 ( breal_max(X0) < XaMax, XaMax < XMax -> 403 wrap_point_lower((XaMax, YMin), 404 LHullPoints0, LHullPoints1) 405 ; 406 LHullPoints1 = LHullPoints0 407 ) 408 ; 409 LHullPoints1 = LHullPoints0 410 ), 411 412 % If (the upper bound of) this segment slopes downwards and 413 % intersects YMax, we need to add this intersection point to 414 % the upper hull. 415 ( breal_max(Y0) > YMax, YMax > breal_max(YRightMin) -> 416 interpolate(Y0, X0, Y1, X1, YMax, Xb), 417 breal_max(Xb, XbMax), 418 419 ( breal_max(X0) < XbMax, XbMax < XMax -> 420 wrap_point_upper((XbMax, YMax), 421 UHullPoints0, UHullPoints1) 422 ; 423 UHullPoints1 = UHullPoints0 424 ) 425 ; 426 UHullPoints1 = UHullPoints0 427 ), 428 429 % Add XMax to the hull 430 wrap_point_upper((XMax, YRightMax), UHullPoints1, UHullPoints), 431 wrap_point_lower((XMax, YRightMin), LHullPoints1, LHullPoints) 432 ). 433 434 435 % Wrap all the points in the chunks from ChunkNum to MaxChunk. 436wrap_chunks(ChunkNum, MaxChunk, Chunks, YMin, YMax, 437 UHullPoints0, UHullPoints, LHullPoints0, LHullPoints) :- 438 ( ChunkNum =< MaxChunk -> 439 Chunk is Chunks[ChunkNum], 440 Chunk = chunk with points:Points, 441 functor(Points, _, MaxPoint), 442 wrap_points(1, MaxPoint, Points, YMin, YMax, between, 443 UHullPoints0, UHullPoints1, 444 LHullPoints0, LHullPoints1), 445 ChunkNum1 is ChunkNum + 1, 446 wrap_chunks(ChunkNum1, MaxChunk, Chunks, YMin, YMax, 447 UHullPoints1, UHullPoints, 448 LHullPoints1, LHullPoints) 449 ; 450 UHullPoints = UHullPoints0, 451 LHullPoints = LHullPoints0 452 ). 453 454 % Wrap all the points from index I to index N. 455 % Prev indicates whether the previous point was `above', `below' or 456 % `between' the Y bounds (for detecting crossing of a bound). 457 % 458 % Note: any point at YMin (YMax) cannot be on the upper (lower) 459 % convex hull, unless it is also at XMin or XMax (in which case it 460 % is handled elsewhere). 461wrap_points(I, N, Points, YMin, YMax, Prev, UHullPoints0, UHullPoints, LHullPoints0, LHullPoints) :- 462 ( I =< N -> 463 Point1 is Points[I], 464 Point1 = (X1, Y1), 465 ( breal_max(Y1) > YMax -> 466 ( I > 1, Prev \== above -> 467 (X0, Y0) is Points[I - 1], 468 ( Prev = below -> 469 % Find intersection with YMin, wrap that point 470 interpolate(Y0, X0, Y1, X1, YMin, X2), 471 breal_max(X2, X2Max), 472 wrap_point_lower((X2Max, YMin), 473 LHullPoints0, LHullPoints1) 474 ; 475 LHullPoints1 = LHullPoints0 476 ), 477 % Find intersection with YMax, wrap that point 478 interpolate(Y0, X0, Y1, X1, YMax, X3), 479 breal_min(X3, X3Min), 480 ( X3Min > breal_min(X0) -> 481 wrap_point_upper((X3Min, YMax), 482 UHullPoints0, UHullPoints1) 483 ; 484 % Don't include the intersection if the fuzz puts it 485 % before the endpoint of this segment. 486 UHullPoints1 = UHullPoints0 487 ) 488 ; 489 % Do nothing 490 UHullPoints1 = UHullPoints0, 491 LHullPoints1 = LHullPoints0 492 ), 493 Curr = above 494 ; breal_min(Y1) < YMin -> 495 ( I > 1, Prev \== below -> 496 (X0, Y0) is Points[I - 1], 497 ( Prev = above -> 498 % Find intersection with YMax, wrap that point 499 interpolate(Y0, X0, Y1, X1, YMax, X2), 500 breal_max(X2, X2Max), 501 wrap_point_upper((X2Max, YMax), 502 UHullPoints0, UHullPoints1) 503 ; 504 UHullPoints1 = UHullPoints0 505 ), 506 % Find intersection with YMin, wrap that point 507 interpolate(Y0, X0, Y1, X1, YMin, X3), 508 breal_min(X3, X3Min), 509 ( X3Min > breal_min(X0) -> 510 wrap_point_lower((X3Min, YMin), 511 LHullPoints0, LHullPoints1) 512 ; 513 % Don't include the intersection if the fuzz puts it 514 % before the endpoint of this segment. 515 LHullPoints1 = LHullPoints0 516 ) 517 ; 518 % Do nothing 519 UHullPoints1 = UHullPoints0, 520 LHullPoints1 = LHullPoints0 521 ), 522 Curr = below 523 ; 524 % Point within current Y bounds 525 ( Prev = above -> 526 % Find intersection with YMax, wrap that point 527 (X0, Y0) is Points[I - 1], 528 interpolate(Y0, X0, Y1, X1, YMax, X2), 529 breal_max(X2, X2Max), 530 ( X2Max < breal_max(X1) -> 531 wrap_point_upper((X2Max, YMax), 532 UHullPoints0, UHullPoints0a) 533 ; 534 % Don't include the intersection if the fuzz puts it 535 % before the endpoint of this segment. 536 UHullPoints0a = UHullPoints0 537 ), 538 LHullPoints0a = LHullPoints0 539 ; Prev = below -> 540 % Find intersection with YMin, wrap that point 541 (X0, Y0) is Points[I - 1], 542 interpolate(Y0, X0, Y1, X1, YMin, X2), 543 breal_max(X2, X2Max), 544 ( X2Max < breal_max(X1) -> 545 wrap_point_lower((X2Max, YMin), 546 LHullPoints0, LHullPoints0a) 547 ; 548 % Don't include the intersection if the fuzz puts it 549 % before the endpoint of this segment. 550 LHullPoints0a = LHullPoints0 551 ), 552 UHullPoints0a = UHullPoints0 553 ; 554 UHullPoints0a = UHullPoints0, 555 LHullPoints0a = LHullPoints0 556 ), 557 558 wrap_point_upper(Point1, UHullPoints0a, UHullPoints1), 559 wrap_point_lower(Point1, LHullPoints0a, LHullPoints1), 560 561 Curr = between 562 ), 563 564 I1 is I + 1, 565 wrap_points(I1, N, Points, YMin, YMax, Curr, 566 UHullPoints1, UHullPoints, LHullPoints1, LHullPoints) 567 ; 568 UHullPoints = UHullPoints0, 569 LHullPoints = LHullPoints0 570 ). 571 572 573 % Add a point to the upper convex hull, dropping any points which 574 % would lie "inside" (under) the hull after the addition. 575wrap_point_upper(P, HullPoints0, HullPoints) :- 576 ( HullPoints0 = [P0 | HullPoints1] -> 577 P = (X, Y), 578 Yu is breal_max(Y), 579 P0 = (X0, Y0), 580 Y0u is breal_max(Y0), 581 ( HullPoints1 = [P1 | _Tail] -> 582 P1 = (X1, Y1), 583 Y1u is breal_max(Y1), 584 ( 585 % Figure out whether P0 should be retained. 586 ( Y0u > Y1u -> 587 ( Yu > Y0u -> 588 % Both segments going up. 589 % Use top left corners to compare gradients. 590 % XXX - assumes min(X1) < min(X0) < min(X) - 591 % must be updated for non-zero-width intervals. 592 % (Y0 - Y1) / (X0 - X1) > (Y - Y0) / (X - X0) 593 % (Y0 - Y1) * (X - X0) > (Y - Y0) * (X0 - X1) 594 % (Y0 - Y1) * (X - X0) - (Y - Y0) * (X0 - X1) > 0 595 Y0u1 is breal(Y0u), 596 X0l1 is breal(breal_min(X0)), 597 LHS is (Y0u1 - Y1u) * (breal_min(X) - X0l1) - 598 (Yu - Y0u1) * (X0l1 - breal_min(X1)), 599 % Assume spanning zero not good enough. 600 breal_min(LHS) > 0 601 ; 602 % First segment going up, second going down. 603 true 604 ) 605 ; 606 % If first segment going down, P0 can only be 607 % retained if second segment going down too. 608 Yu =< Y0u, 609 % Use top right corners to compare gradients. 610 % XXX - assumes max(X1) < max(X0) < max(X) - 611 % must be updated for non-zero-width intervals. 612 % (Y0 - Y1) / (X0 - X1) > (Y - Y0) / (X - X0) 613 % (Y0 - Y1) * (X - X0) > (Y - Y0) * (X0 - X1) 614 % (Y0 - Y1) * (X - X0) - (Y - Y0) * (X0 - X1) > 0 615 Y0u1 is breal(Y0u), 616 X0u1 is breal(breal_max(X0)), 617 LHS is (Y0u1 - Y1u) * (breal_max(X) - X0u1) - 618 (Yu - Y0u1) * (X0u1 - breal_max(X1)), 619 % Assume spanning zero not good enough. 620 breal_min(LHS) > 0 621 ) 622 -> 623 % New hull segment 624 HullPoints = [P | HullPoints0] 625 ; 626 % P0 is (at least approximately) inside the 627 % convex hull, so drop it and recurse. 628 wrap_point_upper(P, HullPoints1, HullPoints) 629 ) 630 ; 631 % P0---P is the initial segment: drop the lower of P0 and P 632 % if this segment is vertical (the corresponding constraint 633 % is implied by the X bound). 634 ( breal_min(X0) == breal_min(X) -> 635 ( Y0u < Yu -> 636 HullPoints = [P] 637 ; 638 HullPoints = HullPoints0 639 ) 640 ; 641 HullPoints = [P | HullPoints0] 642 ) 643 ) 644 ; 645 % There aren't any points on the hull yet, so just add 646 % this one. 647 HullPoints = [P] 648 ). 649 650 % Add a point to the lower convex hull, dropping any points which 651 % would lie "inside" (above) the hull after the addition. 652wrap_point_lower(P, HullPoints0, HullPoints) :- 653 ( HullPoints0 = [P0 | HullPoints1] -> 654 P = (X, Y), 655 Yl is breal_min(Y), 656 P0 = (X0, Y0), 657 Y0l is breal_min(Y0), 658 ( HullPoints1 = [P1 | _Tail] -> 659 P1 = (X1, Y1), 660 Y1l is breal_min(Y1), 661 ( 662 % Figure out whether P0 should be retained. 663 ( Y0l < Y1l -> 664 ( Yl < Y0l -> 665 % Both segments going down. 666 % Use bottom left corners to compare gradients. 667 % XXX - assumes min(X1) < min(X0) < min(X) - 668 % must be updated for non-zero-width intervals. 669 % (Y0 - Y1) / (X0 - X1) < (Y - Y0) / (X - X0) 670 % (Y0 - Y1) * (X - X0) < (Y - Y0) * (X0 - X1) 671 % (Y0 - Y1) * (X - X0) - (Y - Y0) * (X0 - X1) < 0 672 Y0l1 is breal(Y0l), 673 X0l1 is breal(breal_min(X0)), 674 LHS is (Y0l1 - Y1l) * (breal_min(X) - X0l1) - 675 (Yl - Y0l1) * (X0l1 - breal_min(X1)), 676 % Assume spanning zero not good enough. 677 breal_max(LHS) < 0 678 ; 679 % First segment going down, second going up. 680 true 681 ) 682 ; 683 % If first segment going up, P0 can only be 684 % retained if second segment going up too. 685 Yl >= Y0l, 686 % Use bottom right corners to compare gradients. 687 % XXX - assumes max(X1) < max(X0) < max(X) - 688 % must be updated for non-zero-width intervals. 689 % (Y0 - Y1) / (X0 - X1) < (Y - Y0) / (X - X0) 690 % (Y0 - Y1) * (X - X0) < (Y - Y0) * (X0 - X1) 691 % (Y0 - Y1) * (X - X0) - (Y - Y0) * (X0 - X1) < 0 692 Y0l1 is breal(Y0l), 693 X0u1 is breal(breal_max(X0)), 694 LHS is (Y0l1 - Y1l) * (breal_max(X) - X0u1) - 695 (Yl - Y0l1) * (X0u1 - breal_max(X1)), 696 % Assume spanning zero not good enough. 697 breal_max(LHS) < 0 698 ) 699 -> 700 % New hull segment 701 HullPoints = [P | HullPoints0] 702 ; 703 % P0 is (at least approximately) inside the 704 % convex hull, so drop it and recurse. 705 wrap_point_lower(P, HullPoints1, HullPoints) 706 ) 707 ; 708 % P0---P is the initial segment: drop the higher of P0 and P 709 % if this segment is vertical (the corresponding constraint 710 % is implied by the X bound). 711 ( breal_min(X0) == breal_min(X) -> 712 ( Y0l > Yl -> 713 HullPoints = [P] 714 ; 715 HullPoints = HullPoints0 716 ) 717 ; 718 HullPoints = [P | HullPoints0] 719 ) 720 ) 721 ; 722 % There aren't any points on the hull yet, so just add 723 % this one. 724 HullPoints = [P] 725 ). 726 727 728 % Update the hull constraints, based on the new points found for the 729 % convex hull. 730 % UL indicates whether the `upper' or `lower' hull is being updated, 731 % so that the hull-defining constraints can be given the correct 732 % sense. 733update_hull_constraints(UL, Pool, X, Y, OldPoints, NewPoints, 734 OldConstraints, NewConstraints) :- 735 % Find the first point in common between the old and the new hulls, 736 % if it exists. 737 ( find_first_common_point(OldPoints, NewPoints, OldN, NewN, 738 OldPoints1, NewPoints1) -> 739 % Find the last point in common. 740 find_last_common_point(OldPoints1, NewPoints1, N, 741 _OldPoints2, NewHullPoints2), 742 743 % Split the old constraints into those which have been 744 % subsumed, and those which are in common with the new hull. 745 take_drop_n(OldN, OldConstraints, SubsumedConstraints1, 746 OldConstraints1), 747 take_drop_n(N, OldConstraints1, CommonConstraints, 748 SubsumedConstraints2), 749 750 % Extract the points corresonding to the new constraints 751 % which need to be added at the start of the hull. 752 NewN1 is NewN + 1, 753 take_n(NewN1, NewPoints, NewHullPoints1), 754 755 % Subsume the constraints no longer needed. 756 subsumed_constraints(SubsumedConstraints1), 757 subsumed_constraints(SubsumedConstraints2), 758 759 % Add the new constraints required. 760 hull_points_to_constraints(UL, NewHullPoints1, X, Y, Pool, 761 NewConstraints1), 762 hull_points_to_constraints(UL, NewHullPoints2, X, Y, Pool, 763 NewConstraints2), 764 765 % Construct the new list of constraints. 766 append(CommonConstraints, NewConstraints2, NewConstraints3), 767 append(NewConstraints1, NewConstraints3, NewConstraints) 768 ; 769 % No common points, so replace the hull wholesale. 770 subsumed_constraints(OldConstraints), 771 hull_points_to_constraints(UL, NewPoints, X, Y, Pool, 772 NewConstraints) 773 ). 774 775 % Returns the number of mismatching points in each list, as well the 776 % tails of these lists with the mismatching points stripped off. 777find_first_common_point(OldPoints, NewPoints, OldN, NewN, OldTail, NewTail) :- 778 OldPoints = [OldPoint | OldPoints1], 779 NewPoints = [NewPoint | NewPoints1], 780 % XXX - assumes X coordinates are always bounded reals, and 781 % both the lower and upper bounds are non-decreasing (as enforced by 782 % unpack_point/3 and check_sorted_and_trim/2, respectively). 783 compare(R, OldPoint, NewPoint), 784 % Remember: points are in reverse X order. 785 ( R = = -> 786 OldN = 0, 787 NewN = 0, 788 OldTail = OldPoints, 789 NewTail = NewPoints 790 ; R = < -> 791 find_first_common_point(OldPoints, NewPoints1, 792 OldN, NewN0, OldTail, NewTail), 793 NewN is NewN0 + 1 794 ; 795 find_first_common_point(OldPoints1, NewPoints, 796 OldN0, NewN, OldTail, NewTail), 797 OldN is OldN0 + 1 798 ). 799 800 % Strips off all matching points but the last from each list, 801 % returning the number of points removed. 802find_last_common_point(OldPoints, NewPoints, N, OldTail, NewTail) :- 803 OldPoints = [OldPoint | OldPoints1], 804 NewPoints = [NewPoint | NewPoints1], 805 find_last_common_point_2(OldPoint, OldPoints1, NewPoint, NewPoints1, N, 806 OldTail, NewTail). 807 808find_last_common_point_2(PrevOldPoint, OldPoints, PrevNewPoint, NewPoints, N, 809 OldTail, NewTail) :- 810 ( 811 OldPoints = [OldPoint | OldPoints1], 812 NewPoints = [NewPoint | NewPoints1], 813 OldPoint = NewPoint 814 -> 815 find_last_common_point_2(OldPoint, OldPoints1, 816 NewPoint, NewPoints1, N0, OldTail, NewTail), 817 N is N0 + 1 818 ; 819 N = 0, 820 OldTail = [PrevOldPoint | OldPoints], 821 NewTail = [PrevNewPoint | NewPoints] 822 ). 823 824 825 % Converts a list of hull points to constraints on X and Y. 826hull_points_to_constraints(upper, Points, X, Y, Pool, Constraints) :- 827 upper_points_to_constraints(Points, X, Y, Pool, Constraints). 828hull_points_to_constraints(lower, Points, X, Y, Pool, Constraints) :- 829 lower_points_to_constraints(Points, X, Y, Pool, Constraints). 830 831upper_points_to_constraints([], _, _, _, []). % For infinite hulls 832upper_points_to_constraints([Point | Points], X, Y, Pool, Constraints) :- 833 upper_points_to_constraints_2(Point, Points, X, Y, Pool, Constraints). 834 835upper_points_to_constraints_2(_, [], _, _, _, []). 836upper_points_to_constraints_2(PrevP, [P | Ps], X, Y, Pool, [C | Cs]) :- 837 % Remember: points are in reverse order. 838 upper_points_to_constraint(P, PrevP, X, Y, Pool, C), 839 upper_points_to_constraints_2(P, Ps, X, Y, Pool, Cs). 840 841 % Assumes X0 =< X1 842upper_points_to_constraint((X0, Y0), (X1, Y1), X, Y, Pool, C) :- 843 % Y =< Y0 + (Y1 - Y0) / (X1 - X0) * (X - X0) 844 % (X1 - X0) * Y - (Y1 - Y0) * X =< Y0 * (X1 - X0) - X0 * (Y1 - Y0) 845 Y0Max is breal_max(Y0), 846 Y1Max is breal_max(Y1), 847 XCoef0 is breal(Y1Max) - breal(Y0Max), 848 breal_bounds(XCoef0, XCoefMin, XCoefMax), 849 ( XCoefMin > 0 -> 850 % Segment going up. Use top left points. 851 XCoef = XCoef0, 852 X0a is breal_min(X0), 853 X1a is breal_min(X1) 854 ; XCoefMax < 0 -> 855 % Segment going down. Use top right points. 856 XCoef = XCoef0, 857 X0a is breal_max(X0), 858 X1a is breal_max(X1) 859 ; 860 % Segment is horizontal (or close enough). 861 % Constraint should be implied by variable bounds, so just 862 % return something which yields a tautology. 863 XCoef = 0.0__0.0, 864 X0a = 0.0, 865 X1a = 0.0 866 ), 867 YCoef is breal(X1a) - breal(X0a), 868 RHSMax is breal_max(Y0Max * YCoef - X0a * XCoef), 869 870 % Which bound of the X and Y coefficients to use (for safety) 871 % depends on the sign of the X and Y values its going to be used 872 % with... Probably we don't have to be too paranoid about this 873 % because: 874 % a) the coefficient is almost certainly a very narrow interval 875 % (it's a simple difference between two floats), 876 % b) we're passing them to an external solver which is going to 877 % use approximating arithmetic a lot anyway, and presumably 878 % includes its own safeguards, 879 % c) we get some fuzz protection from the RHS bound. 880 % So for segments which span a zero coordinate (which would 881 % otherwise require special processing) we just choose the bound 882 % which is appropriate for the majority of the segment... 883 ( X0a + X1a >= 0.0 -> 884 XCoef1 is breal_max(XCoef) 885 ; 886 XCoef1 is breal_min(XCoef) 887 ), 888 ( Y0Max + Y1Max >= 0.0 -> 889 YCoef1 is breal_min(YCoef) 890 ; 891 YCoef1 is breal_max(YCoef) 892 ), 893 Rep = (YCoef1 * Y - XCoef1 * X =< RHSMax), 894 895 C = constraint with [dummy:Rep], 896 %printf("Adding constraint: %DQw%n", [Rep]), 897 add_pool_constraint(Rep, Pool). 898 899lower_points_to_constraints([], _, _, _, []). % For infinite hulls 900lower_points_to_constraints([Point | Points], X, Y, Pool, Constraints) :- 901 lower_points_to_constraints_2(Point, Points, X, Y, Pool, Constraints). 902 903lower_points_to_constraints_2(_, [], _, _, _, []). 904lower_points_to_constraints_2(PrevP, [P | Ps], X, Y, Pool, [C | Cs]) :- 905 % Remember: points are in reverse order. 906 lower_points_to_constraint(P, PrevP, X, Y, Pool, C), 907 lower_points_to_constraints_2(P, Ps, X, Y, Pool, Cs). 908 909 % Assumes X0 =< X1 910lower_points_to_constraint((X0, Y0), (X1, Y1), X, Y, Pool, C) :- 911 % Y >= Y0 + (Y1 - Y0) / (X1 - X0) * (X - X0) 912 % (X1 - X0) * Y - (Y1 - Y0) * X >= Y0 * (X1 - X0) - X0 * (Y1 - Y0) 913 Y0Min is breal_min(Y0), 914 Y1Min is breal_min(Y1), 915 XCoef0 is breal(Y1Min) - breal(Y0Min), 916 breal_bounds(XCoef0, XCoefMin, XCoefMax), 917 ( XCoefMax < 0 -> 918 % Segment going down. Use bottom left points. 919 XCoef = XCoef0, 920 X0a is breal_min(X0), 921 X1a is breal_min(X1) 922 ; XCoefMin > 0 -> 923 % Segment going up. Use bottom right points. 924 XCoef = XCoef0, 925 X0a is breal_max(X0), 926 X1a is breal_max(X1) 927 ; 928 % Segment is horizontal (or close enough). 929 % Constraint should be implied by variable bounds, so just 930 % return something which yields a tautology. 931 XCoef = 0.0__0.0, 932 X0a = 0.0, 933 X1a = 0.0 934 ), 935 YCoef is breal(X1a) - breal(X0a), 936 RHSMin is breal_min(Y0Min * YCoef - X0a * XCoef), 937 938 % Which bound of the X and Y coefficients to use (for safety) 939 % depends on the sign of the X and Y values its going to be used 940 % with... Probably we don't have to be too paranoid about this 941 % because: 942 % a) the coefficient is almost certainly a very narrow interval 943 % (it's a simple difference between two floats), 944 % b) we're passing them to an external solver which is going to 945 % use approximating arithmetic a lot anyway, and presumably 946 % includes its own safeguards, 947 % c) we get some fuzz protection from the RHS bound. 948 % So for segments which span a zero coordinate (which would 949 % otherwise require special processing) we just choose the bound 950 % which is appropriate for the majority of the segment... 951 ( X0a + X1a >= 0.0 -> 952 XCoef1 is breal_min(XCoef) 953 ; 954 XCoef1 is breal_max(XCoef) 955 ), 956 ( Y0Min + Y1Min >= 0.0 -> 957 YCoef1 is breal_max(YCoef) 958 ; 959 YCoef1 is breal_min(YCoef) 960 ), 961 Rep = (YCoef1 * Y - XCoef1 * X >= RHSMin), 962 963 C = constraint with [dummy:Rep], 964 %printf("Adding constraint: %DQw%n", [Rep]), 965 add_pool_constraint(Rep, Pool). 966 967 968 % Subsume the given list of constraints. 969 % XXX - Since eplex doesn't support this yet and it's safe to leave 970 % them active, for now we just do nothing. 971subsumed_constraints(_SubsumedConstraints) :- 972/* 973 ( 974 foreach(Constraint, SubsumedConstraints) 975 do 976 Constraint = constraint with dummy:Rep, 977 write("Subsuming constraint "), 978 writeln(Rep) 979 ). 980*/ 981 true. 982 983 984 % Evaluate the piecewise linear function at the given point. 985 % ChunkLo and ChunkHi give upper and lower bounds on the chunk 986 % number in which to find X. ChunkNum and SegNum return the chunk 987 % and segment in which X was found, and Y is the function result. 988fx(Chunks, ChunkLo, ChunkHi, X, Y, LeftRight, ChunkNum, SegNum) :- 989 ( ChunkLo = ChunkHi -> 990 Chunk is Chunks[ChunkLo], 991 Chunk = chunk with points:Points, 992 functor(Points, _, N), 993 N_1 is N - 1, % Number of _segments_ 994 fx(Points, 1, N_1, X, Y, LeftSegNum, RightSegNum), 995 ( LeftRight = left -> 996 SegNum = LeftSegNum 997 ; 998 SegNum = RightSegNum 999 ), 1000 ChunkNum = ChunkLo 1001 ; 1002 % Do a binary chop. 1003 ChunkMid is (ChunkLo + ChunkHi + 1) // 2, 1004 Chunk is Chunks[ChunkMid], 1005 Chunk = chunk with [points:Points, left:Left], 1006 (Xp, _) is Points[1], 1007 select_bound(LeftRight, Xp, Xp0), 1008 ( go_left(Left, X, Xp0) -> 1009 NewChunkHi is ChunkMid - 1, 1010 fx(Chunks, ChunkLo, NewChunkHi, X, Y, LeftRight, ChunkNum, SegNum) 1011 ; 1012 fx(Chunks, ChunkMid, ChunkHi, X, Y, LeftRight, ChunkNum, SegNum) 1013 ) 1014 ). 1015 1016 % Select the bound of Xp to use as the cut-off point, depending on 1017 % whether we're looking for the leftmost X or rightmost X. 1018select_bound(left, Xp, Xp0) :- 1019 Xp0 is breal_max(Xp). 1020select_bound(right, Xp, Xp0) :- 1021 Xp0 is breal_min(Xp). 1022 1023 % Decide whether to go left or right in the binary chop, depending on 1024 % whether the discontinuity at the left end of the chunk is open or 1025 % closed. 1026go_left(open, X, Xp0) :- 1027 X =< Xp0. 1028go_left(closed, X, Xp0) :- 1029 X < Xp0. 1030 1031 1032take_n(N, List, Head) :- 1033 take_drop_n(N, List, Head, _). 1034 1035take_drop_n(N, List, Head, Tail) :- 1036 ( 1037 count(I, 1, N), 1038 fromto(List, [H | ListT], ListT, Tail), 1039 fromto(Head, [H | TakenT], TakenT, []) 1040 do 1041 true 1042 ). 1043 1044% ---------------------------------------------------------------------- 1045% Bounds transfer 1046% ---------------------------------------------------------------------- 1047 1048:- demon ic_to_eplex_bounds/3. 1049ic_to_eplex_bounds(V, Pool, S) :- 1050 var(V), 1051 ic: get_bounds(V, Min, Max), 1052 Pool: (V:: Min..Max). 1053ic_to_eplex_bounds(V, _, S) :- 1054 nonvar(V), 1055 kill_suspension(S). 1056