1% 2% CLP(BNR) == Constraints On Boolean, Integer, and Real Intervals 3% 4/* The MIT License (MIT) 5 * 6 * Copyright (c) 2019-2024 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 interval/1, % filter for clpBNR constrained var 33 interval_degree/2, % number of constraints on clpBNR constrained var 34 list/1, % O(1) list filter (also for compatibility) 35 domain/2, range/2, % get type and bounds (domain) 36 delta/2, % width (span) of an interval or numeric (also arithmetic function) 37 midpoint/2, % midpoint of an interval (or numeric) (also arithmetic function) 38 median/2, % median of an interval (or numeric) (also arithmetic function) 39 lower_bound/1, % narrow interval to point equal to lower bound 40 upper_bound/1, % narrow interval to point equal to upper bound 41 42 % additional constraint operators 43 op(200, fy, ~), % boolean 'not' 44 op(500, yfx, and), % boolean 'and' 45 op(500, yfx, or), % boolean 'or' 46 op(500, yfx, nand), % boolean 'nand' 47 op(500, yfx, nor), % boolean 'nor' 48 % op(500, yfx, xor), % boolean 'xor', but operator already defined (P=400) for arithmetic 49 op(700, xfx, <>), % integer not equal 50 op(700, xfx, <=), % included (one way narrowing) 51 52 % utilities 53 print_interval/1, print_interval/2, % pretty print interval with optional stream 54 small/1, small/2, % defines small interval width based on precision value 55 solve/1, solve/2, % solve (list of) intervals using split to find point solutions 56 splitsolve/1, splitsolve/2, % solve (list of) intervals using split 57 absolve/1, absolve/2, % absolve (list of) intervals, narrows by nibbling bounds 58 enumerate/1, % "enumerate" integers 59 global_minimum/2, % find interval containing global minimum(s) for an expression 60 global_minimum/3, % global_minimum/2 with definable precision 61 global_maximum/2, % find interval containing global maximum(s) for an expression 62 global_maximum/3, % global_maximum/2 with definable precision 63 global_minimize/2, % global_minimum/2 plus narrow vars to found minimizers 64 global_minimize/3, % global_minimum/3 plus narrow vars to found minimizers 65 global_maximize/2, % global_maximum/2 plus narrow vars to found maximizers 66 global_maximize/3, % global_maximum/3 plus narrow vars to found maximizers 67 nb_setbounds/2, % non-backtracking set bounds (use with branch and bound) 68 partial_derivative/3, % differentiate Exp wrt. X and simplify 69 clpStatistics/0, % reset 70 clpStatistic/1, % get selected 71 clpStatistics/1, % get all defined in a list 72 watch/2, % enable monitoring of changes for interval or (nested) list of intervals 73 trace_clpBNR/1 % enable/disable tracing of clpBNR narrowing ops 74 ]). 75 76/* missing(?) functionality from original CLP(BNR): utility accumulate/2. */ 77 78/* supported interval relations: 79 80+ - * / %% arithmetic 81** %% includes real exponent, odd/even integer 82abs %% absolute value 83sqrt %% positive square root 84min max %% binary min/max 85== is <> =< >= < > %% comparison (`is` synonym for `==`) 86<= %% included (one way narrowing) 87and or nand nor xor -> , %% boolean (`,` synonym for `and`) 88- ~ %% unary negate and not 89exp log %% exp/ln 90sin asin cos acos tan atan %% trig functions 91integer %% must be an integer value 92 93*/
114version("0.11.7b"). 115 116% debug feature control and messaging 117:- if(exists_source(swish(lib/swish_debug))). 118 :- create_prolog_flag(clpBNR_swish, true, [access(read_only)]). 119 :- use_module(swish(lib/swish_debug)). 120 :- use_module(library(http/html_write)). 121:- else. 122 :- use_module(library(debug)). 123:- endif. 124 125:- use_module(library(prolog_versions)). % SWIP dependency enforcement 126 127:- require_prolog_version('9.1.22', % properly exported arithmetic functions 128 [ rational % require rational number support, implies bounded=false 129 ]). 130 131% 132% Optimize arithmetic, but not debug. Only called via directives. 133% 134set_optimize_flags_ :- % called at start of load 135 set_prolog_flag(optimise,true), % scoped to file/module 136 current_prolog_flag(optimise_debug,ODflag), % save; restore in initialization 137 nb_linkval('$optimise_debug_save',ODflag), 138 set_prolog_flag(optimise_debug,false). % so debug/3, debugging/1 don't get "optimized" 139 140restore_optimize_flags_ :- % called at module initialization (after load) 141 nb_getval('$optimise_debug_save',ODflag), nb_delete('$optimise_debug_save'), 142 set_prolog_flag(optimise_debug,ODflag). 143 144:- set_optimize_flags_. 145 146% local debug and trace message support 147debug_clpBNR_(FString,Args) :- debug(clpBNR,FString,Args). 148 149% sandboxing for SWISH 150:- multifile(sandbox:safe_prolog_flag/1). 151:- multifile(sandbox:safe_global_variable/1). 152:- multifile(sandbox:safe_primitive/1). 153:- multifile(sandbox:safe_meta/2). 154 155current_node_(Node) :- % look back to find current Op being executed for debug messages 156 prolog_current_frame(Frame), % this is a little grungy, but necessary to get intervals 157 prolog_frame_attribute(Frame,parent_goal,doNode_(Args,Op,_,_,_,_,_)), 158 map_constraint_op_(Op,Args,Node), 159 !. 160 161sandbox:safe_primitive(clpBNR:current_node_(_Node)). 162 163% 164% statistics 165% 166 167% assign,increment/read global counter (assumed to be ground value so use _linkval) 168g_assign(G,V) :- nb_linkval(G,V). 169g_inc(G) :- nb_getval(G,N), N1 is N+1, nb_linkval(G,N1). 170g_incb(G) :- nb_getval(G,N), N1 is N+1, b_setval(G,N1). % undone on backtrack 171g_read(G,V) :- nb_getval(G,V). 172 173sandbox:safe_global_variable('clpBNR:thread_init_done'). 174sandbox:safe_global_variable('clpBNR:userTime'). 175sandbox:safe_global_variable('clpBNR:inferences'). 176sandbox:safe_global_variable('clpBNR:gc_time'). 177 178% 179% Global var statistics are per thread and therefore must be "lazily" initialized 180% Also ensures that thread copies of flags are properly set. 181% This exception handler will be invoked the first time 'clpBNR:thread_init_done' is read 182% Public predicates ::, {}, clpStatistics/1 and range read this global before proceeding. 183% 184user:exception(undefined_global_variable,'clpBNR:thread_init_done',retry) :- !, 185 set_prolog_flags, % initialize/check environment flags 186 clpStatistics, % defines remaining globals 187 g_assign('clpBNR:thread_init_done',1). 188 189% 190% Define custom clpBNR flags when module loaded 191% 192:- create_prolog_flag(clpBNR_iteration_limit,3000,[type(integer),keep(true)]). 193:- create_prolog_flag(clpBNR_default_precision,6,[type(integer),keep(true)]). 194:- create_prolog_flag(clpBNR_verbose,false,[type(boolean),keep(true)]). 195 196sandbox:safe_prolog_flag(clpBNR_iteration_limit,_). 197sandbox:safe_prolog_flag(clpBNR_default_precision,_). 198sandbox:safe_prolog_flag(clpBNR_verbose,_). 199% 200% Set public flags (see module/thread initialization) 201% 202set_prolog_flags :- 203 set_prolog_flag(prefer_rationals, true), % enable rational arithmetic 204 (current_prolog_flag(max_rational_size,_) 205 -> true % already defined, leave as is 206 ; set_prolog_flag(max_rational_size, 16) % rational size in bytes before .. 207 ), 208 set_prolog_flag(max_rational_size_action, float), % conversion to float 209 210 set_prolog_flag(float_overflow,infinity), % enable IEEE continuation values 211 set_prolog_flag(float_zero_div,infinity), 212 set_prolog_flag(float_undefined,nan), 213 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).
/
279:- discontiguous clpBNR:clpStatistics/0, clpBNR:clpStatistic/1. 280 281clpStatistics :- 282 % garbage_collect, % ? do gc before time snapshots 283 statistics(cputime,T), g_assign('clpBNR:userTime',T), % thread based 284 statistics(inferences,I), g_assign('clpBNR:inferences',I), 285 statistics(garbage_collection,[_,_,G,_]), g_assign('clpBNR:gc_time',G), 286 fail. % backtrack to reset other statistics. 287 288clpStatistic(_) :- g_read('clpBNR:thread_init_done',0). % ensures per-thread initialization then fails 289 290clpStatistic(userTime(T)) :- statistics(cputime,T1), g_read('clpBNR:userTime',T0), T is T1-T0. 291 292clpStatistic(gcTime(G)) :- statistics(garbage_collection,[_,_,G1,_]), g_read('clpBNR:gc_time',G0), G is (G1-G0)/1000.0. 293 294clpStatistic(globalStack(U/T)) :- statistics(globalused,U), statistics(global,T). 295 296clpStatistic(trailStack(U/T)) :- statistics(trailused,U), statistics(trail,T). 297 298clpStatistic(localStack(U/T)) :- statistics(localused,U), statistics(local,T). 299 300clpStatistic(inferences(I)) :- statistics(inferences,I1), g_read('clpBNR:inferences',I0), I is I1-I0.
308list(X) :- compound(X) -> X=[_|_] ; X==[]. 309 310:- include(clpBNR/ia_primitives). % interval arithmetic relations via evalNode/4. 311 312%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 313% 314% SWI-Prolog implementation of IA 315% 316%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 317% 318% Intervals are constrained (attributed) variables. 319% 320% Current interval bounds updates via setarg(Val) which is backtrackable
clpBNR attribute; otherwise fails.
/
326interval(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.
/
342interval_degree(X, N) :- 343 interval_object(X, _, _, Nodelist) 344 -> len_nodelist(Nodelist,0,N) % current number of elements in (indefinite) Nodelist 345 ; number(X), N = 0. % number -> no constraints ; fail 346 347len_nodelist([T],N,N) :- var(T), !. % end of indefinite list 348len_nodelist([_|T],Nin,N) :- 349 Nout is Nin+1, 350 len_nodelist(T,Nout,N). 351 352% internal abstraction 353interval_object(Int, Type, Val, Nodelist) :- 354 get_attr(Int, clpBNR, interval(Type, Val, Nodelist, _)). 355 356% removes constraints in Nodelist 357reset_interval_nodelist_(Int) :- 358 get_attr(Int, clpBNR, Def) -> setarg(3,Def,_) ; true. 359 360% flags (list) abstraction 361get_interval_flags_(Int, Flags) :- 362 get_attr(Int, clpBNR, interval(_, _, _, Flags)). 363 364set_interval_flags_(Int, Flags) :- % flags assumed to be ground so no copy required 365 interval_object(Int, Type, Val, Nodelist), 366 put_attr(Int, clpBNR, interval(Type, Val, Nodelist, Flags)). 367 368% 369% Interval value constants 370% 371universal_interval((-1.0Inf,1.0Inf)). 372 373empty_interval((1.0Inf,-1.0Inf)). 374 375% Finite intervals - 64 bit IEEE reals, 376finite_interval(real, (-1.0e+16,1.0e+16)). 377finite_interval(integer, (L,H)) :-
378 current_prolog_flag(min_tagged_integer,L), 379 current_prolog_flag(max_tagged_integer,H)
379. 380finite_interval(boolean, (0,1)). 381 382% legal integer bounds values 383integerBnd(1.0Inf). 384integerBnd(-1.0Inf). 385integerBnd(B) :- integer(B). 386 387% precise bounds values 388preciseBnd(1.0Inf). 389preciseBnd(-1.0Inf). 390preciseBnd(1.5NaN) :- !, fail. 391preciseBnd(B) :- 392 rational(B) -> true 393 ; 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). /
402% 403% non-backtracking set bounds for use with branch and bound 404% 405nb_setbounds(Int, [L,U]) :- 406 number(L), number(U), 407 get_attr(Int, clpBNR, Def), 408 arg(2, Def, Val), % WAM code 409 ^(Val,(L,U),NewVal), % new range is intersection (from ia_primitives) 410 nb_setarg(2, Def, NewVal). 411 412% 413% get current value 414% 415getValue(Int, Val) :- 416 number(Int) 417 -> Val=(Int,Int) % numeric point value 418 ; get_attr(Int, clpBNR, interval(_, Val, _, _)). % interval, optimized for SWIP 419 420% 421% set interval value (assumes bounds check has already been done) 422% Note: putValue_ cannot modify the new value since the narrowing op is still in the Agenda 423% and may not be run again. Insert `integral` op for integer bounds adjustment. 424% 425putValue_(New, Int, NodeList) :- 426 get_attr(Int, clpBNR, Def) % still an interval ? 427 -> (debugging(clpBNR) -> check_monitor_(Int, New, Def) ; true), 428 Def = interval(Type,_,Nodes,_), % get Type and Nodes for queueing 429 New = (L,H), 430 ( 0 is cmpr(L,H) % point interval -> 431 -> setarg(3,Def,_NL), % clear node list (so nothing done in unify) 432 pointValue_(L,H,Int), % unify, avoiding -0.0 and preferring rational (invokes hook) 433 NodeList = Nodes % still add Nodes to Agenda 434 ; setarg(2,Def,New), % update value in interval (backtrackable) 435 % if integer has non-integral bounds, schedule `integral` op to fix it 436 (Type == integer 437 -> ( integerBnd(L), integerBnd(H) -> NodeList = Nodes ; NodeList = [node(integral,_,0,$(Int))|_] ) 438 ; NodeList = Nodes 439 ) 440 ) 441 ; true. % no longer an interval, NodeList is empty 442 443pointValue_(-0.0,_,0.0) :-!. 444pointValue_(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).
/
462% 463% range(Int, Bounds) for compatibility 464% 465range(Int, [L,H]) :- getValue(Int, (IL,IH)), !, % existing interval or number => number(IL) 466 (var(L) -> L=IL ; non_empty(L,IL)), % range check (L=<IL,IH=<H), no narrowing 467 (var(H) -> H=IH ; non_empty(IH,H)). 468range(Int, [L,H]) :- var(Int), % for other var(Int), declare it to a real interval with specified bounds 469 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.
489domain(Int, Dom) :- 490 interval_object(Int, Type, Val, _), 491 interval_domain_(Type, Val, Dom). 492 493interval_domain_(integer,(0,1),boolean) :- !. % integer(0,1) is synonymous with boolean 494interval_domain_(T,(L,H),Dom) :- Dom=..[T,L,H]. 495 496:- 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).
/
517% 518% delta(Int, Wid) width/span of an interval or numeric value, can be infinite 519% 520:- arithmetic_function(delta/1). 521 522delta(Int, Wid) :- 523 getValue(Int,(L,H)), 524 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).
/
545% 546% midpoint(Int, Wid) midpoint of an interval or numeric value 547% based on: 548% Frédéric Goualard. How do you compute the midpoint of an interval?. 549% ACM Transactions on Mathematical Software, Association for Computing Machinery, 2014, 550% 40 (2), 10.1145/2493882. hal-00576641v1 551% Exception, single infinite bound treated as largest finite FP value 552% 553:- arithmetic_function(midpoint/1). 554 555midpoint(Int, Mid) :- 556 getValue(Int,(L,H)), 557 midpoint_(L,H,Mid). 558 559midpoint_(L,H,M) :- L =:= -H, !, M=0. % symmetric including (-inf,inf) 560midpoint_(-1.0Inf,H,M) :- !, M is nexttoward(-1.0Inf,0)/2 + H/2. 561midpoint_(L,1.0Inf,M) :- !, M is L/2 + nexttoward(1.0Inf,0)/2. 562midpoint_(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).
/
583% 584% median(Int,Med) from CLP(RI) 585% Med = 0 if Int contains 0, else a number which divides Int into equal 586% numbers of FP values. Med is always a float 587% 588:- arithmetic_function(median/1). 589 590median(Int, Med) :- 591 getValue(Int,(L,H)), 592 median_bound_(lo,L,FL), 593 median_bound_(hi,H,FH), 594 median_(FL,FH,Med), !. 595 596median_bound_(lo,B,FB) :- B=:=0, FB is nexttoward(B,1.0). 597median_bound_(lo,-1.0Inf,FB) :- FB is nexttoward(-1.0Inf,0.0). 598median_bound_(lo,B,FB) :- FB is roundtoward(float(B), to_negative). 599 600median_bound_(hi,B,FB) :- B=:=0, !, FB is nexttoward(B,-1.0). 601median_bound_(hi,1.0Inf,FB) :- FB is nexttoward(1.0Inf,0.0). 602median_bound_(hi,B,FB) :- FB is roundtoward(float(B), to_positive). 603 604median_(B,B,B). % point interval 605median_(L,H,0.0) :- L < 0.0, H > 0.0. % contains 0: handles (-inf,inf) 606median_(L,H,M) :- M is copysign(sqrt(abs(L))*sqrt(abs(H)),L). % L and H have same sign 607 608% 609% lower_bound and upper_bound 610%
?- 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. /
624lower_bound(Int) :-
625 getValue(Int,(L,_H)),
626 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. /
641upper_bound(Int) :-
642 getValue(Int,(_L,H)),
643 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).
/
683Rs::Dom :- list(Rs),!, % list of vars 684 intervals_(Rs,Dom). 685 686R::Dom :- % single var 687 g_read('clpBNR:thread_init_done',_), % ensures per-thread initialization 688 (var(Dom) % domain undefined 689 -> (var(R) -> int_decl_(real,_,R) ; true), % default or domain query (if interval(R) else fail) 690 domain(R,Dom) % extract domain 691 ; Dom=..[Type|Bounds], % domain defined 692 Val=..[','|Bounds], 693 int_decl_(Type,Val,R) 694 ). 695 696intervals_([],_Def). 697intervals_([Int|Ints],Def) :- 698 Int::Def, !, 699 intervals_(Ints,Def). 700 701int_decl_(boolean,_,R) :- !, % boolean is integer; 0=false, 1=true, ignore any bounds. 702 int_decl_(integer,(0,1),R). 703int_decl_(Type,(','),R) :- !, % no bounds, fill with vars 704 int_decl_(Type,(_,_),R). 705int_decl_(Type,Val,R) :- interval_object(R,CType,CVal,_NL), !, % already interval 706 (Type = CType, Val = CVal % query, no change 707 -> true 708 ; Val = (L,H), % changing type,bounds? 709 lower_bound_val_(Type,L,IL), 710 upper_bound_val_(Type,H,IH), 711 applyType_(Type, R, T/T, Agenda), % coerce reals to integers (or no-op). 712 ^(CVal,(IL,IH),New), % low level functional primitive 713 updateValue_(CVal, New, R, 1, Agenda, NewAgenda), % update value (Agenda empty if no value change) 714 stable_(NewAgenda) % then execute Agenda 715 ). 716int_decl_(Type,(L,H),R) :- var(R), !, % new interval (R can't be existing interval) 717 lower_bound_val_(Type,L,IL), 718 upper_bound_val_(Type,H,IH), 719 C is cmpr(IL,IH), % compare bounds 720 (C == 0 721 -> (rational(IL) -> R=IL ; R = IH) % point range, can unify now 722 ; C == -1, % necessary condition: IL < IH 723 put_attr(R, clpBNR, interval(Type, (IL,IH), _NL, [])) % attach clpBNR attribute 724 ). 725int_decl_(Type,(L,H),R) :- constant_type_(Type,R), % R already a point value, check range 726 lower_bound_val_(Type,L,IL), non_empty(IL,R), % IL=<R, 727 upper_bound_val_(Type,H,IH), non_empty(R,IH). % R=<IH. 728 729lower_bound_val_(Type,L,L) :- var(L), !, % unspecified bound, make it finite 730 finite_interval(Type,(L,_)). 731lower_bound_val_(real,L,IL) :- % real: evaluate and round outward (if float) 732 ((L == pi ; L == e) 733 -> IL is roundtoward(L,to_negative) 734 ; Lv is L, 735 (preciseBnd(Lv) -> IL=Lv ; IL is nexttoward(Lv,-1.0Inf)) 736 ). 737lower_bound_val_(integer,L,IL) :- % integer: make integer 738 IL is ceiling(L). 739lower_bound_val_(boolean,L,IL) :- % boolean: narrow to L=0 740 IL is max(0,ceiling(L)). 741 742upper_bound_val_(Type,H,H) :- var(H), !, % unspecified bound, make it finite 743 finite_interval(Type,(_,H)). 744upper_bound_val_(real,H,IH) :- % real: evaluate and round outward (if float) 745 ((H == pi ; H == e) 746 -> IH is roundtoward(H,to_positive) 747 ; Hv is H, 748 (preciseBnd(Hv) -> IH=Hv ; IH is nexttoward(Hv,1.0Inf)) 749 ). 750upper_bound_val_(integer,H,IH) :- % integer: make integer 751 IH is floor(H). 752upper_bound_val_(boolean,H,IH) :- % boolean: narrow to H=1 753 IH is min(1,floor(H)). 754 755constant_type_(real,C) :- number(C). 756constant_type_(integer,C) :- integer(C), !. 757constant_type_(integer,C) :- float(C), float_class(C,infinite). 758 759applyType_(NewType, Int, Agenda, NewAgenda) :- % narrow Int to Type 760 get_attr(Int,clpBNR,interval(Type,Val,NodeList,Flags)), 761 (NewType @< Type % standard order of atoms: boolean @< integer @< real 762 -> (debugging(clpBNR) -> check_monitor_(Int, NewType, interval(Type,Val,NodeList,Flags)) ; true), 763 Val = (L,H), 764 lower_bound_val_(NewType,L,IL), 765 upper_bound_val_(NewType,H,IH), 766 (IL=IH 767 -> Int=IL % narrowed to point 768 ; (put_attr(Int,clpBNR,interval(integer,(IL,IH),NodeList,Flags)), % set Type (only case allowed) 769 linkNodeList_(NodeList, Agenda, NewAgenda) 770 ) 771 ) 772 ; NewAgenda = Agenda 773 ). 774 775% 776% this goal gets triggered whenever an interval is unified, valid for a numeric value or another interval 777% 778attr_unify_hook(IntDef, Num) :- % unify an interval with a numeric 779 number(Num), 780 IntDef = interval(Type,(L,H),Nodelist,_Flags), 781 constant_type_(Type,Num), % numeric consistent with type 782 % L=<Num, Num=<H, assumes L < H 783 cmpr(L,Num) + cmpr(Num,H) < 0, !, % and in range (not NaN) 784 (debugging(clpBNR) -> monitor_unify_(IntDef, Num) ; true), 785 (var(Nodelist) 786 -> true % nothing to do 787 ; linkNodeList_(Nodelist, T/T, Agenda), 788 stable_(Agenda) % broadcast change 789 ). 790 791attr_unify_hook(interval(Type1,V1,Nodelist1,Flags1), Int) :- % unifying two intervals 792 get_attr(Int, clpBNR, interval(Type2,V2,Nodelist2,Flags2)), %%SWIP attribute def. 793 ^(V1,V2,R), % unified range is intersection (from ia_primitives), 794 mergeValues_(Type1, Type2, NewType, R, NewR), !, 795 mergeNodes_(Nodelist2,Nodelist1,Newlist), % unified node list is a merge of two lists 796 mergeFlags_(Flags1,Flags2,Flags), 797 (debugging(clpBNR) -> monitor_unify_(interval(Type1,V1,_,Flags), Int) ; true), 798 % update new type, value and constraint list, undone on backtracking 799 put_attr(Int,clpBNR,interval(NewType,NewR,Newlist,Flags)), 800 (var(Newlist) 801 -> true % nothing to do 802 ; linkNodeList_(Newlist, T/T, Agenda), 803 stable_(Agenda) % broadcast change 804 ). 805 806attr_unify_hook(interval(Type,Val,_Nodelist,_Flags), V) :- % new value out of range 807 g_inc('clpBNR:evalNodeFail'), % count of primitive call failures 808 debugging(clpBNR), % fail immediately unless debugging 809 debug_clpBNR_('Failed to unify ~w::(~w) with ~w',[Type,Val,V]), 810 fail. 811 812% awkward monitor case because original interval gone 813monitor_unify_(IntDef, Update) :- % debugging, manufacture temporary interval 814 put_attr(Temp,clpBNR,IntDef), 815 check_monitor_(Temp, Update, IntDef). 816 817% if types identical, result type and bounds unchanged, 818% else one is integer so result type integer, and integer bounds applied 819mergeValues_(T, T, T, R, R) :- !. 820mergeValues_(_, _, integer, (L,H), (IL,IH)) :- 821 lower_bound_val_(integer,L,IL), % type check bounds 822 upper_bound_val_(integer,H,IH), 823 non_empty(IL,IH). % still non-empty 824 825% optimize for one or both lists (dominant case) 826mergeFlags_([],Flags2,Flags2) :- !. 827mergeFlags_(Flags1,[],Flags1) :- !. 828mergeFlags_([F1|Flags1],Flags2,Flags) :- % discard if F in Flags2 829 functor(F1,N,1), % ignore value 830 functor(F2,N,1), 831 memberchk(F2,Flags2), !, 832 mergeFlags_(Flags1,Flags2,Flags). 833mergeFlags_([F1|Flags1],Flags2,[F1|Flags]) :- % keep F, not in Flags2 834 mergeFlags_(Flags1,Flags2,Flags). 835 836% merge two node lists removing duplicates based on operation and arguments 837mergeNodes_([N],NodeList,NodeList) :- var(N),!. % end of list 838mergeNodes_([node(Op,_,_,Ops)|Ns],NodeList,NewList) :- % if same Op and Ops, discard 839 matching_node_(NodeList,Op,Ops), !, 840 mergeNodes_(Ns,NodeList,NewList). 841mergeNodes_([N|Ns],NodeList,[N|NewList]) :- % not a duplicate, retain 842 mergeNodes_(Ns,NodeList,NewList). 843 844matching_node_([node(Op,_,_,NOps)|_Ns],Op,Ops) :- 845 NOps==Ops, !. % identical args 846matching_node_([N|Ns],Op,Ops) :- 847 nonvar(N), % not end of list 848 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 |
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.
/
898{Cons} :- 899 g_read('clpBNR:thread_init_done',_), % ensures per-thread initialization 900 term_variables(Cons, CVars), 901 declare_vars_(CVars), % convert any plain vars to intervals 902 addConstraints_(Cons,T/T,Agenda), % add constraints 903 stable_(Agenda). % then execute Agenda 904 905declare_vars_([]). 906declare_vars_([CV|CVars]) :- 907 (interval(CV) -> true ; new_interval_(CV,real)), 908 declare_vars_(CVars). 909 910new_interval_(V,Type) :- 911 universal_interval(UI), 912 int_decl_(Type,UI,V). 913 914addConstraints_([],Agenda,Agenda) :- !. 915addConstraints_([C|Cs],Agenda,NewAgenda) :- 916 nonvar(C), 917 addConstraints_(C,Agenda,NextAgenda), !, 918 addConstraints_(Cs,NextAgenda,NewAgenda). 919addConstraints_((C,Cs),Agenda,NewAgenda) :- % Note: comma as operator 920 nonvar(C), 921 addConstraints_(C,Agenda,NextAgenda), !, 922 addConstraints_(Cs,NextAgenda,NewAgenda). 923addConstraints_(C,Agenda,NewAgenda) :- 924 constraint_(C), % a constraint is a boolean expression that evaluates to true 925 simplify(C,CS), % optional 926 buildConstraint_(CS, Agenda, NewAgenda). 927 928% low overhead version for internal use (also used by arithmetic_types pack) 929constrain_(C) :- 930 buildConstraint_(C,T/T,Agenda), 931 stable_(Agenda). 932 933buildConstraint_(C,Agenda,NewAgenda) :- 934 debug_clpBNR_('Add ~p',{C}), 935 % need to catch errors from ground expression evaluation 936 catch(build_(C, 1, boolean, Agenda, NewAgenda),_Err,fail), !. 937buildConstraint_(C,_Agenda,_NewAgenda) :- 938 debug_clpBNR_('{} failure due to bad or inconsistent constraint: ~p',{C}), 939 fail. 940 941:- include(clpBNR/ia_simplify). % simplifies constraints to a hopefully more efficient equivalent 942 943% 944% build a node from an expression 945% 946build_(Int, Int, VarType, Agenda, NewAgenda) :- 947 interval(Int), !, % existing interval object 948 applyType_(VarType, Int, Agenda, NewAgenda). % coerces exsiting intervals to required type 949build_(Var, Var, VarType, Agenda, Agenda) :- % implicit interval creation. 950 var(Var), !, 951 new_interval_(Var,VarType). 952build_(Num, Int, VarType, Agenda, Agenda) :- % floating point constant, may not be precise 953 float(Num), !, 954 ( 955 float_class(Num,infinite) 956 -> Int=Num % infinities are point values 957 ; int_decl_(VarType,(Num,Num),Int) % may be fuzzed, so not a point 958 ). 959build_(::(L,H), Int, VarType, Agenda, Agenda) :- % hidden :: feature: interval bounds literal (without fuzzing) 960 number(L), number(H), !, 961 C is cmpr(L,H), % compare bounds 962 (C == 0 963 -> (rational(L) -> Int=L ; Int=H) % point value, if either bound rational (precise), use it 964 ; C == -1, % necessary condition: L < H 965 once(VarType = real ; true), % if undefined Type, use 'real' 966 put_attr(Int, clpBNR, interval(VarType, (L,H), _NL, [])) % create clpBNR attribute 967 ). 968build_(Num, Int, VarType, Agenda, Agenda) :- % pre-compile constants pi and e 969 (Num == pi ; Num == e), !, 970 int_decl_(VarType,(Num,Num),Int). 971build_(Exp, Num, _, Agenda, Agenda) :- % pre-compile any ground precise Exp 972 ground(Exp), 973 safe_(Exp), % safe to evaluate using is/2 974 Num is Exp, 975 preciseBnd(Num), % precise result 976 !. 977build_(Exp, Z, _, Agenda, NewAgenda) :- % deconstruct to primitives 978 Exp =.. [F|Args], 979 fmap_(F,Op,[Z|Args],ArgsR,Types), !, % supported arithmetic op 980 build_args_(ArgsR,Objs,Types,Agenda,ObjAgenda), 981 newNode_(Op,Objs,ObjAgenda,NewAgenda). 982build_(Exp, Z, _, Agenda, NewAgenda) :- % user defined 983 Exp =.. [Prim|Args], 984 chk_primitive_(Prim), 985 build_args_([Z|Args],Objs,_Types,Agenda,ObjAgenda), 986 newNode_(user(Prim),Objs,ObjAgenda,NewAgenda). 987 988build_args_([],[],_,Agenda,Agenda). 989build_args_([Arg|ArgsR],[Obj|Objs],[Type|Types],Agenda,NewAgenda) :- 990 (var(Type) -> Type=real ; true), % var if user defined 991 build_(Arg,Obj,Type,Agenda,NxtAgenda), 992 build_args_(ArgsR,Objs,Types,NxtAgenda,NewAgenda). 993 994chk_primitive_(Prim) :- % wraps safe usage of unsafe current_predicate/2 995 UsrHead =..[Prim,'$op',_,_,_], 996 current_predicate(_,clpBNR:UsrHead). 997 998sandbox:safe_primitive(clpBNR:chk_primitive_(_Prim)). 999 1000% to invoke user defined primitive 1001call_user_primitive(Prim, P, InArgs, OutArgs) :- % wraps unsafe meta call/N 1002 call(clpBNR:Prim, '$op', InArgs, OutArgs, P). 1003 1004% really unsafe, but in a pengine a user can't define any predicates in another module, so this is safe 1005sandbox:safe_meta(clpBNR:call_user_primitive(_Prim, _P, _InArgs, _OutArgs), []). 1006 1007% only called when argument is ground 1008safe_(E) :- atomic(E), !. % all atomics, including [] - allows possibility of user defined arithmetic types 1009safe_([A|As]) :- !, 1010 safe_(A), 1011 safe_(As). 1012safe_(_ xor _) :- !, % clpBNR xor incompatible with `is` xor 1013 fail. 1014safe_(integer(_)) :- !, % clpBNR integer incompatible with `is` integer 1015 fail. 1016safe_(_/Z) :- 0.0 is abs(Z), !, % division by 0.0 (or -0.0) requires fuzzed 0.0 1017 fail. 1018safe_(F) :- 1019 current_arithmetic_function(F), % evaluable by is/2 1020 F =.. [_Op|Args], 1021 safe_(Args). 1022 1023% a constraint must evaluate to a boolean 1024constraint_(C) :- nonvar(C), C =..[Op|_], fmap_(Op,_,_,_,[boolean|_]), !. 1025 1026% map constraint operator to primitive/arity/types 1027fmap_(+, add, ZXY, ZXY, [real,real,real]). 1028fmap_(-, add, [Z,X,Y], [X,Z,Y], [real,real,real]). % note subtract before minus 1029fmap_(*, mul, ZXY, ZXY, [real,real,real]). 1030fmap_(/, mul, [Z,X,Y], [X,Z,Y], [real,real,real]). 1031fmap_(**, pow, ZXY, ZXY, [real,real,real]). 1032fmap_(min, min, ZXY, ZXY, [real,real,real]). 1033fmap_(max, max, ZXY, ZXY, [real,real,real]). 1034fmap_(==, eq, ZXY, ZXY, [boolean,real,real]). % strict equality 1035fmap_(=:=, eq, ZXY, ZXY, [boolean,real,real]). % Prolog compatible, strict equality 1036fmap_(is, eq, ZXY, ZXY, [boolean,real,real]). 1037fmap_(<>, ne, ZXY, ZXY, [boolean,integer,integer]). 1038fmap_(=\=, ne, ZXY, ZXY, [boolean,integer,integer]). % Prolog compatible 1039fmap_(=<, le, ZXY, ZXY, [boolean,real,real]). 1040fmap_(>=, le, [Z,X,Y], [Z,Y,X], [boolean,real,real]). 1041fmap_(<, lt, ZXY, ZXY, [boolean,real,real]). 1042fmap_(>, lt, [Z,X,Y], [Z,Y,X], [boolean,real,real]). 1043fmap_(<=, in, ZXY, ZXY, [boolean,real,real]). % inclusion/subinterval 1044 1045fmap_(and, and, ZXY, ZXY, [boolean,boolean,boolean]). 1046fmap_(',', and, ZXY, ZXY, [boolean,boolean,boolean]). % for usability 1047fmap_(or, or, ZXY, ZXY, [boolean,boolean,boolean]). 1048fmap_(nand, nand, ZXY, ZXY, [boolean,boolean,boolean]). 1049fmap_(nor, nor, ZXY, ZXY, [boolean,boolean,boolean]). 1050fmap_(xor, xor, ZXY, ZXY, [boolean,boolean,boolean]). 1051fmap_(->, imB, ZXY, ZXY, [boolean,boolean,boolean]). 1052 1053fmap_(sqrt, sqrt, ZX, ZX, [real,real]). % pos. root version vs. **(1/2) 1054fmap_(-, minus, ZX, ZX, [real,real]). 1055fmap_(~, not, ZX, ZX, [boolean,boolean]). 1056fmap_(integer,int, ZX, ZX, [boolean,real]). 1057fmap_(exp, exp, ZX, ZX, [real,real]). 1058fmap_(log, exp, [Z,X], [X,Z], [real,real]). 1059fmap_(abs, abs, ZX, ZX, [real,real]). 1060fmap_(sin, sin, ZX, ZX, [real,real]). 1061fmap_(asin, sin, [Z,X], [X,Z], [real,real]). 1062fmap_(cos, cos, ZX, ZX, [real,real]). 1063fmap_(acos, cos, [Z,X], [X,Z], [real,real]). 1064fmap_(tan, tan, ZX, ZX, [real,real]). 1065fmap_(atan, tan, [Z,X], [X,Z], [real,real]). 1066 1067% reverse map node info to a facsimile of the original constraint 1068map_constraint_op_(integral,$(V),integral(V)) :- !. 1069map_constraint_op_(user(Func),Args,C) :- !, 1070 remap_(Func,Args,C). 1071map_constraint_op_(Op,Args,C) :- 1072 fmap_(COp,Op,_,_,_), 1073 remap_(COp,Args,C), 1074 !. 1075 1076remap_(Op,$(Z,X,Y),C) :- constraint_(Op), Z==1, !, % simplification for binary constraints 1077 C=..[Op,X,Y]. 1078remap_(Op,$(Z,X),C) :- constraint_(Op), Z==1, !, % simplification for unary constraints (~) 1079 C=..[Op,X]. 1080remap_(Op,$(Z,X,Y),Z==C) :- !, 1081 C=..[Op,X,Y]. 1082remap_(Op,$(Z,X),Z==C) :- 1083 C=..[Op,X]. 1084 1085% 1086% Node constructor 1087% 1088% First clause is an optimization for E1 == E2 1089% In that case just unify E1 and E2; heavy lifting done by attr_unify_hook. 1090newNode_(eq, [Z,X,Y], Agenda, Agenda) :- Z==1, !, X=Y. 1091newNode_(Op, Objs, Agenda, NewAgenda) :- 1092 Args =.. [$|Objs], % store arguments as $/N where N=1..3 1093 NewNode = node(Op, _P, 0, Args), % L=0 1094 addNode_(Objs,NewNode), 1095 % increment count of added nodes, will be decremented on backtracking/failure 1096 g_incb('clpBNR:node_count'), 1097 linkNode_(Agenda, NewNode, NewAgenda). 1098 1099addNode_([],_Node). 1100addNode_([Arg|Args],Node) :- 1101 (interval_object(Arg, _Type, _Val, Nodelist) -> newmember(Nodelist, Node) ; true), 1102 addNode_(Args,Node). 1103 1104sandbox:safe_global_variable('clpBNR:node_count'). 1105 1106clpStatistics :- 1107 g_assign('clpBNR:node_count',0), % reset/initialize node count to 0 1108 fail. % backtrack to reset other statistics. 1109 1110clpStatistic(node_count(C)) :- 1111 g_read('clpBNR:node_count',C). 1112 1113% extend list with N 1114newmember([X|Xs],N) :- 1115 (nonvar(X) 1116 -> newmember(Xs,N) % not end of (indefinite) list. 1117 ; X = N % end of list, insert N and new tail Xs 1118 ). 1119 1120% 1121% Process Agenda to narrow intervals (fixed point iteration) 1122% 1123stable_([]/[]) :- !. % small optimization when agenda is empty 1124stable_(Agenda) :- 1125 current_prolog_flag(clpBNR_iteration_limit,Ops), % budget for current operation 1126 stableLoop_(Agenda,Ops), 1127 !. % achieved stable state with empty Agenda -> commit. 1128 1129stableLoop_([]/[], OpsLeft) :- !, % terminate successfully when agenda comes to an end 1130 g_read('clpBNR:iterations',Cur), % maintain "low" water mark (can be negative) 1131 (OpsLeft<Cur -> g_assign('clpBNR:iterations',OpsLeft) ; true), 1132 (OpsLeft<0 -> E is -OpsLeft, debug_clpBNR_('Iteration throttle limit exceeded by ~w ops.',E) ; true). 1133stableLoop_([Node|Agenda]/T, OpsLeft) :- 1134 Node = node(Op,P,_,Args), % if node on queue ignore link bit (was: Node = node(Op,P,1,Args)) 1135 doNode_(Args, Op, P, OpsLeft, DoOver, Agenda/T, NxtAgenda), % undoable on backtrack 1136 nb_setarg(3,Node,0), % reset linked bit 1137 % if doNode_ requested DoOver and above Ops threshold, put Node back at the end of the queue 1138 (atom(DoOver), OpsLeft > 0 -> linkNode_(NxtAgenda,Node,NewAgenda) ; NewAgenda = NxtAgenda), 1139 RemainingOps is OpsLeft-1, % decrement OpsLeft (can go negative) 1140 stableLoop_(NewAgenda,RemainingOps). 1141 1142% support for max_iterations statistic 1143sandbox:safe_global_variable('clpBNR:iterations'). 1144 1145clpStatistics :- 1146 current_prolog_flag(clpBNR_iteration_limit,L), 1147 g_assign('clpBNR:iterations',L), % set "low" water mark to upper limit 1148 fail. % backtrack to reset other statistics. 1149 1150clpStatistic(max_iterations(O/L)) :- 1151 g_read('clpBNR:iterations',Ops), 1152 current_prolog_flag(clpBNR_iteration_limit,L), 1153 O is L-Ops. % convert "low" water mark to high water mark 1154 1155% 1156% doNode_/7 : Evaluate a node and add new nodes to end of queue. `evalNode` primitives can 1157% fail, resulting in eventual failure of `stable_`, i.e., inconsistent constraint set. 1158% 1159% Note: primitives operate on interval values (sets) only, unaware of common arguments, 1160% so additional consistency checks required on update. 1161% 1162doNode_($(ZArg,XArg,YArg), Op, P, OpsLeft, DoOver, Agenda, NewAgenda) :- % Arity 3 Op 1163 (var(P) % check persistent bit 1164 -> getValue(ZArg,ZVal), 1165 getValue(XArg,XVal), 1166 getValue(YArg,YVal), 1167 evalNode(Op, P, $(ZVal,XVal,YVal), $(NZVal,NXVal,NYVal)), % can fail 1168 % enforce consistency for common arguments by intersecting and redoing as required. 1169 (var(ZArg) % if Z == X and/or Y 1170 -> (ZArg==XArg -> consistent_value_(NZVal,NXVal,NZ1,DoOver) ; NZ1 = NZVal), 1171 (ZArg==YArg -> consistent_value_(NZ1, NYVal,NZ2,DoOver) ; NZ2 = NZ1), 1172 updateValue_(ZVal, NZ2, ZArg, OpsLeft, Agenda, AgendaZ) 1173 ; AgendaZ = Agenda 1174 ), 1175 (var(XArg), XArg==YArg % if X==Y 1176 -> consistent_value_(NXVal,NYVal,NVal,DoOver), 1177 updateValue_(XVal, NVal, XArg, OpsLeft, AgendaZ, NewAgenda) % only one update needed 1178 ; updateValue_(XVal, NXVal, XArg, OpsLeft, AgendaZ, AgendaZX), 1179 updateValue_(YVal, NYVal, YArg, OpsLeft, AgendaZX, NewAgenda) 1180 ) 1181 ; % P = p, trim persistent nodes from all arguments 1182 trim_op_(ZArg), trim_op_(XArg), trim_op_(YArg), 1183 NewAgenda = Agenda 1184 ). 1185 1186doNode_($(ZArg,XArg), Op, P, OpsLeft, DoOver, Agenda, NewAgenda) :- % Arity 2 Op 1187 (var(P) % check persistent bit 1188 -> getValue(ZArg,ZVal), 1189 getValue(XArg,XVal), 1190 evalNode(Op, P, $(ZVal,XVal), $(NZVal,NXVal)), % can fail 1191 % enforce consistency for common arguments by intersecting and redoing as required. 1192 (var(ZArg), ZArg==XArg % if Z==X 1193 -> consistent_value_(NZVal,NXVal,NVal,DoOver), % consistent value, DoOver if needed 1194 updateValue_(ZVal, NVal, ZArg, OpsLeft, Agenda, NewAgenda) % only one update needed 1195 ; updateValue_(ZVal, NZVal, ZArg, OpsLeft, Agenda, AgendaZ), 1196 updateValue_(XVal, NXVal, XArg, OpsLeft, AgendaZ, NewAgenda) 1197 ) 1198 ; % P = p, trim persistent nodes from all arguments 1199 trim_op_(ZArg), trim_op_(XArg), 1200 NewAgenda = Agenda 1201 ). 1202 1203doNode_($(Arg), Op, P, _OpsLeft, _, Agenda, NewAgenda) :- % Arity 1 Op 1204 (var(P) % check persistent bit 1205 -> getValue(Arg,Val), 1206 evalNode(Op, P, $(Val), $(NVal)), % can fail causing stable_ to fail => backtracking 1207 updateValue_(Val, NVal, Arg, 1, Agenda,NewAgenda) % always update value regardless of OpsLeft limiter 1208 ; % P = p, trim persistent nodes from argument 1209 trim_op_(Arg), 1210 NewAgenda = Agenda 1211 ). 1212 1213consistent_value_(Val,Val,Val,_) :- !. % same value 1214consistent_value_(Val1,Val2,Val,true) :- ^(Val1,Val2,Val). % different values, intersect 1215 1216% 1217% remove any persistent nodes from Arg 1218% called whenever a persistent node is encountered in FP iteration 1219% 1220trim_op_(Arg) :- 1221 ( get_attr(Arg, clpBNR, Def) % an interval ? 1222 -> arg(3,Def,NList), %%Def = interval(_, _, NList, _), 1223 trim_persistent_(NList,TrimList), 1224 % if trimmed list empty, set to a new unshared var to avoid cycles(?) on backtracking 1225 (var(TrimList) -> setarg(3,Def,_) ; setarg(3,Def,TrimList)) % write trimmed node list 1226 ; true % assumed to be a number, nothing to trim 1227 ). 1228 1229trim_persistent_(T,T) :- var(T), !. % end of indefinite node list 1230trim_persistent_([node(_,P,_,_)|Ns],TNs) :- nonvar(P), !, trim_persistent_(Ns,TNs). 1231trim_persistent_([N|Ns],[N|TNs]) :- trim_persistent_(Ns,TNs). 1232 1233% 1234% Any changes in interval values should come through here. 1235% Note: This captures all updated state for undoing on backtracking 1236% 1237updateValue_(Old, Old, _, _, Agenda, Agenda) :- !. % no change in value (constant?) 1238 1239updateValue_(Old, New, Int, OpsLeft, Agenda, NewAgenda) :- % set interval value to New 1240 (OpsLeft>0 -> true ; propagate_if_(Old, New)), !, % if OpsLeft >0 or narrowing sufficent 1241 putValue_(New, Int, Nodelist), % update value (may fail) 1242 linkNodeList_(Nodelist, Agenda, NewAgenda). % then propagate change 1243 1244updateValue_(_, _, _, _, Agenda, Agenda). % otherwise just continue with Agenda 1245 1246% propgate if sufficient narrowing (> 10%) 1247propagate_if_((OL,OH), (NL,NH)) :- (NH-NL)/(OH-OL) < 0.9. % any overflow in calculation will propagate 1248 1249linkNodeList_([X|Xs], List, NewList) :- 1250 (var(X) 1251 -> List = NewList % end of indefinite list 1252 ; (arg(3,X,Linked), Linked == 1 % test linked flag (only SWIP VM codes) 1253 -> linkNodeList_(Xs, List, NewList) % add rest of the node list 1254 ; linkNode_(List, X, NextList), % not linked add it to list 1255 linkNodeList_(Xs, NextList, NewList) % add rest of the node list 1256 ) 1257 ). 1258 1259linkNode_(List/[X|NextTail], X, List/NextTail) :- % add to list 1260 setarg(3,X,1). % set linked bit 1261 1262:- 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.
/1269watch(Int,Action) :- 1270 atom(Action), 1271 current_module(clpBNR), % this just marks watch/2 as unsafe regarding body 1272 get_interval_flags_(Int,Flags), !, 1273 remove_(Flags,watch(_),Flags1), 1274 (Action = none -> Flags2=Flags1 ; Flags2=[watch(Action)|Flags1]), 1275 set_interval_flags_(Int,Flags2). 1276watch(Ints,Action) :- 1277 list(Ints), 1278 watch_list_(Ints,Action). 1279 1280remove_([],_,[]). 1281remove_([X|Xs],X,Xs) :- !. 1282remove_([X|Xs],X,[X|Ys]) :- 1283 remove_(Xs,X,Ys). 1284 1285watch_list_([],_Action). 1286watch_list_([Int|Ints],Action) :- 1287 watch(Int,Action), 1288 watch_list_(Ints,Action). 1289 1290% check if watch enabled on this interval 1291check_monitor_(Int, Update, interval(_Type,_Val,_Nodelist,Flags)) :- 1292 (memberchk(watch(Action), Flags) 1293 -> once(monitor_action_(Action, Update, Int)) 1294 ; true 1295 ). 1296 1297% 1298% for monitoring changes - all actions defined here 1299% 1300monitor_action_(trace, Update, Int) :- !, % log changes to console and enter debugger 1301 monitor_action_(log, Update, Int), 1302 trace. 1303monitor_action_(log, Update, Int) :- var(Update), !, % log interval unify 1304 debug_clpBNR_('Unify ~p with ~p',[Int,Update]). 1305monitor_action_(log, Update, Int) :- number(Update), !, % log interval unify with number 1306 domain(Int,Dom), 1307 debug_clpBNR_('Unify _?{~p} with ~p',[Dom,Update]). 1308monitor_action_(log, integer, Int) :- !, % change type to integer 1309 debug_clpBNR_('Set type of ~p to ~p',[Int,integer]). 1310monitor_action_(log, Val, Int) :- !, % narrow range 1311 debug_clpBNR_('Set value of ~p to (~p)',[Int,Val]). 1312monitor_action_(_, _, _). % default to noop (as in 'none') 1313 1314sandbox:safe_primitive(clpBNR:watch(_Int,Action)) :- % watch(X,trace) is unsafe. 1315 Action \= trace. 1316% only safe because watch(X,trace) is unsafe. 1317sandbox: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.
/1324% 1325% tracing doNode_ support - using wrap_predicate(:Head, +Name, -Wrapped, +Body)/4 1326% trace_clpBNR/1 is unsafe (wrapping is global) 1327% 1328:- use_module(library(prolog_wrap)). 1329 1330trace_clpBNR(Bool) :- % query or already in defined state 1331 ( current_predicate_wrapper(clpBNR:doNode_(_Args, _Op, _P, _OpsLeft, _DoOver, _Agenda, _NewAgenda), 1332 'clpBNR:doNode_', _Wrapped, _Body) 1333 -> Bool = true ; Bool = false 1334 ), 1335 !. 1336trace_clpBNR(true) :- % turn on wrapper 1337 wrap_predicate(clpBNR:doNode_(Args, Op, _P, _OpsLeft, _DoOver, _Agenda, _NewAgenda), 1338 'clpBNR:doNode_', 1339 Wrapped, 1340 doNode_wrap_(Wrapped, Args,Op)). 1341trace_clpBNR(false) :- % turn off wrapper 1342 unwrap_predicate(clpBNR:doNode_/7, %(Args, Op, P, OpsLeft, DoOver, Agenda, NewAgenda), 1343 'clpBNR:doNode_'). 1344 1345doNode_wrap_(Wrapped, Args,Op) :- 1346 map_constraint_op_(Op,Args,C), 1347 , % execute wrapped doNode_ 1348 debug_clpBNR_("~p.",C). % print trace message , fail messages from evalNode_, attr_unify_hook 1349 1350% 1351% Get all defined statistics 1352% 1353clpStatistics(Ss) :- findall(S, clpStatistic(S), Ss). 1354 1355% end of reset chain succeeds. 1356clpStatistics. 1357 1358% 1359% module initialization 1360% 1361init_clpBNR :- 1362 restore_optimize_flags_, 1363 print_message(informational, clpBNR(versionInfo)), 1364 print_message(informational, clpBNR(arithmeticFlags)). % cautionary, set on first use 1365 1366:- if(false). % test code used to test sandbox worthiness of hooks 1367check_hooks_safety :- % Note: calls must have no side effects 1368 ignore(attr_portray_hook([],_)), % fails 1369 ignore(user:exception(undefined_global_variable,'clpBNR:thread_init_done',[])), % fails 1370% ignore(term_html:portray('$clpBNR...'(_),_,_,_)), % fails 1371 ignore(user:portray('$clpBNR...'(_))). % fails 1372:- endif. 1373 1374:- multifile prolog:message//1. 1375 1376prologmessage(clpBNR(versionInfo)) --> 1377 { version(Version) }, 1378 [ '*** clpBNR v~w ***.'-[Version] ]. 1379 1380prologmessage(clpBNR(arithmeticFlags)) --> 1381 [ ' Arithmetic global flags will be set to prefer rationals and IEEE continuation values.'-[] ]. 1382 1383:- 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