1/*	The MIT License (MIT)
    2 *
    3 *	Copyright (c) 2022-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:- module(clpBNR_toolkit,       % SWI module declaration
   24	[
   25	iterate_until/3,            % general purpose iterator
   26	mid_split_one/1,            % contractor to split largest interval at midpoint
   27	mid_split/1,                % contractor to split an interval at midpoint
   28	taylor_contractor/2,        % build cf_contractor based on Taylor expansion
   29	taylor_merged_contractor/2, % build merged Taylor cf_contractor from list of equations
   30	cf_contractor/2,            % execute cf_contractor
   31	cf_solve/1, cf_solve/2,     % a solve predicate for centre form contractors
   32
   33	integrate/3, integrate/4,   % simple numerical integration
   34	boundary_values/2, boundary_values/3, boundary_values/4,  % solve boundary value problems
   35
   36	lin_minimum/3,              % find minimum of linear problem using library(simplex)
   37	lin_maximum/3,              % find maximum of linear problem using library(simplex)
   38	lin_minimize/3,             % lin_minimum/3 plus bind vars to solution minimizers
   39	lin_maximize/3,             % lin_maximum/3 plus bind vars to solution maximizers
   40
   41	local_minima/1,             % apply KT constraints for objective function expression (OFE)
   42	local_maxima/1,             % semantically equivalent to local_minima/1
   43	local_minima/2,             % apply KT constraints for minima with constraints
   44	local_maxima/2              % apply KT constraints for maxima with constraints
   45	]).

clpBNR_toolkit: Toolkit of various utilities used for solving problems with clpBNR

CLP(BNR) (library(clpBNR)) is a CLP over the domain of real numbers extended with ±∞. This module contains a number of useful utilities for specific problem domains like the optimization of linear systems, enforcing local optima conditions, and constructing centre form contractors to improve performance (e.g., Taylor extensions of constraints). For more detailed discussion, see A Guide to CLP(BNR) (HTML version included with this pack in directory docs/).

Documentation for exported predicates follows. The "custom" types include:

   56:- use_module(library(apply),[maplist/3]).   57:- use_module(library(apply_macros)).  % compiler support for `maplist`, helps a bit
   58:- use_module(library(clpBNR)).   59:- use_module(library(simplex)).   60
   61% sandboxing for SWISH
   62:- multifile(sandbox:safe_primitive/1).   63
   64% messages for noisy failure
   65fail_msg_(FString,Args) :-
   66	debug(clpBNR,FString,Args), fail.
   67	
   68:- set_prolog_flag(optimise,true).  % for arithmetic, this module only
 iterate_until(+Count:integer, Test, Goal) is nondet
Succeeds when Goal succeeds; otherwise fails. Goal will be called recursively up to Count times as long as Test succeeds Example using this predicate for simple labelling using Test small/2 and Goal midsplit/1 :
?- X::real(-1,1),iterate_until(10,small(X,0),mid_split(X)),format("X = ~w\n",X),fail;true.
X = _6288{real(-1,-1r2)}
X = _6288{real(-1r2,0)}
X = _6288{real(0,1r2)}
X = _6288{real(1r2,1)}
true.

The specific intended use case is to provide an iterator for meta-contractors such as the centre-form contractor such as midsplit/1 (example above) or as constructed by taylor_contractor/2 as in:

?- X::real,taylor_contractor({X**4-4*X**3+4*X**2-4*X+3==0},T),
iterate_until(50,small(X),(T,mid_split_one([X]))),format("X = ~w\n",X),fail;true.
X = _150{real(0.999999999926943,1.00000000007306)}
X = _150{real(2.999999999484828,3.0000000005152105)}
true.

(Aside: For some problems, solving with Taylor contractors can be a faster and more precise alternative to clpBNR:solve/1.) */

   92%
   93% General purpose iterator: execute Goal a maximum of N times or until Test succeeds
   94%
   95iterate_until(N,Test,Goal) :- N>0, !,
   96	Goal,
   97	N1 is N-1,
   98	(Test
   99	 -> true
  100	  ; iterate_until(N1,Test,Goal)
  101	).
  102iterate_until(_N,_,_).  % non-positive N --> exit
  103
  104sandbox:safe_meta(clpBNR_toolkit:iterate_until(_N,Test,Goal), [Test, Goal]).
 mid_split_one(+Xs:numeric_list) is nondet
Succeeds splitting the widest interval in Xs, a list of intervals; fails if Xs is not a list of intervals. See mid_split for details of interval splitting for this predicate.
See also
- mid_split/1 */
  113mid_split_one(Xs) :-
  114	select_split(Xs,X),  % select largest interval with largest width
  115	mid_split(X).        % split it
 mid_split(X:numeric) is nondet
Succeeds if X is an interval that can be split at its midpoint narrowing X to it's lower "half"; on backtracking X is constrained to the upper half; fails if X is not a numeric. If X is an interval, defined as:
mid_split(X) :-
        M is midpoint(X),
        ({X=<M} ; {M=<X}).

Note that mid_split succeeds if X is a number, but doesn't do anything.

Use clpBNR:small as a pre-test to avoid splitting intervals which are already small enough.

See also
- clpBNR:small/1 */
  131mid_split(X) :-
  132	(number(X)                    % optimise number case
  133	 -> true
  134	 ;  (small(X)
  135	     -> true
  136	     ;  midpoint(X,M),       % fails if not an interval
  137	        ({X=<M} ; {M=<X})    % possible choicepoint, Note: can split on solution leaving CP
  138	    )
  139	).
  140%
  141% select interval with largest width
  142%
  143select_split([X],X) :- !.       % select last remaining element
  144select_split([X1,X2|Xs],X) :-   % compare widths and discard one interval
  145	delta(X1,D1),
  146	delta(X2,D2),
  147	(D1 >= D2
  148	 -> select_split([X1|Xs],X)
  149	 ;  select_split([X2|Xs],X)
  150	).
 cf_contractor(Xs:interval_list, As:interval_list) is semidet
Succeeds if each interval in As can be unified with the midpoints of the respective interval in Xs; otherwise fails. This predicate executes one narrowing step of the centre form contractor such as that generated by taylor_contractor. In normal usage, a direct call to cf_contractor does appear; instead use cf_contractor or in a Goal for iterate_until/3.
See also
- taylor_contractor/2, cf_solve/1, iterate_until/3 */
  159%
  160% centred form contractor
  161%
  162% Bind the values of As to the midpoints of Xs. To support repetitive application 
  163% of the contractor (required by the iterator), the contractor should not permanently 
  164% bind anything so findall/3 will be used to achieve this "forward checking" 
  165% (as suggested in [CLIP]). After the call to findall, the bounds of the resulting list
  166% of narrowed domains (XDs) are then applied to Xs.
  167%
  168% This contractor can be used with any "centred form", e.g., Newton or Krawczyk, since it
  169% only depends on intervals and their midpoints, hence its name `cf_contractor`. The 
  170% details which distinguish the variety of centred form are built into the variables' 
  171% constraints.
  172%
  173cf_contractor(Xs,As) :-
  174	findall(Ds,(maplist(bind_to_midpoint,Xs,As),maplist(cf_domain,Xs,Ds)),[XDs]),
  175	maplist(set_domain,Xs,XDs).
  176	
  177bind_to_midpoint(X,A) :- A is float(midpoint(X)).
  178
  179cf_domain(X,D) :- 
  180	number(X) -> D = X ; domain(X,D).  % in case X narrowed to a point
  181	
  182set_domain(X,D) :- 
  183	number(D) -> X = D ; X::D.
 cf_solve(+Contractor) is nondet
Equivalent to cf_solve/2 using default precision. */
 cf_solve(+Contractor, +Precision:integer) is nondet
Succeeds if a solution can be found for all variables in the centre form contractor, Contractor, where the resultant domain of any variable is narrower than the limit specified by Precision (for cf_solve/1, default Precision is number of digits as defined by the environment flag clpBNR_default_precision); otherwise fails.

This is done by using iterate_until/3 limited to a count determined by the flag clpBNR_iteration_limit. Examples:

?- X::real, taylor_contractor({X**4-4*X**3+4*X**2-4*X+3==0},T), cf_solve(T).
T = cf_contractor([X], [_A]),
X:: 1.000000000...,
_A::real(-1.0Inf, 1.0Inf) ;
T = cf_contractor([X], [_A]),
X:: 3.00000000...,
_A::real(-1.0Inf, 1.0Inf) ;
false.

?- taylor_contractor({2*X1+5*X1**3+1==X2*(1+X2), 2*X2+5*X2**3+1==X1*(1+X1)},T), cf_solve(T).
T = cf_contractor([X2, X1], [_A, _B]),
X1:: -0.42730462...,
X2:: -0.42730462...,
_B::real(-1.0Inf, 1.0Inf),
_A::real(-1.0Inf, 1.0Inf) ;
false.
See also
- taylor_contractor/2, cf_contractor/2 */
  217cf_solve(T) :-
  218    current_prolog_flag(clpBNR_default_precision,P),
  219    cf_solve(T,P).
  220cf_solve(cf_contractor(Xs,As),P) :-
  221    current_prolog_flag(clpBNR_iteration_limit,L),
  222    Count is L div 10,          % heuristic - primitive iteration limit/10    
  223    cf_iterate_(Count,Xs,As,P).
  224
  225cf_iterate_(Count,Xs,As,P) :-
  226	Count > 0,
  227	\+ small(Xs,P),             % at least one var not narrow enough
  228	!, 
  229	cf_contractor(Xs,As),       % execute contractor
  230	select_split(Xs,X),         % select widest
  231	(small(X,P)                 % still wide enough to split?
  232	 -> true                    % no, we're done 
  233	  ; cf_split(X,Pt),         % yes, split it
  234	    ({X=<Pt} ; {Pt=<X}),
  235	    Count1 is Count-1,
  236	    cf_iterate_(Count1,Xs,As,P)  % and iterate
  237	).
  238cf_iterate_(_,_,_,_).           % done (Count=<0 or all small Xs)
  239
  240cf_split(X,Pt) :-
  241	range(X,[L,H]),
  242	cf_split_point(L,H,X,M),    % modest attempt to find non-solution
  243	!,
  244	(M = 1.5NaN, L<0,0<H -> Pt = 0 ; Pt = M).  % (-inf,inf) case: if NaN and spans 0, use 0
  245
  246cf_split_point(L,H,X,M) :-	M is L/2 + H/2, \+ X = M.          % midpoint, not a solution
  247cf_split_point(L,H,X,M) :-	M is 0.625*L + 0.375*H, \+ X = M.  % 0.375*width, not a solution
  248cf_split_point(L,H,_,M) :-	M is 0.375*L + 0.625*H.            % 0.625*width, may be a solution
 taylor_contractor(+Constraints, -Contractor) is semidet
Succeeds if a centre form contractor Contractor can be generated from one or more multivariate equality (== or =:=) constraints Constraints; otherwise fails. Example:
?- taylor_contractor({X**4-4*X**3+4*X**2-4*X+3==0},T).
T = cf_contractor([X], [_A]),
X::real(-1.509169756145379, 4.18727500493995),
_A::real(-1.0Inf, 1.0Inf).

Use the contractor with cf_solve to search for solutions, as in:

?- X::real,taylor_contractor({X**4-4*X**3+4*X**2-4*X+3==0},T), cf_solve(T).
T = cf_contractor([X], [_A]),
X:: 1.000000000...,
_A::real(-1.0Inf, 1.0Inf) ;
T = cf_contractor([X], [_A]),
X:: 3.00000000...,
_A::real(-1.0Inf, 1.0Inf) ;
false.

Multiple equality constraints are supported, as in this example of the Broyden banded problem (N=2):

?- taylor_contractor({2*X1+5*X1**3+1==X2*(1+X2), 2*X2+5*X2**3+1==X1*(1+X1)},T), cf_solve(T).
T = cf_contractor([X2, X1], [_A, _B]),
X1:: -0.42730462...,
X2:: -0.42730462...,
_B::real(-1.0Inf, 1.0Inf),
_A::real(-1.0Inf, 1.0Inf) ;
false.

Centre form contractors can converge faster than the general purpose builtin fixed point iteration provided by solve/1.

See also
- cf_solve/1 */
  285%
  286% build a cf_contractor for a multivariate expression based on Taylor expansion
  287%
  288taylor_contractor({E1=:=E2},CF) :-
  289	taylor_contractor({E1==E2},CF).
  290taylor_contractor({E1==E2},cf_contractor(Xs,As)) :-
  291	Exp=E1-E2,
  292	term_variables(Exp,Xs),              % original arguments, bound to TXs on call
  293	make_EQ_(Exp,TEQ),                   % original constraint with arguments
  294	% build constraint list starting with Z's and ending with TEQ and DEQ ()
  295	T::real(0,1),
  296	make_As_and_Zs_(Xs,T,As,Zs,Cs,[TEQ,DEQ]),  % T per Z
  297	% now build Taylor constraint, DEQ
  298	copy_term_nat(Exp,AExp),              % copy of original constraint with  As
  299	term_variables(AExp,As),
  300	sum_diffs(Xs, As, Zs, Zs, Exp, AExp, DEQ),  % add on D(Z)'s'
  301	% make any vars in original equation and contractor arguments finite real intervals
  302	!,
  303	Xs::real,  % all vars are intervals
  304	{Cs}.      % apply constraints
  305taylor_contractor({Es},CF) :-
  306	taylor_merged_contractor({Es},CF),  % list or sequence	
  307	!.
  308taylor_contractor(Eq,_) :-
  309	fail_msg_('Invalid constraint for Taylor contractor: ~w',[Eq]).
  310
  311make_As_and_Zs_([],_,[],[],Tail,Tail).
  312make_As_and_Zs_([X|Xs],T,[A|As],[Z|Zs],[Z==A+T*(X-A)|CZs],Tail) :-
  313	make_As_and_Zs_(Xs,T,As,Zs,CZs,Tail).
  314
  315sum_diffs([], [], [], _AllZs, _Exp, ExpIn, EQ) :- make_EQ_(ExpIn,EQ).
  316sum_diffs([X|Xs], [A|As], [Z|Zs], AllZs, Exp, AExp, DEQ) :-
  317	copy_term_nat(Exp,NExp),        % copy expression and replace Xs by Zs
  318	term_variables(NExp,AllZs),
  319	partial_derivative(NExp,Z,DZ),  % differentiate wrt. Z and add to generated expression
  320	sum_diffs(Xs, As, Zs, AllZs, Exp, AExp+DZ*(X-A), DEQ).
  321
  322% map expression Exp to an equation equivalent to Exp==0 with numeric RHS
  323make_EQ_(Exp,LHS==RHS) :-    % turn expression into equation equivalent to Exp==0.
  324	make_EQ_(Exp,LHS,RHS).
  325	
  326make_EQ_(E,E,0)         :- var(E), !.
  327make_EQ_(X+Y,X,SY)      :- number(Y), !, SY is -Y.
  328make_EQ_(X-Y,X,Y)       :- number(Y), !.
  329make_EQ_(X+Y,Y,SX)      :- number(X), !, SX is -X.
  330make_EQ_(X-Y,SY,SX)     :- number(X), !, SX is -X, negate_sum_(Y,SY).
  331make_EQ_(X+Y,LHS+Y,RHS) :- !, make_EQ_(X,LHS,RHS).
  332make_EQ_(X-Y,LHS-Y,RHS) :- !, make_EQ_(X,LHS,RHS).
  333make_EQ_(E,E,0).        % default (non +/- subexpression)
  334
  335negate_sum_(Y,-Y) :- var(Y), !.
  336negate_sum_(X+Y,NX-Y) :- !, negate_sum_(X,NX).
  337negate_sum_(X-Y,NX+Y) :- !, negate_sum_(X,NX).
  338negate_sum_(E,-E).
 taylor_merged_contractor(+Constraints, -Contractor) is semidet
Succeeds if a merged centre form contractor Contractor can be generated from each equality (== or =:=) constraint in Constraints; otherwise fails.
deprecated
-
  • use taylor_contractor/2 instead */
  347%
  348% build a cf_contractor by merging a list of cf_contractor's
  349%
  350taylor_merged_contractor({Es},T) :-
  351	cf_list(Es,Ts),
  352	cf_merge(Ts,T).
  353
  354cf_list([],[]) :- !.
  355cf_list([C|Cs],[CF|CFs]) :- !,
  356	cf_list(C, CF),
  357	cf_list(Cs,CFs).
  358cf_list((C,Cs),[CF|CFs]) :-  !,
  359	cf_list(C, CF),
  360	cf_list(Cs,CFs).
  361cf_list(C,CF) :-
  362	taylor_contractor({C},CF).
  363
  364cf_merge(CFs,CF) :- cf_merge(CFs,cf_contractor([],[]),CF).
  365
  366cf_merge([],CF,CF).
  367cf_merge([CF|CFs],CFIn,CFOut) :-
  368	cf_merge(CF,CFIn,CFNxt),
  369	cf_merge(CFs,CFNxt,CFOut).	
  370cf_merge(cf_contractor(Xs,As),cf_contractor(XsIn,AsIn),cf_contractor(XsOut,AsOut)) :-
  371	cf_add(Xs,As,XsIn,AsIn,XsOut,AsOut).
  372
  373cf_add([],[],Xs,As,Xs,As).
  374cf_add([X|Xs],[A|As],XsIn,AsIn,XsOut,AsOut) :-
  375	var_existing(XsIn,AsIn,X,A), !,
  376	cf_add(Xs,As,XsIn,AsIn,XsOut,AsOut).
  377cf_add([X|Xs],[A|As],XsIn,AsIn,XsOut,AsOut) :-
  378	cf_add(Xs,As,[X|XsIn],[A|AsIn],XsOut,AsOut).
  379
  380var_existing([Xex|Xs],[Aex|As], X,A) :- Xex==X -> Aex=A ; var_existing(Xs,As,X,A).
 integrate(+F, -X:interval, -R) is semidet
Equivalent to integrate/4 with default precision. */
 integrate(+F, -X:interval, -R, +P:integer) is semidet
Succeeds if R is the numerical integration of F over the range of interval X. F must be continuously differentiable over this range or it fails.

The number of integration steps (= 2**P) is determined by the precision parameter P (default is value of environment flag clpBNR_default_precision). Example of use with increasing precision values:

?- X::real(0.0,1.0), F=X**2, between(2,10,P),integrate(F,X,RV,P), range(RV,R), format('~w:~w\n',[P,R]), fail.
2:[0.328125,0.359375]
3:[0.33203125,0.33984375]
4:[0.3330078125,0.3349609375]
5:[0.333251953125,0.333740234375]
6:[0.33331298828125,0.33343505859375]
7:[0.3333282470703125,0.3333587646484375]
8:[0.3333320617675781,0.3333396911621094]
9:[0.33333301544189453,0.33333492279052734]
10:[0.33333325386047363,0.33333373069763184]
false.

*/

  408% integrate(F,X,R) where X is an interval over which to integrate and F = f(X)
  409% Note that integrate(F,X,R) is equivalent to boundary_values(X,[dV(_,F,0,R])
  410integrate(F,X,R) :-
  411	current_prolog_flag(clpBNR_default_precision,P),
  412	integrate(F,X,R,P,_).
  413integrate(F,X,R,P) :-
  414	integrate(F,X,R,P,_).
  415integrate(F,X,R,P,Steps) :-                    % internal arity 5 for development
  416	compound(F),                               % F must be an expression in X
  417    interval(X),                               % X must be an interval 
  418    integer(P), P>0,                           % P must be positive integer
  419	!,                                         % args OK, commit 
  420	boundary_values(X,[dV(_,F,0,R)],P,Steps). % use integration in `boundary_values`
  421integrate(F,X,R,P,_Steps) :-
  422	fail_msg_('Invalid argument(s): ~w',[integrate(F,X,R,P)]).
 boundary_values(-X:interval, +Dvars:list) is semidet
Equivalent to boundary_values/4 with default precision and discarding steps list. */
 boundary_values(-X:interval, +Dvars:list, +P:integer) is semidet
Equivalent to boundary_values/4 discarding steps list. */
 boundary_values(-X:interval, +Dvars:list, +P:integer, -Steps:list) is semidet
Succeeds if a solution can be found to the boundary value problem specified by the independent variable X and the list of dependent variables defined by Dvars. Each dependent variable definition is a compound term of the form dV(Y, Fxy, Yi, Yf). The initial and final values of the independent variable for the purposes of the boundary value problem are specified by the lower and upper values of the domain of X. The arguments of for dvar/4 are:
  472boundary_values(X,YDefs) :-
  473	current_prolog_flag(clpBNR_default_precision,P),
  474	boundary_values_(X,YDefs,P,_).
  475	
  476boundary_values(X,YDefs,P) :-
  477	boundary_values_(X,YDefs,P,_).
  478
  479boundary_values(X,YDefs,P,Steps) :-
  480	boundary_values_(X,YDefs,P,Steps/_).
  481
  482boundary_values_(X,YDefs,P,Out/[(Xf,Yfs)]) :-
  483    integer(P), P>0,                             % P must be positive integer
  484    domain(X,Xdom), Xdom =.. [_Type,Xi,Xf],      % X must be an interval 
  485	eval_dvars(YDefs,Ys,Yis,Yfs,Ydoms),          % too many args for maplist	
  486	maplist(fXY_lambda(X,Ys),YDefs,Fxys),        % list of partial derivative lambda args
  487	(maplist(total_derivative_(Fxys),Fxys,DFxys) % list of connective derivative lambda args
  488	 -> true
  489	 ;  DFxys = none                             % F non-differentiable(?), use euler step
  490	),
  491	!,                                           % args all good, commit
  492    integrate_(P,Fxys,DFxys,(Xi,Yis),(Xf,Yfs),Ydoms,Out/[(Xf,Yfs)]).
  493boundary_values_(X,YDefs,P,_) :-
  494	fail_msg_('Invalid argument(s): ~w',[boundary_values(X,YDefs,P)]).
  495
  496eval_dvars([],[],[],[],[]).
  497eval_dvars([dV(Y, _PD, Lexp, Uexp)|YDefs],[Y|Ys],[Yi|Yis],[Yf|Yfs],[Ydom|Ydoms]) :-
  498	(domain(Y,Ydom) -> true ; Ydom = real),  % Ydom defaults to real 
  499	(var(Lexp) -> Yi=Lexp, Yi::Ydom ; Yi is Lexp),
  500	(var(Uexp) -> Yf=Uexp, Yf::Ydom ; Yf is Uexp),
  501	eval_dvars(YDefs,Ys,Yis,Yfs,Ydoms).
  502
  503% construct Lambda args for Fxy
  504fXY_lambda(X,Ys,dV(_Y,Fxy,_,_),FxyArgs) :- 
  505	lambda_for_(Fxy,X,Ys,FxyArgs).
  506
  507% construct Lambda args for derivative function of Fxy from Lambda of Fxy
  508
  509total_derivative_(Fxys,_Free/Ps,DxyArgs) :- !,  % ignore free variables
  510	total_derivative_(Fxys,Ps,DxyArgs).
  511total_derivative_(Fxys,[X,Ys,Fxy],DxyArgs) :-
  512	partial_derivative(Fxy,X,DFDX),          % clpBNR built-in
  513	sumYpartials(Fxys,Ys,Fxy,0,DYsum),
  514	simplify_sum_(DFDX, DYsum, DExp), !,
  515	lambda_for_(DExp,X,Ys,DxyArgs).
  516	
  517sumYpartials([],[],_Fxy,Acc,Acc).
  518sumYpartials([_Free/FxyI|FxyIs],YIs,Fxy,Acc,Sum) :- !,
  519	sumYpartials([FxyI|FxyIs],YIs,Fxy,Acc,Sum).
  520sumYpartials([[_X,_Ys,FxyI]|FxyIs],[YI|YIs],Fxy,Acc,Sum) :-
  521	partial_derivative(Fxy,YI,DFDYI),
  522	(number(DFDYI), DFDYI =:= 0 -> NxtAcc = Acc ; simplify_sum_(Acc,FxyI*DFDYI,NxtAcc)),
  523	!,
  524	sumYpartials(FxyIs,YIs,Fxy,NxtAcc,Sum).
  525
  526simplify_sum_(X,Y,Y) :- number(X),X=:=0.
  527simplify_sum_(X,Y,X) :- number(Y),Y=:=0.
  528simplify_sum_(X,Y,X+Y).
  529
  530% construct args for Lambda expression
  531lambda_for_(Fxy,X,Ys,Args) :-
  532	Lambda_parms = [X,Ys,Fxy],
  533	term_variables(Fxy,FVs),
  534	exclude(var_member_([X|Ys]),FVs,EVs),    % EVs = free variables
  535	(comma_op_(EVs,EV) -> Args = {EV}/Lambda_parms ; Args = Lambda_parms).
  536
  537var_member_([L|Ls],E) :- L==E -> true ; var_member_(Ls,E).
  538
  539comma_op_([X],X).  % assumes use in if-then
  540comma_op_([X|Xs],(X,Y)) :- comma_op_(Xs,Y).	
  541
  542% integration loop
  543integrate_(0, Fxys, Dxys, Initial, Final, Ydomains, [Initial|Ps]/Ps) :- !,
  544	% select integration step
  545	( Dxys == none
  546	 -> step_euler(Fxys, Initial, Final, Ydomains)
  547	 ;  step_trap(Fxys, Dxys, Initial, Final, Ydomains)
  548	).
  549integrate_(P, Fxys, Dxys, Initial, Final, Ydomains, L/E) :-
  550	% create interpolation point and integrate two halves
  551	interpolate_(Initial, Final, Ydomains, Middle),
  552	Pn is P - 1,
  553	integrate_(Pn, Fxys, Dxys, Initial, Middle, Ydomains, L/M),
  554	integrate_(Pn, Fxys, Dxys, Middle,  Final,  Ydomains, M/E).
  555
  556interpolate_((X0,_Y0s), (X1,_Y1s), Ydomains, (XM,YMs)) :-
  557	XM is (X0 + X1)/2,            % XM is midpoint of (X0,X1)
  558	maplist(::,YMs,Ydomains).     % corresponding YMs
  559
  560step_euler(Fxys, (X0,Y0), (X1,Y1), Ydoms) :-
  561	X01:: real(X0,X1),            % range of X in step
  562	maplist(lambda_constrain_(X01,Y01),Fxys,Fs),   % approx f' over X0 
  563	Dx is X1 - X0,                % assumed (X1>X0)
  564	DX :: real(0,Dx),             % range for estimate
  565	euler_constraints(Y0,Y1,Y01,Ydoms,Dx,DX,Fs,In/In,Cs/[]),  % flatten with diff list
  566	{Cs}.
  567	
  568step_trap(Fxys, Dxys, (X0,Y0), (X1,Y1), Ydoms) :-
  569	X01:: real(X0,X1),            % range of X in step
  570	maplist(lambda_constrain_(X0,Y0),Fxys,F0s),    % F0s = slopes at X0
  571	maplist(lambda_constrain_(X1,Y1),Fxys,F1s),    % F1s = slopes at X1 
  572	maplist(lambda_constrain_(X01,Y01),Dxys,Ds),   % approx f' over X0 
  573	Dx is X1 - X0,                % assumed (X1>X0)
  574	DX :: real(0,Dx),             % range for estimate
  575	trap_constraints(Y0,Y1,Y01,Ydoms,Dx,DX,F0s,F1s,Ds,In/In,Cs/[]),  % flatten with diff list
  576	{Cs}.  %%, absolve(Y1,2).
  577
  578lambda_constrain_(X,Ys,Args,F) :-  % reorder args for yall: >>
  579	yall: >>(Args,true,X,Ys,F).    % avoid meta-call (basically just makes copy)
  580%  known safe since lambda Goal=true
  581sandbox:safe_primitive(clpBNR_toolkit:lambda_constrain_(_X,_Ys,_Args,_F)). 
  582
  583/* see Carleton notes:
  584https://www.softwarepreservation.org/projects/prolog/bnr/doc/Older-Introduction_to_CLP%28BNR%29-1995.pdf/view
  585*/
  586euler_constraints([],[],[],[],_Dx,_DX,[],In,In).
  587euler_constraints([Y0|Y0s],[Y1|Y1s],[Y01|Y01s],[Ydom|Ydoms],Dx,DX,[F|Fs],
  588	In/[FM <= F,                      % FM = slope inclusion 
  589	    Y01 - Y0 is DX*FM,
  590	    Y1  - Y0 is Dx*FM
  591	   |Cs],
  592	Out) :-
  593	Y01:: Ydom,
  594	euler_constraints(Y0s,Y1s,Y01s,Ydoms,Dx,DX,Fs,In/Cs,Out).
  595
  596trap_constraints([],[],[],[],_Dx,_DX,[],[],[],In,In).
  597trap_constraints([Y0|Y0s],[Y1|Y1s],[Y01|Y01s],[Ydom|Ydoms],Dx,DX,[F0|F0s],[F1|F1s],[D|Ds],
  598	% use `is` to circumvent `simplify`
  599	In/[FM <= (F0+F1)/2,              % FM = average slope using step endpoints (one-way)	    % Note that the following must not be symbolically simplified to eliminate D	    8*DR is D - D,                % 4*deltaR == (D-D)/2 (for error term)	    Y01 - Y0 is DX*(FM + DR*DX),	    Y1  - Y0 is Dx*(FM + DR*Dx)	   |Cs],	Out) :-	Y01:: Ydom, DR::real,	trap_constraints(Y0s,Y1s,Y01s,Ydoms,Dx,DX,F0s,F1s,Ds,In/Cs,Out).
 lin_minimum(+ObjF, Constraints:{}, ?Min:numeric) is semidet
Succeeds if the global minimum value of the objective function ObjF subject to Constraints can be unified with Min; otherwise fails. Both the objective function and Constraints must be linear, i.e., only subexpressions summing of the form X*C (or C*X) are permitted since the actual computation is done using library(simplex). Narrowing of minimizers (variables in ObjF) is limited to that constrained by the Min result to accomodate multiple sets of minimizers. (See lin_minimize/3 to use minimizers used to derive Min.) A solution generator, e.g., clpBNR:solve/1 can be used to search for alternative sets of minimizers. "Universal Mines" example from the User Guide:
?- [M_Idays,M_IIdays,M_IIIdays]::integer(0,7),
lin_minimum(20*M_Idays+22*M_IIdays+18*M_IIIdays,
{4*M_Idays+6*M_IIdays+M_IIIdays>=54,4*M_Idays+4*M_IIdays+6*M_IIIdays>=65}, Min).
Min = 284,
M_Idays::integer(2, 7),
M_IIdays::integer(4, 7),
M_IIIdays::integer(2, 7).

?- [M_Idays,M_IIdays,M_IIIdays]::integer(0,7),
lin_minimum(20*M_Idays+22*M_IIdays+18*M_IIIdays,
{4*M_Idays+6*M_IIdays+M_IIIdays>=54,4*M_Idays+4*M_IIdays+6*M_IIIdays>=65}, Min),
solve([M_Idays,M_IIdays,M_IIIdays]).
M_Idays = 2,
M_IIdays = 7,
M_IIIdays = 5,
Min = 284 ;
false.

For linear systems, lin_minimum/3, lin_maximum/3 can be significantly faster than using the more general purpose clpBNR:global_minimum/3, clpBNR:global_maximum/3

See also
- lin_minimize/3, library(simplex) */
 lin_maximum(+ObjF, Constraints:{}, ?Max:numeric) is semidet
Succeeds if the global minimum value of the objective function ObjF subject to Constraints can be unified with Max; otherwise fails. This is the corresponding predicate to lin_minimum/3 for finding global maxima.
See also
- lin_minimum/3, lin_maximize/3 */
  644lin_minimum(ObjF,{Constraints},MinValue) :-
  645	lin_minimum_(ObjF,{Constraints},MinValue,false).
  646
  647lin_maximum(ObjF,{Constraints},MinValue) :-
  648	lin_maximum_(ObjF,{Constraints},MinValue,false).
 lin_minimize(+ObjF, Constraints:{}, ?Min:numeric) is semidet
Succeeds if the global minimum value of the objective function ObjF subject to Constraints can be unified with Min; otherwise fails. This behaves identically to lin_minimum/3 except variables in ObjF will be narrowed to the values used in calculating the final value of Min. Any other sets of minimizers corresponding to Min are removed from the solution space. "Universal Mines" example from the User Guide:
?- [M_Idays,M_IIdays,M_IIIdays]::integer(0,7),
lin_minimize(20*M_Idays+22*M_IIdays+18*M_IIIdays,
{4*M_Idays+6*M_IIdays+M_IIIdays>=54,4*M_Idays+4*M_IIdays+6*M_IIIdays>=65}, Min).
M_Idays = 2,
M_IIdays = 7,
M_IIIdays = 5,
Min = 284.
See also
- lin_minimum/3 */
 lin_maximize(+ObjF, Constraints:{}, ?Max:numeric) is semidet
Succeeds if the global maximum value of the objective function ObjF subject to Constraints can be unified with Max; otherwise fails. This behaves identically to lin_maximum/3 except variables in ObjF will be narrowed to the values used in calculating the final value of Max. Any other sets of minimizers corresponding to Min are removed from the solution space.
See also
- lin_maximum/3 */
  671lin_minimize(ObjF,{Constraints},MinValue) :-
  672	lin_minimum_(ObjF,{Constraints},MinValue,true).
  673
  674lin_maximize(ObjF,{Constraints},MinValue) :-
  675	lin_maximum_(ObjF,{Constraints},MinValue,true).
  676	
  677
  678lin_minimum_(ObjF,{Constraints},MinValue,BindVars) :-
  679	init_simplex_(ObjF,Constraints,Objective,S0,Vs),
  680	(minimize(Objective,S0,S)
  681	 -> objective(S,MinValue), {ObjF == MinValue},
  682	    (BindVars == true
  683	     -> bind_vars_(Vs,S)
  684	     ;  remove_names_(Vs),
  685	        {Constraints}  % apply constraints
  686	    )
  687	 ;  fail_msg_('Failed to minimize: ~w',[ObjF])
  688	).
  689	
  690lin_maximum_(ObjF,{Constraints},MaxValue,BindVars) :-
  691	init_simplex_(ObjF,Constraints,Objective,S0,Vs),
  692	(maximize(Objective,S0,S)
  693	 -> objective(S,MaxValue), {ObjF == MaxValue},
  694	    (BindVars == true
  695	     -> bind_vars_(Vs,S)
  696	     ;  remove_names_(Vs),
  697	        {Constraints}  % apply constraints
  698	    )
  699	 ;  fail_msg_('Failed to maximize: ~w',[ObjF])
  700	).
  701
  702init_simplex_(ObjF,Constraints,Objective,S,Vs) :-
  703	gen_state(S0),
  704	term_variables((ObjF,Constraints),Vs),
  705	(Vs = []
  706	 -> fail_msg_('No variables in expression to optimize: ~w',[ObjF])
  707	 ;  sim_constraints_(Constraints,S0,S1),
  708	    _::real(_,Max),  % max value to constrain for simplex
  709	    constrain_ints_(Vs,Max,S1,S),
  710	    (map_simplex_(ObjF,T/T,Objective/[])
  711	     -> true
  712	     ;  fail_msg_('Illegal linear objective: ~w',[ObjF])
  713	    )
  714	).
  715
  716% use an attribute to associate a var with a unique simplex variable name
  717simplex_var_(V,SV) :- var(V),
  718	(get_attr(V,clpBNR_toolkit,SV)
  719	 -> true
  720	 ;  statistics(inferences,VID), SV = var(VID), put_attr(V,clpBNR_toolkit,SV)
  721	).
  722
  723% Name attribute removed before exit.
  724remove_names_([]).
  725remove_names_([V|Vs]) :-
  726	del_attr(V,clpBNR_toolkit),
  727	remove_names_(Vs).
  728		 
  729attr_unify_hook(var(_),_).     % unification always does nothing and succeeds
  730
  731constrain_ints_([],_,S,S).
  732constrain_ints_([V|Vs],Max,Sin,Sout) :- 
  733	% Note: library(simplex) currently constrains all variables to be non-negative
  734	simplex_var_(V,SV),               % corresponding simplex variable name
  735	% if not already an interval, make it one with finite non-negative value
  736	(domain(V,D) -> true ; V::real(0,_), domain(V,D)),
  737	(D == boolean -> Dom = integer(0,1); Dom = D),
  738	Dom =.. [Type,L,H],
  739	(Type == integer -> constraint(integral(SV),Sin,S1) ; S1 = Sin),
  740	(L < 0
  741	 -> % apply non-negativity condition    
  742	    ({V >= 0} -> L1 = 0 ; fail_msg_('Negative vars not supported by \'simplex\': ~w',[V]))
  743	 ;  L1 = L
  744	),
  745	% explicitly constrain any vars not (0,Max-eps)
  746	(L1 > 0   -> constraint([SV] >= L1,S1,S2)   ; S2 = S1),    % L1 not negative
  747	(H  < Max -> constraint([SV] =< H,S2,SNxt)  ; SNxt = S2),
  748	constrain_ints_(Vs,Max,SNxt,Sout).
  749
  750bind_vars_([],_).
  751bind_vars_([V|Vs],S) :-  
  752	% Note: skip anything nonvar (already bound due to active constraints)
  753	(simplex_var_(V,SV) -> variable_value(S,SV,V) ; true),
  754	bind_vars_(Vs,S).
  755
  756% clpBNR constraints have already been applied so worst errors have been detected
  757sim_constraints_([],S,S) :- !.
  758sim_constraints_([C|Cs],Sin,Sout) :- !,
  759	sim_constraints_(C, Sin,Snxt),
  760	sim_constraints_(Cs,Snxt,Sout).
  761sim_constraints_((C,Cs),Sin,Sout) :-  !,
  762	sim_constraints_(C, Sin,Snxt),
  763	sim_constraints_(Cs,Snxt,Sout).
  764sim_constraints_(C,Sin,Sout) :- 
  765	sim_constraint_(C,SC),
  766	constraint(SC,Sin,Sout).  % from simplex
  767
  768sim_constraint_(C,SC) :-
  769	C=..[Op,LHS,RHS],              % decompose
  770	constraint_op(Op,COp),         % acceptable to simplex
  771	number(RHS), RHS >= 0,         % requirement of simplex
  772	map_simplex_(LHS,T/T,Sim/[]),  % map to simplex list of terms
  773	!,
  774	SC=..[COp,Sim,RHS].            % recompose
  775sim_constraint_(C,_) :-
  776	fail_msg_('Illegal linear constraint: ~w',[C]).	
  777
  778map_simplex_(T,CT/[S|Tail],CT/Tail) :-
  779	map_simplex_term_(T,S),
  780	!.
  781map_simplex_(A+B,Cin,Cout) :- !,
  782	map_simplex_(A,Cin,Cnxt),
  783	map_simplex_(B,Cnxt,Cout).
  784map_simplex_(A-B,Cin,Cout) :- !,
  785	map_simplex_(A,Cin,Cnxt),
  786	map_simplex_(-B,Cnxt,Cout).
  787			
  788map_simplex_term_(V,1*SV)   :- simplex_var_(V,SV), !.
  789map_simplex_term_(-T,NN*V)  :- !,
  790	map_simplex_term_(T,N*V),
  791	NN is -N.
  792map_simplex_term_(N*V,N*SV) :- number(N), simplex_var_(V,SV), !.
  793map_simplex_term_(V*N,N*SV) :- number(N), simplex_var_(V,SV).
  794
  795constraint_op(==,=).
  796constraint_op(=<,=<).
  797constraint_op(>=,>=).
 local_minima(+ObjF) is semidet
Succeeds if the value of objective function ObjF can be constrained to be a local minimum, i.e, it's "slope" is 0 in every dimension; otherwise fails. This requires that a partial derivative of ObjF exists for each variable. local_minima should be executed prior to a call to clpBNR:global_minimum using the same objective function, e.g.,
?- X::real(0,10), OF=X**3-6*X**2+9*X+6, local_minima(OF), global_minimum(OF,Z).
OF = X**3-6*X**2+9*X+6,
X:: 3.00000000000000...,
Z:: 6.000000000000... .

Using any local optima predicate can significantly improve performance compared to searching for global optima (clpBNR:global_*) without local constraints.

See also
- clpBNR:local_minima/2 */
 local_maxima(+ObjF) is semidet
Succeeds if the value of objective function ObjF can be constrained to be a local maximum, i.e, it's "slope" is 0 in every dimension; otherwise fails. This requires that a partial derivative of ObjF exists for each variable. local_maxima should be executed prior to a call to clpBNR:global_maximum using the same objective function, e.g.,
?- X::real(0,10), OF=X**3-6*X**2+9*X+6, local_maxima(OF), global_maximum(OF,Z).
OF = X**3-6*X**2+9*X+6,
X:: 1.000000000000000...,
Z:: 10.0000000000000... .
See also
- clpBNR:local_maxima/2 */
  825%
  826%	local_minima/1,       % apply KT constraints for objective function expression (OFE)
  827%	local_maxima/1,       % semantically equivalent to local_minima/1
  828%
  829local_minima(ObjExp) :-
  830	local_optima_(min,ObjExp,[]).
  831
  832local_maxima(ObjExp) :-
  833	local_optima_(max,ObjExp,[]).
 local_minima(+ObjF, +Constraints:{}) is semidet
Succeeds if the value of objective function ObjF can be constrained to be a local minimum, i.e, it's "slope" is 0 in every dimension, subject to Constraints; otherwise fails. This requires that a partial derivative of ObjF, and any subexpression in Constraints, exists for each variable. local_minima should be executed prior to a call to clpBNR:global_minimum using the same objective function, e.g.,
?- [X1,X2]::real, OF=X1**4*exp(-0.01*(X1*X2)**2),
local_minima(OF,{2*X1**2+X2**2==10}), global_minimum(OF,Z), solve([X1,X2]).
OF = X1**4*exp(-0.01*(X1*X2)**2),
X1::real(-1.703183936003284e-108, 1.703183936003284e-108),
X2:: -3.16227766016838...,
Z:: 0.0000000000000000... ;
OF = X1**4*exp(-0.01*(X1*X2)**2),
X1::real(-1.703183936003284e-108, 1.703183936003284e-108),
X2:: 3.16227766016838...,
Z:: 0.0000000000000000... .
See also
- clpBNR:local_minima/1 */
 local_maxima(+ObjF, +Constraints:{}) is semidet
Succeeds if the value of objective function ObjF can be constrained to be a local maximum, i.e, it's "slope" is 0 in every dimension; otherwise fails. This requires that a partial derivative of ObjF, and any subexpression in Constraints, exists for each variable. local_maxima should be executed prior to a call to clpBNR:global_maximum using the same objective function, e.g.,
?- [X1,X2]::real,OF=X1**4*exp(-0.01*(X1*X2)**2),
local_maxima(OF,{2*X1**2+X2**2==10}), global_maximum(OF,Z),solve([X1,X2]).
OF = X1**4*exp(-0.01*(X1*X2)**2),
X1:: -2.23606797749979...,
X2:: 0.0000000000000000...,
Z:: 25.0000000000000... ;
OF = X1**4*exp(-0.01*(X1*X2)**2),
X1:: 2.23606797749979...,
X2:: 0.0000000000000000...,
Z:: 25.0000000000000... .
See also
- clpBNR:local_maxima/1 */
  872%
  873%	local_minima/2,       % apply KT constraints for minima with constraints
  874%	local_maxima/2        % apply KT constraints for maxima with constraints
  875%
  876local_minima(ObjExp,{Constraints}) :-
  877	local_optima_(min,ObjExp,Constraints).
  878	
  879local_maxima(ObjExp,{Constraints}) :-
  880	local_optima_(max,ObjExp,Constraints).
  881	
  882	
  883local_optima_(MinMax,ObjF,Constraints) :-
  884	local_optima_(MinMax,ObjF,Constraints,Cs),  % generate constraints
  885	{Cs}.                                       % then apply
  886
  887local_optima_(MinMax,ObjF,Constraints,[Constraints,GCs,LCs]) :-
  888	lagrangian_(Constraints,MinMax,ObjF,LObjF,LCs),
  889	term_variables((Constraints,ObjF),Vs),
  890	gradient_constraints_(Vs,GCs,LObjF).
  891
  892gradient_constraints_([],[],_Exp).
  893gradient_constraints_([X|Xs],[C|Cs],Exp) :-
  894	partial_derivative(Exp,X,D),
  895	(number(D) -> C=[] ; C=(D==0)),  % no constraint if PD is a constant
  896	gradient_constraints_(Xs,Cs,Exp).
  897
  898% for each constraint add to Lagrangian expression with auxiliary KKT constraints
  899lagrangian_(C,MinMax,Exp,LExp, LC) :- nonvar(C),
  900	kt_constraint_(C,CExp, LC), % generate langrange term with multiplier
  901	lexp(MinMax,Exp,CExp,LExp),
  902	!.
  903lagrangian_([],_,L,L,[]).
  904lagrangian_([C|Cs],MinMax,Exp,LExp,[LC|LCs]) :-
  905	lagrangian_(C, MinMax,Exp,NExp,LC),
  906	lagrangian_(Cs,MinMax,NExp,LExp,LCs).
  907lagrangian_((C,Cs),MinMax,Exp,LExp,[LC|LCs]) :-
  908	lagrangian_(C,MinMax,Exp,NExp,LC),
  909	lagrangian_(Cs,MinMax,NExp,LExp,LCs).
  910
  911lexp(min,Exp,CExp,Exp+CExp).
  912lexp(max,Exp,CExp,Exp-CExp).
  913
  914kt_constraint_(LHS == RHS, M*(LHS-RHS), []) :-
  915	M::real.                          % finite multiplier only
  916kt_constraint_(LHS =< RHS, MGx,     MGx==0) :-
  917	MGx = M*(LHS-RHS), M::real(0,_).  % positive multiplier and KKT non-negativity condition
  918kt_constraint_(LHS >= RHS, Exp,         LC) :-
  919	kt_constraint_(RHS =< LHS, Exp, LC).  % >= convert to =<