1% 2% CLP(BNR) == Constraints On Boolean, Integer, and Real Intervals 3% 4/* The MIT License (MIT) 5 * 6 * Copyright (c) 2019-2025 Rick Workman 7 * 8 * Permission is hereby granted, free of charge, to any person obtaining a copy 9 * of this software and associated documentation files (the "Software"), to deal 10 * in the Software without restriction, including without limitation the rights 11 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 12 * copies of the Software, and to permit persons to whom the Software is 13 * furnished to do so, subject to the following conditions: 14 * 15 * The above copyright notice and this permission notice shall be included in all 16 * copies or substantial portions of the Software. 17 * 18 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 19 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 20 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 21 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 22 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 23 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 24 * SOFTWARE. 25 */ 26 27:- module(clpBNR, % SWI module declaration 28 [ 29 op(700, xfx, ::), 30 (::)/2, % declare interval 31 {}/1, % define constraint 32 add_constraint/1, % more primitive define single constraint, bypass simplify 33 interval/1, % filter for clpBNR constrained var 34 interval_degree/2, % number of constraints on clpBNR constrained var 35 list/1, % O(1) list filter (also for compatibility) 36 domain/2, range/2, % get type and bounds (domain) 37 delta/2, % width (span) of an interval or numeric (also arithmetic function) 38 midpoint/2, % midpoint of an interval (or numeric) (also arithmetic function) 39 median/2, % median of an interval (or numeric) (also arithmetic function) 40 lower_bound/1, % narrow interval to point equal to lower bound 41 upper_bound/1, % narrow interval to point equal to upper bound 42 43 % additional constraint operators 44 op(200, fy, ~), % boolean 'not' 45 op(500, yfx, and), % boolean 'and' 46 op(500, yfx, or), % boolean 'or' 47 op(500, yfx, nand), % boolean 'nand' 48 op(500, yfx, nor), % boolean 'nor' 49 % op(500, yfx, xor), % boolean 'xor', but operator already defined (P=400) for arithmetic 50 op(700, xfx, <>), % integer not equal 51 op(700, xfx, <=), % included (one way narrowing) 52 53 % utilities 54 print_interval/1, print_interval/2, % pretty print interval with optional stream 55 small/1, small/2, % defines small interval width based on precision value 56 solve/1, solve/2, % solve (list of) intervals using split to find point solutions 57 splitsolve/1, splitsolve/2, % solve (list of) intervals using split 58 absolve/1, absolve/2, % absolve (list of) intervals, narrows by nibbling bounds 59 enumerate/1, % "enumerate" integers 60 global_minimum/2, % find interval containing global minimum(s) for an expression 61 global_minimum/3, % global_minimum/2 with definable precision 62 global_maximum/2, % find interval containing global maximum(s) for an expression 63 global_maximum/3, % global_maximum/2 with definable precision 64 global_minimize/2, % global_minimum/2 plus narrow vars to found minimizers 65 global_minimize/3, % global_minimum/3 plus narrow vars to found minimizers 66 global_maximize/2, % global_maximum/2 plus narrow vars to found maximizers 67 global_maximize/3, % global_maximum/3 plus narrow vars to found maximizers 68 nb_setbounds/2, % non-backtracking set bounds (use with branch and bound) 69 partial_derivative/3, % differentiate Exp wrt. X and simplify 70 clpStatistics/0, % reset 71 clpStatistic/1, % get selected 72 clpStatistics/1, % get all defined in a list 73 watch/2, % enable monitoring of changes for interval or (nested) list of intervals 74 trace_clpBNR/1 % enable/disable tracing of clpBNR narrowing ops 75 ]). 76 77/* missing(?) functionality from original CLP(BNR): utility accumulate/2. */ 78 79/* supported interval relations: 80 81+ - * / %% arithmetic 82** %% includes real exponent, odd/even integer 83abs %% absolute value 84sqrt %% positive square root 85min max %% binary min/max 86== is <> =< >= < > %% comparison (`is` synonym for `==`) 87<= %% included (one way narrowing) 88and or nand nor xor -> , %% boolean (`,` synonym for `and`) 89- ~ %% unary negate and not 90exp log %% exp/ln 91sin asin cos acos tan atan %% trig functions 92integer %% must be an integral value, e.g., 1 and 1.0 are both integral 93sig %% signum of real, (-1,0,+1) 94 95*/
116version("0.12.1"). 117 118% support various optimizations via goal expansion 119:- discontiguous clpBNR:goal_expansion/2. 120 121% debug feature control and messaging 122:- if(exists_source(swish(lib/swish_debug))). 123 :- create_prolog_flag(clpBNR_swish, true, [access(read_only)]). 124 :- use_module(swish(lib/swish_debug)). 125 :- use_module(library(http/html_write)). 126:- else. 127 :- use_module(library(debug)). 128:- endif. 129 130:- use_module(library(prolog_versions)). % SWIP dependency enforcement 131 132:- require_prolog_version('9.1.22', % properly exported arithmetic functions 133 [ rational % require rational number support, implies bounded=false 134 ]). 135 136% 137% Optimize arithmetic, but not debug. Only called via directives. 138% 139set_optimize_flags_ :- % called at start of load 140 set_prolog_flag(optimise,true), % scoped to file/module 141 current_prolog_flag(optimise_debug,ODflag), % save; restore in initialization 142 nb_linkval('$optimise_debug_save',ODflag), 143 set_prolog_flag(optimise_debug,false). % so debug/3, debugging/1 don't get "optimized" 144 145restore_optimize_flags_ :- % called at module initialization (after load) 146 nb_getval('$optimise_debug_save',ODflag), nb_delete('$optimise_debug_save'), 147 set_prolog_flag(optimise_debug,ODflag). 148 149:- set_optimize_flags_. 150 151% local debug and trace message support 152debug_clpBNR_(FString,Args) :- debug(clpBNR,FString,Args). 153 154% sandboxing for SWISH 155:- multifile(sandbox:safe_prolog_flag/1). 156:- multifile(sandbox:safe_global_variable/1). 157:- multifile(sandbox:safe_primitive/1). 158:- multifile(sandbox:safe_meta/2). 159 160current_node_(Node) :- % look back to find current Op being executed for debug messages 161 prolog_current_frame(Frame), % this is a little grungy, but necessary to get intervals 162 prolog_frame_attribute(Frame,parent_goal,doNode_(Args,Op,_,_,_,_,_)), 163 map_constraint_op_(Op,Args,Node), 164 !. 165 166sandbox:safe_primitive(clpBNR:current_node_(_Node)). 167 168% 169% statistics 170% 171 172% assign,increment/read global counter (assumed to be ground value so use _linkval) 173g_assign(G,V) :- nb_linkval(G,V). 174g_inc(G) :- nb_getval(G,N), N1 is N+1, nb_linkval(G,N1). 175g_incb(G) :- nb_getval(G,N), N1 is N+1, b_setval(G,N1). % undone on backtrack 176g_read(G,V) :- nb_getval(G,V). 177 178sandbox:safe_global_variable('$clpBNR:thread_init_done'). 179sandbox:safe_global_variable('$clpBNR:userTime'). 180sandbox:safe_global_variable('$clpBNR:inferences'). 181sandbox:safe_global_variable('$clpBNR:gc_time'). 182 183% 184% Global var statistics are per thread and therefore must be "lazily" initialized 185% Also ensures that thread copies of flags are properly set. 186% This exception handler will be invoked the first time '$clpBNR:thread_init_done' is read 187% Public predicates ::, {}, clpStatistics/1 and range read this global before proceeding. 188% 189user:exception(undefined_global_variable,'$clpBNR:thread_init_done',retry) :- !, 190 set_prolog_flags, % initialize/check environment flags 191 clpStatistics, % defines remaining globals 192 g_assign('$clpBNR:thread_init_done',1). 193 194% 195% Define custom clpBNR flags when module loaded 196% 197:- create_prolog_flag(clpBNR_iteration_limit,3000,[type(integer),keep(true)]). 198:- create_prolog_flag(clpBNR_default_precision,6,[type(integer),keep(true)]). 199:- create_prolog_flag(clpBNR_verbose,false,[type(boolean),keep(true)]). 200 201sandbox:safe_prolog_flag(clpBNR_iteration_limit,_). 202sandbox:safe_prolog_flag(clpBNR_default_precision,_). 203sandbox:safe_prolog_flag(clpBNR_verbose,_). 204% 205% Set public flags (see module/thread initialization) 206% 207set_prolog_flags :- 208 set_prolog_flag(prefer_rationals, true), % enable rational arithmetic 209 (current_prolog_flag(max_rational_size,_) 210 -> true % already defined, leave as is 211 ; set_prolog_flag(max_rational_size, 16) % rational size in bytes before .. 212 ), 213 set_prolog_flag(max_rational_size_action, float), % conversion to float 214 215 set_prolog_flag(float_overflow,infinity), % enable IEEE continuation values 216 set_prolog_flag(float_zero_div,infinity), 217 set_prolog_flag(float_undefined,nan), 218 set_prolog_flag(write_attributes,portray). % thread-local, init to 'portray'
clpBNR statistics - always succeeds.
clpBNR collects a number of "operational measurements" on a per-thread basis and combines them with some system statistics for subsequent querying. clpBNR measurements include:
narrowingOps | number of interval primitives called |
narrowingFails | number of interval primitive failures |
node_count | number of nodes in clpBNR constraint network |
max_iterations | maximum number of iterations before throttling occurs (max/limit |
System statistics included in clpStatistics:
userTime | from statistics:cputime |
gcTime | from statistics:garbage_collection.Time |
globalStack | from statistics:globalused/statistics:global |
trailStack | from statistics:trailused/statistics:trail |
localStack | from statistics:localused/statistics:local |
inferences | from statistics:inferences |
/
clpStatistic value; otherwise fails. On backtracking all values that unify with S will be generated. Examples:
?- clpStatistics, X::real, {X**4-4*X**3+4*X**2-4*X+3==0}, clpStatistic(narrowingOps(Ops)).
Ops = 2245,
X::real(-1.509169756145379, 4.18727500493995).
?- clpStatistics, X::real, {X**4-4*X**3+4*X**2-4*X+3==0}, clpStatistic(S).
S = userTime(0.02277600000000035),
X::real(-1.509169756145379, 4.18727500493995) ;
S = gcTime(0.0),
X::real(-1.509169756145379, 4.18727500493995) ;
S = globalStack(43696/524256),
X::real(-1.509169756145379, 4.18727500493995) ;
S = trailStack(664/133096),
X::real(-1.509169756145379, 4.18727500493995) ;
S = localStack(1864/118648),
X::real(-1.509169756145379, 4.18727500493995) ;
S = inferences(86215),
X::real(-1.509169756145379, 4.18727500493995) ;
S = narrowingOps(2245),
X::real(-1.509169756145379, 4.18727500493995) ;
S = narrowingFails(0),
X::real(-1.509169756145379, 4.18727500493995) ;
S = node_count(9),
X::real(-1.509169756145379, 4.18727500493995) ;
S = max_iterations(2245/3000),
X::real(-1.509169756145379, 4.18727500493995).
*/
clpStatistic's values; otherwise fails. Example:
?- clpStatistics, X::real, {X**4-4*X**3+4*X**2-4*X+3==0}, clpStatistics(Ss).
Ss = [userTime(0.023398999999999504), gcTime(0.001), globalStack(19216/131040), trailStack(1296/133096), localStack(2184/118648), inferences(82961), narrowingOps(2245), narrowingFails(0), node_count(9), max_iterations(2245/3000)],
X::real(-1.509169756145379, 4.18727500493995).
/
284:- discontiguous clpBNR:clpStatistics/0, clpBNR:clpStatistic/1. 285 286clpStatistics :- 287 % garbage_collect, % ? do gc before time snapshots 288 statistics(cputime,T), g_assign('$clpBNR:userTime',T), % thread based 289 statistics(inferences,I), g_assign('$clpBNR:inferences',I), 290 statistics(garbage_collection,[_,_,G,_]), g_assign('$clpBNR:gc_time',G), 291 fail. % backtrack to reset other statistics. 292 293clpStatistic(_) :- g_read('$clpBNR:thread_init_done',0). % ensures per-thread initialization then fails 294 295clpStatistic(userTime(T)) :- statistics(cputime,T1), g_read('$clpBNR:userTime',T0), T is T1-T0. 296 297clpStatistic(gcTime(G)) :- statistics(garbage_collection,[_,_,G1,_]), g_read('$clpBNR:gc_time',G0), G is (G1-G0)/1000.0. 298 299clpStatistic(globalStack(U/T)) :- statistics(globalused,U), statistics(global,T). 300 301clpStatistic(trailStack(U/T)) :- statistics(trailused,U), statistics(trail,T). 302 303clpStatistic(localStack(U/T)) :- statistics(localused,U), statistics(local,T). 304 305clpStatistic(inferences(I)) :- statistics(inferences,I1), g_read('$clpBNR:inferences',I0), I is I1-I0.
313list(X) :- compound(X) -> X=[_|_] ; X==[]. 314 315:- include(clpBNR/ia_primitives). % interval arithmetic relations via evalNode/4. 316 317%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 318% 319% SWI-Prolog implementation of IA 320% 321%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 322% 323% Intervals are constrained (attributed) variables. 324% 325% Current interval bounds updates via setarg(Val) which is backtrackable
clpBNR attribute; otherwise fails.
/
331interval(Int) :- get_attr(Int, clpBNR, _).clpBNR constraints on X; otherwise fails. If X is a number, N = 0. Examples:
?- {X==Y+1}, interval_degree(X,N).
N = 1,
X::real(-1.0Inf, 1.0Inf),
Y::real(-1.0Inf, 1.0Inf).
?- interval_degree(42,N).
N = 0.
/
347interval_degree(X, N) :- 348 number(X) 349 -> N = 0 350 ; interval_object(X, _, _, Nodelist), % fail if not a number or interval 351 len_nodelist(Nodelist,0,N). % current number of elements in (indefinite) Nodelist 352 353len_nodelist(T,N,N) :- var(T), !. % end of indefinite list 354len_nodelist([_|T],Nin,N) :- 355 Nout is Nin+1, 356 len_nodelist(T,Nout,N). 357 358% internal abstraction 359interval_object(Int, Type, Val, Nodelist) :- 360 get_attr(Int, clpBNR, interval(Type, Val, Nodelist, _)). 361 362% removes constraints in Nodelist 363reset_interval_nodelist_(Int) :- 364 get_attr(Int, clpBNR, Def) -> setarg(3,Def,_) ; true. 365 366% flags (list) abstraction 367get_interval_flags_(Int, Flags) :- 368 get_attr(Int, clpBNR, interval(_, _, _, Flags)). 369 370set_interval_flags_(Int, Flags) :- % flags assumed to be ground so no copy required 371 interval_object(Int, Type, Val, Nodelist), 372 put_attr(Int, clpBNR, interval(Type, Val, Nodelist, Flags)). 373 374% 375% Interval value constants 376% 377universal_interval((-1.0Inf,1.0Inf)). 378 379empty_interval((1.0Inf,-1.0Inf)). 380 381new_universal_interval(boolean,Int) :- !, 382 put_attr(Int, clpBNR, interval(integer, (0,1), _NL, [])). 383new_universal_interval(Type,Int) :- 384 universal_interval(UI), 385 put_attr(Int, clpBNR, interval(Type, UI, _NL, [])). 386 387% Finite intervals - 64 bit IEEE reals, 388finite_interval(real, (-1.0e+16,1.0e+16)). 389finite_interval(integer, (L,H)) :-
390 current_prolog_flag(min_tagged_integer,L), 391 current_prolog_flag(max_tagged_integer,H)
391. 392finite_interval(boolean, (0,1)). 393 394% legal integer bounds values 395integerBnd(1.0Inf). 396integerBnd(-1.0Inf). 397integerBnd(B) :- integer(B). 398 399% precise bounds values - Note: assumes external commit so cuts not req'd in first two clauses 400preciseBnd(1.0Inf). 401preciseBnd(-1.0Inf). 402preciseBnd(1.5NaN) :- !, fail. 403preciseBnd(B) :- preciseNumber(B). 404 405preciseNumber(B) :- 406 rational(B) -> true 407 ; number(B), 0 is cmpr(B,rationalize(B)). % rational(B)=:=rationalize(B), fails if float not precise
Bs = [L,U]; otherwise fails. On backtracking, this value is not undone.
Caution: this predicate is non-logical and intended for specialized use case, e.g., some branch-and-bound algorithms (narrow to current solution, then backtrack to next solution). /
416% 417% non-backtracking set bounds for use with branch and bound 418% 419nb_setbounds(Int, [L,U]) :- 420 number(L), number(U), 421 get_attr(Int, clpBNR, Def), 422 arg(2, Def, Val), % WAM code 423 ^(Val,(L,U),NewVal), % new range is intersection (from ia_primitives) 424 nb_setarg(2, Def, NewVal). 425 426% 427% get current value 428% 429getValue(Int, Val) :- 430 number(Int) 431 -> Val=(Int,Int) % numeric point value 432 ; get_attr(Int, clpBNR, interval(_, Val, _, _)). % interval, optimized for SWIP 433 434% 435% set interval value (assumes bounds check has already been done) 436% Note: putValue_ cannot modify the new value since the narrowing op is still in the Agenda 437% and may not be run again. Insert `integral` op for integer bounds adjustment. 438% 439putValue_(New, Int, NodeList) :- 440 get_attr(Int, clpBNR, Def) % still an interval ? 441 -> (debugging(clpBNR) -> check_monitor_(Int, New, Def) ; true), 442 Def = interval(Type,_,Nodes,_), % get Type and Nodes for queueing 443 New = (L,H), 444 ( 0 is cmpr(L,H) % point interval -> 445 -> setarg(3,Def,_NL), % clear node list (so nothing done in unify) 446 pointValue_(L,H,Int), % unify, avoiding -0.0 and preferring rational (invokes hook) 447 NodeList = Nodes % still add Nodes to Agenda 448 ; setarg(2,Def,New), % update value in interval (backtrackable) 449 % if integer has non-integral bounds, schedule `integral` op to fix it 450 (Type == integer 451 -> ( integerBnd(L), integerBnd(H) -> NodeList = Nodes ; NodeList = [node(integral,_,0,$(Int))|_] ) 452 ; NodeList = Nodes 453 ) 454 ) 455 ; true. % no longer an interval, NodeList is empty 456 457pointValue_(-0.0,_,0.0) :-!. 458pointValue_(L,H,Int) :- (rational(L) -> Int = L ; Int = H).
range(X,[2,3]) is equivalent to X::real(2,3). If X is a number the lower and upper bounds are the same. Examples:
?- X::integer(1,10), range(X,Bs). Bs = [1, 10], X::integer(1, 10). ?- range(42,Bs). Bs = [42, 42]. ?- range(X,[2,3]). X::real(2, 3).
/
476% 477% range(Int, Bounds) for compatibility 478% 479range(Int, [L,H]) :- getValue(Int, (IL,IH)), !, % existing interval or number => number(IL) 480 (var(L) -> L=IL ; non_empty(L,IL)), % range check (L=<IL,IH=<H), no narrowing 481 (var(H) -> H=IH ; non_empty(IH,H)). 482range(Int, [L,H]) :- var(Int), % for other var(Int), declare it to a real interval with specified bounds 483 Int::real(L,H).
integer(0,1), Dom will be boolean. Examples:
?- range(X,[2,3]), domain(X,Dom). Dom = real(2, 3), X::real(2, 3). ?- X::integer(0,1),domain(X,Dom). Dom = boolean, X::boolean. ?- domain(X,Dom). false.
503domain(Int, Dom) :- 504 interval_object(Int, Type, Val, _), 505 interval_domain_(Type, Val, Dom). 506 507interval_domain_(integer,(0,1),boolean) :- !. % integer(0,1) is synonymous with boolean 508interval_domain_(T,(L,H),Dom) :- Dom=..[T,L,H]. 509 510:- use_module(library(arithmetic), [arithmetic_function/1]). % extended arithmetic functions
upper bound-lowerbound); otherwise fails. Examples:
?- X:: real(1r2,5r3),delta(X,D). D = 7r6, X::real(0.5, 1.6666666666666667). ?- delta(42,W). W = 0.
delta is also available as an arithmetic function:
?- X::real(1r2,pi), W is delta(X). W = 2.6415926535897936, X::real(0.5, 3.1415926535897936).
/
531% 532% delta(Int, Wid) width/span of an interval or numeric value, can be infinite 533% 534:- arithmetic_function(delta/1). 535 536delta(Int, Wid) :- 537 getValue(Int,(L,H)), 538 Wid is roundtoward(H-L,to_positive).
?- X:: real(1r2,5r3), midpoint(X,M). M = 13r12, X::real(0.5, 1.6666666666666667). ?- midpoint(42,M). M = 42.
midpoint is also available as an arithmetic function:
?- X::real(1r2,pi), M is midpoint(X). M = 1.8207963267948968, X::real(0.5, 3.1415926535897936).
/
559% 560% midpoint(Int, Wid) midpoint of an interval or numeric value 561% based on: 562% Frédéric Goualard. How do you compute the midpoint of an interval?. 563% ACM Transactions on Mathematical Software, Association for Computing Machinery, 2014, 564% 40 (2), 10.1145/2493882. hal-00576641v1 565% Exception, single infinite bound treated as largest finite FP value 566% 567:- arithmetic_function(midpoint/1). 568 569midpoint(Int, Mid) :- 570 getValue(Int,(L,H)), 571 midpoint_(L,H,Mid). 572 573midpoint_(L,H,M) :- L =:= -H, !, M=0. % symmetric including (-inf,inf) 574midpoint_(-1.0Inf,H,M) :- !, M is nexttoward(-1.0Inf,0)/2 + H/2. 575midpoint_(L,1.0Inf,M) :- !, M is L/2 + nexttoward(1.0Inf,0)/2. 576midpoint_(L,H,M) :- M1 is L/2 + H/2, M=M1. % general case
?- X:: real(1r2,5r3), median(X,M). M = 0.9128709291752769, X::real(0.5, 1.6666666666666667). ?- median(42,M). M = 42.0.
median is also available as an arithmetic function:
?- X::real(1r2,pi), M is median(X). M = 1.2533141373155003, X::real(0.5, 3.1415926535897936).
/
597% 598% median(Int,Med) from CLP(RI) 599% Med = 0 if Int contains 0, else a number which divides Int into equal 600% numbers of FP values. Med is always a float 601% 602:- arithmetic_function(median/1). 603 604median(Int, Med) :- 605 getValue(Int,(L,H)), 606 median_bound_(lo,L,FL), 607 median_bound_(hi,H,FH), 608 median_(FL,FH,Med), !. 609 610median_bound_(lo,B,FB) :- B=:=0, FB is nexttoward(B,1.0). 611median_bound_(lo,-1.0Inf,FB) :- FB is nexttoward(-1.0Inf,0.0). 612median_bound_(lo,B,FB) :- FB is roundtoward(float(B), to_negative). 613 614median_bound_(hi,B,FB) :- B=:=0, !, FB is nexttoward(B,-1.0). 615median_bound_(hi,1.0Inf,FB) :- FB is nexttoward(1.0Inf,0.0). 616median_bound_(hi,B,FB) :- FB is roundtoward(float(B), to_positive). 617 618median_(B,B,B). % point interval 619median_(L,H,0.0) :- L < 0.0, H > 0.0. % contains 0: handles (-inf,inf) 620median_(L,H,M) :- M is copysign(sqrt(abs(L))*sqrt(abs(H)),L). % L and H have same sign 621 622% 623% lower_bound and upper_bound 624%
?- X::integer(1,10),lower_bound(X). X = 1. ?- X = 42, lower_bound(X). X = 42.
Note that lower_bound will unify X with a number on success, but it may fail if this value is inconsistent with current constraints. /
638lower_bound(Int) :-
639 getValue(Int,(L,_H)),
640 Int=L.?- X::integer(1,10),upper_bound(X). X = 10. ?- X = 42, upper_bound(X). X = 42.
Note that upper_bound will unify X with a number on success, but it may fail if this value is inconsistent with current constraints. /
655upper_bound(Int) :-
656 getValue(Int,(_L,H)),
657 Int=H.real(-1.0e+16, 1.0e+16) and integer(-72057594037927936, 72057594037927935). Examples:
?- X::real(-pi/2,pi/2). X::real(-1.5707963267948968, 1.5707963267948968). ?- X::real, Y::integer. X::real(-1.0e+16, 1.0e+16), Y::integer(-72057594037927936, 72057594037927935). ?- Y::integer(1,_), Y::Dom. Dom = integer(1, 72057594037927935), Y::integer(1, 72057594037927935). ?- B::boolean. B::boolean. ?- 42::Dom. false.
Note that bounds can be defined using arithmetic expressions.
Alternatively, the first argument may be a list of variables:
?- [B1,B2,B3]::boolean. B1::boolean, B2::boolean, B3::boolean. ?- length(Vs,3), Vs::real(-1,1). Vs = [_A, _B, _C], _A::real(-1, 1), _B::real(-1, 1), _C::real(-1, 1).
/
697Rs::Dom :- list(Rs),!, % list of vars 698 intervals_(Rs,Dom). 699 700R::Dom :- % single var 701 g_read('$clpBNR:thread_init_done',_), % ensures per-thread initialization 702 (var(Dom) % domain undefined 703 -> (var(R) -> int_decl_(real,_,R) ; true), % default or domain query (if interval(R) else fail) 704 domain(R,Dom) % extract domain 705 ; Dom=..[Type|Bounds], % domain defined 706 Val=..[','|Bounds], 707 int_decl_(Type,Val,R) 708 ). 709 710intervals_([],_Dom). 711intervals_([Int|Ints],Dom) :- 712 copy_term(Dom,CDom), % need a fresh copy of domain for each interval 713 Int::CDom, !, 714 intervals_(Ints,Dom). 715 716int_decl_(boolean,_,R) :- !, % boolean is integer; 0=false, 1=true, ignore any bounds. 717 int_decl_(integer,(0,1),R). 718int_decl_(Type,(','),R) :- !, % no bounds, fill with vars 719 int_decl_(Type,(_,_),R). 720int_decl_(Type,Val,R) :- interval_object(R,CType,CVal,_NL), !, % already interval 721 (Type = CType, Val = CVal % query, no change 722 -> true 723 ; Val = (L,H), % changing type,bounds? 724 lower_bound_val_(Type,L,IL), 725 upper_bound_val_(Type,H,IH), 726 applyType_(Type, R, T/T, Agenda), % coerce reals to integers (or no-op). 727 ^(CVal,(IL,IH),New), % low level functional primitive 728 updateValue_(CVal, New, R, 1, Agenda, NewAgenda), % update value (Agenda empty if no value change) 729 stable_(NewAgenda) % then execute Agenda 730 ). 731int_decl_(Type,(L,H),R) :- var(R), !, % new interval (R can't be existing interval) 732 lower_bound_val_(Type,L,IL), 733 upper_bound_val_(Type,H,IH), 734 C is cmpr(IL,IH), % compare bounds 735 (C == 0 736 -> (rational(IL) -> R=IL ; R = IH) % point range, can unify now 737 ; C == -1, % necessary condition: IL < IH 738 put_attr(R, clpBNR, interval(Type, (IL,IH), _NL, [])) % attach clpBNR attribute 739 ). 740int_decl_(Type,(L,H),R) :- constant_type_(Type,R), % R already a point value, check range 741 lower_bound_val_(Type,L,IL), non_empty(IL,R), % IL=<R, 742 upper_bound_val_(Type,H,IH), non_empty(R,IH). % R=<IH. 743 744lower_bound_val_(Type,L,L) :- var(L), !, % unspecified bound, make it finite 745 finite_interval(Type,(L,_)). 746lower_bound_val_(real,L,IL) :- % real: evaluate and round outward (if float) 747 ((L == pi ; L == e) 748 -> IL is roundtoward(L,to_negative) 749 ; Lv is L, 750 (preciseBnd(Lv) -> IL=Lv ; IL is nexttoward(Lv,-1.0Inf)) 751 ). 752lower_bound_val_(integer,L,IL) :- % integer: make integer 753 IL is ceiling(L). 754lower_bound_val_(boolean,L,IL) :- % boolean: narrow to L=0 755 IL is max(0,ceiling(L)). 756 757upper_bound_val_(Type,H,H) :- var(H), !, % unspecified bound, make it finite 758 finite_interval(Type,(_,H)). 759upper_bound_val_(real,H,IH) :- % real: evaluate and round outward (if float) 760 ((H == pi ; H == e) 761 -> IH is roundtoward(H,to_positive) 762 ; Hv is H, 763 (preciseBnd(Hv) -> IH=Hv ; IH is nexttoward(Hv,1.0Inf)) 764 ). 765upper_bound_val_(integer,H,IH) :- % integer: make integer 766 IH is floor(H). 767upper_bound_val_(boolean,H,IH) :- % boolean: narrow to H=1 768 IH is min(1,floor(H)). 769 770constant_type_(real,C) :- number(C). 771constant_type_(integer,C) :- integer(C), !. 772constant_type_(integer,C) :- float(C), float_class(C,infinite). 773 774applyType_(NewType, Int, Agenda, NewAgenda) :- % narrow Int to Type 775 get_attr(Int,clpBNR,interval(Type,Val,NodeList,Flags)), 776 (NewType @< Type % standard order of atoms: boolean @< integer @< real 777 -> (debugging(clpBNR) -> check_monitor_(Int, NewType, interval(Type,Val,NodeList,Flags)) ; true), 778 Val = (L,H), 779 lower_bound_val_(NewType,L,IL), 780 upper_bound_val_(NewType,H,IH), 781 (IL == IH 782 -> Int=IL % narrowed to point 783 ; (put_attr(Int,clpBNR,interval(integer,(IL,IH),NodeList,Flags)), % set Type (only case allowed) 784 linkNodeList_(NodeList, Agenda, NewAgenda) 785 ) 786 ) 787 ; NewAgenda = Agenda 788 ). 789 790% 791% this goal gets triggered whenever an interval is unified, valid for a numeric value or another interval 792% 793attr_unify_hook(IntDef, Num) :- % unify an interval with a numeric 794 number(Num), 795 IntDef = interval(Type,(L,H),Nodelist,_Flags), 796 constant_type_(Type,Num), % numeric consistent with type 797 % L=<Num, Num=<H, assumes L < H 798 cmpr(L,Num) + cmpr(Num,H) < 0, !, % and in range (not NaN) 799 (debugging(clpBNR) -> monitor_unify_(IntDef, Num) ; true), 800 (var(Nodelist) 801 -> true % nothing to do 802 ; linkNodeList_(Nodelist, T/T, Agenda), 803 stable_(Agenda) % broadcast change 804 ). 805 806attr_unify_hook(interval(Type1,V1,Nodelist1,Flags1), Int) :- % unifying two intervals 807 get_attr(Int, clpBNR, interval(Type2,V2,Nodelist2,Flags2)), %%SWIP attribute def. 808 ^(V1,V2,R), % unified range is intersection (from ia_primitives), 809 mergeValues_(Type1, Type2, NewType, R, NewR), !, 810 mergeNodes_(Nodelist2,Nodelist1,Newlist), % unified node list is a merge of two lists 811 mergeFlags_(Flags1,Flags2,Flags), 812 (debugging(clpBNR) -> monitor_unify_(interval(Type1,V1,_,Flags), Int) ; true), 813 % update new type, value and constraint list, undone on backtracking 814 put_attr(Int,clpBNR,interval(NewType,NewR,Newlist,Flags)), 815 (var(Newlist) 816 -> true % nothing to do 817 ; linkNodeList_(Newlist, T/T, Agenda), 818 stable_(Agenda) % broadcast change 819 ). 820 821attr_unify_hook(interval(Type,Val,_Nodelist,_Flags), V) :- % new value out of range 822 g_inc('$clpBNR:evalNodeFail'), % count of primitive call failures 823 debugging(clpBNR), % fail immediately unless debugging 824 debug_clpBNR_('Failed to unify ~w(~w) with ~w',[Type,Val,V]), 825 fail. 826 827% awkward monitor case because original interval gone 828monitor_unify_(IntDef, Update) :- % debugging, manufacture temporary interval 829 put_attr(Temp,clpBNR,IntDef), 830 check_monitor_(Temp, Update, IntDef). 831 832% if types identical, result type and bounds unchanged, 833% else one is integer so result type integer, and integer bounds applied 834mergeValues_(T, T, T, R, R) :- !. 835mergeValues_(_, _, integer, (L,H), (IL,IH)) :- 836 lower_bound_val_(integer,L,IL), % type check bounds 837 upper_bound_val_(integer,H,IH), 838 non_empty(IL,IH). % still non-empty 839 840% optimize for one or both lists (dominant case) 841mergeFlags_([],Flags2,Flags2) :- !. 842mergeFlags_(Flags1,[],Flags1) :- !. 843mergeFlags_([F1|Flags1],Flags2,Flags) :- % discard if F in Flags2 844 functor(F1,N,1), % ignore value 845 functor(F2,N,1), 846 memberchk(F2,Flags2), !, 847 mergeFlags_(Flags1,Flags2,Flags). 848mergeFlags_([F1|Flags1],Flags2,[F1|Flags]) :- % keep F, not in Flags2 849 mergeFlags_(Flags1,Flags2,Flags). 850 851% merge two node lists removing duplicates based on operation and arguments 852mergeNodes_([N],NodeList,NodeList) :- var(N),!. % end of list 853mergeNodes_([node(Op,_,_,Ops)|Ns],NodeList,NewList) :- % if same Op and Ops, discard 854 matching_node_(NodeList,Op,Ops), !, 855 mergeNodes_(Ns,NodeList,NewList). 856mergeNodes_([N|Ns],NodeList,[N|NewList]) :- % not a duplicate, retain 857 mergeNodes_(Ns,NodeList,NewList). 858 859matching_node_([node(Op,_,_,NOps)|_Ns],Op,Ops) :- 860 NOps==Ops, !. % identical args 861matching_node_([N|Ns],Op,Ops) :- 862 nonvar(N), % not end of list 863 matching_node_(Ns,Op,Ops).
Table of supported interval relations:
+ - * / | arithmetic |
** | includes real exponent, odd/even integer |
abs | absolute value |
sqrt | positive square root |
min max | binary min/max |
== is <> =\= =< >= < > | comparison (is and =\= synonyms for == and <>) |
<= | included (one way narrowing) |
and or nand nor xor -> , | boolean (`,` synonym for and) |
- ~ | unary negate and not (boolean) |
exp log | exp/ln |
sin asin cos acos tan atan | trig functions |
integer | must be an integer value |
sig | signum of real (-1,0,+1) |
clpBNR defines the following additional operators for use in constraint expressions:
op(200, fy, ~) | boolean 'not' |
op(500, yfx, and) | boolean 'and' |
op(500, yfx, or) | boolean 'or' |
op(500, yfx, nand) | boolean 'nand' |
op(500, yfx, nor) | boolean 'nor' |
Note that the comparison operators <>, =\=, '<' and '>' are unsound (due to incompleteness) over the real domain but sound over the integer domain. Strict inequality (<> and =\=) is disallowed for type real (will be converted to type integer) but < and > are supported for reals since they may be useful for things like branch and bound searches (with caution). The boolean functions are restricted to type 'boolean', constraining their argument types accordingly. Some examples:
?- {X == Y+1, Y >= 1}.
X::real(2, 1.0Inf),
Y::real(1, 1.0Inf).
?- {X == cos(X)}.
X:: 0.73908513321516... .
?- X::real, {X**4-4*X**3+4*X**2-4*X+3==0}.
X::real(-1.509169756145379, 4.18727500493995).
?- {A or B, C and D}.
C = D, D = 1,
A::boolean,
B::boolean.
Note that any variable in a constraint expression with no domain will be assigned the most general value consistent with the operator types, e.g., real(-1.0Inf,1.0Inf), boolean, etc.
/
914{Cons} :- 915 g_read('$clpBNR:thread_init_done',_), % ensures per-thread initialization 916 term_variables(Cons, CVars), % predeclare variables to maintain ordering 917 declare_vars_(CVars), % convert any plain vars to intervals 918 addConstraints_(Cons,T/T,Agenda), % add constraints 919 stable_(Agenda). % then execute Agenda
{}/1. (Conjunction of contraints is supported using `,` operator).
add_constraint/1 is a primitive version of {}/1, without expression simplification, reducing overhead in some scenarios.
*/
928add_constraint(Con) :- 929 g_read('$clpBNR:thread_init_done',_), % ensures per-thread initialization 930 term_variables(Con, CVars), % predeclare variables to maintain ordering 931 declare_vars_(CVars), 932 constraint_(Con), % a constraint is a boolean expression that evaluates to true 933 buildConstraint_(Con,T/T,Agenda), 934 stable_(Agenda). 935 936% **Obsolete**, use public `add_constraint` 937% low overhead version for internal use (also used by arithmetic_types pack) 938constrain_(C) :- 939 add_constraint(C). 940 941declare_vars_([]). 942declare_vars_([CV|CVars]) :- 943 % if not already an interval, constrain to real(-inf,inf) 944 (interval(CV) -> true ; new_universal_interval(real,CV)), 945 declare_vars_(CVars). 946 947addConstraints_([],Agenda,Agenda) :- !. 948addConstraints_([C|Cs],Agenda,NewAgenda) :- 949 nonvar(C), 950 addConstraints_(C,Agenda,NextAgenda), !, 951 addConstraints_(Cs,NextAgenda,NewAgenda). 952addConstraints_((C,Cs),Agenda,NewAgenda) :- % Note: comma as operator 953 nonvar(C), 954 addConstraints_(C,Agenda,NextAgenda), !, 955 addConstraints_(Cs,NextAgenda,NewAgenda). 956addConstraints_(C,Agenda,NewAgenda) :- 957 constraint_(C), % a constraint is a boolean expression that evaluates to true 958 simplify(C,CS), % possible rewrite 959 buildConstraint_(CS, Agenda, NewAgenda). 960 961buildConstraint_(C,Agenda,NewAgenda) :- 962 debug_clpBNR_('Add ~p',{C}), 963 % need to catch errors from ground expression evaluation 964 catch(build_(C, 1, boolean, Agenda, NewAgenda),_Err,fail), !. 965buildConstraint_(C,_Agenda,_NewAgenda) :- 966 debug_clpBNR_('{} failure due to bad or inconsistent constraint: ~p',{C}), 967 fail. 968 969:- include(clpBNR/ia_simplify). % simplifies constraints to a hopefully more efficient equivalent 970 971% 972% build a node from an expression 973% 974build_(Var, Var, VarType, Agenda, NewAgenda) :- % already an interval or new variable 975 var(Var), !, 976 (interval_object(Var,Type,_,_) % already an interval? 977 -> (Type \== VarType % type narrowing? 978 -> applyType_(VarType, Var, Agenda, NewAgenda) % coerces exisiting intervals to required type 979 ; NewAgenda = Agenda % nothing to do 980 ) 981 ; new_universal_interval(VarType,Var), % implicit interval creation 982 NewAgenda = Agenda % nothing to schedule 983 ). 984build_(::(L,H), Int, VarType, Agenda, Agenda) :- % hidden :: feature: interval bounds literal (without fuzzing) 985 number(L), number(H), !, 986 C is cmpr(L,H), % compare bounds 987 (C == 0 988 -> (rational(L) -> Int=L ; Int=H) % point value, if either bound rational (precise), use it 989 ; C == -1, % necessary condition: L < H 990 (VarType = real -> true ; true), % if undefined Type, use 'real' (efficient ignore/1) 991 put_attr(Int, clpBNR, interval(VarType, (L,H), _NL, [])) % create clpBNR attribute 992 ). 993build_(Num, Int, VarType, Agenda, Agenda) :- % atomic value representing a numeric 994 atomic(Num), !, % includes inf, nan, pi, e 995 % imprecise numbers will be fuzzed 996 (preciseBnd(Num) -> Int = Num ; int_decl_(VarType,(Num,Num),Int)). 997build_(Exp, Num, _, Agenda, Agenda) :- % pre-compile any ground precise Exp 998 ground(Exp), 999 safe_(Exp), % safe to evaluate using is/2 1000 Num is Exp, 1001 preciseBnd(Num), % precise result 1002 !. 1003build_(Exp, Z, _, Agenda, NewAgenda) :- % deconstruct to primitives 1004 Exp =.. [F|Args], 1005 fmap_(F,Op,[Z|Args],ArgsR,Types), !, % supported arithmetic op 1006 build_args_(ArgsR,Objs,Types,Agenda,ObjAgenda), 1007 newNode_(Op,Objs,ObjAgenda,NewAgenda). 1008build_(Exp, Z, _, Agenda, NewAgenda) :- % user defined 1009 Exp =.. [Prim|Args], 1010 chk_primitive_(Prim), 1011 build_args_([Z|Args],Objs,_Types,Agenda,ObjAgenda), 1012 newNode_(user(Prim),Objs,ObjAgenda,NewAgenda). 1013 1014build_args_([],[],_,Agenda,Agenda). 1015build_args_([Arg|ArgsR],[Obj|Objs],[Type|Types],Agenda,NewAgenda) :- 1016 (var(Type) -> Type=real ; true), % var if user defined 1017 build_(Arg,Obj,Type,Agenda,NxtAgenda), 1018 build_args_(ArgsR,Objs,Types,NxtAgenda,NewAgenda). 1019 1020chk_primitive_(Prim) :- % wraps safe usage of unsafe current_predicate/2 1021 UsrHead =..[Prim,'$op',_,_,_], 1022 current_predicate(_,clpBNR:UsrHead). 1023 1024sandbox:safe_primitive(clpBNR:chk_primitive_(_Prim)). 1025 1026% to invoke user defined primitive 1027call_user_primitive(Prim, P, InArgs, OutArgs) :- % wraps unsafe meta call/N 1028 call(clpBNR:Prim, '$op', InArgs, OutArgs, P). 1029 1030% really unsafe, but in a pengine a user can't define any predicates in another module, so this is safe 1031sandbox:safe_meta(clpBNR:call_user_primitive(_Prim, _P, _InArgs, _OutArgs), []). 1032 1033% only called when argument(s) ground 1034safe_(_ xor _) :- !, % clpBNR xor incompatible with `is` xor 1035 fail. 1036safe_(integer(_)) :- !, % clpBNR integer incompatible with `is` integer 1037 fail. 1038safe_(asin(_)) :- !, % asin multi-valued, can't use `is` 1039 fail. 1040safe_(acos(_)) :- !, % acos multi-valued, can't use `is` 1041 fail. 1042safe_(atan(_)) :- !, % atan multi-valued, can't use `is` 1043 fail. 1044safe_(_ ** E) :- !, % ** multi-valued for rational even denominator 1045 rational(E,_N,D), 1046 1 is D mod 2. 1047safe_(F) :- 1048 current_arithmetic_function(F), % evaluable by is/2 1049 F =.. [_Op|Args], 1050 safe_args_(Args). 1051 1052safe_args_([]). 1053safe_args_([A|Args]) :- 1054 (atomic(A) -> true ; safe_(A)), 1055 safe_args_(Args). 1056 1057% a constraint must evaluate to a boolean 1058constraint_(C) :- nonvar(C), C =..[Op|_], fmap_(Op,_,_,_,[boolean|_]), !. 1059 1060% map constraint operator to primitive/arity/types 1061fmap_(+, add, ZXY, ZXY, [real,real,real]). 1062fmap_(-, add, [Z,X,Y], [X,Z,Y], [real,real,real]). % note subtract before minus 1063fmap_(*, mul, ZXY, ZXY, [real,real,real]). 1064fmap_(/, mul, [Z,X,Y], [X,Z,Y], [real,real,real]). 1065fmap_(**, pow, ZXY, ZXY, [real,real,real]). 1066fmap_(min, min, ZXY, ZXY, [real,real,real]). 1067fmap_(max, max, ZXY, ZXY, [real,real,real]). 1068fmap_(==, eq, ZXY, ZXY, [boolean,real,real]). % strict equality 1069fmap_(=:=, eq, ZXY, ZXY, [boolean,real,real]). % Prolog compatible, strict equality 1070fmap_(is, eq, ZXY, ZXY, [boolean,real,real]). 1071fmap_(<>, ne, ZXY, ZXY, [boolean,integer,integer]). 1072fmap_(=\=, ne, ZXY, ZXY, [boolean,integer,integer]). % Prolog compatible 1073fmap_(=<, le, ZXY, ZXY, [boolean,real,real]). 1074fmap_(>=, le, [Z,X,Y], [Z,Y,X], [boolean,real,real]). 1075fmap_(<, lt, ZXY, ZXY, [boolean,real,real]). 1076fmap_(>, lt, [Z,X,Y], [Z,Y,X], [boolean,real,real]). 1077fmap_(<=, in, ZXY, ZXY, [boolean,real,real]). % inclusion/subinterval 1078 1079fmap_(and, and, ZXY, ZXY, [boolean,boolean,boolean]). 1080fmap_(',', and, ZXY, ZXY, [boolean,boolean,boolean]). % for usability 1081fmap_(or, or, ZXY, ZXY, [boolean,boolean,boolean]). 1082fmap_(nand, nand, ZXY, ZXY, [boolean,boolean,boolean]). 1083fmap_(nor, nor, ZXY, ZXY, [boolean,boolean,boolean]). 1084fmap_(xor, xor, ZXY, ZXY, [boolean,boolean,boolean]). 1085fmap_(->, imB, ZXY, ZXY, [boolean,boolean,boolean]). 1086 1087fmap_(sqrt, sqrt, ZX, ZX, [real,real]). % pos. root version vs. **(1/2) 1088fmap_(-, minus, ZX, ZX, [real,real]). 1089fmap_(~, not, ZX, ZX, [boolean,boolean]). 1090fmap_(integer,int, ZX, ZX, [boolean,real]). 1091fmap_(sgn, sgn, ZX, ZX, [integer,real]). % signum, Z::integer(-1,1) 1092fmap_(exp, exp, ZX, ZX, [real,real]). 1093fmap_(log, exp, [Z,X], [X,Z], [real,real]). 1094fmap_(abs, abs, ZX, ZX, [real,real]). 1095fmap_(sin, sin, ZX, ZX, [real,real]). 1096fmap_(asin, sin, [Z,X], [X,Z], [real,real]). 1097fmap_(cos, cos, ZX, ZX, [real,real]). 1098fmap_(acos, cos, [Z,X], [X,Z], [real,real]). 1099fmap_(tan, tan, ZX, ZX, [real,real]). 1100fmap_(atan, tan, [Z,X], [X,Z], [real,real]). 1101 1102% reverse map node info to a facsimile of the original constraint 1103map_constraint_op_(integral,$(V),integral(V)) :- !. 1104map_constraint_op_(user(Func),Args,C) :- !, 1105 remap_(Func,Args,C). 1106map_constraint_op_(Op,Args,C) :- 1107 fmap_(COp,Op,_,_,_), 1108 remap_(COp,Args,C), 1109 !. 1110 1111remap_(Op,$(Z,X,Y),C) :- constraint_(Op), Z==1, !, % simplification for binary constraints 1112 C=..[Op,X,Y]. 1113remap_(Op,$(Z,X),C) :- constraint_(Op), Z==1, !, % simplification for unary constraints (~) 1114 C=..[Op,X]. 1115remap_(Op,$(Z,X,Y),Z==C) :- !, 1116 C=..[Op,X,Y]. 1117remap_(Op,$(Z,X),Z==C) :- 1118 C=..[Op,X]. 1119 1120% 1121% Node constructor 1122% 1123% First clause is an optimization for E1 == E2 1124% In that case just unify E1 and E2; heavy lifting done by attr_unify_hook. 1125newNode_(eq, [Z,X,Y], Agenda, Agenda) :- Z==1, !, % no node required 1126 (number(X), number(Y) -> 0 is cmpr(X,Y) ; X=Y). 1127newNode_(Op, Objs, Agenda, NewAgenda) :- 1128 Args =.. [$|Objs], % store arguments as $/N where N=1..3 1129 NewNode = node(Op, _P, 0, Args), % L=0 1130 addNode_(Objs,NewNode), 1131 % increment count of added nodes, will be decremented on backtracking/failure 1132 g_incb('$clpBNR:node_count'), 1133 linkNode_(Agenda, NewNode, NewAgenda). 1134 1135addNode_([],_Node). 1136addNode_([Arg|Args],Node) :- 1137 (number(Arg) 1138 -> true % if Arg a number, nothing to do 1139 ; interval_object(Arg, _Type, _Val, Nodelist), % else add Node to Arg's Nodelist 1140 newmember(Nodelist, Node) 1141 ), 1142 addNode_(Args,Node). 1143 1144sandbox:safe_global_variable('$clpBNR:node_count'). 1145 1146clpStatistics :- 1147 g_assign('$clpBNR:node_count',0), % reset/initialize node count to 0 1148 fail. % backtrack to reset other statistics. 1149 1150clpStatistic(node_count(C)) :- 1151 g_read('$clpBNR:node_count',C). 1152 1153% extend list with N 1154newmember([X|Xs],N) :- 1155 (nonvar(X) 1156 -> newmember(Xs,N) % not end of (indefinite) list. 1157 ; X = N % end of list, insert N and new tail Xs 1158 ). 1159 1160% 1161% Process Agenda to narrow intervals (fixed point iteration) 1162% 1163stable_([]/[]) :- !. % small optimization when agenda is empty 1164stable_(Agenda) :- 1165 current_prolog_flag(clpBNR_iteration_limit,Ops), % budget for current operation 1166 stableLoop_(Agenda,Ops), 1167 !. % achieved stable state with empty Agenda -> commit. 1168 1169stableLoop_([]/[], OpsLeft) :- !, % terminate successfully when agenda comes to an end 1170 g_read('$clpBNR:iterations',Cur), % maintain "low" water mark (can be negative) 1171 (OpsLeft<Cur -> g_assign('$clpBNR:iterations',OpsLeft) ; true), 1172 (OpsLeft<0 -> E is -OpsLeft, debug_clpBNR_('Iteration throttle limit exceeded by ~w ops.',E) ; true). 1173stableLoop_([Node|Agenda]/T, OpsLeft) :- 1174 Node = node(Op,P,_,Args), % if node on queue ignore link bit (was: Node = node(Op,P,1,Args)) 1175 doNode_(Args, Op, P, OpsLeft, DoOver, Agenda/T, NxtAgenda), % undoable on backtrack 1176 nb_setarg(3,Node,0), % reset linked bit 1177 % if doNode_ requested DoOver and above Ops threshold, put Node back at the end of the queue 1178 (atom(DoOver), OpsLeft > 0 -> linkNode_(NxtAgenda,Node,NewAgenda) ; NewAgenda = NxtAgenda), 1179 RemainingOps is OpsLeft-1, % decrement OpsLeft (can go negative) 1180 stableLoop_(NewAgenda,RemainingOps). 1181 1182% support for max_iterations statistic 1183sandbox:safe_global_variable('$clpBNR:iterations'). 1184 1185clpStatistics :- 1186 current_prolog_flag(clpBNR_iteration_limit,L), 1187 g_assign('$clpBNR:iterations',L), % set "low" water mark to upper limit 1188 fail. % backtrack to reset other statistics. 1189 1190clpStatistic(max_iterations(O/L)) :- 1191 g_read('$clpBNR:iterations',Ops), 1192 current_prolog_flag(clpBNR_iteration_limit,L), 1193 O is L-Ops. % convert "low" water mark to high water mark 1194 1195% 1196% doNode_/7 : Evaluate a node and add new nodes to end of queue. `evalNode` primitives can 1197% fail, resulting in eventual failure of `stable_`, i.e., inconsistent constraint set. 1198% 1199% Note: primitives operate on interval values (sets) only, unaware of common arguments, 1200% so additional consistency checks required on update. 1201% 1202doNode_($(ZArg,XArg,YArg), Op, P, OpsLeft, DoOver, Agenda, NewAgenda) :- % Arity 3 Op 1203 (var(P) % check persistent bit 1204 -> getValue(ZArg,ZVal), 1205 getValue(XArg,XVal), 1206 getValue(YArg,YVal), 1207 evalNode(Op, P, $(ZVal,XVal,YVal), $(NZVal,NXVal,NYVal)), % can fail 1208 % enforce consistency for common arguments by intersecting and redoing as required. 1209 (var(ZArg) % if Z == X and/or Y 1210 -> (ZArg==XArg -> consistent_value_(NZVal,NXVal,NZ1,DoOver) ; NZ1 = NZVal), 1211 (ZArg==YArg -> consistent_value_(NZ1, NYVal,NZ2,DoOver) ; NZ2 = NZ1), 1212 updateValue_(ZVal, NZ2, ZArg, OpsLeft, Agenda, AgendaZ) 1213 ; AgendaZ = Agenda 1214 ), 1215 (var(XArg), XArg==YArg % if X==Y 1216 -> consistent_value_(NXVal,NYVal,NVal,DoOver), 1217 updateValue_(XVal, NVal, XArg, OpsLeft, AgendaZ, NewAgenda) % only one update needed 1218 ; updateValue_(XVal, NXVal, XArg, OpsLeft, AgendaZ, AgendaZX), 1219 updateValue_(YVal, NYVal, YArg, OpsLeft, AgendaZX, NewAgenda) 1220 ) 1221 ; % P = p, trim persistent nodes from all arguments 1222 trim_op_(ZArg), trim_op_(XArg), trim_op_(YArg), 1223 NewAgenda = Agenda 1224 ). 1225 1226doNode_($(ZArg,XArg), Op, P, OpsLeft, DoOver, Agenda, NewAgenda) :- % Arity 2 Op 1227 (var(P) % check persistent bit 1228 -> getValue(ZArg,ZVal), 1229 getValue(XArg,XVal), 1230 evalNode(Op, P, $(ZVal,XVal), $(NZVal,NXVal)), % can fail 1231 % enforce consistency for common arguments by intersecting and redoing as required. 1232 (var(ZArg), ZArg==XArg % if Z==X 1233 -> consistent_value_(NZVal,NXVal,NVal,DoOver), % consistent value, DoOver if needed 1234 updateValue_(ZVal, NVal, ZArg, OpsLeft, Agenda, NewAgenda) % only one update needed 1235 ; updateValue_(ZVal, NZVal, ZArg, OpsLeft, Agenda, AgendaZ), 1236 updateValue_(XVal, NXVal, XArg, OpsLeft, AgendaZ, NewAgenda) 1237 ) 1238 ; % P = p, trim persistent nodes from all arguments 1239 trim_op_(ZArg), trim_op_(XArg), 1240 NewAgenda = Agenda 1241 ). 1242 1243doNode_($(Arg), Op, P, _OpsLeft, _, Agenda, NewAgenda) :- % Arity 1 Op 1244 (var(P) % check persistent bit 1245 -> getValue(Arg,Val), 1246 evalNode(Op, P, $(Val), $(NVal)), % can fail causing stable_ to fail => backtracking 1247 updateValue_(Val, NVal, Arg, 1, Agenda,NewAgenda) % always update value regardless of OpsLeft limiter 1248 ; % P = p, trim persistent nodes from argument 1249 trim_op_(Arg), 1250 NewAgenda = Agenda 1251 ). 1252 1253consistent_value_(Val,Val,Val,_) :- !. % same value 1254consistent_value_(Val1,Val2,Val,true) :- ^(Val1,Val2,Val). % different values, intersect 1255 1256% 1257% remove any persistent nodes from Arg 1258% called whenever a persistent node is encountered in FP iteration 1259% 1260trim_op_(Arg) :- 1261 number(Arg) 1262 -> true % if a number, nothing to trim 1263 ; get_attr(Arg, clpBNR, Def), % an interval 1264 arg(3,Def,NList), % Def = interval(_, _, NList, _), 1265 trim_persistent_(NList,TrimList), 1266 % if trimmed list empty, set to a new unshared var to avoid cycles(?) on backtracking 1267 (var(TrimList) -> setarg(3,Def,_) ; setarg(3,Def,TrimList)). % write trimmed node list 1268 1269trim_persistent_(T,T) :- var(T), !. % end of indefinite node list 1270trim_persistent_([node(_,P,_,_)|Ns],TNs) :- nonvar(P), !, trim_persistent_(Ns,TNs). 1271trim_persistent_([N|Ns],[N|TNs]) :- trim_persistent_(Ns,TNs). 1272 1273% 1274% Any changes in interval values should come through here. 1275% Note: This captures all updated state for undoing on backtracking 1276% 1277updateValue_(Old, New, Int, OpsLeft, Agenda, NewAgenda) :- % set interval value to New 1278 Old \== New, % if value changes 1279 (OpsLeft>0 ; propagate_if_(Old, New)), !, % if OpsLeft >0 or narrowing sufficent 1280 putValue_(New, Int, Nodelist), % update value (may fail) 1281 linkNodeList_(Nodelist, Agenda, NewAgenda). % then propagate change 1282updateValue_(_, _, _, _, Agenda, Agenda). % otherwise just continue with Agenda 1283 1284% propgate if sufficient narrowing (> 10%) 1285propagate_if_((OL,OH), (NL,NH)) :- (NH-NL)/(OH-OL) < 0.9. % any overflow in calculation will propagate 1286 1287linkNodeList_(X, List, List) :- var(X), !. % end of indefinite list 1288linkNodeList_([X|Xs], List, NewList) :- 1289 (arg(3,X,Linked), Linked == 1 % test linked flag (only SWIP VM codes) 1290 -> linkNodeList_(Xs, List, NewList) % add rest of the node list 1291 ; linkNode_(List, X, NextList), % not linked add it to list 1292 linkNodeList_(Xs, NextList, NewList) % add rest of the node list 1293 ). 1294 1295linkNode_(List/[X|NextTail], X, List/NextTail) :- % add to list 1296 setarg(3,X,1). % set linked bit 1297 1298:- include(clpBNR/ia_utilities). % print,solve, etc.
none, a watchpoint is placed on X. Watchpoints are only "actioned" when the debug topic clpBNR is enabled. If Action = log a debug message is printed when the interval doman narrows. If Action = trace the debugger is invoked. If Action = none the watchpoint is removed.
/1305watch(Int,Action) :- 1306 atom(Action), 1307 current_module(clpBNR), % this just marks watch/2 as unsafe regarding body 1308 get_interval_flags_(Int,Flags), !, 1309 remove_(Flags,watch(_),Flags1), 1310 (Action = none -> Flags2=Flags1 ; Flags2=[watch(Action)|Flags1]), 1311 set_interval_flags_(Int,Flags2). 1312watch(Ints,Action) :- 1313 list(Ints), 1314 watch_list_(Ints,Action). 1315 1316remove_([],_,[]). 1317remove_([X|Xs],X,Xs) :- !. 1318remove_([X|Xs],X,[X|Ys]) :- 1319 remove_(Xs,X,Ys). 1320 1321watch_list_([],_Action). 1322watch_list_([Int|Ints],Action) :- 1323 watch(Int,Action), 1324 watch_list_(Ints,Action). 1325 1326% check if watch enabled on this interval 1327check_monitor_(Int, Update, interval(_Type,_Val,_Nodelist,Flags)) :- 1328 (memberchk(watch(Action), Flags) 1329 -> once(monitor_action_(Action, Update, Int)) 1330 ; true 1331 ). 1332 1333% 1334% for monitoring changes - all actions defined here 1335% 1336monitor_action_(trace, Update, Int) :- !, % log changes to console and enter debugger 1337 monitor_action_(log, Update, Int), 1338 trace. 1339monitor_action_(log, Update, Int) :- var(Update), !, % log interval unify 1340 debug_clpBNR_('Unify ~p with ~p',[Int,Update]). 1341monitor_action_(log, Update, Int) :- number(Update), !, % log interval unify with number 1342 domain(Int,Dom), 1343 debug_clpBNR_('Unify _?{~p} with ~p',[Dom,Update]). 1344monitor_action_(log, integer, Int) :- !, % change type to integer 1345 debug_clpBNR_('Set type of ~p to ~p',[Int,integer]). 1346monitor_action_(log, Val, Int) :- !, % narrow range 1347 debug_clpBNR_('Set value of ~p to (~p)',[Int,Val]). 1348monitor_action_(_, _, _). % default to noop (as in 'none') 1349 1350sandbox:safe_primitive(clpBNR:watch(_Int,Action)) :- % watch(X,trace) is unsafe. 1351 Action \= trace. 1352% only safe because watch(X,trace) is unsafe. 1353sandbox:safe_primitive(clpBNR:monitor_action_(_Action, _Update, _Int)).
clpBNR trace flag or if the trace flag can be set to B (true or false); otherwise fails. If the trace flag is true and the clpBNR debug topic is enabled, a trace of the fixed point iteration is displayed.
/1360% 1361% tracing doNode_ support - using wrap_predicate(:Head, +Name, -Wrapped, +Body)/4 1362% trace_clpBNR/1 is unsafe (wrapping is global) 1363% 1364:- use_module(library(prolog_wrap)). 1365 1366trace_clpBNR(Bool) :- % query or already in defined state 1367 ( current_predicate_wrapper(clpBNR:doNode_(_Args, _Op, _P, _OpsLeft, _DoOver, _Agenda, _NewAgenda), 1368 'clpBNR:doNode_', _Wrapped, _Body) 1369 -> Bool = true ; Bool = false 1370 ), 1371 !. 1372trace_clpBNR(true) :- % turn on wrapper 1373 wrap_predicate(clpBNR:doNode_(Args, Op, _P, _OpsLeft, _DoOver, _Agenda, _NewAgenda), 1374 'clpBNR:doNode_', 1375 Wrapped, 1376 doNode_wrap_(Wrapped, Args,Op)). 1377trace_clpBNR(false) :- % turn off wrapper 1378 unwrap_predicate(clpBNR:doNode_/7, %(Args, Op, P, OpsLeft, DoOver, Agenda, NewAgenda), 1379 'clpBNR:doNode_'). 1380 1381doNode_wrap_(Wrapped, Args,Op) :- 1382 map_constraint_op_(Op,Args,C), 1383 , % execute wrapped doNode_ 1384 debug_clpBNR_("~p.",C). % print trace message , fail messages from evalNode_, attr_unify_hook 1385 1386% 1387% Get all defined statistics 1388% 1389clpStatistics(Ss) :- findall(S, clpStatistic(S), Ss). 1390 1391% end of reset chain succeeds. 1392clpStatistics. 1393 1394% 1395% module initialization 1396% 1397init_clpBNR :- 1398 restore_optimize_flags_, 1399 print_message(informational, clpBNR(versionInfo)), 1400 print_message(informational, clpBNR(arithmeticFlags)). % cautionary, set on first use 1401 1402:- if(false). % test code used to test sandbox worthiness of hooks 1403check_hooks_safety :- % Note: calls must have no side effects 1404 ignore(attr_portray_hook([],_)), % fails 1405 ignore(user:exception(undefined_global_variable,'$clpBNR:thread_init_done',[])), % fails 1406% ignore(term_html:portray('$clpBNR...'(_),_,_,_)), % fails 1407 ignore(user:portray('$clpBNR...'(_))). % fails 1408:- endif. 1409 1410:- version(clpBNR(versionInfo)). % add message to system version 1411 1412:- multifile prolog:message//1. 1413 1414prologmessage(clpBNR(versionInfo)) --> 1415 { version(Version) }, 1416 [ '*** clpBNR v~w ***.'-[Version] ]. 1417 1418prologmessage(clpBNR(arithmeticFlags)) --> 1419 [ ' Arithmetic global flags will be set to prefer rationals and IEEE continuation values.'-[] ]. 1420 1421:- initialization(init_clpBNR, now). % Most initialization deferred to "first use" - see user:exception/3
clpBNR: Constraint Logic Programming over Continuous Domain of Reals
CLP(BNR) (
library(clpbnr), henceforth justclpBNR) is a CLP over the domain of real numbers extended with ±∞. Since integers are a proper subset of reals, and booleans (0 or 1) a subset of integers, these "sub-domains" are also supported.Since the set of real numbers is continuous it's not possible to represent an aribitray real number, e.g., π in the finite resources of a computer. So
clpBNRuses intervals to represent the domain of a numeric variable. A real variable X has a domain of (L,U) if L ≤ X ≤ U where L and U are numeric values which can be finitely represented, e.g., floats, integers or rationals.The use of intervals (and interval arithmetic) provides guarantees of completeness and correctness - unlike floating point arithmetic - by sacrificing some precision since calulations using floating point domain bounds will be outward rounded.
Finiteness is guaranteed since intervals can only get narrower over the course of a computation. Certainty is only guaranteed if there are no solutions (i.e., query fails) - final interval values may contain 0, 1, or many solutions. When this occurs, the application can further constrain the solution, e.g., by testing specific (point) values in the domain, or by making use of some external knowledge of the problem being solved.
More extensive documentation and many examples are provided in A Guide to CLP(BNR) (HTML version included with this pack in directory
docs/).Documentation for exported predicates follows. The "custom" types include:
clpBNRattribute