1/*	The MIT License (MIT)
    2 *
    3 *	Copyright (c) 2019-2024 Rick Workman
    4 *
    5 *	Permission is hereby granted, free of charge, to any person obtaining a copy
    6 *	of this software and associated documentation files (the "Software"), to deal
    7 *	in the Software without restriction, including without limitation the rights
    8 *	to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
    9 *	copies of the Software, and to permit persons to whom the Software is
   10 *	furnished to do so, subject to the following conditions:
   11 *
   12 *	The above copyright notice and this permission notice shall be included in all
   13 *	copies or substantial portions of the Software.
   14 *
   15 *	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
   16 *	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   17 *	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
   18 *	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   19 *	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
   20 *	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
   21 *	SOFTWARE.
   22 */
   23
   24%
   25% compatibility predicate
   26%
   27% print_interval(Term), print_interval(Stream,Term)
   28%	prints Term to optional Stream with intervals expanded to domains
   29%	uses format/3 so extended Stream options, e.g., atom(A), are supported
   30%
   31print_interval(Term) :- 
   32	printed_term_(Term,Out),
   33	format('~w',[Out]).          % safe
   34
   35print_interval(Stream,Term) :-
   36	printed_term_(Term,Out),
   37	format(Stream,'~w',[Out]).   % unsafe
   38
   39printed_term_(Term,Out) :-
   40	copy_term_nat(Term,Out),     % copy of Term without attributes for output
   41	term_variables(Term,TVars),
   42	term_variables(Out,OVars),
   43	name_vars_(TVars,OVars,0).   % name variables, intervals replaced by declarations
   44
   45name_vars_([],[],_).
   46name_vars_([TVar|TVars],[OVar|OVars],N) :-
   47	(domain(TVar,Dom)
   48	 -> OVar = (V::Dom)
   49	  ; OVar = V
   50	),
   51	number_chars(N,C),atom_chars(V,['V'|C]),
   52	N1 is N+1,
   53	name_vars_(TVars,OVars,N1).
   54
   55%
   56% SWIP attribute portray hook - used if write_attributes=portray
   57%
   58attr_portray_hook(interval(Type,Val,_Node,_Flags),_Int) :-
   59	interval_domain_(Type,Val,Dom),
   60	format('~w',Dom).                       % safe
   61
   62%
   63% Support ellipsis format in answers
   64%
   65user:portray('$clpBNR...'(Out)) :-          % remove quotes on stringified number in ellipsis format
   66	string(Out),
   67	format('~W...',[Out,[quoted(false)]]).  % safe
   68
   69%
   70% SWIP hooks for answer formatting - two modes : 
   71%	Verbose=true,  display all interval domains and associated constraints
   72%	Verbose=false, display var domains with no constraints - internal vars omitted
   73%
   74user:expand_answer(Bindings,Bindings) :-          % for toplevel answer control
   75	current_prolog_flag(clpBNR_verbose,Verbose),  % save verbose mode in attributes
   76	add_names_(Bindings,Verbose).
   77
   78% annotate interval variables in Bindings
   79add_names_([],_).
   80add_names_([Var|Bindings],Verbose) :- var(Var), !,
   81	add_names_(['' = Var|Bindings],Verbose).
   82add_names_([Name = Var|Bindings],Verbose) :-
   83	(get_interval_flags_(Var,Flags)
   84	 -> set_interval_flags_(Var,[name(Name,Verbose)|Flags]),        % mainly to attach Verbose
   85	    (Verbose == false -> reset_interval_nodelist_(Var) ; true)  % Nodes restored on backtrack
   86	 ;  (compound(Var)                 % Var not a clpBNR interval, if compound...
   87	     -> term_variables(Var,Vars),
   88	        add_names_(Vars,Verbose)   % mark any internal intervals  
   89	     ; true                        % else nothing to mark
   90	    )
   91	),
   92	add_names_(Bindings,Verbose).
   93
   94%
   95%  SWISH versions:
   96%
   97:- if(current_prolog_flag(clpBNR_swish, true)).   98
   99% portray (HTML)
  100:- use_module(library(http/html_write)).  101:- multifile(term_html:portray//2).  102
  103term_html:portray('$clpBNR...'(Out),_) -->   % avoid quotes on stringified number in ellipsis format
  104	{ string(Out) },
  105	html(span(class('pl-float'), [Out, '...'])).
  106
  107% for answer formatting
  108:- multifile(swish_trace:post_context/1).  109
  110swish_trace:post_context(Dict) :-
  111	expand_answer(Dict.bindings,_).
  112
  113:- endif.  % current_prolog_flag(clpBNR_swish, true)
  114
  115%
  116% Constructs goals to build interval(X)
  117%
  118attribute_goals(X) -->
  119	{(interval_object(X, Type, Val, Nodes)
  120 	  -> (get_interval_flags_(X,Flags), memberchk(name(_,false),Flags), var(Nodes) 
  121 	      -> intValue_(Val,Type,Dom),           % non-verbose and empty nodelist => "pretty" value
  122	         Goals = [X::Dom]
  123	      ;  interval_domain_(Type, Val, Dom),  % full copy 
  124	         constraints_(Nodes,X,Cs),          % reverse map to list of original constraints
  125	         to_comma_exp_(Cs,X::Dom,Goals)
  126	     )
  127	  ;  Goals = []
  128	 )
  129	},
  130	list_dcg_(Goals).  % avoids "unsafe" call to phrase/3
  131
  132list_dcg_([]) --> [].
  133list_dcg_([H|T]) --> [H], list_dcg_(T).
  134
  135constraints_([Node],_,[]) :- var(Node), !.  % end of indefinite list
  136constraints_([node(Op,P,_,Args)|Nodes],X,[C|Cs]) :-
  137	var(P),                       % skip persistent nodes (already in value)
  138	term_variables(Args,[V|_]),   % this ensures constraints only get output once
  139	V==X,
  140	map_constraint_op_(Op,Args,C),
  141	constraints_(Nodes,X,Cs).
  142constraints_([_|Nodes],X,Cs) :-
  143	constraints_(Nodes,X,Cs).
  144		
  145to_comma_exp_([],Decl,[Decl]).
  146to_comma_exp_([N|NCs],Decl,[(Decl,{Cs})]) :-
  147	to_comma_cexp_(NCs,N,Cs).
  148
  149to_comma_cexp_([],In,In).
  150to_comma_cexp_([C|Cs],In,Out) :-
  151	to_comma_cexp_(Cs,(In,C),Out).
  152
  153%
  154%  Simplified output format for interval ranges
  155%
  156% 1. Any 'real' rational bounds output as floats (conservatively rounded).
  157% 2. Any subnormal bounds converted to zero.
  158% 3. Sufficiently narrow reals output in ellipsis format.
  159% 4. Any real intervals with bounds zero or subnormal output as 0.0...
  160% 5. Query interval variables output as V = V::Domain
  161% 6. Subterm intervals as sub-terms printed as domain (Type(L,H)).
  162%
  163intValue_((0,1),integer,boolean).                  % boolean
  164intValue_((L,H),real,'$clpBNR...'(Out)) :-         % virtual zero (zero or subnormal) 
  165	zero_float_(L),
  166	zero_float_(H), !,
  167	format(chars(Zero),"~16f",0.0),
  168	string_chars(Out,[' '|Zero]).                  % space prefix 
  169intValue_((L,H),real,'$clpBNR...'(Out)) :-         % two floats with minimal leading match
  170	float_chars_(L,LC),
  171	float_chars_(H,HC),
  172	sign_chars_(LC,HC,ULC,UHC,Codes/Match),
  173	leading_zero(ULC,UHC,ZLC),
  174	matching_(ZLC,UHC,Match,0,Dec,MLen),
  175	MLen>Dec+1, !,	% minimum of one matching digit after decimal point
  176	string_codes(Out,Codes).
  177intValue_((L,H),Type,Dom) :-                       % default, just convert rationals
  178	rational_fraction_(L,FL),
  179	rational_fraction_(H,FH),
  180	interval_domain_(Type, (FL,FH), Dom).
  181
  182zero_float_(B) :- float(B), float_class(B,C), (C=zero ; C=subnormal).
  183
  184float_chars_(B,Cs) :-
  185	rational_fraction_(B,F),
  186	float(F), float_class(F, normal),
  187	format(chars(Cs),'~16f',F),  % same length after '.' (pads with trailing 0's)
  188	length(Cs,Len), Len=<32.     % reverts to precise format if too long
  189
  190rational_fraction_(B,FB) :-
  191	rational(B,_,D), D \= 1, !,  % non-integer rational
  192	FB is float(B).
  193rational_fraction_(F,0.0) :- 
  194	float(F),
  195	float_class(F,subnormal), !.
  196rational_fraction_(B,B).
  197
  198sign_chars_(['-'|LC],['-'|HC],HC,LC,[' ','-'|T]/T) :- !.  % negate, save sign and reverse bounds
  199sign_chars_(LC,HC,LC,HC,[' '|T]/T).  % need a character that doesn't print: 2=STX (Start of Text)
  200
  201leading_zero(['9'|ULC],['1','0'|_],['0','9'|ULC]) :- !.
  202leading_zero(ULC,_,ULC).
  203
  204% Note: match should never be exact (L=H) so [] case not required
  205% matching_([],[],[],N,Dec,N).
  206matching_([C,'.'|LCs],[C,'.'|HCs],[C,'.'|Cs],N,Dec,Nout) :-   !, % absorbing "."
  207	digit_(C,_),
  208	Dec is N+1,
  209	N1 is N+2,
  210	matching_(LCs,HCs,Cs,N1,Dec,Nout).
  211matching_([LC,'.',LC1|LCs],[HC,'.',HC1|HCs],[HC,'.'|Cs],N,Dec,Nout) :-  !,  % absorbing "."
  212	digit_match_(LC,LC1,HC,HC1),
  213	Dec is N+1,
  214	N1 is N+2,
  215	matching_([LC1|LCs],[HC1|HCs],Cs,N1,Dec,Nout).
  216matching_([C|LCs],[C|HCs],[C|Cs],N,Dec,Nout) :- !,  % matching digit
  217	digit_(C,_),
  218	N1 is N+1,
  219	matching_(LCs,HCs,Cs,N1,Dec,Nout).
  220matching_([LC,LC1|LCs],[HC,HC1|HCs],[HC|Cs],N,Dec,Nout) :- % match after rounding
  221	digit_match_(LC,LC1,HC,HC1), !,
  222	N1 is N+1,
  223	matching_([LC1|LCs],[HC1|HCs],Cs,N1,Dec,Nout).
  224matching_(_LCs,_HCs,[],N,Dec,N) :-                  % non-matching after '.'
  225	nonvar(Dec). 
  226
  227digit_match_(LC,LC1,HC,HC1) :-  % rounding test if first digits different
  228	digit_(LC,L), digit_(LC1,L1), digit_(HC,H), digit_(HC1,H1),
  229	H is (L+1) mod 10,
  230	L1 >= 5, H1 < 5.
  231
  232% char test for properly formatted number 
  233digit_(DC,D) :- atom_number(DC,D), integer(D), 0=<D,D=<9.
  234
  235%
  236%  enumerate integer and boolean intervals
  237%
  238enumerate(X) :-
  239	interval_object(X, integer, (L,H), _), !,
  240	between(L,H,X).             % gen values, constraints run on unification
  241enumerate(X) :- list(X), !,     % list of ..
  242	enumerate_list_(X).
  243enumerate(_X).                  % X not enumerable, skip it
  244
  245enumerate_list_([]).
  246enumerate_list_([X|Xs]) :-
  247	(integer(X) -> true ; enumerate(X)),  % optimization: X already an integer, skip it
  248	enumerate_list_(Xs).
  249
  250%
  251% Definition of "small" interval based on width and precision value
  252%   (defaults to clpBNR_default_precision)
  253%
  254small(V) :-
  255	current_prolog_flag(clpBNR_default_precision,P),
  256	small(V,P).
  257
  258small(V,P) :- 
  259	Err is 10.0**(-P),
  260	small_(V,Err).
  261
  262small_(V,Err) :- 
  263	getValue(V,(L,H)), !,
  264	chk_small(L,H,Err).
  265small_(L,Err) :-
  266	list(L),
  267	small_list_(L,Err).
  268
  269small_list_([],_).
  270small_list_([X|Xs],Err) :-
  271	small_(X,Err),
  272	small_list_(Xs,Err).
  273
  274chk_small(L,H,Err) :-                     % from CLIP?
  275	W is H-L,
  276	( W < Err                             % absolute error test  
  277	 -> true
  278	 ;  ErrH is Err*sqrt((L**2+H**2)/2),  % guaranteed to be a float
  279	    W < ErrH,                         % relative error test                 
  280	    ErrH \= 1.0Inf                    % overflow check
  281	).
  282
  283%
  284%  global_minimum(Exp,Z) - Z is an interval containing one or more global minimums of Exp.
  285%  global_maximum(Exp,Z) - as above but global maximums
  286%
  287%  Uses Moore-Skelboe algorithm documented in:
  288%  Ratschek & Rokne (June, 2007) "New Computer Methods for Global Optimization", 
  289%
  290%  The main complication is that evaluation of Exp may yield a false positive after 
  291%  fixed point iteration so any result must be "validated" either by a subsequent MS
  292%  iteration of by a final absolve (see continue/5).
  293%
  294global_maximum(Exp,Z) :-
  295	global_minimum(-Exp,NZ),
  296	constrain_(Z== -NZ).
  297global_maximum(Exp,Z,P) :-
  298	global_minimum(-Exp,NZ,P),
  299	constrain_(Z== -NZ).
  300
  301global_minimum(Exp,Z) :-
  302	current_prolog_flag(clpBNR_default_precision,P),
  303	global_minimum(Exp,Z,P).
  304global_minimum(Exp,Z,_P) :- ground(Exp), !,
  305	Z is Exp.
  306global_minimum(Exp,Z,P) :-
  307	global_optimum_(Exp,Z,P,false).
  308
  309% Copy of global_maximum & global_minimum with BindVars=true
  310global_maximize(Exp,Z) :-
  311	global_minimize(-Exp,NZ),
  312	constrain_(Z== -NZ).
  313global_maximize(Exp,Z,P) :-
  314	global_minimize(-Exp,NZ,P),
  315	constrain_(Z== -NZ).
  316
  317global_minimize(Exp,Z) :-
  318	current_prolog_flag(clpBNR_default_precision,P),
  319	global_minimize(Exp,Z,P).
  320global_minimize(Exp,Z,_P) :- ground(Exp), !,
  321	Z is Exp.
  322global_minimize(Exp,Z,P) :-
  323	global_optimum_(Exp,Z,P,true).
  324
  325global_optimum_(Exp,Z,P,BindVars) :-
  326	term_variables(Exp,Xs),                           % vars to search on
  327	{Z==Exp},                                         % Steps 1. - 4.
  328	box_([Z|Xs],[(Zl,Zh)|XVs]),                       % construct initial box
  329	iterate_MS(Z,Xs,P,Zl-(Zh,XVs),_ZTree,BindVars).   % and start iteration
  330
  331iterate_MS(Z,Xs,P,Zl-(Zh,XVs),ZTree,BindVars) :-
  332	continue_MS(Zl,Zh,P,Xs,XVs,False),                % Step 12., check termination condition
  333	widest_MS(Xs,XVs,Xf,XfMid), !,                    % Step 5., get midpoint of widest variable
  334	eval_MS(False,Z,Xs,XVs,Xf=< ::(XfMid,XfMid),V1),  % Step 6., 7. & 8. for V1
  335	tree_insert(ZTree,V1,ZTree1),                     % Step 10. for V1
  336	eval_MS(False,Z,Xs,XVs,Xf>= ::(XfMid,XfMid),V2),  % Step 6., 7. & 8. for V2
  337	tree_insert(ZTree1,V2,ZTree2),                    % Step 10. for V2
  338	select_min(ZTree2,NxtY,ZTreeY),                   % Step 9. and 11., remove Y from tree
  339	iterate_MS(Z,Xs,P,NxtY,ZTreeY,BindVars).          % Step 13.
  340iterate_MS(Z,Xs,_P,Zl-(Zh,XVs),_ZTree,BindVars) :-    % termination criteria (Step 12.) met
  341	Z::real(Zl,Zh),
  342	(BindVars -> optimize_vars_(Xs,XVs) ; true).      % optional minimizer narrowing
  343
  344optimize_vars_([],[]).
  345optimize_vars_([X|Xs],[(Xl,Xh)|XVs]) :-
  346	X::real(Xl,Xh),
  347	optimize_vars_(Xs,XVs).
  348
  349continue_MS(Zl,Zh,P,_,_,_) :-                    % w(Y) termination criteria
  350	Err is 10.0**(-P),
  351	\+chk_small(Zl,Zh,Err),                      % fail continue_MS if narrow enough
  352	!.
  353continue_MS(_,_,P,Xs,XVs,_) :-                   % test for false positive
  354	build_box_MS(Xs,XVs,T/T),
  355	SErr is 10.0**(-min(6,P+2)), % ?? heuristic
  356	simplesolveall_(Xs,SErr),                    % validate solution
  357	!, fail.                                     % no narrowing escapes
  358continue_MS(_Zl,_Zh,_,_,_XVs,false).             % continue with false positive
  359
  360% calculate resultant box and return it, original box left unchanged (uses global var) 
  361eval_MS(False,_,_,_,_,[]) :- False == false, !.  % if false positive return "fail" result
  362eval_MS(_,Z,Xs,XVs,XConstraint,_FV) :-           % Step 7., calculate F(V) and save
  363	nb_setval('clpBNR:eval_MS',[]),                      % [] means failure to evaluate
  364	buildConstraint_(XConstraint,T/T,Agenda),
  365	build_box_MS(Xs,XVs,Agenda),
  366	box_([Z|Xs],[(Zl,Zh)|NXVs]),                 % copy Z and Xs solution bounds
  367	nb_setval('clpBNR:eval_MS',Zl-(Zh,NXVs)),            % save solution in format for Z tree
  368	fail.                                        % backtack to unwind 
  369eval_MS(_,_,_,_,_,FV) :-
  370	nb_getval('clpBNR:eval_MS',FV).                      % retrieve solution ([] on failure)
  371
  372sandbox:safe_global_variable('clpBNR:eval_MS').
  373
  374build_box_MS([],[],Agenda) :-
  375	stable_(Agenda).
  376build_box_MS([X|Xs],[XNew|XVs],Agenda) :-
  377	getValue(X,XCurrent),
  378	updateValue_(XCurrent,XNew,X,1,Agenda,NextAgenda), !, %%% unnecessary cut??
  379	build_box_MS(Xs,XVs,NextAgenda).
  380
  381tree_insert(Tree,[],Tree) :-!.
  382tree_insert(MT,Data,Ntree) :- var(MT), !,
  383	Ntree = n(_L,_R,Data).
  384tree_insert(n(L,R,K-Data), NK-NData, Ntree) :- -1 is cmpr(NK,K), !,  % NK<K
  385	tree_insert(L, NK-NData, NewL),
  386	Ntree = n(NewL,R,K-Data).
  387tree_insert(n(L,R,K-Data), NK-NData, Ntree) :-  1 is cmpr(NK,K), !,  % NK>K
  388	tree_insert(R, NK-NData, NewR),
  389	Ntree = n(L,NewR,K-Data).
  390tree_insert(n(L,R,K-(U,Data)), K-(NU,NData), Ntree) :- -1 is cmpr(NU,U), !,  % K=NK, NU < U
  391	tree_insert(L, K-(NU,NData), NewL),
  392	Ntree = n(NewL,R,K-(U,Data)).
  393tree_insert(n(L,R,Data), NData, n(L,NewR,Data)) :-                   % K=NK, NU>=U
  394	Data = NData
  395	 -> NewR = R                      % if duplicate node, discard
  396	 ;  tree_insert(R, NData, NewR).  % else add to right branch
  397	
  398select_min(n(L,R,Min),Min,R) :- var(L), !, 
  399	nonvar(Min).  % fail if tree was empty
  400select_min(n(L,R,KData),Min,n(NewL,R,KData)) :-
  401	select_min(L,Min,NewL).
  402
  403% construct box from interval bounds
  404box_([],[]).
  405box_([X|Xs],[R|Rs]) :-
  406	getValue(X,R),
  407	box_(Xs,Rs).
  408
  409% find widest interval and its midpoint
  410widest_MS([X|Xs],[XV|XVs],Xf,XfMid) :-
  411	widest1_MS(Xs,XVs,XV,X,Xf,XfMid).
  412	
  413widest1_MS([],[],(L,H),Xf,Xf,XfMid) :-
  414	midpoint_(L,H,XfMid), !,              % must be "splittable" so
  415	-2 is cmpr(L,XfMid) + cmpr(XfMid,H).  % L < XfMid < H
  416widest1_MS([X|Xs],[(L,H)|XVs],(L0,H0),_X0,Xf,XfMid) :-
  417	1 is cmpr((H-L),(H0-L0)), !,  % (H-L) > (H0-L0)
  418	widest1_MS(Xs,XVs,(L,H),X,Xf,XfMid).
  419widest1_MS([_X|Xs],[_XV|XVs],W0,X0,Xf,XfMid) :-
  420	widest1_MS(Xs,XVs,W0,X0,Xf,XfMid).
  421
  422% Midpoint test, Step 11+:
  423% Discard all pairs (Z,z) from the list that satisfy F(C)<z where c = mid Y.
  424%midpoint_MS(L,H,M) :-  % L and H finite, non-zero ==> geometric/arithmetic mean
  425%	M is min(sqrt(abs(L)),abs(L)/2)*sign(L)+min(sqrt(abs(H)),abs(H)/2)*sign(H).
  426
  427%
  428%  splitsolve(Int) - joint search on list of intervals
  429%  simple split, e.g., no filtering on splutions, etc.
  430%
  431splitsolve(X) :-
  432	current_prolog_flag(clpBNR_default_precision,P),
  433	splitsolve(X,P).
  434
  435splitsolve(X,_P) :-
  436	number(X), !.                 % already a point value
  437splitsolve(X,P) :-
  438	interval(X), !,               % if single interval, put it into a list
  439	Err is 10.0**(-P),
  440	simplesolveall_([X],Err).
  441splitsolve(X,P) :-                % assumed to be a list
  442	flatten_list(X,[],XF),        % flatten before iteration
  443	Err is 10.0**(-P),
  444	simplesolveall_(XF,Err).
  445
  446% flatten list(s) using difference lists
  447flatten_list(V,Tail,[V|Tail])   :- var(V), !.
  448flatten_list([],Tail,Tail)      :- !.
  449flatten_list([H|T], Tail, List) :- !,
  450	flatten_list(H, HTail, List),
  451	flatten_list(T, Tail, HTail).
  452flatten_list(N,Tail,[N|Tail]).
  453
  454simplesolveall_(Xs,Err) :-
  455	select_wide_(Xs,_,X),
  456	interval_object(X, Type, V, _),
  457	split_choices_(Type,X,V,Err,Choices),
  458	!,
  459	xpsolve_choice(Choices),  % from solve/1, generate alternatives
  460	simplesolveall_(Xs,Err).
  461simplesolveall_(_Xs,_Err).  % terminated
  462
  463select_wide_([X],_,X) :- !.        % select last remaining element
  464select_wide_([X1,X2|Xs],D1,X) :-   % compare widths and discard one interval
  465	(var(D1) -> delta(X1,D1) ; true),
  466	delta(X2,D2),
  467	(D1 >= D2
  468	 -> select_wide_([X1|Xs],D1,X)
  469	 ;  select_wide_([X2|Xs],D2,X)
  470	).
  471
  472split_choices_(real,X,V,Err,split(X,::(SPt,SPt))) :- !,
  473	splitinterval_real_(V,SPt,Err).            % from solve/1, but splits anywhere
  474split_choices_(integer,X,V,_Err,Cons) :-
  475	splitinterval_(integer,X,V,_,Cons).        % from solve/1
  476
  477%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  478%
  479%					absolve( X )
  480%					absolve( X, Precision)
  481%
  482%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  483%
  484%  absolve( X ), for X an interval or list of intervals, is intended 
  485%  solely to trim up the boundaries of what is essentially a single
  486%  (non-point)  solution to a problem. Hence absolve- unlike the rest
  487%  of the solve family -  is deterministic!
  488%	 The strategy used in absolve( derived from the old V 3 solve) is to 
  489%  work in from the edges of the interval ("nibbling away") until you
  490%  cannot go nay farther, then reduce step size and resume nibbling.
  491%  Control parameters limit the number of successive halvings of the 
  492%  step size, which is initially half the interval width. 
  493%  		Note that absolve and solve each abstract a different part of
  494%  of the strategy used in the solve in BNRP V. 3. In this sense, the
  495%  combination: " solve(X), absolve(X) "  { in that order } does something
  496%  like what "solve(X)"did under the old system.
  497
  498absolve( X ):-
  499	current_prolog_flag(clpBNR_default_precision,P),
  500	absolve(X,P),!.
  501
  502absolve(X, _ ):- number(X), !.
  503absolve(X, Limit ):- interval_object(X, Type, Val, _), !,  % interval
  504	delta_(Type,Val,Delta),
  505	% if bound is already a solution avoid the work
  506	(not(not(lower_bound(X))) -> true ; absolve_l(X,Type,Delta,1,Limit)),
  507	(not(not(upper_bound(X))) -> true ; absolve_r(X,Type,Delta,1,Limit)).
  508
  509absolve([],_).		% extends to lists
  510absolve([X|Xs],Lim):- absolve(X,Lim),!, absolve(Xs,Lim).
  511
  512delta_(integer,(L,U),D) :- D is U div 2 - L div 2.
  513delta_(real,   (L,U),D) :- D is U/2 - L/2.
  514
  515absolve_l(X, Type, DL, NL, Limit):- NL<Limit, % work on left side
  516	getValue(X,(LB1,UB1)), 
  517	trim_point_(NL,NL1,Type,Limit,DL,DL1),    % generates trim points
  518	Split is LB1 + DL1,
  519	LB1 < Split, Split < UB1,                 % in range, not endpoint
  520	not(constrain_(X=< ::(Split,Split))),!,
  521	constrain_(X>= ::(Split,Split)),          % so X must be >
  522	absolve_l(X,Type, DL1, NL1, Limit).
  523absolve_l(_,_,_,_,_).                         % final result
  524         
  525absolve_r(X, Type, DU, NU, Limit):- NU<Limit, % work on right side
  526	getValue(X,(LB1,UB1)), 
  527	trim_point_(NU,NU1,Type,Limit,DU,DU1),    % generates trim points
  528	Split is UB1 - DU1,
  529	LB1 < Split, Split < UB1,                 % in range, not endpoint
  530	not(constrain_(X>= ::(Split,Split))),!,
  531	constrain_(X=< ::(Split,Split)),          % so X must be <
  532	absolve_r(X,Type, DU1, NU1,Limit).
  533absolve_r(_,_,_,_,_).                         % final result
  534
  535trim_point_( N,N, _Type, _Limit, Delta, Delta).
  536trim_point_( N,M, integer, Limit, Delta, Result):- N<Limit,N1 is N + 1,
  537       D is  Delta div 2,
  538       trim_point_(N1,M, integer, Limit,D, Result).
  539trim_point_( N,M, real, Limit, Delta, Result):- N<Limit,N1 is N + 1,
  540       D is  Delta/2,
  541       trim_point_(N1,M,real, Limit,D, Result).
  542
  543%
  544%  solve(Int) - joint search on list of intervals
  545%
  546solve(X) :-
  547	current_prolog_flag(clpBNR_default_precision,P),
  548	solve(X,P).
  549
  550solve(X,P) :-
  551	interval(X), !,               % if single interval, put it into a list
  552	solve([X],P).
  553solve(X,_P) :-
  554	number(X), !.                 % already a point value
  555solve(X,P) :-                     % assume list
  556	Err is 10.0**(-(P+1)),        % convert digits of precision to normalized error value 
  557	xpsolveall_(X,Err).
  558
  559xpsolveall_([],_Err) :- !.
  560xpsolveall_(Xs,Err) :-
  561	xpsolve_each_(Xs,Us,Err),     % split once and collect successes
  562	xpsolveall_(Us,Err).          % continue to split remaining
  563
  564xpsolve_each_([],[],_Err).
  565xpsolve_each_([X|Xs],[X|Us],Err) :-
  566	interval_object(X,Type,V,_),          % get interval type and value
  567	splitinterval_(Type,X,V,Err,Choices), % split interval
  568	!,
  569	xpsolve_choice(Choices),              % create choice(point)
  570	xpsolve_each_(Xs,Us,Err).             % split others in list
  571xpsolve_each_([X|Xs],Us,Err) :-           % avoid unfreeze overhead if [] unified in head
  572	X==[], !,                             % end of nested listed, discard
  573	xpsolve_each_(Xs,Us,Err).             % split remaining
  574xpsolve_each_([X|Xs],[U|Us],Err) :-
  575	list(X), !,                           % nested list
  576	xpsolve_each_(X,U,Err),               % split nested list
  577	xpsolve_each_(Xs,Us,Err).             % then others in main list
  578xpsolve_each_([_X|Xs],Us,Err) :-
  579	xpsolve_each_(Xs,Us,Err).             % split failed or already a number, drop interval from list, and keep going
  580
  581xpsolve_choice(split(X,SPt)) :- constrain_(X =< SPt).  % avoid meta call
  582xpsolve_choice(split(X,SPt)) :- constrain_(SPt =< X).
  583xpsolve_choice(split_integer(X,SPt)) :- constrain_(X =< SPt).
  584xpsolve_choice(split_integer(X,SPt)) :- constrain_(SPt < X).
  585xpsolve_choice(enumerate(X)) :- enumerate(X).
  586
  587%
  588% try to split interval - fails if unsplittable (too narrow)
  589%
  590splitinterval_(real,X,V,Err,split(X,::(SPt,SPt))) :-  % split on point value
  591	splitinterval_real_(V,Pt,Err),          % initial guess
  592	split_real_(X,V,Pt,Err,SPt).
  593splitinterval_(integer,X,V,_,Cons) :- 
  594	V = (L,H),
  595	( H-L < 15
  596	 -> Cons = enumerate(X)                              % small enough to enumerate
  597	 ;  splitinterval_integer_(V,Pt),                    % splittable at Pt
  598	    Cons = split_integer(X,Pt)                       % success
  599	).
  600%splitinterval_(boolean,X,Err,Choices) :-
  601%	splitinterval_(integer,X,Err,Choices).
  602
  603%  split a real interval
  604split_real_(X,_,Pt,_,Pt) :-                % Pt not in solution space, split here
  605	X = Pt -> fail ; !.  % not({X==Pt}).
  606split_real_(X,(L,_H),Pt,Err,NPt) :-        % Pt in current solution space, try lower range
  607	split_real_lo(X,L,Pt,NPt,Err), !.
  608split_real_(X,(_L,H),Pt,Err,NPt) :-        % Pt in current solution space, try upper range
  609	split_real_hi(X,Pt,H,NPt,Err).
  610
  611split_real_lo(X,L,Pt,NPt,Err) :-           % search lower range for a split point 
  612	splitinterval_real_((L,Pt),SPt,Err),
  613	(X\=SPt -> NPt=SPt ; split_real_lo(X,L,SPt,NPt,Err)).
  614
  615split_real_hi(X,Pt,H,NPt,Err) :-           % search upper range for a split point 
  616	splitinterval_real_((Pt,H),SPt,Err),
  617	(X\=SPt -> NPt=SPt ; split_real_hi(X,SPt,H,NPt,Err)).
  618
  619%
  620% splitinterval_integer_ and splitinterval_real_
  621%
  622splitinterval_integer_((L,H),0) :-
  623	L < 0, H > 0, !.                  % contains 0 but not bounded by 0 
  624splitinterval_integer_((-1.0Inf,H),Pt) :-  !,  % infinite lower bound, integers unbounded
  625	finite_interval(integer, (Pt,_)), % split at finite integer min if in range
  626	H > Pt.
  627splitinterval_integer_((L,1.0Inf), Pt) :-  !,  % infinite upper bound, integers unbounded
  628	finite_interval(integer, (_,Pt)), % split at finite integer max if in range
  629	L < Pt.
  630splitinterval_integer_((L,H),Pt) :-   % all positive or negative (including zero)  
  631	Pt is (L div 2) + (H div 2).      % avoid overflow
  632
  633splitinterval_real_((L,H),0,E) :-     % if interval contains 0, split on (precisely) 0.
  634	L < 0, H > 0, !,                  % cut in case error criteria fails
  635	(H-L) > E.                        % fail if width is less than error criteria, overflow is OK
  636splitinterval_real_((-1.0Inf,H),Pt,_) :-  % neg. infinity to H=<0
  637	!,  % if following overflows, split failed
  638	Pt is float(H*10-1),               % subtract 1 in case H is 0. (-1, -11, -101, -1001, ...)
  639	Pt > -1.0Inf.                      % Pt must be finite
  640splitinterval_real_((L,1.0Inf),Pt,_) :-   % L>=0 to pos. infinity
  641	!,  % if following overflows, split failed
  642	Pt is float(L*10+1),               % add 1 in case L is 0. (1, 11, 101, 1001, ...)
  643	Pt < 1.0Inf.                       % Pt as float must be finite
  644splitinterval_real_((L,H),Pt,E) :-     % finite L,H, positive or negative but not split, Pt\=0.
  645	\+ chk_small(L,H,E),               % only split if not small 
  646	splitMean_(L,H,Pt), !,
  647	L < Pt,Pt < H.                     % split point must be between L and H
  648
  649%splitinterval_real_([L,H],Pt,E) :-
  650%	writeln('FAIL'([L,H],Pt,E)),fail.
  651
  652% (approx.) geometric mean of L and H (same sign)
  653splitMean_(L,H,M) :-  L >=0, !,  % positive range
  654	(L=:=0 -> M is min(H/2, sqrt( H)) ; M is sqrt(L)*sqrt(H)).
  655splitMean_(L,H,M) :-             % negative range
  656	(H=:=0 -> M is max(L/2,-sqrt(-L)) ; M is -sqrt(-L)*sqrt(-H)). 
  657
  658
  659%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  660%
  661%					partial_derivative(E, X, D)
  662%
  663%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  664%
  665%  Perform symbolic differentiation followed by simplify/2:  D = dE/dX
  666%
  667%  located in utilities as it's a common requirement for applying series expansion
  668%  to improve convergance. Avoids necessity of exposing simplify/2.
  669%  Note that the set of rules only covers functions supported by {}.
  670%
  671partial_derivative(F,X,DFX) :-
  672	pd(F,X,DD),
  673	simplify(DD,DFX).
  674	%	simplify(DD,DX),
  675	%	simplify(DX,DFX).  % post re-simplify for insurance??
  676
  677pd(F,X,DD) :-
  678	( var(F)
  679	 -> ( F==X -> DD=1 ; DD=0)
  680	 ;  ( atomic(F)
  681	     -> DD=0
  682	     ; pd_f(F,X,DD)
  683	    )
  684	).
  685
  686:- style_check(-singleton).        % for pd_f
  687
  688pd_f(-U,X,DX) :-                   % DX = -DU
  689	pd(U,X,DU),
  690	pd_minus(DU,DX).
  691
  692pd_f(U+V,X,DX) :-                  % DX = DU+DV
  693	pd(U,X,DU), pd(V,X,DV),
  694	pd_add(DU,DV,DX).
  695
  696pd_f(U-V,X,DX) :-                  % DX = DU-DV
  697	pd(U,X,DU), pd(V,X,DV),
  698	pd_sub(DU,DV,DX).
  699
  700pd_f(U*V,X,DX) :-                  % DX = V*DU+U*DV
  701	pd(U,X,DU), pd(V,X,DV), 
  702	pd_mul(V,DU,VDU),
  703	pd_mul(U,DV,UDV),
  704	pd_add(VDU,UDV,DX).
  705
  706pd_f(U/V,X,DX) :-                  % DX = (V*DU-U*DV)/(V**2)
  707	pd(U,X,DU), pd(V,X,DV),
  708	pd_mul(V,DU,VDU),
  709	pd_mul(U,DV,UDV),
  710	pd_sub(VDU,UDV,DXN),
  711	pd_pow(V,2,DXD),
  712	pd_div(DXN,DXD,DX).
  713
  714pd_f(U**N,X,DX) :-                 % DX = (N*U**(N-1))*DU
  715	pd(U,X,DU),
  716	pd_sub(N,1,N1),  %N1 is N-1,
  717	pd_pow(U,N1,UN1),
  718	pd_mul(N,UN1,DX1),
  719	pd_mul(DX1,DU,DX).
  720
  721pd_f(log(U),X,DX) :-               % DX = DU/U
  722	pd(U,X,DU),
  723	pd_div(DU,U,DX).
  724
  725pd_f(exp(U),X,DX) :-               % DX = exp(U)*DU
  726	pd(U,X,DU),
  727	pd_exp(U,ExpU),
  728	pd_mul(ExpU,DU,DX).
  729
  730pd_f(sqrt(U),X,DX) :-              % DX = DU/(2*sqrt(U))
  731	pd(U,X,DU),
  732	pd_sqrt(U,SqrtU),
  733	pd_mul(2,SqrtU,DXD),
  734	pd_div(DU,DXD,DX).
  735
  736pd_f(sin(U),X,DX) :-               % DX = cos(U)*DU
  737	pd(U,X,DU),
  738	pd_cos(U,CosU),
  739	pd_mul(CosU,DU,DX).
  740
  741pd_f(cos(U),X,DX) :-               % DX = -sin(U)*DU
  742	pd(U,X,DU),
  743	pd_sin(U,SinU),
  744	pd_mul(SinU,DU,DSU),
  745	pd_minus(DSU,DX).
  746
  747pd_f(tan(U),X,DX) :-               % DX = DU/cos(U)**2
  748	pd(U,X,DU),
  749	pd_cos(U,CosU),
  750	pd_pow(CosU,2,DXD),
  751	pd_div(DU,DXD,DX).
  752
  753pd_f(asin(U),X,DX) :-              % DX = DU/sqrt(1-U**2)
  754	pd(U,X,DU),
  755	pd_pow(U,2,USQ),
  756	pd_sub(1,USQ,USQ1),
  757	pd_sqrt(USQ1,DXD),
  758	pd_div(DU,DXD,DX).
  759
  760pd_f(acos(U),X,DX) :-              % DX = -DU/sqrt(1-U**2) = -D(asin(U))
  761	pd(asin(U),X,DasinX),
  762	pd_minus(DasinX,DX).
  763
  764pd_f(atan(U),X,DX) :-              % DX = DU/(1+U**2)
  765	pd(U,X,DU),
  766	pd_pow(U,2,USQ),
  767	pd_add(1,USQ,USQ1),
  768	pd_div(DU,USQ1,DX).
  769
  770
  771% optimizations
  772% also facilitates compiled arithmetic, avoids using catch/3 for instantiation errors
  773pd_minus(DU,DX)  :- ground(DU) -> DX is -DU  ; DX = -DU.
  774
  775pd_add(DU,DV,DV) :- DU==0, !.
  776pd_add(DU,DV,DU) :- DV==0, !.
  777pd_add(DU,DV,DX) :- ground((DU,DV)) -> DX is DU+DV ; DX = DU+DV.
  778
  779pd_sub(DU,DV,DX) :- DU==0, !, pd_minus(DV,DX).
  780pd_sub(DU,DV,DU) :- DV==0, !.
  781pd_sub(DU,DV,DX) :- ground((DU,DV)) -> DX is DU-DV ; DX = DU-DV.
  782
  783pd_mul(DU,DV,0)  :- DU==0, !.
  784pd_mul(DU,DV,0)  :- DV==0, !.
  785pd_mul(DU,DV,DV) :- DU==1, !.
  786pd_mul(DU,DV,DU) :- DV==1, !.
  787pd_mul(DU,DV,DX) :- DU== -1, !, pd_minus(DV,DX).
  788pd_mul(DU,DV,DX) :- DV== -1, !, pd_minus(DU,DX).
  789pd_mul(DU,DV,DX) :- ground((DU,DV)) -> DX is DU*DV ; DX = DU*DV.
  790
  791pd_div(DU,DV,0)  :- DU==0, !.
  792pd_div(DU,DV,0)  :- DV==0, !, fail.
  793pd_div(DU,DV,DU) :- DV==1, !.
  794pd_div(DU,DV,DU) :- DV== -1, !, pd_minus(DU,DX).
  795pd_div(DU,DV,DX) :- ground((DU,DV)) -> DX is DU/DV ; DX = DU/DV.
  796
  797pd_pow(DU,DV,0)  :- DU==0, !.
  798pd_pow(DU,DV,1)  :- DV==0, !.
  799pd_pow(DU,DV,1)  :- DU==1, !.
  800pd_pow(DU,DV,DU) :- DV==1, !.
  801pd_pow(DU,DV,DX) :- ground((DU,DV)) -> DX is DU**DV ; DX = DU**DV.
  802
  803pd_exp(DU,DX)    :- ground(DU) -> DX is exp(DU)  ; DX = exp(DU).
  804
  805pd_sqrt(DU,DX)   :- ground(DU) -> DX is sqrt(DU) ; DX = sqrt(DU).
  806
  807pd_sin(DU,DX)    :- ground(DU) -> DX is sin(DU)  ; DX = sin(DU).
  808
  809pd_cos(DU,DX)    :- ground(DU) -> DX is cos(DU)  ; DX = cos(DU).
  810
  811:- style_check(+singleton).  % exit pd_f