1/* The MIT License (MIT) 2 * 3 * Copyright (c) 2022-2026 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 taylor_solve/1, taylor_solve/2, % combines taylor_contractor/2 and cf_solve/1 33 34 integrate/3, integrate/4, % simple numerical integration 35 boundary_values/2, boundary_values/3, boundary_values/4, % solve boundary value problems 36 37 lin_minimum/3, % find minimum of linear problem using library(simplex) 38 lin_maximum/3, % find maximum of linear problem using library(simplex) 39 lin_minimize/3, % lin_minimum/3 plus bind vars to solution minimizers 40 lin_maximize/3, % lin_maximum/3 plus bind vars to solution maximizers 41 42 local_minima/1, % apply KT constraints for objective function expression (OFE) 43 local_maxima/1, % semantically equivalent to local_minima/1 44 local_minima/2, % apply KT constraints for minima with constraints 45 local_maxima/2 % apply KT constraints for maxima with constraints 46 ]).
58% module dependencies: 59:- use_module(library(lists),[member/2]). 60:- use_module(library(apply),[maplist/3]). 61:- use_module(library(apply_macros)). % compiler support for `maplist`, helps a bit 62:- use_module(library(clpBNR)). 63:- use_module(library(simplex)). % finding optima of linear systems 64 65% sandboxing for SWISH 66:- multifile(sandbox:safe_primitive/1). 67 68% message support before optimise 69debug_msg_(FString,Args) :- debug(clpBNR,FString,Args). 70 71fail_msg_(FString,Args) :- debug_msg_(FString,Args), fail. 72 73:- 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.)
*/
97% 98% General purpose iterator: execute Goal a maximum of N times or until Test succeeds 99% 100iterate_until(N,Test,Goal) :- N>0, !, 101 , 102 N1 is N-1, 103 ( 104 -> true 105 ; iterate_until(N1,Test,Goal) 106 ). 107iterate_until(_N,_,_). % non-positive N --> exit 108 109sandbox:safe_meta(clpBNR_toolkit:iterate_until(_N,Test,Goal), [Test, Goal]).
mid_split for details of interval splitting for this predicate.
118mid_split_one(Xs) :- 119 predsort(delta_order_,Xs,[X|_]), % select largest interval with largest width 120 mid_split(X). % split it 121 122% used with predsort/2 to order intervals in decreasing width 123delta_order_(COp,X1,X2) :- 124 delta(X1,D1), delta(X2,D2), 125 % equal delta isn't duplicate, note descending order 126 (D2 > D1 -> COp = > ; COp = <). % never '='
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.
143mid_split(X) :-
144 (number(X) % optimise number case
145 -> true
146 ; (small(X)
147 -> true
148 ; midpoint(X,M), % fails if not an interval
149 ({X=<M} ; {M=<X}) % possible choicepoint, Note: can split on solution leaving CP
150 )
151 ).160taylor_solve(C) :- 161 current_prolog_flag(clpBNR_default_precision,P), 162 taylor_solve(C,P). 163 164taylor_solve(C,P) :- 165 taylor_contractor(C,T), 166 cf_solve(T,P). % leaves CP's for multiple solutions
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.
175% 176% centred form contractor 177% 178% Bind the values of As to the midpoints of Xs. To support repetitive application 179% of the contractor (required by the iterator), the contractor should not permanently 180% bind anything so findall/3 will be used to achieve this "forward checking" 181% (as suggested in [CLIP]). After the call to findall, the bounds of the resulting list 182% of narrowed domains (XDs) are then applied to Xs. 183% 184% This contractor can be used with any "centred form", e.g., Newton or Krawczyk, since it 185% only depends on intervals and their midpoints, hence its name `cf_contractor`. The 186% details which distinguish the variety of centred form are built into the variables' 187% constraints. 188% 189cf_contractor(Xs,As) :- 190 findall(Ds,(maplist(bind_to_midpoint,Xs,As),maplist(cf_domain,Xs,Ds)),[XDs]), 191 maplist(set_domain,Xs,XDs). 192 193bind_to_midpoint(X,A) :- A is float(midpoint(X)). 194 195cf_domain(X,D) :- 196 number(X) -> D = X ; domain(X,D). % in case X narrowed to a point 197 198set_domain(X,D) :- 199 number(D) -> X = D 200 ; domain(X,D) -> true
201 ; X::D
201. 202
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.
235cf_solve(T) :- 236 current_prolog_flag(clpBNR_default_precision,P), 237 cf_solve(T,P). 238cf_solve(cf_contractor(Xs,As),P) :- 239 cf_iterate_(Xs,As,P). 240 241cf_iterate_(Xs, _,P) :- small(Xs,P), !. % all vars small, normal exit 242cf_iterate_(Xs,As,P) :- 243 cf_contractor(Xs,As), % execute contractor, must succeed 244 predsort(delta_order_,Xs,SXs), % sort intervals by width 245 (cf_split(SXs,P,X,Pt) % find a non-small interval with a split point and 246 -> ({X=<Pt} ; {Pt=<X}), % split, then 247 cf_iterate_(Xs,As,P) % iterate 248 ; true % nothing splittable, just exit 249 ). 250 251cf_split(SXs,P,X,Pt) :- % semi deterministic 252 member(X,SXs), % generate intervals from widest to narrowest 253 \+ small(X,P), % large enough 254 real_split_point(X,P,Pt), 255 !. % cut generator
== 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.
292% 293% build a cf_contractor for a multivariate expression based on Taylor expansion 294% 295taylor_contractor({E1=:=E2},CF) :- !, 296 taylor_contractor({E1==E2},CF). 297taylor_contractor({E1==E2},cf_contractor(Xs,As)) :- 298 Exp=E1-E2, 299 term_variables(Exp,Xs), % original arguments, bound to TXs on call 300 make_EQ_(Exp,TEQ), % original constraint with arguments 301 % build constraint list starting with Z's and ending with TEQ and DEQ () 302 T::real(0,1), 303 make_As_and_Zs_(Xs,T,As,Zs,Cs,[TEQ,DEQ]), % T per Z 304 % now build Taylor constraint, DEQ 305 copy_term_nat(Exp,AExp), % copy of original constraint with As 306 term_variables(AExp,As), 307 sum_diffs(Xs, As, Zs, Zs, Exp, AExp, DEQ), % add on D(Z)'s' 308 % make any vars in original equation and contractor arguments finite real intervals 309 !, 310 Xs::real, % all vars are intervals 311 {Cs}. % apply constraints 312taylor_contractor({Es},CF) :- (list(Es) ; Es = (_,_)), !, % list or sequence 313 taylor_merged_contractor({Es},CF). 314taylor_contractor({C},cf_contractor([],[])) :- 315 C=..[Op,_,_], memberchk(Op,[<, =<, >=, >]), !, % inequalities ... 316 {C}. % are just constraints 317taylor_contractor(Eq,_) :- 318 fail_msg_('Invalid constraint for Taylor contractor: ~w',[Eq]). 319 320make_As_and_Zs_([],_,[],[],Tail,Tail). 321make_As_and_Zs_([X|Xs],T,[A|As],[Z|Zs],[Z==A+T*(X-A)|CZs],Tail) :- 322 make_As_and_Zs_(Xs,T,As,Zs,CZs,Tail). 323 324sum_diffs([], [], [], _AllZs, _Exp, ExpIn, EQ) :- make_EQ_(ExpIn,EQ). 325sum_diffs([X|Xs], [A|As], [Z|Zs], AllZs, Exp, AExp, DEQ) :- 326 copy_term_nat(Exp,NExp), % copy expression and replace Xs by Zs 327 term_variables(NExp,AllZs), 328 partial_derivative(NExp,Z,DZ), % differentiate wrt. Z and add to generated expression 329 sum_diffs(Xs, As, Zs, AllZs, Exp, AExp+DZ*(X-A), DEQ). 330 331% map expression Exp to an equation equivalent to Exp==0 with numeric RHS 332make_EQ_(Exp,LHS==RHS) :- % turn expression into equation equivalent to Exp==0. 333 make_EQ_(Exp,LHS,RHS). 334 335make_EQ_(E,E,0) :- var(E), !. 336make_EQ_(X+Y,X,SY) :- number(Y), !, SY is -Y. 337make_EQ_(X-Y,X,Y) :- number(Y), !. 338make_EQ_(X+Y,Y,SX) :- number(X), !, SX is -X. 339make_EQ_(X-Y,SY,SX) :- number(X), !, SX is -X, negate_sum_(Y,SY). 340make_EQ_(X+Y,LHS+Y,RHS) :- !, make_EQ_(X,LHS,RHS). 341make_EQ_(X-Y,LHS-Y,RHS) :- !, make_EQ_(X,LHS,RHS). 342make_EQ_(E,E,0). % default (non +/- subexpression) 343 344negate_sum_(Y,-Y) :- var(Y), !. 345negate_sum_(X+Y,NX-Y) :- !, negate_sum_(X,NX). 346negate_sum_(X-Y,NX+Y) :- !, negate_sum_(X,NX). 347negate_sum_(E,-E).
== or =:=) constraint in Constraints; otherwise fails.
356% 357% build a cf_contractor by merging a list of cf_contractor's 358% 359taylor_merged_contractor({Es},T) :- 360 cf_list(Es,Ts), 361 cf_merge(Ts,T). 362 363cf_list([],[]) :- !. 364cf_list([C|Cs],[CF|CFs]) :- !, 365 cf_list(C, CF), 366 cf_list(Cs,CFs). 367cf_list((C,Cs),[CF|CFs]) :- !, 368 cf_list(C, CF), 369 cf_list(Cs,CFs). 370cf_list(C,CF) :- 371 taylor_contractor({C},CF). 372 373cf_merge(CFs,CF) :- cf_merge(CFs,cf_contractor([],[]),CF). 374 375cf_merge([],CF,CF). 376cf_merge([CF|CFs],CFIn,CFOut) :- 377 cf_merge(CF,CFIn,CFNxt), 378 cf_merge(CFs,CFNxt,CFOut). 379cf_merge(cf_contractor(Xs,As),cf_contractor(XsIn,AsIn),cf_contractor(XsOut,AsOut)) :- 380 cf_add(Xs,As,XsIn,AsIn,XsOut,AsOut). 381 382cf_add([],[],Xs,As,Xs,As). 383cf_add([X|Xs],[A|As],XsIn,AsIn,XsOut,AsOut) :- 384 var_existing(XsIn,AsIn,X,A), !, 385 cf_add(Xs,As,XsIn,AsIn,XsOut,AsOut). 386cf_add([X|Xs],[A|As],XsIn,AsIn,XsOut,AsOut) :- 387 cf_add(Xs,As,[X|XsIn],[A|AsIn],XsOut,AsOut). 388 389var_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.
*/
417% integrate(F,X,R) where X is an interval over which to integrate and F = f(X) 418% Note that integrate(F,X,R) is equivalent to boundary_values(X,[dV(_,F,0,R]) 419integrate(F,X,R) :- 420 current_prolog_flag(clpBNR_default_precision,P), 421 integrate(F,X,R,P,_). 422integrate(F,X,R,P) :- 423 integrate(F,X,R,P,_). 424integrate(F,X,R,P,Steps) :- % internal arity 5 for development 425 interval(X), % X must be an interval 426 integer(P), P>0, % P must be positive integer 427 !, % args OK, commit 428 boundary_values(X,[dV(_,F,0,R)],P,Steps). % use integration in `boundary_values` 429integrate(F,X,R,P,_Steps) :- 430 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). */
480boundary_values(X,YDefs) :- 481 current_prolog_flag(clpBNR_default_precision,P), 482 boundary_values_(X,YDefs,P,_). 483 484boundary_values(X,YDefs,P) :- 485 boundary_values_(X,YDefs,P,_). 486 487boundary_values(X,YDefs,P,Steps) :- 488 boundary_values_(X,YDefs,P,Steps/_). 489 490boundary_values_(X,YDefs,P,Out/[(Xf,Yfs)]) :- 491 integer(P), P>0, % P must be positive integer 492 domain(X,Xdom), Xdom =.. [_Type,Xi,Xf], % X must be an interval 493 eval_dvars(YDefs,Ys,Yis,Yfs,Ydoms), % too many args for maplist 494 maplist(fXY_lambda(X,Ys),YDefs,Fxys), % list of partial derivative lambda args 495 (maplist(total_derivative_(Fxys),Fxys,DFxys) % list of connective derivative lambda args 496 -> true 497 ; DFxys = none % F non-differentiable(?), use euler step 498 ), 499 !, % args all good, commit 500 integrate_(P,Fxys,DFxys,(Xi,Yis),(Xf,Yfs),Ydoms,Out/[(Xf,Yfs)]). 501boundary_values_(X,YDefs,P,_) :- 502 fail_msg_('Invalid argument(s): ~w',[boundary_values(X,YDefs,P)]). 503 504eval_dvars([],[],[],[],[]). 505eval_dvars([dV(Y, _PD, Lexp, Uexp)|YDefs],[Y|Ys],[Yi|Yis],[Yf|Yfs],[Ydom|Ydoms]) :- 506 (domain(Y,Ydom) -> true ; Ydom = real), % Ydom defaults to real 507 (var(Lexp) -> Yi=Lexp, Yi::Ydom ; Yi is Lexp), 508 (var(Uexp) -> Yf=Uexp, Yf::Ydom ; Yf is Uexp), 509 eval_dvars(YDefs,Ys,Yis,Yfs,Ydoms). 510 511% construct Lambda args for Fxy 512fXY_lambda(X,Ys,dV(_Y,Fxy,_,_),FxyArgs) :- 513 lambda_for_(Fxy,X,Ys,FxyArgs). 514 515% construct Lambda args for derivative function of Fxy from Lambda of Fxy 516 517total_derivative_(Fxys,_Free/Ps,DxyArgs) :- !, % ignore free variables 518 total_derivative_(Fxys,Ps,DxyArgs). 519total_derivative_(Fxys,[X,Ys,Fxy],DxyArgs) :- 520 partial_derivative(Fxy,X,DFDX), % clpBNR built-in 521 sumYpartials(Fxys,Ys,Fxy,0,DYsum), 522 simplify_sum_(DFDX, DYsum, DExp), !, 523 lambda_for_(DExp,X,Ys,DxyArgs). 524 525sumYpartials([],[],_Fxy,Acc,Acc). 526sumYpartials([_Free/FxyI|FxyIs],YIs,Fxy,Acc,Sum) :- !, 527 sumYpartials([FxyI|FxyIs],YIs,Fxy,Acc,Sum). 528sumYpartials([[_X,_Ys,FxyI]|FxyIs],[YI|YIs],Fxy,Acc,Sum) :- 529 partial_derivative(Fxy,YI,DFDYI), 530 (number(DFDYI), DFDYI =:= 0 -> NxtAcc = Acc ; simplify_sum_(Acc,FxyI*DFDYI,NxtAcc)), 531 !, 532 sumYpartials(FxyIs,YIs,Fxy,NxtAcc,Sum). 533 534simplify_sum_(X,Y,Y) :- number(X),X=:=0. 535simplify_sum_(X,Y,X) :- number(Y),Y=:=0. 536simplify_sum_(X,Y,X+Y). 537 538% construct args for Lambda expression 539lambda_for_(Fxy,X,Ys,Args) :- 540 Lambda_parms = [X,Ys,Fxy], 541 term_variables(Fxy,FVs), 542 exclude(var_member_([X|Ys]),FVs,EVs), % EVs = free variables 543 (comma_op_(EVs,EV) -> Args = {EV}/Lambda_parms ; Args = Lambda_parms). 544 545var_member_([L|Ls],E) :- L==E -> true ; var_member_(Ls,E). 546 547comma_op_([X],X). % assumes use in if-then 548comma_op_([X|Xs],(X,Y)) :- comma_op_(Xs,Y). 549 550% integration loop 551integrate_(0, Fxys, Dxys, Initial, Final, Ydomains, [Initial|Ps]/Ps) :- !, 552 % select integration step 553 ( Dxys == none 554 -> step_euler(Fxys, Initial, Final, Ydomains) 555 ; step_trap(Fxys, Dxys, Initial, Final, Ydomains) 556 ). 557integrate_(P, Fxys, Dxys, Initial, Final, Ydomains, L/E) :- 558 % create interpolation point and integrate two halves 559 interpolate_(Initial, Final, Ydomains, Middle), 560 Pn is P - 1, 561 integrate_(Pn, Fxys, Dxys, Initial, Middle, Ydomains, L/M), 562 integrate_(Pn, Fxys, Dxys, Middle, Final, Ydomains, M/E). 563 564interpolate_((X0,_Y0s), (X1,_Y1s), Ydomains, (XM,YMs)) :- 565 XM is (X0 + X1)/2, % XM is midpoint of (X0,X1) 566 maplist(::,YMs,Ydomains). % corresponding YMs 567 568step_euler(Fxys, (X0,Y0), (X1,Y1), Ydoms) :- 569 X01:: real(X0,X1), % range of X in step 570 maplist(lambda_constrain_(X01,Y01),Fxys,Fs), % approx f' over X0 571 Dx is X1 - X0, % assumed (X1>X0) 572 DX :: real(0,Dx), % range for estimate 573 euler_constraints(Y0,Y1,Y01,Ydoms,Dx,DX,Fs,In/In,Cs/[]), % flatten with diff list 574 {Cs}. 575 576step_trap(Fxys, Dxys, (X0,Y0), (X1,Y1), Ydoms) :- 577 X01:: real(X0,X1), % range of X in step 578 maplist(lambda_constrain_(X0,Y0),Fxys,F0s), % F0s = slopes at X0 579 maplist(lambda_constrain_(X1,Y1),Fxys,F1s), % F1s = slopes at X1 580 maplist(lambda_constrain_(X01,Y01),Dxys,Ds), % approx f' over X0 581 Dx is X1 - X0, % assumed (X1>X0) 582 DX :: real(0,Dx), % range for estimate 583 trap_constraints(Y0,Y1,Y01,Ydoms,Dx,DX,F0s,F1s,Ds,In/In,Cs/[]), % flatten with diff list 584 {Cs}. %%, absolve(Y1,2). 585 586lambda_constrain_(X,Ys,Args,F) :- % reorder args for yall: >> 587 yall: >>(Args,true,X,Ys,F). % avoid meta-call (basically just makes copy) 588% known safe since lambda Goal=true 589sandbox:safe_primitive(clpBNR_toolkit:lambda_constrain_(_X,_Ys,_Args,_F)). 590 591/* see Carleton notes: 592https://www.softwarepreservation.org/projects/prolog/bnr/doc/Older-Introduction_to_CLP%28BNR%29-1995.pdf/view 593*/ 594euler_constraints([],[],[],[],_Dx,_DX,[],In,In). 595euler_constraints([Y0|Y0s],[Y1|Y1s],[Y01|Y01s],[Ydom|Ydoms],Dx,DX,[F|Fs], 596 In/[FM <= F, % FM = slope inclusion 597 Y01 - Y0 is DX*FM, 598 Y1 - Y0 is Dx*FM 599 |Cs], 600 Out) :- 601 Y01:: Ydom, 602 euler_constraints(Y0s,Y1s,Y01s,Ydoms,Dx,DX,Fs,In/Cs,Out). 603 604trap_constraints([],[],[],[],_Dx,_DX,[],[],[],In,In). 605trap_constraints([Y0|Y0s],[Y1|Y1s],[Y01|Y01s],[Ydom|Ydoms],Dx,DX,[F0|F0s],[F1|F1s],[D|Ds], 606 % use `is` to circumvent `simplify` 607 In/[FM <= (F0+F1)/2, % FM = average slope using step endpoints (one-way) 608 % Note that the following must not be symbolically simplified to eliminate D 609 8*DR is D - D, % 4*deltaR == (D-D)/2 (for error term) 610 Y01 - Y0 is DX*(FM + DR*DX), 611 Y1 - Y0 is Dx*(FM + DR*Dx) 612 |Cs], 613 Out) :- 614 Y01:: Ydom, DR::real, 615 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.
652lin_minimum(ObjF,{Constraints},MinValue) :- 653 lin_minimum_(ObjF,{Constraints},MinValue,false). 654 655lin_maximum(ObjF,{Constraints},MinValue) :- 656 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.
679lin_minimize(ObjF,{Constraints},MinValue) :- 680 lin_minimum_(ObjF,{Constraints},MinValue,true). 681 682lin_maximize(ObjF,{Constraints},MinValue) :- 683 lin_maximum_(ObjF,{Constraints},MinValue,true). 684 685 686lin_minimum_(ObjF,{Constraints},MinValue,BindVars) :- 687 init_simplex_(ObjF,Constraints,Objective,S0,Vs), 688 (minimize(Objective,S0,S) 689 -> objective(S,MinValue), {ObjF == MinValue}, 690 (BindVars == true 691 -> bind_vars_(Vs,S) 692 ; remove_names_(Vs), 693 {Constraints} % apply constraints 694 ) 695 ; fail_msg_('Failed to minimize: ~w',[ObjF]) 696 ). 697 698lin_maximum_(ObjF,{Constraints},MaxValue,BindVars) :- 699 init_simplex_(ObjF,Constraints,Objective,S0,Vs), 700 (maximize(Objective,S0,S) 701 -> objective(S,MaxValue), {ObjF == MaxValue}, 702 (BindVars == true 703 -> bind_vars_(Vs,S) 704 ; remove_names_(Vs), 705 {Constraints} % apply constraints 706 ) 707 ; fail_msg_('Failed to maximize: ~w',[ObjF]) 708 ). 709 710init_simplex_(ObjF,Constraints,Objective,S,Vs) :- 711 gen_state(S0), 712 term_variables((ObjF,Constraints),Vs), 713 (Vs = [] 714 -> fail_msg_('No variables in expression to optimize: ~w',[ObjF]) 715 ; sim_constraints_(Constraints,S0,S1), 716 _::real(_,Max), % max value to constrain for simplex 717 constrain_ints_(Vs,Max,S1,S), 718 (map_simplex_(ObjF,T/T,Objective/[]) 719 -> true 720 ; fail_msg_('Illegal linear objective: ~w',[ObjF]) 721 ) 722 ). 723 724% use an attribute to associate a var with a unique simplex variable name 725simplex_var_(V,SV) :- var(V), 726 (get_attr(V,clpBNR_toolkit,SV) 727 -> true 728 ; statistics(inferences,VID), SV = var(VID), put_attr(V,clpBNR_toolkit,SV) 729 ). 730 731% Name attribute removed before exit. 732remove_names_([]). 733remove_names_([V|Vs]) :- 734 del_attr(V,clpBNR_toolkit), 735 remove_names_(Vs). 736 737attr_unify_hook(var(_),_). % unification always does nothing and succeeds 738 739constrain_ints_([],_,S,S). 740constrain_ints_([V|Vs],Max,Sin,Sout) :- 741 % Note: library(simplex) currently constrains all variables to be non-negative 742 simplex_var_(V,SV), % corresponding simplex variable name 743 % if not already an interval, make it one with finite non-negative value 744 (domain(V,D) -> true ; V::real(0,_), domain(V,D)), 745 (D == boolean -> Dom = integer(0,1); Dom = D), 746 Dom =.. [Type,L,H], 747 (Type == integer -> constraint(integral(SV),Sin,S1) ; S1 = Sin), 748 (L < 0 749 -> % apply non-negativity condition 750 ({V >= 0} -> L1 = 0 ; fail_msg_('Negative vars not supported by \'simplex\': ~w',[V])) 751 ; L1 = L 752 ), 753 % explicitly constrain any vars not (0,Max-eps) 754 (L1 > 0 -> constraint([SV] >= L1,S1,S2) ; S2 = S1), % L1 not negative 755 (H < Max -> constraint([SV] =< H,S2,SNxt) ; SNxt = S2), 756 constrain_ints_(Vs,Max,SNxt,Sout). 757 758bind_vars_([],_). 759bind_vars_([V|Vs],S) :- 760 % Note: skip anything nonvar (already bound due to active constraints) 761 (simplex_var_(V,SV) -> variable_value(S,SV,V) ; true), 762 bind_vars_(Vs,S). 763 764% clpBNR constraints have already been applied so worst errors have been detected 765sim_constraints_([],S,S) :- !. 766sim_constraints_([C|Cs],Sin,Sout) :- !, 767 sim_constraints_(C, Sin,Snxt), 768 sim_constraints_(Cs,Snxt,Sout). 769sim_constraints_((C,Cs),Sin,Sout) :- !, 770 sim_constraints_(C, Sin,Snxt), 771 sim_constraints_(Cs,Snxt,Sout). 772sim_constraints_(C,Sin,Sout) :- 773 sim_constraint_(C,SC), 774 constraint(SC,Sin,Sout). % from simplex 775 776sim_constraint_(C,SC) :- 777 C=..[Op,LHS,RHS], % decompose 778 constraint_op(Op,COp), % acceptable to simplex 779 number(RHS), RHS >= 0, % requirement of simplex 780 map_simplex_(LHS,T/T,Sim/[]), % map to simplex list of terms 781 !, 782 SC=..[COp,Sim,RHS]. % recompose 783sim_constraint_(C,_) :- 784 fail_msg_('Illegal linear constraint: ~w',[C]). 785 786map_simplex_(T,CT/[S|Tail],CT/Tail) :- 787 map_simplex_term_(T,S), 788 !. 789map_simplex_(A+B,Cin,Cout) :- !, 790 map_simplex_(A,Cin,Cnxt), 791 map_simplex_(B,Cnxt,Cout). 792map_simplex_(A-B,Cin,Cout) :- !, 793 map_simplex_(A,Cin,Cnxt), 794 map_simplex_(-B,Cnxt,Cout). 795 796map_simplex_term_(V,1*SV) :- simplex_var_(V,SV), !. 797map_simplex_term_(-T,NN*V) :- !, 798 map_simplex_term_(T,N*V), 799 NN is -N. 800map_simplex_term_(N*V,N*SV) :- number(N), simplex_var_(V,SV), !. 801map_simplex_term_(V*N,N*SV) :- number(N), simplex_var_(V,SV). 802 803constraint_op(==,=). 804constraint_op(=<,=<). 805constraint_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... .
833% 834% local_minima/1, % apply KT constraints for objective function expression (OFE) 835% local_maxima/1, % semantically equivalent to local_minima/1 836% 837local_minima(ObjExp) :- 838 local_optima_(min,ObjExp,[]). 839 840local_maxima(ObjExp) :- 841 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... .
880% 881% local_minima/2, % apply KT constraints for minima with constraints 882% local_maxima/2 % apply KT constraints for maxima with constraints 883% 884local_minima(ObjExp,{Constraints}) :- 885 local_optima_(min,ObjExp,Constraints). 886 887local_maxima(ObjExp,{Constraints}) :- 888 local_optima_(max,ObjExp,Constraints). 889 890 891local_optima_(MinMax,ObjF,Constraints) :- 892 local_optima_(MinMax,ObjF,Constraints,Cs), % generate constraints 893 {Cs}. % then apply 894 895local_optima_(MinMax,ObjF,Constraints,[Constraints,GCs,LCs]) :- 896 lagrangian_(Constraints,MinMax,ObjF,LObjF,LCs), 897 term_variables((Constraints,ObjF),Vs), 898 gradient_constraints_(Vs,GCs,LObjF). 899 900gradient_constraints_([],[],_Exp). 901gradient_constraints_([X|Xs],[C|Cs],Exp) :- 902 partial_derivative(Exp,X,D), 903 (number(D) -> C=[] ; C=(D==0)), % no constraint if PD is a constant 904 gradient_constraints_(Xs,Cs,Exp). 905 906% for each constraint add to Lagrangian expression with auxiliary KKT constraints 907lagrangian_(C,MinMax,Exp,LExp, LC) :- nonvar(C), 908 kt_constraint_(C,CExp, LC), % generate langrange term with multiplier 909 lexp(MinMax,Exp,CExp,LExp), 910 !. 911lagrangian_([],_,L,L,[]). 912lagrangian_([C|Cs],MinMax,Exp,LExp,[LC|LCs]) :- 913 lagrangian_(C, MinMax,Exp,NExp,LC), 914 lagrangian_(Cs,MinMax,NExp,LExp,LCs). 915lagrangian_((C,Cs),MinMax,Exp,LExp,[LC|LCs]) :- 916 lagrangian_(C,MinMax,Exp,NExp,LC), 917 lagrangian_(Cs,MinMax,NExp,LExp,LCs). 918 919lexp(min,Exp,CExp,Exp+CExp). 920lexp(max,Exp,CExp,Exp-CExp). 921 922kt_constraint_(LHS == RHS, M*(LHS-RHS), []) :- 923 M::real. % finite multiplier only 924kt_constraint_(LHS =< RHS, MGx, MGx==0) :- 925 MGx = M*(LHS-RHS), M::real(0,_). % positive multiplier and KKT non-negativity condition 926kt_constraint_(LHS >= RHS, Exp, LC) :- 927 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