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*/
   94
   95version("0.11.5").
   96
   97% debug feature control and messaging
   98:- if(exists_source(swish(lib/swish_debug))).   99	:- create_prolog_flag(clpBNR_swish, true, [access(read_only)]).  100	:- use_module(swish(lib/swish_debug)).  101:- else.  102	:- use_module(library(debug)).  103:- endif.  104
  105:- use_module(library(prolog_versions)).  % SWIP dependency enforcement
  106
  107:- require_prolog_version('9.1.22',       % properly exported arithmetic functions
  108	          [ rational   % require rational number support, implies bounded=false
  109	          ]).  110
  111%
  112% Optimize arithmetic, but not debug. Only called via directives.
  113%
  114set_optimize_flags_ :-      % called at start of load
  115	set_prolog_flag(optimise,true),              % scoped to file/module
  116	current_prolog_flag(optimise_debug,ODflag),  % save; restore in initialization
  117	nb_linkval('$optimise_debug_save',ODflag),
  118	set_prolog_flag(optimise_debug,false).       % so debug/3, debugging/1 don't get "optimized"
  119
  120restore_optimize_flags_ :-  % called at module initialization (after load)
  121	nb_getval('$optimise_debug_save',ODflag), nb_delete('$optimise_debug_save'),
  122	set_prolog_flag(optimise_debug,ODflag).
  123
  124:- set_optimize_flags_.  125
  126% local debug and trace message support
  127debug_clpBNR_(FString,Args) :- debug(clpBNR,FString,Args).
  128
  129% sandboxing for SWISH
  130:- multifile(sandbox:safe_prolog_flag/1).  131:- multifile(sandbox:safe_global_variable/1).  132:- multifile(sandbox:safe_primitive/1).  133:- multifile(sandbox:safe_meta/2).  134
  135current_node_(Node) :-  % look back to find current Op being executed for debug messages
  136	prolog_current_frame(Frame),  % this is a little grungy, but necessary to get intervals
  137	prolog_frame_attribute(Frame,parent_goal,doNode_(Args,Op,_,_,_,_,_)),
  138	map_constraint_op_(Op,Args,Node),
  139	!.
  140
  141sandbox:safe_primitive(clpBNR:current_node_(_Node)). 
  142
  143%
  144% statistics
  145%
  146
  147% assign,increment/read global counter (assumed to be ground value so use _linkval)
  148g_assign(G,V)  :- nb_linkval(G,V).
  149g_inc(G)       :- nb_getval(G,N), N1 is N+1, nb_linkval(G,N1).
  150g_incb(G)      :- nb_getval(G,N), N1 is N+1, b_setval(G,N1).    % undone on backtrack
  151g_read(G,V)    :- nb_getval(G,V).
  152
  153sandbox:safe_global_variable('clpBNR:thread_init_done').
  154sandbox:safe_global_variable('clpBNR:userTime').
  155sandbox:safe_global_variable('clpBNR:inferences').
  156sandbox:safe_global_variable('clpBNR:gc_time').
  157
  158%  
  159% Global var statistics are per thread and therefore must be "lazily" initialized
  160% Also ensures that thread copies of flags are properly set.
  161% This exception handler will be invoked the first time 'clpBNR:thread_init_done' is read
  162% Public predicates ::, {}, clpStatistics/1 and range read this global before proceeding. 
  163%
  164user:exception(undefined_global_variable,'clpBNR:thread_init_done',retry) :- !,
  165	set_prolog_flags,     % initialize/check environment flags  
  166	clpStatistics,        % defines remaining globals 
  167	g_assign('clpBNR:thread_init_done',1).
  168
  169%
  170% Define custom clpBNR flags when module loaded
  171%
  172:- create_prolog_flag(clpBNR_iteration_limit,3000,[type(integer),keep(true)]).  173:- create_prolog_flag(clpBNR_default_precision,6,[type(integer),keep(true)]).  174:- create_prolog_flag(clpBNR_verbose,false,[type(boolean),keep(true)]).  175
  176sandbox:safe_prolog_flag(clpBNR_iteration_limit,_).
  177sandbox:safe_prolog_flag(clpBNR_default_precision,_).
  178sandbox:safe_prolog_flag(clpBNR_verbose,_).
  179%
  180% Set public flags (see module/thread initialization)
  181%
  182set_prolog_flags :-
  183	set_prolog_flag(prefer_rationals, true),           % enable rational arithmetic
  184	(current_prolog_flag(max_rational_size,_)
  185	 -> true                                           % already defined, leave as is
  186	 ;  set_prolog_flag(max_rational_size, 16)         % rational size in bytes before ..
  187	),
  188	set_prolog_flag(max_rational_size_action, float),  % conversion to float
  189
  190	set_prolog_flag(float_overflow,infinity),          % enable IEEE continuation values
  191	set_prolog_flag(float_zero_div,infinity),
  192	set_prolog_flag(float_undefined,nan),
  193	set_prolog_flag(write_attributes,portray).         % thread-local, init to 'portray'
  194
  195:- discontiguous clpBNR:clpStatistics/0, clpBNR:clpStatistic/1.  196
  197clpStatistics :-
  198	% garbage_collect,  % ? do gc before time snapshots
  199	statistics(cputime,T), g_assign('clpBNR:userTime',T),   % thread based
  200	statistics(inferences,I), g_assign('clpBNR:inferences',I),
  201	statistics(garbage_collection,[_,_,G,_]), g_assign('clpBNR:gc_time',G),
  202	fail.  % backtrack to reset other statistics.
  203
  204clpStatistic(_) :- g_read('clpBNR:thread_init_done',0).  % ensures per-thread initialization then fails
  205
  206clpStatistic(userTime(T)) :- statistics(cputime,T1), g_read('clpBNR:userTime',T0), T is T1-T0.
  207
  208clpStatistic(gcTime(G)) :- statistics(garbage_collection,[_,_,G1,_]), g_read('clpBNR:gc_time',G0), G is (G1-G0)/1000.0.
  209
  210clpStatistic(globalStack(U/T)) :- statistics(globalused,U), statistics(global,T).
  211
  212clpStatistic(trailStack(U/T)) :- statistics(trailused,U), statistics(trail,T).
  213
  214clpStatistic(localStack(U/T)) :- statistics(localused,U), statistics(local,T).
  215
  216clpStatistic(inferences(I)) :- statistics(inferences,I1), g_read('clpBNR:inferences',I0), I is I1-I0.
  217
  218%
  219% list filter for compatibility
  220%	Note: not equivalent to is_list/1 but executes in O(1) time.
  221%
  222list(X) :- compound(X) ->  X=[_|_] ; X==[].
  223
  224:- include(clpBNR/ia_primitives).  % interval arithmetic relations via evalNode/4.
  225
  226%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  227%
  228%  SWI-Prolog implementation of IA
  229%
  230%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  231%
  232% Intervals are constrained (attributed) variables.
  233%
  234% Current interval bounds updates via setarg(Val) which is backtrackable
  235%
  236%  interval(Int)  - filter
  237%
  238interval(Int) :- get_attr(Int, clpBNR, _).
  239
  240%
  241%  interval_degree/2   % N = number of constraints on interval(Int), 0 for number(Int)
  242%
  243interval_degree(X, N) :- 
  244	interval_object(X, _, _, Nodelist)
  245	-> system:'$skip_list'(N, Nodelist, _) % optimized for SWIP, handles indefinite lists
  246	;  number(X), N = 0.                   % number -> no constraints ; fail
  247
  248% internal abstraction
  249interval_object(Int, Type, Val, Nodelist) :-
  250	get_attr(Int, clpBNR, interval(Type, Val, Nodelist, _)).
  251
  252% flags (list)  abstraction
  253get_interval_flags_(Int, Flags) :-
  254	get_attr(Int, clpBNR, interval(_, _, _, Flags)).
  255
  256set_interval_flags_(Int, Flags) :-  % flags assumed to be ground so no copy required
  257	interval_object(Int, Type, Val, Nodelist),
  258	put_attr(Int, clpBNR, interval(Type, Val, Nodelist, Flags)).
  259
  260reset_interval_nodelist_(Int) :-
  261	get_attr(Int, clpBNR, Def) -> setarg(3,Def,_) ; true.
  262
  263%
  264% Interval value constants
  265%
  266universal_interval((-1.0Inf,1.0Inf)).
  267
  268empty_interval((1.0Inf,-1.0Inf)).
  269
  270% Finite intervals - 64 bit IEEE reals, 
  271finite_interval(real,    (-1.0e+16,1.0e+16)).
  272finite_interval(integer, (L,H)) :-

current_prolog_flag(bounded,false), % required - unbounded integers

  274	current_prolog_flag(min_tagged_integer,L),
  275	current_prolog_flag(max_tagged_integer,H)

  275.
  276finite_interval(boolean, (0,1)).
  277
  278% legal integer bounds values
  279integerBnd(1.0Inf).
  280integerBnd(-1.0Inf).
  281integerBnd(B) :- integer(B).
  282
  283% precise bounds values
  284preciseBnd(1.0Inf).
  285preciseBnd(-1.0Inf).
  286preciseBnd(B) :- rational(B).
  287
  288%
  289%  non-backtracking set bounds for use with branch and bound
  290%
  291nb_setbounds(Int, [L,U]) :- 
  292	get_attr(Int, clpBNR, Def),
  293	arg(2, Def, Val),             % WAM code
  294	^(Val,(L,U),NewVal),          % new range is intersection (from ia_primitives)
  295	nb_setarg(2, Def, NewVal).
  296
  297%
  298% get current value
  299%
  300getValue(Int, Val) :- 
  301	number(Int)
  302	 -> Val=(Int,Int)                                   % numeric point value
  303	 ;  get_attr(Int, clpBNR, interval(_, Val, _, _)).  % interval, optimized for SWIP
  304
  305%
  306% set interval value (assumes bounds check has already been done)
  307% Note: putValue_ cannot modify the new value since the narrowing op is still in the Agenda
  308%	and may not be run again. Insert `integral` op for integer bounds adjustment.
  309%
  310putValue_(New, Int, NodeList) :- 
  311	get_attr(Int, clpBNR, Def)             % still an interval ?
  312	 -> (debugging(clpBNR) -> check_monitor_(Int, New, Def) ; true),
  313	    Def = interval(Type,_,Nodes,_),    % get Type and Nodes for queueing
  314	    New = (L,H),
  315	    ( 0 is cmpr(L,H)                   % point interval ->
  316	     -> setarg(3,Def,_NL),             % clear node list (so nothing done in unify)
  317	        pointValue_(L,H,Int),          % unify, avoiding -0.0 and preferring rational (invokes hook)
  318	        NodeList = Nodes               % still add Nodes to Agenda
  319	     ;  setarg(2,Def,New),             % update value in interval (backtrackable)
  320	        % if integer has non-integral bounds, schedule `integral` op to fix it
  321	        (Type == integer
  322	         -> ( integerBnd(L), integerBnd(H) -> NodeList = Nodes ; NodeList = [node(integral,_,0,$(Int))|_] )
  323	         ;  NodeList = Nodes
  324	        )
  325	    )
  326	 ; true.  % no longer an interval, NodeList is empty
  327
  328pointValue_(-0.0,_,0.0) :-!.
  329pointValue_(L,H,Int) :- (rational(L) -> Int = L ; Int = H).
  330
  331%
  332%  range(Int, Bounds) for compatibility 
  333%
  334range(Int, [L,H]) :- getValue(Int, (IL,IH)), !,  % existing interval or number =>  number(IL)
  335	(var(L) -> L=IL ; non_empty(L,IL)),          % range check (L=<IL,IH=<H), no narrowing
  336	(var(H) -> H=IH ; non_empty(IH,H)).
  337range(Int, [L,H]) :- var(Int),  % for other var(Int), declare it to a real interval with specified bounds
  338	Int::real(L,H).
  339
  340%
  341%  domain(Int, Dom) for interval(Int)
  342%
  343domain(Int, Dom) :-
  344	interval_object(Int, Type, Val, _),
  345	interval_domain_(Type, Val, Dom).
  346
  347interval_domain_(integer,(0,1),boolean) :- !.  % integer(0,1) is synonymous with boolean
  348interval_domain_(T,(L,H),Dom) :- Dom=..[T,L,H].
  349
  350:- use_module(library(arithmetic), []).   % extended arithmetic functions
  351%
  352%  delta(Int, Wid) width/span of an interval or numeric value, can be infinite
  353%
  354:- arithmetic_function(delta/1).  355
  356delta(Int, Wid) :-
  357	getValue(Int,(L,H)),
  358	Wid is roundtoward(H-L,to_positive).
  359
  360%
  361%  midpoint(Int, Wid) midpoint of an interval or numeric value
  362% based on:
  363%	Frédéric Goualard. How do you compute the midpoint of an interval?. 
  364%	ACM Transactions on Mathematical Software, Association for Computing Machinery, 2014,
  365%	40 (2), 10.1145/2493882. hal-00576641v1
  366% Exception, single infinite bound treated as largest finite FP value
  367%
  368:- arithmetic_function(midpoint/1).  369
  370midpoint(Int, Mid) :-
  371	getValue(Int,(L,H)),
  372	midpoint_(L,H,Mid).
  373
  374midpoint_(L,H,M)       :- L =:= -H, !, M=0.              % symmetric including (-inf,inf)
  375midpoint_(-1.0Inf,H,M) :- !, M is nexttoward(-1.0Inf,0)/2 + H/2.
  376midpoint_(L,1.0Inf,M)  :- !, M is L/2 + nexttoward(1.0Inf,0)/2.
  377midpoint_(L,H,M)       :- M1 is L/2 + H/2, M=M1.        % general case
  378
  379%
  380% median(Int,Med) from CLP(RI)
  381% Med = 0 if Int contains 0, else a number which divides Int into equal
  382% numbers of FP values. Med is always a float
  383%
  384:- arithmetic_function(median/1).  385
  386median(Int, Med) :-
  387	getValue(Int,(L,H)),
  388	median_bound_(lo,L,FL),
  389	median_bound_(hi,H,FH),
  390	median_(FL,FH,Med), !.
  391	
  392median_bound_(lo,B,FB) :- B=:=0, FB is nexttoward(B,1.0).
  393median_bound_(lo,-1.0Inf,FB) :- FB is nexttoward(-1.0Inf,0.0).
  394median_bound_(lo,B,FB) :- FB is roundtoward(float(B), to_negative).
  395
  396median_bound_(hi,B,FB) :- B=:=0, !, FB is nexttoward(B,-1.0).
  397median_bound_(hi,1.0Inf,FB) :- FB is nexttoward(1.0Inf,0.0).
  398median_bound_(hi,B,FB) :- FB is roundtoward(float(B), to_positive).
  399
  400median_(B,B,B).                          % point interval
  401median_(L,H,0.0) :- L < 0.0, H > 0.0.    % contains 0: handles (-inf,inf)
  402median_(L,H,M)   :- M is copysign(sqrt(abs(L))*sqrt(abs(H)),L).      % L and H have same sign
  403
  404%
  405%  lower_bound and upper_bound
  406%
  407lower_bound(Int) :-
  408	getValue(Int,(L,_H)),
  409	Int=L.
  410
  411upper_bound(Int) :-
  412	getValue(Int,(_L,H)),
  413	Int=H.
  414
  415%
  416% Interval declarations
  417%
  418Rs::Dom :- list(Rs),!,                    % list of vars
  419	intervals_(Rs,Dom).
  420
  421R::Dom :-                                 % single var
  422	g_read('clpBNR:thread_init_done',_),  % ensures per-thread initialization 
  423	(var(Dom)                             % domain undefined
  424	 -> (var(R) -> int_decl_(real,_,R) ; true),  % default or domain query (if interval(R) else fail)
  425	    domain(R,Dom)                     % extract domain
  426	 ;  Dom=..[Type|Bounds],              % domain defined
  427	    Val=..[','|Bounds],
  428	    int_decl_(Type,Val,R)
  429	).
  430
  431intervals_([],_Def).
  432intervals_([Int|Ints],Def) :-
  433	Int::Def, !,
  434	intervals_(Ints,Def).
  435
  436int_decl_(boolean,_,R) :- !,          % boolean is integer; 0=false, 1=true, ignore any bounds.
  437	int_decl_(integer,(0,1),R).
  438int_decl_(Type,(','),R) :- !,         % no bounds, fill with vars
  439	int_decl_(Type,(_,_),R).
  440int_decl_(Type,Val,R) :- interval_object(R,CType,CVal,_NL), !,  % already interval
  441	(Type = CType, Val = CVal         % query, no change
  442	 -> true
  443	 ;	Val = (L,H),                  % changing type,bounds?
  444	    lower_bound_val_(Type,L,IL),
  445	    upper_bound_val_(Type,H,IH),
  446	    applyType_(Type, R, T/T, Agenda),           % coerce reals to integers (or no-op).
  447	    ^(CVal,(IL,IH),New),          % low level functional primitive
  448	    updateValue_(CVal, New, R, 1, Agenda, NewAgenda),  % update value (Agenda empty if no value change)
  449	    stable_(NewAgenda)            % then execute Agenda
  450	).
  451int_decl_(Type,(L,H),R) :- var(R), !,    % new interval (R can't be existing interval)
  452	lower_bound_val_(Type,L,IL),
  453	upper_bound_val_(Type,H,IH),
  454	C is cmpr(IL,IH),  % compare bounds
  455	(C == 0
  456	 -> (rational(IL) -> R=IL ; R = IH)  % point range, can unify now
  457	 ;  C == -1,                         % necessary condition: IL < IH
  458	    put_attr(R, clpBNR, interval(Type, (IL,IH), _NL, []))  % attach clpBNR attribute
  459	).
  460int_decl_(Type,(L,H),R) :- constant_type_(Type,R), % R already a point value, check range
  461	lower_bound_val_(Type,L,IL), non_empty(IL,R),  % IL=<R,
  462	upper_bound_val_(Type,H,IH), non_empty(R,IH).  % R=<IH.
  463
  464lower_bound_val_(Type,L,L) :- var(L), !,   % unspecified bound, make it finite
  465	finite_interval(Type,(L,_)).
  466lower_bound_val_(real,L,IL) :-             % real: evaluate and round outward (if float)
  467	((L == pi ; L == e)
  468	 -> IL is roundtoward(L,to_negative)
  469	 ;  Lv is L,
  470	    (preciseBnd(Lv) -> IL=Lv ; IL is nexttoward(Lv,-1.0Inf)) 
  471	).
  472lower_bound_val_(integer,L,IL) :-          % integer: make integer
  473	IL is ceiling(L).
  474lower_bound_val_(boolean,L,IL) :-          % boolean: narrow to L=0
  475	IL is max(0,ceiling(L)).
  476
  477upper_bound_val_(Type,H,H) :- var(H), !,   % unspecified bound, make it finite
  478	finite_interval(Type,(_,H)).
  479upper_bound_val_(real,H,IH) :-             % real: evaluate and round outward (if float)
  480	((H == pi ; H == e)
  481	 -> IH is roundtoward(H,to_positive)
  482	 ;  Hv is H,
  483	    (preciseBnd(Hv) -> IH=Hv ; IH is nexttoward(Hv,1.0Inf)) 
  484	).
  485upper_bound_val_(integer,H,IH) :-          % integer: make integer
  486	IH is floor(H).
  487upper_bound_val_(boolean,H,IH) :-          % boolean: narrow to H=1
  488	IH is min(1,floor(H)).
  489
  490constant_type_(real,C) :- number(C).
  491constant_type_(integer,C) :- integer(C), !.
  492constant_type_(integer,C) :- float(C), float_class(C,infinite).
  493
  494applyType_(NewType, Int, Agenda, NewAgenda) :-      % narrow Int to Type
  495	get_attr(Int,clpBNR,interval(Type,Val,NodeList,Flags)),
  496	(NewType @< Type    % standard order of atoms:  boolean @< integer @< real
  497	 -> (debugging(clpBNR) -> check_monitor_(Int, NewType, interval(Type,Val,NodeList,Flags)) ; true),
  498	    Val = (L,H),
  499	    lower_bound_val_(NewType,L,IL),
  500	    upper_bound_val_(NewType,H,IH),
  501	    (IL=IH
  502	     -> Int=IL  % narrowed to point
  503	     ; 	(put_attr(Int,clpBNR,interval(integer,(IL,IH),NodeList,Flags)),  % set Type (only case allowed)
  504		     linkNodeList_(NodeList, Agenda, NewAgenda)
  505		    )
  506	    )
  507	 ; NewAgenda = Agenda
  508	).
  509
  510%
  511% this goal gets triggered whenever an interval is unified, valid for a numeric value or another interval
  512%
  513attr_unify_hook(IntDef, Num) :-         % unify an interval with a numeric
  514	number(Num),
  515	IntDef = interval(Type,(L,H),Nodelist,_Flags),
  516	constant_type_(Type,Num),                % numeric consistent with type
  517	% L=<Num, Num=<H, assumes L < H
  518	cmpr(L,Num) + cmpr(Num,H) < 0, !,        % and in range (not NaN)
  519	(debugging(clpBNR) -> monitor_unify_(IntDef, Num) ; true),
  520	(var(Nodelist)
  521	 -> true                                 % nothing to do
  522	 ;  linkNodeList_(Nodelist, T/T, Agenda),
  523	    stable_(Agenda)                      % broadcast change
  524	).
  525
  526attr_unify_hook(interval(Type1,V1,Nodelist1,Flags1), Int) :-   % unifying two intervals
  527	get_attr(Int, clpBNR, interval(Type2,V2,Nodelist2,Flags2)),  %%SWIP attribute def.
  528	^(V1,V2,R),                     % unified range is intersection (from ia_primitives),
  529	mergeValues_(Type1, Type2, NewType, R, NewR), !,
  530	mergeNodes_(Nodelist2,Nodelist1,Newlist),  % unified node list is a merge of two lists
  531	mergeFlags_(Flags1,Flags2,Flags),
  532	(debugging(clpBNR) -> monitor_unify_(interval(Type1,V1,_,Flags), Int) ; true),
  533	% update new type, value and constraint list, undone on backtracking
  534	put_attr(Int,clpBNR,interval(NewType,NewR,Newlist,Flags)),
  535	(var(Newlist)
  536	 -> true                                 % nothing to do
  537	 ;  linkNodeList_(Newlist, T/T, Agenda),
  538	    stable_(Agenda)                      % broadcast change
  539	).
  540
  541attr_unify_hook(interval(Type,Val,_Nodelist,_Flags), V) :-   % new value out of range
  542	g_inc('clpBNR:evalNodeFail'),  % count of primitive call failures
  543	debugging(clpBNR),             % fail immediately unless debugging
  544	debug_clpBNR_('Failed to unify ~w::(~w) with ~w',[Type,Val,V]),
  545	fail.
  546
  547% awkward monitor case because original interval gone
  548monitor_unify_(IntDef, Update) :-  % debugging, manufacture temporary interval
  549	put_attr(Temp,clpBNR,IntDef),
  550	check_monitor_(Temp, Update, IntDef).
  551
  552% if types identical, result type and bounds unchanged,
  553% else one is integer so result type integer, and integer bounds applied
  554mergeValues_(T, T, T, R, R) :- !.
  555mergeValues_(_, _, integer, (L,H), (IL,IH)) :-
  556	lower_bound_val_(integer,L,IL),     % type check bounds
  557	upper_bound_val_(integer,H,IH),
  558	non_empty(IL,IH).                   % still non-empty
  559
  560% optimize for one or both lists (dominant case)
  561mergeFlags_([],Flags2,Flags2) :- !.
  562mergeFlags_(Flags1,[],Flags1) :- !.
  563mergeFlags_([F1|Flags1],Flags2,Flags) :-   % discard if F in Flags2 
  564	functor(F1,N,1),                       % ignore value
  565	functor(F2,N,1),
  566	memberchk(F2,Flags2), !,
  567	mergeFlags_(Flags1,Flags2,Flags).
  568mergeFlags_([F1|Flags1],Flags2,[F1|Flags]) :-  % keep F, not in Flags2 
  569	mergeFlags_(Flags1,Flags2,Flags).
  570
  571% merge two node lists removing duplicates based on operation and arguments
  572mergeNodes_([N],NodeList,NodeList) :- var(N),!.         % end of list
  573mergeNodes_([node(Op,_,_,Ops)|Ns],NodeList,NewList) :-  % if same Op and Ops, discard
  574	matching_node_(NodeList,Op,Ops), !,
  575	mergeNodes_(Ns,NodeList,NewList).
  576mergeNodes_([N|Ns],NodeList,[N|NewList]) :-             % not a duplicate, retain
  577	mergeNodes_(Ns,NodeList,NewList).
  578
  579matching_node_([node(Op,_,_,NOps)|_Ns],Op,Ops) :-
  580	NOps==Ops, !.  % identical args
  581matching_node_([N|Ns],Op,Ops) :-
  582	nonvar(N),     % not end of list
  583	matching_node_(Ns,Op,Ops).
  584
  585%
  586% New Constraints use { ... } syntax.
  587%
  588{Cons} :-
  589	g_read('clpBNR:thread_init_done',_),      % ensures per-thread initialization
  590	term_variables(Cons, CVars),
  591	declare_vars_(CVars),                     % convert any plain vars to intervals
  592	addConstraints_(Cons,T/T,Agenda),         % add constraints
  593	stable_(Agenda).                          % then execute Agenda
  594
  595declare_vars_([]).
  596declare_vars_([CV|CVars]) :-
  597	(interval(CV) -> true ; new_interval_(CV,real)),
  598	declare_vars_(CVars).
  599
  600new_interval_(V,Type) :-
  601	universal_interval(UI),
  602	int_decl_(Type,UI,V).
  603
  604addConstraints_([],Agenda,Agenda) :- !.
  605addConstraints_([C|Cs],Agenda,NewAgenda) :-
  606	nonvar(C),
  607	addConstraints_(C,Agenda,NextAgenda), !,
  608	addConstraints_(Cs,NextAgenda,NewAgenda).
  609addConstraints_((C,Cs),Agenda,NewAgenda) :-  % Note: comma as operator
  610	nonvar(C),
  611	addConstraints_(C,Agenda,NextAgenda), !,
  612	addConstraints_(Cs,NextAgenda,NewAgenda).
  613addConstraints_(C,Agenda,NewAgenda) :-
  614	constraint_(C),    % a constraint is a boolean expression that evaluates to true
  615	simplify(C,CS),    % optional
  616	buildConstraint_(CS, Agenda, NewAgenda).
  617
  618% low overhead version for internal use (also used by arithmetic_types pack)
  619constrain_(C) :- 
  620	buildConstraint_(C,T/T,Agenda),
  621	stable_(Agenda).
  622	
  623buildConstraint_(C,Agenda,NewAgenda) :-
  624	debug_clpBNR_('Add ~p',{C}),
  625	% need to catch errors from ground expression evaluation
  626	catch(build_(C, 1, boolean, Agenda, NewAgenda),_Err,fail), !.
  627buildConstraint_(C,_Agenda,_NewAgenda) :-
  628	debug_clpBNR_('{} failure due to bad or inconsistent constraint: ~p',{C}),
  629	fail.
  630
  631:- include(clpBNR/ia_simplify).  % simplifies constraints to a hopefully more efficient equivalent
  632
  633%
  634% build a node from an expression
  635%
  636build_(Int, Int, VarType, Agenda, NewAgenda) :-
  637	interval(Int), !,                                   % existing interval object
  638	applyType_(VarType, Int, Agenda, NewAgenda).        % coerces exsiting intervals to required type
  639build_(Var, Var, VarType, Agenda, Agenda) :-            % implicit interval creation.
  640	var(Var), !,
  641	new_interval_(Var,VarType).
  642build_(Num, Int, VarType, Agenda, Agenda) :-            % floating point constant, may not be precise
  643	float(Num), !,
  644	(
  645	float_class(Num,infinite)
  646	 -> Int=Num                                         % infinities are point values
  647	 ;	int_decl_(VarType,(Num,Num),Int)                % may be fuzzed, so not a point
  648	).
  649build_(::(L,H), Int, VarType, Agenda, Agenda) :-        % hidden :: feature: interval bounds literal
  650	number(L), number(H), !,
  651	C is cmpr(L,H),  % compare bounds
  652	(C == 0
  653	 -> (rational(L) -> Int=L ; Int=H)                  % point value, if either bound rational (precise), use it
  654	 ;	C == -1,                                        % necessary condition: L < H
  655	    once(VarType == real ; true),                   % if undefined Type, use 'real'
  656	    put_attr(Int, clpBNR, interval(VarType, (L,H), _NL, []))  % create clpBNR attribute
  657	).
  658build_(Num, Int, VarType, Agenda, Agenda) :-            % pre-compile constants pi and e
  659	(Num == pi ; Num == e), !,
  660	int_decl_(VarType,(Num,Num),Int).  
  661build_(Exp, Num, _, Agenda, Agenda) :-                  % pre-compile any ground precise Exp
  662	ground(Exp),
  663	safe_(Exp),                                         % safe to evaluate using is/2
  664	Num is Exp,
  665	% (Num == 1.5NaN -> !, fail ; preciseBnd(Num))      % fail build_ if nan
  666	preciseBnd(Num),                                    % precise result
  667	!.
  668build_(Exp, Z, _, Agenda, NewAgenda) :-                 % deconstruct to primitives
  669	Exp =.. [F|Args],
  670	fmap_(F,Op,[Z|Args],ArgsR,Types), !,                % supported arithmetic op
  671	build_args_(ArgsR,Objs,Types,Agenda,ObjAgenda),
  672	newNode_(Op,Objs,ObjAgenda,NewAgenda).
  673build_(Exp, Z, _, Agenda, NewAgenda) :-                 % user defined
  674	Exp =.. [Prim|Args],
  675	chk_primitive_(Prim),
  676	build_args_([Z|Args],Objs,_Types,Agenda,ObjAgenda),
  677	newNode_(user(Prim),Objs,ObjAgenda,NewAgenda).
  678
  679build_args_([],[],_,Agenda,Agenda).
  680build_args_([Arg|ArgsR],[Obj|Objs],[Type|Types],Agenda,NewAgenda) :-
  681	(var(Type) -> Type=real ; true),                    % var if user defined
  682	build_(Arg,Obj,Type,Agenda,NxtAgenda),
  683	build_args_(ArgsR,Objs,Types,NxtAgenda,NewAgenda).
  684
  685chk_primitive_(Prim) :-  % wraps safe usage of unsafe current_predicate/2
  686	UsrHead =..[Prim,'$op',_,_,_],
  687	current_predicate(_,clpBNR:UsrHead).
  688
  689sandbox:safe_primitive(clpBNR:chk_primitive_(_Prim)).
  690
  691% to invoke user defined primitive
  692call_user_primitive(Prim, P, InArgs, OutArgs) :-  % wraps unsafe meta call/N
  693	call(clpBNR:Prim, '$op', InArgs, OutArgs, P).
  694
  695% really unsafe, but in a pengine a user can't define any predicates in another module, so this is safe
  696sandbox:safe_meta(clpBNR:call_user_primitive(_Prim, _P, _InArgs, _OutArgs), []).
  697
  698% only called when argument is ground
  699safe_(E) :- atomic(E), !.  % all atomics, including [] - allows possibility of user defined arithmetic types
  700safe_([A|As]) :- !,
  701	safe_(A),
  702	safe_(As).
  703safe_(_ xor _) :- !,                                    % clpBNR xor incompatible with `is` xor
  704	fail.
  705safe_(integer(_)) :- !,                                 % clpBNR integer incompatible with `is` integer
  706	fail.
  707safe_(_/Z) :- 0.0 is abs(Z), !,                         % division by 0.0 (or -0.0) requires fuzzed 0.0
  708	fail.
  709safe_(F) :- 
  710	current_arithmetic_function(F),                     % evaluable by is/2
  711	F =.. [_Op|Args],
  712	safe_(Args).
  713
  714%  a constraint must evaluate to a boolean 
  715constraint_(C) :- nonvar(C), C =..[Op|_], fmap_(Op,_,_,_,[boolean|_]), !.
  716
  717%  map constraint operator to primitive/arity/types
  718fmap_(+,    add,   ZXY,     ZXY,     [real,real,real]).
  719fmap_(-,    add,   [Z,X,Y], [X,Z,Y], [real,real,real]).     % note subtract before minus 
  720fmap_(*,    mul,   ZXY,     ZXY,     [real,real,real]).
  721fmap_(/,    mul,   [Z,X,Y], [X,Z,Y], [real,real,real]).
  722fmap_(**,   pow,   ZXY,     ZXY,     [real,real,real]).
  723fmap_(min,  min,   ZXY,     ZXY,     [real,real,real]).
  724fmap_(max,  max,   ZXY,     ZXY,     [real,real,real]).
  725fmap_(==,   eq,    ZXY,     ZXY,     [boolean,real,real]).  % strict equality
  726fmap_(=:=,  eq,    ZXY,     ZXY,     [boolean,real,real]).  % Prolog compatible, strict equality
  727fmap_(is,   eq,    ZXY,     ZXY,     [boolean,real,real]).
  728fmap_(<>,   ne,    ZXY,     ZXY,     [boolean,integer,integer]).
  729fmap_(=\=,  ne,    ZXY,     ZXY,     [boolean,integer,integer]).  % Prolog compatible
  730fmap_(=<,   le,    ZXY,     ZXY,     [boolean,real,real]).
  731fmap_(>=,   le,    [Z,X,Y], [Z,Y,X], [boolean,real,real]).
  732fmap_(<,    lt,    ZXY,     ZXY,     [boolean,real,real]).
  733fmap_(>,    lt,    [Z,X,Y], [Z,Y,X], [boolean,real,real]).
  734fmap_(<=,   in,    ZXY,     ZXY,     [boolean,real,real]).  % inclusion/subinterval
  735
  736fmap_(and,  and,   ZXY,     ZXY,     [boolean,boolean,boolean]).
  737fmap_(',',  and,   ZXY,     ZXY,     [boolean,boolean,boolean]).  % for usability
  738fmap_(or,   or,    ZXY,     ZXY,     [boolean,boolean,boolean]).
  739fmap_(nand, nand,  ZXY,     ZXY,     [boolean,boolean,boolean]).
  740fmap_(nor,  nor,   ZXY,     ZXY,     [boolean,boolean,boolean]).
  741fmap_(xor,  xor,   ZXY,     ZXY,     [boolean,boolean,boolean]).
  742fmap_(->,   imB,   ZXY,     ZXY,     [boolean,boolean,boolean]).
  743
  744fmap_(sqrt, sqrt,  ZX,      ZX,      [real,real]).          % pos. root version vs. **(1/2)
  745fmap_(-,    minus, ZX,      ZX,      [real,real]).
  746fmap_(~,    not,   ZX,      ZX,      [boolean,boolean]).
  747fmap_(integer,int, ZX,      ZX,      [boolean,real]).
  748fmap_(exp,  exp,   ZX,      ZX,      [real,real]).
  749fmap_(log,  exp,   [Z,X],   [X,Z],   [real,real]).
  750fmap_(abs,  abs,   ZX,      ZX,      [real,real]).
  751fmap_(sin,  sin,   ZX,      ZX,      [real,real]).
  752fmap_(asin, sin,   [Z,X],   [X,Z],   [real,real]).
  753fmap_(cos,  cos,   ZX,      ZX,      [real,real]).
  754fmap_(acos, cos,   [Z,X],   [X,Z],   [real,real]).
  755fmap_(tan,  tan,   ZX,      ZX,      [real,real]).
  756fmap_(atan, tan,   [Z,X],   [X,Z],   [real,real]).
  757
  758% reverse map node info to a facsimile of the original constraint
  759map_constraint_op_(integral,$(V),integral(V)) :- !.
  760map_constraint_op_(user(Func),Args,C) :- !,
  761	remap_(Func,Args,C).
  762map_constraint_op_(Op,Args,C) :-
  763	fmap_(COp,Op,_,_,_),
  764	remap_(COp,Args,C),
  765	!.
  766
  767remap_(Op,$(Z,X,Y),C) :- constraint_(Op), Z==1, !,  % simplification for binary constraints
  768	C=..[Op,X,Y]. 
  769remap_(Op,$(Z,X),C) :- constraint_(Op), Z==1, !,    % simplification for unary constraints (~)
  770	C=..[Op,X].
  771remap_(Op,$(Z,X,Y),Z==C) :- !,
  772	C=..[Op,X,Y].
  773remap_(Op,$(Z,X),Z==C) :-
  774	C=..[Op,X].
  775
  776%
  777% Node constructor
  778%
  779% First clause is an optimization for E1 == E2
  780%	In that case just unify E1 and E2; heavy lifting done by attr_unify_hook.
  781newNode_(eq, [Z,X,Y], Agenda, Agenda) :- Z==1, !, X=Y.
  782newNode_(Op, Objs, Agenda, NewAgenda) :-
  783	Args =.. [$|Objs],  % store arguments as $/N where N=1..3
  784	NewNode = node(Op, _P, 0, Args),  % L=0
  785	addNode_(Objs,NewNode),
  786	% increment count of added nodes, will be decremented on backtracking/failure
  787	g_incb('clpBNR:node_count'),
  788	linkNode_(Agenda, NewNode, NewAgenda).
  789
  790addNode_([],_Node).
  791addNode_([Arg|Args],Node) :-
  792	(interval_object(Arg, _Type, _Val, Nodelist) -> newmember(Nodelist, Node) ; true),
  793	addNode_(Args,Node).
  794
  795sandbox:safe_global_variable('clpBNR:node_count').
  796
  797clpStatistics :-
  798	g_assign('clpBNR:node_count',0),  % reset/initialize node count to 0
  799	fail.  % backtrack to reset other statistics.
  800
  801clpStatistic(node_count(C)) :-
  802	g_read('clpBNR:node_count',C).
  803
  804% extend list with N
  805newmember([X|Xs],N) :- 
  806	(nonvar(X)
  807	 -> newmember(Xs,N) % not end of (indefinite) list.
  808     ;  X = N           % end of list, insert N and new tail Xs
  809    ).
  810
  811%
  812% Process Agenda to narrow intervals (fixed point iteration)
  813%
  814stable_([]/[]) :- !.  % small optimization when agenda is empty
  815stable_(Agenda) :-
  816	current_prolog_flag(clpBNR_iteration_limit,Ops),  % budget for current operation
  817	stableLoop_(Agenda,Ops),
  818	!.  % achieved stable state with empty Agenda -> commit.
  819
  820stableLoop_([]/[], OpsLeft) :- !,           % terminate successfully when agenda comes to an end
  821	g_read('clpBNR:iterations',Cur),        % maintain "low" water mark (can be negative)
  822	(OpsLeft<Cur -> g_assign('clpBNR:iterations',OpsLeft) ; true),
  823	(OpsLeft<0 -> E is -OpsLeft, debug_clpBNR_('Iteration throttle limit exceeded by ~w ops.',E) ; true).
  824stableLoop_([Node|Agenda]/T, OpsLeft) :-
  825	Node = node(Op,P,_,Args),  % if node on queue ignore link bit (was: Node = node(Op,P,1,Args))
  826	doNode_(Args, Op, P, OpsLeft, DoOver, Agenda/T, NxtAgenda),  % undoable on backtrack
  827	nb_setarg(3,Node,0),                    % reset linked bit
  828	% if doNode_ requested DoOver and above Ops threshold, put Node back at the end of the queue 
  829	(atom(DoOver), OpsLeft > 0 -> linkNode_(NxtAgenda,Node,NewAgenda) ; NewAgenda = NxtAgenda),
  830	RemainingOps is OpsLeft-1,              % decrement OpsLeft (can go negative)
  831	stableLoop_(NewAgenda,RemainingOps).
  832
  833% support for max_iterations statistic
  834sandbox:safe_global_variable('clpBNR:iterations').
  835
  836clpStatistics :-
  837	current_prolog_flag(clpBNR_iteration_limit,L), 
  838	g_assign('clpBNR:iterations',L),  % set "low" water mark to upper limit
  839	fail.  % backtrack to reset other statistics.
  840
  841clpStatistic(max_iterations(O/L)) :-
  842	g_read('clpBNR:iterations',Ops),
  843	current_prolog_flag(clpBNR_iteration_limit,L),
  844	O is L-Ops.  % convert "low" water mark to high water mark
  845
  846%
  847% doNode_/7 : Evaluate a node and add new nodes to end of queue. `evalNode` primitives can
  848%	 fail, resulting in eventual failure of `stable_`, i.e., inconsistent constraint set.
  849%
  850% Note: primitives operate on interval values (sets) only, unaware of common arguments,
  851% so additional consistency checks required on update.
  852%
  853doNode_($(ZArg,XArg,YArg), Op, P, OpsLeft, DoOver, Agenda, NewAgenda) :-  % Arity 3 Op
  854	(var(P)                                          % check persistent bit
  855	 -> getValue(ZArg,ZVal),
  856	    getValue(XArg,XVal),
  857	    getValue(YArg,YVal),
  858	    evalNode(Op, P, $(ZVal,XVal,YVal), $(NZVal,NXVal,NYVal)),  % can fail
  859	    % enforce consistency for common arguments by intersecting and redoing as required.
  860	    (var(ZArg)                % if Z == X and/or Y
  861	     -> (ZArg==XArg -> consistent_value_(NZVal,NXVal,NZ1,DoOver) ; NZ1 = NZVal),
  862	        (ZArg==YArg -> consistent_value_(NZ1,  NYVal,NZ2,DoOver) ; NZ2 = NZ1),
  863	        updateValue_(ZVal, NZ2, ZArg, OpsLeft, Agenda, AgendaZ)
  864	     ;  AgendaZ = Agenda
  865	    ),
  866	    (var(XArg), XArg==YArg    % if X==Y
  867	     -> consistent_value_(NXVal,NYVal,NVal,DoOver),
  868	        updateValue_(XVal, NVal,  XArg, OpsLeft, AgendaZ, NewAgenda)  % only one update needed
  869	     ;  updateValue_(XVal, NXVal, XArg, OpsLeft, AgendaZ, AgendaZX),
  870	        updateValue_(YVal, NYVal, YArg, OpsLeft, AgendaZX, NewAgenda)
  871	    )
  872	 ; % P = p, trim persistent nodes from all arguments
  873	    trim_op_(ZArg), trim_op_(XArg), trim_op_(YArg),  
  874	    NewAgenda = Agenda
  875	).
  876
  877doNode_($(ZArg,XArg), Op, P, OpsLeft, DoOver, Agenda, NewAgenda) :-       % Arity 2 Op
  878	(var(P)                                          % check persistent bit
  879	 -> getValue(ZArg,ZVal),
  880	    getValue(XArg,XVal),
  881	    evalNode(Op, P, $(ZVal,XVal), $(NZVal,NXVal)),      % can fail
  882	    % enforce consistency for common arguments by intersecting and redoing as required.
  883	    (var(ZArg), ZArg==XArg    % if Z==X
  884	     -> consistent_value_(NZVal,NXVal,NVal,DoOver),     % consistent value, DoOver if needed
  885	        updateValue_(ZVal, NVal,  ZArg, OpsLeft, Agenda, NewAgenda)  % only one update needed
  886	     ;  updateValue_(ZVal, NZVal, ZArg, OpsLeft, Agenda, AgendaZ),
  887	        updateValue_(XVal, NXVal, XArg, OpsLeft, AgendaZ, NewAgenda)
  888	    )
  889	 ; % P = p, trim persistent nodes from all arguments
  890	    trim_op_(ZArg), trim_op_(XArg),
  891	    NewAgenda = Agenda
  892	).
  893 
  894doNode_($(Arg), Op, P, _OpsLeft, _, Agenda, NewAgenda) :-                 % Arity 1 Op
  895	(var(P)                                          % check persistent bit
  896	 -> getValue(Arg,Val),
  897	    evalNode(Op, P, $(Val), $(NVal)),                   % can fail causing stable_ to fail => backtracking
  898	    updateValue_(Val, NVal, Arg, 1, Agenda,NewAgenda)   % always update value regardless of OpsLeft limiter	
  899	 ;  % P = p, trim persistent nodes from argument
  900	    trim_op_(Arg),
  901	    NewAgenda = Agenda
  902	).
  903
  904consistent_value_(Val,Val,Val,_) :- !.                       % same value
  905consistent_value_(Val1,Val2,Val,true) :- ^(Val1,Val2,Val).   % different values, intersect
  906
  907%
  908% remove any persistent nodes from Arg
  909%	called whenever a persistent node is encountered in FP iteration
  910%
  911trim_op_(Arg) :- 
  912	( get_attr(Arg, clpBNR, Def)     % an interval ?
  913	 -> arg(3,Def,NList),  %%Def = interval(_, _, NList, _),
  914		trim_persistent_(NList,TrimList),
  915		% if trimmed list empty, set to a new unshared var to avoid cycles(?) on backtracking
  916		(var(TrimList) -> setarg(3,Def,_) ; setarg(3,Def,TrimList))  % write trimmed node list
  917	 ; true  % assumed to be a number, nothing to trim
  918	).
  919
  920trim_persistent_(T,T) :- var(T), !.    % end of indefinite node list
  921trim_persistent_([node(_,P,_,_)|Ns],TNs) :- nonvar(P), !, trim_persistent_(Ns,TNs).
  922trim_persistent_([N|Ns],[N|TNs]) :- trim_persistent_(Ns,TNs).
  923
  924%
  925% Any changes in interval values should come through here.
  926% Note: This captures all updated state for undoing on backtracking
  927%
  928updateValue_(Old, Old, _, _, Agenda, Agenda) :- !.                  % no change in value (constant?)
  929
  930updateValue_(Old, New, Int, OpsLeft, Agenda, NewAgenda) :-          % set interval value to New
  931	(OpsLeft>0 -> true ; propagate_if_(Old, New)), !,  % if OpsLeft >0 or narrowing sufficent
  932	putValue_(New, Int, Nodelist),               % update value (may fail)
  933	linkNodeList_(Nodelist, Agenda, NewAgenda).  % then propagate change
  934
  935updateValue_(_, _, _, _, Agenda, Agenda).        % otherwise just continue with Agenda
  936
  937% propgate if sufficient narrowing (> 10%)
  938propagate_if_((OL,OH), (NL,NH)) :- (NH-NL)/(OH-OL) < 0.9.  % any overflow in calculation will propagate
  939
  940linkNodeList_([X|Xs], List, NewList) :-
  941	(var(X)
  942	 -> List = NewList                               % end of indefinite list
  943	 ;  (arg(3,X,Linked), Linked == 1                % test linked flag (only SWIP VM codes)
  944	     -> linkNodeList_(Xs, List, NewList)         % add rest of the node list
  945	     ;  linkNode_(List, X, NextList),            % not linked add it to list
  946	        linkNodeList_(Xs, NextList, NewList)     % add rest of the node list
  947	    )
  948	).
  949
  950linkNode_(List/[X|NextTail], X, List/NextTail) :-    % add to list
  951	setarg(3,X,1).                                   % set linked bit
  952
  953:- include(clpBNR/ia_utilities).  % print,solve, etc.
  954
  955%
  956% set monitor action on an interval using `watch` flags
  957%
  958watch(Int,Action) :-
  959	atom(Action), 
  960	current_module(clpBNR),  % this just marks watch/2 as unsafe regarding body
  961	get_interval_flags_(Int,Flags), !,
  962	remove_(Flags,watch(_),Flags1),
  963	(Action = none -> Flags2=Flags1 ; Flags2=[watch(Action)|Flags1]),
  964	set_interval_flags_(Int,Flags2).
  965watch(Ints,Action) :-
  966	list(Ints),
  967	watch_list_(Ints,Action).
  968
  969remove_([],_,[]).
  970remove_([X|Xs],X,Xs) :- !.
  971remove_([X|Xs],X,[X|Ys]) :-
  972	remove_(Xs,X,Ys).
  973
  974watch_list_([],_Action).
  975watch_list_([Int|Ints],Action) :-
  976	watch(Int,Action),
  977	watch_list_(Ints,Action).
  978
  979% check if watch enabled on this interval
  980check_monitor_(Int, Update, interval(_Type,_Val,_Nodelist,Flags)) :-
  981	(memberchk(watch(Action), Flags)
  982	 -> once(monitor_action_(Action, Update, Int))
  983	 ; true
  984	).
  985
  986%
  987% for monitoring changes - all actions defined here
  988%
  989monitor_action_(trace, Update, Int) :-  !, % log changes to console and enter debugger
  990	monitor_action_(log, Update, Int),
  991	trace.
  992monitor_action_(log, Update, Int) :-  var(Update), !,  % log interval unify
  993	debug_clpBNR_('Unify ~p with ~p',[Int,Update]).
  994monitor_action_(log, Update, Int) :-  number(Update), !,    % log interval unify with number
  995	domain(Int,Dom),
  996	debug_clpBNR_('Unify _?{~p} with ~p',[Dom,Update]).
  997monitor_action_(log, integer, Int) :-  !,  % change type to integer
  998	debug_clpBNR_('Set type of  ~p to ~p',[Int,integer]).
  999monitor_action_(log, Val, Int) :-  !,  % narrow range
 1000	debug_clpBNR_('Set value of ~p to (~p)',[Int,Val]).
 1001monitor_action_(_, _, _).  % default to noop (as in 'none')
 1002
 1003sandbox:safe_primitive(clpBNR:watch(_Int,Action)) :- % watch(X,trace) is unsafe.
 1004	Action \= trace.
 1005% only safe because watch(X,trace) is unsafe.
 1006sandbox:safe_primitive(clpBNR:monitor_action_(_Action, _Update, _Int)).
 1007
 1008%
 1009% tracing doNode_ support - using wrap_predicate(:Head, +Name, -Wrapped, +Body)/4
 1010% trace_clpBNR/1 is unsafe (wrapping is global)
 1011%
 1012:- use_module(library(prolog_wrap)). 1013
 1014trace_clpBNR(Bool)  :-                  % query or already in defined state
 1015	( current_predicate_wrapper(clpBNR:doNode_(_Args, _Op, _P, _OpsLeft, _DoOver, _Agenda, _NewAgenda), 
 1016	                            'clpBNR:doNode_', _Wrapped, _Body)
 1017	 -> Bool = true ; Bool = false
 1018	),
 1019	!.
 1020trace_clpBNR(true)  :-                  % turn on wrapper
 1021	wrap_predicate(clpBNR:doNode_(Args, Op, _P, _OpsLeft, _DoOver, _Agenda, _NewAgenda),
 1022	                   'clpBNR:doNode_', 
 1023	                   Wrapped, 
 1024	                   doNode_wrap_(Wrapped, Args,Op)).
 1025trace_clpBNR(false) :-                  % turn off wrapper
 1026	unwrap_predicate(clpBNR:doNode_/7,  %(Args, Op, P, OpsLeft, DoOver, Agenda, NewAgenda),
 1027	                   'clpBNR:doNode_').
 1028
 1029doNode_wrap_(Wrapped, Args,Op) :-
 1030	map_constraint_op_(Op,Args,C),
 1031	Wrapped,                 % execute wrapped doNode_
 1032	debug_clpBNR_("~p.",C).  % print trace message , fail messages from evalNode_, attr_unify_hook
 1033
 1034%
 1035% Get all defined statistics
 1036%
 1037clpStatistics(Ss) :- findall(S, clpStatistic(S), Ss).
 1038
 1039% end of reset chain succeeds.
 1040clpStatistics.
 1041
 1042%
 1043% module initialization
 1044%
 1045init_clpBNR :-
 1046	restore_optimize_flags_,
 1047	print_message(informational, clpBNR(versionInfo)),
 1048	print_message(informational, clpBNR(arithmeticFlags)).  % cautionary, set on first use
 1049
 1050check_hooks_safety :-  % safety check on any hooks (test only code)
 1051	% Note: calls must have no side effects
 1052	ignore(attr_portray_hook([],_)),                                            % fails
 1053	ignore(user:exception(undefined_global_variable,'clpBNR:thread_init_done',[])),  % fails
 1054%	ignore(term_html:portray('$clpBNR...'(_),_,_,_)),                           % fails
 1055	ignore(user:portray('$clpBNR...'(_))).                                      % fails
 1056
 1057:- multifile prolog:message//1. 1058
 1059prolog:message(clpBNR(versionInfo)) -->
 1060	{ version(Version) },
 1061	[ '*** clpBNR v~w ***.'-[Version] ].
 1062
 1063prolog:message(clpBNR(arithmeticFlags)) -->
 1064	[ '  Arithmetic global flags will be set to prefer rationals and IEEE continuation values.'-[] ].
 1065
 1066:- initialization(init_clpBNR, now).  % Most initialization deferred to "first use" - see user:exception/3