26:- module(clpBNR_toolkit, 27 [
28 iterate_until/3, 29 mid_split_one/1, 30 mid_split/1, 31 taylor_contractor/2, 32 taylor_merged_contractor/2, 33 cf_contractor/2, 34
35 lin_minimum/3, 36 lin_maximum/3, 37 lin_minimize/3, 38 lin_maximize/3, 39
40 local_minima/1, 41 local_maxima/1, 42 local_minima/2, 43 local_maxima/2 44 ]). 45
46:- use_module(library(apply),[maplist/3]). 47:- use_module(library(clpBNR)). 48:- use_module(library(simplex)). 49
51fail_msg_(FString,Args) :-
52 debug(clpBNR,FString,Args), fail.
53
57iterate_until(N,Test,Goal) :- N>0, !,
58 Goal,
59 N1 is N-1,
60 (Test
61 -> true
62 ; iterate_until(N1,Test,Goal)
63 ).
64iterate_until(_N,_,_). 65
66sandbox:safe_meta(clpBNR_toolkit:iterate_until(_N,Test,Goal), [Test, Goal]).
67
71mid_split_one(Xs) :-
72 select_split(Xs,X), 73 mid_split(X). 77mid_split(X) :- small(X), !. 78mid_split(X) :-
79 M is midpoint(X),
80 ({X=<M} ; {M=<X}).
84select_split([X],X) :- !. 85select_split([X1,X2|Xs],X) :- 86 delta(X1,D1),
87 delta(X2,D2),
88 (D1 >= D2
89 -> select_split([X1|Xs],X)
90 ; select_split([X2|Xs],X)
91 ).
106cf_contractor(Xs,As) :-
107 findall(Ds,(maplist(bind_to_midpoint,Xs,As),maplist(cf_domain,Xs,Ds)),[XDs]),
108 maplist(set_domain,Xs,XDs).
109
110bind_to_midpoint(X,A) :- A is float(midpoint(X)).
111
112cf_domain(X,D) :-
113 number(X) -> D = X ; domain(X,D). % in case X narrowed to a point
114
115set_domain(X,D) :-
116 number(D) -> X = D ; X::D.
120taylor_contractor({E1=:=E2},CF) :-
121 taylor_contractor({E1==E2},CF).
122taylor_contractor({E1==E2},cf_contractor(Xs,As)) :-
123 Exp=E1-E2,
124 term_variables(Exp,Xs), % original arguments, bound to TXs on call
125 make_EQ_(Exp,TEQ), % original constraint with arguments
126 % build constraint list starting with Z's and ending with TEQ and DEQ ()
127 T::real(0,1),
128 make_As_and_Zs_(Xs,T,As,Zs,Cs,[TEQ,DEQ]), % T per Z
129 % now build Taylor constraint, DEQ
130 copy_term_nat(Exp,AExp), % copy of original constraint with As
131 term_variables(AExp,As),
132 sum_diffs(Xs, As, Zs, Zs, Exp, AExp, DEQ), % add on D(Z)'s'
133 % make any vars in original equation and contractor arguments finite real intervals
134 !,
135 Xs::real, % all vars are intervals
136 {Cs}. 137taylor_contractor(Eq,_) :-
138 fail_msg_('Invalid equation for Taylor contractor: ~w',[Eq]).
139
140make_As_and_Zs_([],_,[],[],Tail,Tail).
141make_As_and_Zs_([X|Xs],T,[A|As],[Z|Zs],[Z==A+T*(X-A)|CZs],Tail) :-
142 make_As_and_Zs_(Xs,T,As,Zs,CZs,Tail).
143
144sum_diffs([], [], [], _AllZs, _Exp, ExpIn, EQ) :- make_EQ_(ExpIn,EQ).
145sum_diffs([X|Xs], [A|As], [Z|Zs], AllZs, Exp, AExp, DEQ) :-
146 copy_term_nat(Exp,NExp), 147 term_variables(NExp,AllZs),
148 partial_derivative(NExp,Z,DZ), 149 sum_diffs(Xs, As, Zs, AllZs, Exp, AExp+DZ*(X-A), DEQ).
150
152make_EQ_(Exp,LHS==RHS) :- 153 make_EQ_(Exp,LHS,RHS).
154
155make_EQ_(E,E,0) :- var(E), !.
156make_EQ_(X+Y,X,SY) :- number(Y), !, SY is -Y.
157make_EQ_(X-Y,X,Y) :- number(Y), !.
158make_EQ_(X+Y,Y,SX) :- number(X), !, SX is -X.
159make_EQ_(X-Y,SY,SX) :- number(X), !, SX is -X, negate_sum_(Y,SY).
160make_EQ_(X+Y,LHS+Y,RHS) :- !, make_EQ_(X,LHS,RHS).
161make_EQ_(X-Y,LHS-Y,RHS) :- !, make_EQ_(X,LHS,RHS).
162make_EQ_(E,E,0). 163
164negate_sum_(Y,-Y) :- var(Y), !.
165negate_sum_(X+Y,NX-Y) :- !, negate_sum_(X,NX).
166negate_sum_(X-Y,NX+Y) :- !, negate_sum_(X,NX).
167negate_sum_(E,-E).
171taylor_merged_contractor({Es},T) :-
172 cf_list(Es,Ts),
173 cf_merge(Ts,T).
174
175cf_list([],[]) :- !.
176cf_list([C|Cs],[CF|CFs]) :- !,
177 cf_list(C, CF),
178 cf_list(Cs,CFs).
179cf_list((C,Cs),[CF|CFs]) :- !,
180 cf_list(C, CF),
181 cf_list(Cs,CFs).
182cf_list(C,CF) :-
183 taylor_contractor({C},CF).
184
185cf_merge(CFs,CF) :- cf_merge(CFs,cf_contractor([],[]),CF).
186
187cf_merge([],CF,CF).
188cf_merge([CF|CFs],CFIn,CFOut) :-
189 cf_merge(CF,CFIn,CFNxt),
190 cf_merge(CFs,CFNxt,CFOut).
191cf_merge(cf_contractor(Xs,As),cf_contractor(XsIn,AsIn),cf_contractor(XsOut,AsOut)) :-
192 cf_add(Xs,As,XsIn,AsIn,XsOut,AsOut).
193
194cf_add([],[],Xs,As,Xs,As).
195cf_add([X|Xs],[A|As],XsIn,AsIn,XsOut,AsOut) :-
196 var_existing(XsIn,AsIn,X,A), !,
197 cf_add(Xs,As,XsIn,AsIn,XsOut,AsOut).
198cf_add([X|Xs],[A|As],XsIn,AsIn,XsOut,AsOut) :-
199 cf_add(Xs,As,[X|XsIn],[A|AsIn],XsOut,AsOut).
200
201var_existing([Xex|Xs],[Aex|As], X,A) :- Xex==X -> Aex=A ; var_existing(Xs,As,X,A).
202
207lin_minimum(ObjF,{Constraints},MinValue) :-
208 lin_minimum_(ObjF,{Constraints},MinValue,false).
209
210lin_maximum(ObjF,{Constraints},MinValue) :-
211 lin_maximum_(ObjF,{Constraints},MinValue,false).
216lin_minimize(ObjF,{Constraints},MinValue) :-
217 lin_minimum_(ObjF,{Constraints},MinValue,true).
218
219lin_maximize(ObjF,{Constraints},MinValue) :-
220 lin_maximum_(ObjF,{Constraints},MinValue,true).
221
222
223lin_minimum_(ObjF,{Constraints},MinValue,BindVars) :-
224 init_simplex_(ObjF,Constraints,Objective,S0,Vs),
225 (minimize(Objective,S0,S)
226 -> objective(S,MinValue), {ObjF == MinValue},
227 (BindVars == true
228 -> bind_vars_(Vs,S)
229 ; remove_names_(Vs),
230 {Constraints} 231 )
232 ; fail_msg_('Failed to minimize: ~w',[ObjF])
233 ).
234
235lin_maximum_(ObjF,{Constraints},MaxValue,BindVars) :-
236 init_simplex_(ObjF,Constraints,Objective,S0,Vs),
237 (maximize(Objective,S0,S)
238 -> objective(S,MaxValue), {ObjF == MaxValue},
239 (BindVars == true
240 -> bind_vars_(Vs,S)
241 ; remove_names_(Vs),
242 {Constraints} 243 )
244 ; fail_msg_('Failed to maximize: ~w',[ObjF])
245 ).
246
247init_simplex_(ObjF,Constraints,Objective,S,Vs) :-
248 gen_state(S0),
249 term_variables((ObjF,Constraints),Vs),
250 (Vs = []
251 -> fail_msg_('No variables in expression to optimize: ~w',[ObjF])
252 ; sim_constraints_(Constraints,S0,S1),
253 constrain_ints_(Vs,S1,S),
254 (map_simplex_(ObjF,T/T,Objective/[])
255 -> true
256 ; fail_msg_('Illegal linear objective: ~w',[ObjF])
257 )
258 ).
259
261simplex_var_(V,SV) :- var(V),
262 (get_attr(V,clpBNR_toolkit,SV)
263 -> true
264 ; statistics(inferences,VID), SV = var(VID), put_attr(V,clpBNR_toolkit,SV)
265 ).
266
268remove_names_([]).
269remove_names_([V|Vs]) :-
270 del_attr(V,clpBNR_toolkit),
271 remove_names_(Vs).
272
273attr_unify_hook(var(_),_). 274
275
276constrain_ints_([],S,S).
277constrain_ints_([V|Vs],Sin,Sout) :-
278 % Note: library(simplex) currently constrains all variables to be non-negative
279 simplex_var_(V,SV), % corresponding simplex variable name
280 % if not already an interval, make it one with finite non-negative value
281 (domain(V,D) -> true ; V::real(0,_), domain(V,D)),
282 (D == boolean -> Dom = integer(0,1); Dom = D),
283 Dom =.. [Type,L,H],
284 (Type == integer -> constraint(integral(SV),Sin,S1) ; S1 = Sin),
285 (L < 0
286 -> % apply non-negativity condition
287 ({V >= 0} -> L1 = 0 ; fail_msg_('Negative vars not supported by \'simplex\': ~w',[V]))
288 ; L1 = L
289 ),
290 % simplex requires all vars be explicitly constrained
291 (L1 > -inf -> constraint([SV] >= L1,S1,S1a) ; S1a = S1),
292 (H < inf -> constraint([SV] =< H,S1a,S2) ; S2 = S1a),
293 constrain_ints_(Vs,S2,Sout).
294
295bind_vars_([],_).
296bind_vars_([V|Vs],S) :-
297 298 (simplex_var_(V,SV) -> variable_value(S,SV,V) ; true),
299 bind_vars_(Vs,S).
300
302sim_constraints_([],S,S) :- !.
303sim_constraints_([C|Cs],Sin,Sout) :- !,
304 sim_constraints_(C, Sin,Snxt),
305 sim_constraints_(Cs,Snxt,Sout).
306sim_constraints_((C,Cs),Sin,Sout) :- !,
307 sim_constraints_(C, Sin,Snxt),
308 sim_constraints_(Cs,Snxt,Sout).
309sim_constraints_(C,Sin,Sout) :-
310 sim_constraint_(C,SC),
311 constraint(SC,Sin,Sout). 312
313sim_constraint_(C,SC) :-
314 C=..[Op,LHS,RHS], 315 constraint_op(Op,COp), 316 number(RHS), RHS >= 0, 317 map_simplex_(LHS,T/T,Sim/[]), 318 !,
319 SC=..[COp,Sim,RHS]. 320sim_constraint_(C,_) :-
321 fail_msg_('Illegal linear constraint: ~w',[C]).
322
323map_simplex_(T,CT/[S|Tail],CT/Tail) :-
324 map_simplex_term_(T,S),
325 !.
326map_simplex_(A+B,Cin,Cout) :- !,
327 map_simplex_(A,Cin,Cnxt),
328 map_simplex_(B,Cnxt,Cout).
329map_simplex_(A-B,Cin,Cout) :- !,
330 map_simplex_(A,Cin,Cnxt),
331 map_simplex_(-B,Cnxt,Cout).
332
333map_simplex_term_(V,1*SV) :- simplex_var_(V,SV), !.
334map_simplex_term_(-T,NN*V) :- !,
335 map_simplex_term_(T,N*V),
336 NN is -N.
337map_simplex_term_(N*V,N*SV) :- number(N), simplex_var_(V,SV), !.
338map_simplex_term_(V*N,N*SV) :- number(N), simplex_var_(V,SV).
339
340constraint_op(==,=).
341constraint_op(=<,=<).
342constraint_op(>=,>=).
343
348local_minima(ObjExp) :-
349 local_optima_(min,ObjExp,[]).
350
351local_maxima(ObjExp) :-
352 local_optima_(max,ObjExp,[]).
357local_minima(ObjExp,{Constraints}) :-
358 local_optima_(min,ObjExp,Constraints).
359
360local_maxima(ObjExp,{Constraints}) :-
361 local_optima_(max,ObjExp,Constraints).
362
363
364local_optima_(MinMax,ObjF,Constraints) :-
365 local_optima_(MinMax,ObjF,Constraints,Cs), 366 {Cs}. 367
368local_optima_(MinMax,ObjF,Constraints,[Constraints,GCs,LCs]) :-
369 lagrangian_(Constraints,MinMax,ObjF,LObjF,LCs),
370 term_variables((Constraints,ObjF),Vs),
371 gradient_constraints_(Vs,GCs,LObjF).
372
373gradient_constraints_([],[],_Exp).
374gradient_constraints_([X|Xs],[C|Cs],Exp) :-
375 partial_derivative(Exp,X,D),
376 (number(D) -> C=[] ; C=(D==0)), 377 gradient_constraints_(Xs,Cs,Exp).
378
380lagrangian_(C,MinMax,Exp,LExp, LC) :- nonvar(C),
381 kt_constraint_(C,CExp, LC), 382 lexp(MinMax,Exp,CExp,LExp),
383 !.
384lagrangian_([],_,L,L,[]).
385lagrangian_([C|Cs],MinMax,Exp,LExp,[LC|LCs]) :-
386 lagrangian_(C, MinMax,Exp,NExp,LC),
387 lagrangian_(Cs,MinMax,NExp,LExp,LCs).
388lagrangian_((C,Cs),MinMax,Exp,LExp,[LC|LCs]) :-
389 lagrangian_(C,MinMax,Exp,NExp,LC),
390 lagrangian_(Cs,MinMax,NExp,LExp,LCs).
391
392lexp(min,Exp,CExp,Exp+CExp).
393lexp(max,Exp,CExp,Exp-CExp).
394
395kt_constraint_(LHS == RHS, M*(LHS-RHS), []) :-
396 M::real. % finite multiplier only
397kt_constraint_(LHS =< RHS, MGx, MGx==0) :-
398 MGx = M*(LHS-RHS), M::real(0,_). 399kt_constraint_(LHS >= RHS, Exp, LC) :-
400 kt_constraint_(RHS =< LHS, Exp, LC).