1/* The MIT License (MIT) 2 * 3 * Copyright (c) 2022-2025 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 ]).
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
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 , 97 N1 is N-1, 98 ( 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 for details of interval splitting for this predicate.
113mid_split_one(Xs) :- 114 select_split(Xs,X), % select largest interval with largest width 115 mid_split(X). % split it
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.
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 ).
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.
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/2 using default precision.
*/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.
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
== 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.
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) :- (list(Es) ; Es = (_,_)), !, % list or sequence 306 taylor_merged_contractor({Es},CF). 307taylor_contractor(Eq,_) :- 308 fail_msg_('Invalid constraint for Taylor contractor: ~w',[Eq]). 309 310make_As_and_Zs_([],_,[],[],Tail,Tail). 311make_As_and_Zs_([X|Xs],T,[A|As],[Z|Zs],[Z==A+T*(X-A)|CZs],Tail) :- 312 make_As_and_Zs_(Xs,T,As,Zs,CZs,Tail). 313 314sum_diffs([], [], [], _AllZs, _Exp, ExpIn, EQ) :- make_EQ_(ExpIn,EQ). 315sum_diffs([X|Xs], [A|As], [Z|Zs], AllZs, Exp, AExp, DEQ) :- 316 copy_term_nat(Exp,NExp), % copy expression and replace Xs by Zs 317 term_variables(NExp,AllZs), 318 partial_derivative(NExp,Z,DZ), % differentiate wrt. Z and add to generated expression 319 sum_diffs(Xs, As, Zs, AllZs, Exp, AExp+DZ*(X-A), DEQ). 320 321% map expression Exp to an equation equivalent to Exp==0 with numeric RHS 322make_EQ_(Exp,LHS==RHS) :- % turn expression into equation equivalent to Exp==0. 323 make_EQ_(Exp,LHS,RHS). 324 325make_EQ_(E,E,0) :- var(E), !. 326make_EQ_(X+Y,X,SY) :- number(Y), !, SY is -Y. 327make_EQ_(X-Y,X,Y) :- number(Y), !. 328make_EQ_(X+Y,Y,SX) :- number(X), !, SX is -X. 329make_EQ_(X-Y,SY,SX) :- number(X), !, SX is -X, negate_sum_(Y,SY). 330make_EQ_(X+Y,LHS+Y,RHS) :- !, make_EQ_(X,LHS,RHS). 331make_EQ_(X-Y,LHS-Y,RHS) :- !, make_EQ_(X,LHS,RHS). 332make_EQ_(E,E,0). % default (non +/- subexpression) 333 334negate_sum_(Y,-Y) :- var(Y), !. 335negate_sum_(X+Y,NX-Y) :- !, negate_sum_(X,NX). 336negate_sum_(X-Y,NX+Y) :- !, negate_sum_(X,NX). 337negate_sum_(E,-E).
== or =:=) constraint in Constraints; otherwise fails.
346% 347% build a cf_contractor by merging a list of cf_contractor's 348% 349taylor_merged_contractor({Es},T) :- 350 cf_list(Es,Ts), 351 cf_merge(Ts,T). 352 353cf_list([],[]) :- !. 354cf_list([C|Cs],[CF|CFs]) :- !, 355 cf_list(C, CF), 356 cf_list(Cs,CFs). 357cf_list((C,Cs),[CF|CFs]) :- !, 358 cf_list(C, CF), 359 cf_list(Cs,CFs). 360cf_list(C,CF) :- 361 taylor_contractor({C},CF). 362 363cf_merge(CFs,CF) :- cf_merge(CFs,cf_contractor([],[]),CF). 364 365cf_merge([],CF,CF). 366cf_merge([CF|CFs],CFIn,CFOut) :- 367 cf_merge(CF,CFIn,CFNxt), 368 cf_merge(CFs,CFNxt,CFOut). 369cf_merge(cf_contractor(Xs,As),cf_contractor(XsIn,AsIn),cf_contractor(XsOut,AsOut)) :- 370 cf_add(Xs,As,XsIn,AsIn,XsOut,AsOut). 371 372cf_add([],[],Xs,As,Xs,As). 373cf_add([X|Xs],[A|As],XsIn,AsIn,XsOut,AsOut) :- 374 var_existing(XsIn,AsIn,X,A), !, 375 cf_add(Xs,As,XsIn,AsIn,XsOut,AsOut). 376cf_add([X|Xs],[A|As],XsIn,AsIn,XsOut,AsOut) :- 377 cf_add(Xs,As,[X|XsIn],[A|AsIn],XsOut,AsOut). 378 379var_existing([Xex|Xs],[Aex|As], X,A) :- Xex==X -> Aex=A ; var_existing(Xs,As,X,A).
integrate/4 with default precision.
*/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.
*/
407% integrate(F,X,R) where X is an interval over which to integrate and F = f(X) 408% Note that integrate(F,X,R) is equivalent to boundary_values(X,[dV(_,F,0,R]) 409integrate(F,X,R) :- 410 current_prolog_flag(clpBNR_default_precision,P), 411 integrate(F,X,R,P,_). 412integrate(F,X,R,P) :- 413 integrate(F,X,R,P,_). 414integrate(F,X,R,P,Steps) :- % internal arity 5 for development 415 interval(X), % X must be an interval 416 integer(P), P>0, % P must be positive integer 417 !, % args OK, commit 418 boundary_values(X,[dV(_,F,0,R)],P,Steps). % use integration in `boundary_values` 419integrate(F,X,R,P,_Steps) :- 420 fail_msg_('Invalid argument(s): ~w',[integrate(F,X,R,P)]).
boundary_values/4 with default precision and discarding steps list.
*/boundary_values/4 discarding steps list.
*/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:
The optional third argument P defines a precision value, a positive integer (default = environment flag clpBNR_default_precision), which controls the the numerical integration; larger P means smaller step size.
The arity 4 version has an additional (final) argument which is unified with a list of the step values generated by the integration; each value is a tuple of the form (X,Ys).
This predicate fails if any of the arguments are invalid (generates an error message if clpBNR debug topic is enabled) or if a solution to the boundary value problem cannot be found. Examples for X in the range 0..1 and derivative of Y = -2*X*Y:
?- X::real(0,1), boundary_values(X,[dV(Y, -2*X*Y,1,Yf)]). X::real(0, 1), Yf:: 0.368... . ?- X::real(0,1),boundary_values(X,[dV(Y, -2*X*Y,1,Yf)],9). X::real(0, 1), Yf:: 0.36788... . ?- X::real(0,1), boundary_values(X,[dV(Y, -2*X*Y,Yi,1/e)]). X::real(0, 1), Yi:: 1.00... . ?- debug(clpBNR). true. ?- X=42, boundary_values(X,[dV(Y, -2*X*Y,Yi,1/e)]). % Invalid argument(s): boundary_values(42,[dV(_10836,-2*42*_10836,_10850,1/e)],6) false.
As with any application requiring numerical integration, care must be taken to avoid instability problems (more discussion in A Guide to CLP(BNR). */
470boundary_values(X,YDefs) :- 471 current_prolog_flag(clpBNR_default_precision,P), 472 boundary_values_(X,YDefs,P,_). 473 474boundary_values(X,YDefs,P) :- 475 boundary_values_(X,YDefs,P,_). 476 477boundary_values(X,YDefs,P,Steps) :- 478 boundary_values_(X,YDefs,P,Steps/_). 479 480boundary_values_(X,YDefs,P,Out/[(Xf,Yfs)]) :- 481 integer(P), P>0, % P must be positive integer 482 domain(X,Xdom), Xdom =.. [_Type,Xi,Xf], % X must be an interval 483 eval_dvars(YDefs,Ys,Yis,Yfs,Ydoms), % too many args for maplist 484 maplist(fXY_lambda(X,Ys),YDefs,Fxys), % list of partial derivative lambda args 485 (maplist(total_derivative_(Fxys),Fxys,DFxys) % list of connective derivative lambda args 486 -> true 487 ; DFxys = none % F non-differentiable(?), use euler step 488 ), 489 !, % args all good, commit 490 integrate_(P,Fxys,DFxys,(Xi,Yis),(Xf,Yfs),Ydoms,Out/[(Xf,Yfs)]). 491boundary_values_(X,YDefs,P,_) :- 492 fail_msg_('Invalid argument(s): ~w',[boundary_values(X,YDefs,P)]). 493 494eval_dvars([],[],[],[],[]). 495eval_dvars([dV(Y, _PD, Lexp, Uexp)|YDefs],[Y|Ys],[Yi|Yis],[Yf|Yfs],[Ydom|Ydoms]) :- 496 (domain(Y,Ydom) -> true ; Ydom = real), % Ydom defaults to real 497 (var(Lexp) -> Yi=Lexp, Yi::Ydom ; Yi is Lexp), 498 (var(Uexp) -> Yf=Uexp, Yf::Ydom ; Yf is Uexp), 499 eval_dvars(YDefs,Ys,Yis,Yfs,Ydoms). 500 501% construct Lambda args for Fxy 502fXY_lambda(X,Ys,dV(_Y,Fxy,_,_),FxyArgs) :- 503 lambda_for_(Fxy,X,Ys,FxyArgs). 504 505% construct Lambda args for derivative function of Fxy from Lambda of Fxy 506 507total_derivative_(Fxys,_Free/Ps,DxyArgs) :- !, % ignore free variables 508 total_derivative_(Fxys,Ps,DxyArgs). 509total_derivative_(Fxys,[X,Ys,Fxy],DxyArgs) :- 510 partial_derivative(Fxy,X,DFDX), % clpBNR built-in 511 sumYpartials(Fxys,Ys,Fxy,0,DYsum), 512 simplify_sum_(DFDX, DYsum, DExp), !, 513 lambda_for_(DExp,X,Ys,DxyArgs). 514 515sumYpartials([],[],_Fxy,Acc,Acc). 516sumYpartials([_Free/FxyI|FxyIs],YIs,Fxy,Acc,Sum) :- !, 517 sumYpartials([FxyI|FxyIs],YIs,Fxy,Acc,Sum). 518sumYpartials([[_X,_Ys,FxyI]|FxyIs],[YI|YIs],Fxy,Acc,Sum) :- 519 partial_derivative(Fxy,YI,DFDYI), 520 (number(DFDYI), DFDYI =:= 0 -> NxtAcc = Acc ; simplify_sum_(Acc,FxyI*DFDYI,NxtAcc)), 521 !, 522 sumYpartials(FxyIs,YIs,Fxy,NxtAcc,Sum). 523 524simplify_sum_(X,Y,Y) :- number(X),X=:=0. 525simplify_sum_(X,Y,X) :- number(Y),Y=:=0. 526simplify_sum_(X,Y,X+Y). 527 528% construct args for Lambda expression 529lambda_for_(Fxy,X,Ys,Args) :- 530 Lambda_parms = [X,Ys,Fxy], 531 term_variables(Fxy,FVs), 532 exclude(var_member_([X|Ys]),FVs,EVs), % EVs = free variables 533 (comma_op_(EVs,EV) -> Args = {EV}/Lambda_parms ; Args = Lambda_parms). 534 535var_member_([L|Ls],E) :- L==E -> true ; var_member_(Ls,E). 536 537comma_op_([X],X). % assumes use in if-then 538comma_op_([X|Xs],(X,Y)) :- comma_op_(Xs,Y). 539 540% integration loop 541integrate_(0, Fxys, Dxys, Initial, Final, Ydomains, [Initial|Ps]/Ps) :- !, 542 % select integration step 543 ( Dxys == none 544 -> step_euler(Fxys, Initial, Final, Ydomains) 545 ; step_trap(Fxys, Dxys, Initial, Final, Ydomains) 546 ). 547integrate_(P, Fxys, Dxys, Initial, Final, Ydomains, L/E) :- 548 % create interpolation point and integrate two halves 549 interpolate_(Initial, Final, Ydomains, Middle), 550 Pn is P - 1, 551 integrate_(Pn, Fxys, Dxys, Initial, Middle, Ydomains, L/M), 552 integrate_(Pn, Fxys, Dxys, Middle, Final, Ydomains, M/E). 553 554interpolate_((X0,_Y0s), (X1,_Y1s), Ydomains, (XM,YMs)) :- 555 XM is (X0 + X1)/2, % XM is midpoint of (X0,X1) 556 maplist(::,YMs,Ydomains). % corresponding YMs 557 558step_euler(Fxys, (X0,Y0), (X1,Y1), Ydoms) :- 559 X01:: real(X0,X1), % range of X in step 560 maplist(lambda_constrain_(X01,Y01),Fxys,Fs), % approx f' over X0 561 Dx is X1 - X0, % assumed (X1>X0) 562 DX :: real(0,Dx), % range for estimate 563 euler_constraints(Y0,Y1,Y01,Ydoms,Dx,DX,Fs,In/In,Cs/[]), % flatten with diff list 564 {Cs}. 565 566step_trap(Fxys, Dxys, (X0,Y0), (X1,Y1), Ydoms) :- 567 X01:: real(X0,X1), % range of X in step 568 maplist(lambda_constrain_(X0,Y0),Fxys,F0s), % F0s = slopes at X0 569 maplist(lambda_constrain_(X1,Y1),Fxys,F1s), % F1s = slopes at X1 570 maplist(lambda_constrain_(X01,Y01),Dxys,Ds), % approx f' over X0 571 Dx is X1 - X0, % assumed (X1>X0) 572 DX :: real(0,Dx), % range for estimate 573 trap_constraints(Y0,Y1,Y01,Ydoms,Dx,DX,F0s,F1s,Ds,In/In,Cs/[]), % flatten with diff list 574 {Cs}. %%, absolve(Y1,2). 575 576lambda_constrain_(X,Ys,Args,F) :- % reorder args for yall: >> 577 yall: >>(Args,true,X,Ys,F). % avoid meta-call (basically just makes copy) 578% known safe since lambda Goal=true 579sandbox:safe_primitive(clpBNR_toolkit:lambda_constrain_(_X,_Ys,_Args,_F)). 580 581/* see Carleton notes: 582https://www.softwarepreservation.org/projects/prolog/bnr/doc/Older-Introduction_to_CLP%28BNR%29-1995.pdf/view 583*/ 584euler_constraints([],[],[],[],_Dx,_DX,[],In,In). 585euler_constraints([Y0|Y0s],[Y1|Y1s],[Y01|Y01s],[Ydom|Ydoms],Dx,DX,[F|Fs], 586 In/[FM <= F, % FM = slope inclusion 587 Y01 - Y0 is DX*FM, 588 Y1 - Y0 is Dx*FM 589 |Cs], 590 Out) :- 591 Y01:: Ydom, 592 euler_constraints(Y0s,Y1s,Y01s,Ydoms,Dx,DX,Fs,In/Cs,Out). 593 594trap_constraints([],[],[],[],_Dx,_DX,[],[],[],In,In). 595trap_constraints([Y0|Y0s],[Y1|Y1s],[Y01|Y01s],[Ydom|Ydoms],Dx,DX,[F0|F0s],[F1|F1s],[D|Ds], 596 % use `is` to circumvent `simplify` 597 In/[FM <= (F0+F1)/2, % FM = average slope using step endpoints (one-way) 598 % Note that the following must not be symbolically simplified to eliminate D 599 8*DR is D - D, % 4*deltaR == (D-D)/2 (for error term) 600 Y01 - Y0 is DX*(FM + DR*DX), 601 Y1 - Y0 is Dx*(FM + DR*Dx) 602 |Cs], 603 Out) :- 604 Y01:: Ydom, DR::real, trap_constraints(Y0s,Y1s,Y01s,Ydoms,Dx,DX,F0s,F1s,Ds,In/Cs,Out).
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
lin_minimum/3 for finding global maxima.
642lin_minimum(ObjF,{Constraints},MinValue) :- 643 lin_minimum_(ObjF,{Constraints},MinValue,false). 644 645lin_maximum(ObjF,{Constraints},MinValue) :- 646 lin_maximum_(ObjF,{Constraints},MinValue,false).
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.
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.
669lin_minimize(ObjF,{Constraints},MinValue) :- 670 lin_minimum_(ObjF,{Constraints},MinValue,true). 671 672lin_maximize(ObjF,{Constraints},MinValue) :- 673 lin_maximum_(ObjF,{Constraints},MinValue,true). 674 675 676lin_minimum_(ObjF,{Constraints},MinValue,BindVars) :- 677 init_simplex_(ObjF,Constraints,Objective,S0,Vs), 678 (minimize(Objective,S0,S) 679 -> objective(S,MinValue), {ObjF == MinValue}, 680 (BindVars == true 681 -> bind_vars_(Vs,S) 682 ; remove_names_(Vs), 683 {Constraints} % apply constraints 684 ) 685 ; fail_msg_('Failed to minimize: ~w',[ObjF]) 686 ). 687 688lin_maximum_(ObjF,{Constraints},MaxValue,BindVars) :- 689 init_simplex_(ObjF,Constraints,Objective,S0,Vs), 690 (maximize(Objective,S0,S) 691 -> objective(S,MaxValue), {ObjF == MaxValue}, 692 (BindVars == true 693 -> bind_vars_(Vs,S) 694 ; remove_names_(Vs), 695 {Constraints} % apply constraints 696 ) 697 ; fail_msg_('Failed to maximize: ~w',[ObjF]) 698 ). 699 700init_simplex_(ObjF,Constraints,Objective,S,Vs) :- 701 gen_state(S0), 702 term_variables((ObjF,Constraints),Vs), 703 (Vs = [] 704 -> fail_msg_('No variables in expression to optimize: ~w',[ObjF]) 705 ; sim_constraints_(Constraints,S0,S1), 706 _::real(_,Max), % max value to constrain for simplex 707 constrain_ints_(Vs,Max,S1,S), 708 (map_simplex_(ObjF,T/T,Objective/[]) 709 -> true 710 ; fail_msg_('Illegal linear objective: ~w',[ObjF]) 711 ) 712 ). 713 714% use an attribute to associate a var with a unique simplex variable name 715simplex_var_(V,SV) :- var(V), 716 (get_attr(V,clpBNR_toolkit,SV) 717 -> true 718 ; statistics(inferences,VID), SV = var(VID), put_attr(V,clpBNR_toolkit,SV) 719 ). 720 721% Name attribute removed before exit. 722remove_names_([]). 723remove_names_([V|Vs]) :- 724 del_attr(V,clpBNR_toolkit), 725 remove_names_(Vs). 726 727attr_unify_hook(var(_),_). % unification always does nothing and succeeds 728 729constrain_ints_([],_,S,S). 730constrain_ints_([V|Vs],Max,Sin,Sout) :- 731 % Note: library(simplex) currently constrains all variables to be non-negative 732 simplex_var_(V,SV), % corresponding simplex variable name 733 % if not already an interval, make it one with finite non-negative value 734 (domain(V,D) -> true ; V::real(0,_), domain(V,D)), 735 (D == boolean -> Dom = integer(0,1); Dom = D), 736 Dom =.. [Type,L,H], 737 (Type == integer -> constraint(integral(SV),Sin,S1) ; S1 = Sin), 738 (L < 0 739 -> % apply non-negativity condition 740 ({V >= 0} -> L1 = 0 ; fail_msg_('Negative vars not supported by \'simplex\': ~w',[V])) 741 ; L1 = L 742 ), 743 % explicitly constrain any vars not (0,Max-eps) 744 (L1 > 0 -> constraint([SV] >= L1,S1,S2) ; S2 = S1), % L1 not negative 745 (H < Max -> constraint([SV] =< H,S2,SNxt) ; SNxt = S2), 746 constrain_ints_(Vs,Max,SNxt,Sout). 747 748bind_vars_([],_). 749bind_vars_([V|Vs],S) :- 750 % Note: skip anything nonvar (already bound due to active constraints) 751 (simplex_var_(V,SV) -> variable_value(S,SV,V) ; true), 752 bind_vars_(Vs,S). 753 754% clpBNR constraints have already been applied so worst errors have been detected 755sim_constraints_([],S,S) :- !. 756sim_constraints_([C|Cs],Sin,Sout) :- !, 757 sim_constraints_(C, Sin,Snxt), 758 sim_constraints_(Cs,Snxt,Sout). 759sim_constraints_((C,Cs),Sin,Sout) :- !, 760 sim_constraints_(C, Sin,Snxt), 761 sim_constraints_(Cs,Snxt,Sout). 762sim_constraints_(C,Sin,Sout) :- 763 sim_constraint_(C,SC), 764 constraint(SC,Sin,Sout). % from simplex 765 766sim_constraint_(C,SC) :- 767 C=..[Op,LHS,RHS], % decompose 768 constraint_op(Op,COp), % acceptable to simplex 769 number(RHS), RHS >= 0, % requirement of simplex 770 map_simplex_(LHS,T/T,Sim/[]), % map to simplex list of terms 771 !, 772 SC=..[COp,Sim,RHS]. % recompose 773sim_constraint_(C,_) :- 774 fail_msg_('Illegal linear constraint: ~w',[C]). 775 776map_simplex_(T,CT/[S|Tail],CT/Tail) :- 777 map_simplex_term_(T,S), 778 !. 779map_simplex_(A+B,Cin,Cout) :- !, 780 map_simplex_(A,Cin,Cnxt), 781 map_simplex_(B,Cnxt,Cout). 782map_simplex_(A-B,Cin,Cout) :- !, 783 map_simplex_(A,Cin,Cnxt), 784 map_simplex_(-B,Cnxt,Cout). 785 786map_simplex_term_(V,1*SV) :- simplex_var_(V,SV), !. 787map_simplex_term_(-T,NN*V) :- !, 788 map_simplex_term_(T,N*V), 789 NN is -N. 790map_simplex_term_(N*V,N*SV) :- number(N), simplex_var_(V,SV), !. 791map_simplex_term_(V*N,N*SV) :- number(N), simplex_var_(V,SV). 792 793constraint_op(==,=). 794constraint_op(=<,=<). 795constraint_op(>=,>=).
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.
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... .
823% 824% local_minima/1, % apply KT constraints for objective function expression (OFE) 825% local_maxima/1, % semantically equivalent to local_minima/1 826% 827local_minima(ObjExp) :- 828 local_optima_(min,ObjExp,[]). 829 830local_maxima(ObjExp) :- 831 local_optima_(max,ObjExp,[]).
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... .
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... .
870% 871% local_minima/2, % apply KT constraints for minima with constraints 872% local_maxima/2 % apply KT constraints for maxima with constraints 873% 874local_minima(ObjExp,{Constraints}) :- 875 local_optima_(min,ObjExp,Constraints). 876 877local_maxima(ObjExp,{Constraints}) :- 878 local_optima_(max,ObjExp,Constraints). 879 880 881local_optima_(MinMax,ObjF,Constraints) :- 882 local_optima_(MinMax,ObjF,Constraints,Cs), % generate constraints 883 {Cs}. % then apply 884 885local_optima_(MinMax,ObjF,Constraints,[Constraints,GCs,LCs]) :- 886 lagrangian_(Constraints,MinMax,ObjF,LObjF,LCs), 887 term_variables((Constraints,ObjF),Vs), 888 gradient_constraints_(Vs,GCs,LObjF). 889 890gradient_constraints_([],[],_Exp). 891gradient_constraints_([X|Xs],[C|Cs],Exp) :- 892 partial_derivative(Exp,X,D), 893 (number(D) -> C=[] ; C=(D==0)), % no constraint if PD is a constant 894 gradient_constraints_(Xs,Cs,Exp). 895 896% for each constraint add to Lagrangian expression with auxiliary KKT constraints 897lagrangian_(C,MinMax,Exp,LExp, LC) :- nonvar(C), 898 kt_constraint_(C,CExp, LC), % generate langrange term with multiplier 899 lexp(MinMax,Exp,CExp,LExp), 900 !. 901lagrangian_([],_,L,L,[]). 902lagrangian_([C|Cs],MinMax,Exp,LExp,[LC|LCs]) :- 903 lagrangian_(C, MinMax,Exp,NExp,LC), 904 lagrangian_(Cs,MinMax,NExp,LExp,LCs). 905lagrangian_((C,Cs),MinMax,Exp,LExp,[LC|LCs]) :- 906 lagrangian_(C,MinMax,Exp,NExp,LC), 907 lagrangian_(Cs,MinMax,NExp,LExp,LCs). 908 909lexp(min,Exp,CExp,Exp+CExp). 910lexp(max,Exp,CExp,Exp-CExp). 911 912kt_constraint_(LHS == RHS, M*(LHS-RHS), []) :- 913 M::real. % finite multiplier only 914kt_constraint_(LHS =< RHS, MGx, MGx==0) :- 915 MGx = M*(LHS-RHS), M::real(0,_). % positive multiplier and KKT non-negativity condition 916kt_constraint_(LHS >= RHS, Exp, LC) :- 917 kt_constraint_(RHS =< LHS, Exp, LC). % >= convert to =<
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 directorydocs/).Documentation for exported predicates follows. The "custom" types include:
clpBNRattribute