1%
    2% Toolkit of useful utilities for CLP(BNR)
    3%
    4/*	The MIT License (MIT)
    5 *
    6 *	Copyright (c) 2022-2024 Rick Workman
    7 *
    8 *	Permission is hereby granted, free of charge, to any person obtaining a copy
    9 *	of this software and associated documentation files (the "Software"), to deal
   10 *	in the Software without restriction, including without limitation the rights
   11 *	to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
   12 *	copies of the Software, and to permit persons to whom the Software is
   13 *	furnished to do so, subject to the following conditions:
   14 *
   15 *	The above copyright notice and this permission notice shall be included in all
   16 *	copies or substantial portions of the Software.
   17 *
   18 *	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
   19 *	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   20 *	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
   21 *	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   22 *	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
   23 *	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
   24 *	SOFTWARE.
   25 */
   26:- module(clpBNR_toolkit,       % SWI module declaration
   27	[
   28	iterate_until/3,            % general purpose iterator
   29	mid_split_one/1,            % contractor to split largest interval at midpoint
   30	mid_split/1,                % contractor to split an interval at midpoint
   31	taylor_contractor/2,        % build cf_contractor based on Taylor expansion
   32	taylor_merged_contractor/2, % build merged Taylor cf_contractor from list of equations
   33	cf_contractor/2,            % execute cf_contractor
   34
   35	lin_minimum/3,              % find minimum of linear problem using library(simplex)
   36	lin_maximum/3,              % find maximum of linear problem using library(simplex)
   37	lin_minimize/3,             % lin_minimum/3 plus bind vars to solution minimizers
   38	lin_maximize/3,             % lin_maximum/3 plus bind vars to solution maximizers
   39
   40	local_minima/1,             % apply KT constraints for objective function expression (OFE)
   41	local_maxima/1,             % semantically equivalent to local_minima/1
   42	local_minima/2,             % apply KT constraints for minima with constraints
   43	local_maxima/2              % apply KT constraints for maxima with constraints
   44	]).   45	
   46:- use_module(library(apply),[maplist/3]).   47:- use_module(library(clpBNR)).   48:- use_module(library(simplex)).   49
   50% messages for noisy failure
   51fail_msg_(FString,Args) :-
   52	debug(clpBNR,FString,Args), fail.
   53
   54%
   55% General purpose iterator: execute Goal a maximum of N times or until Test succeeds
   56%
   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,_,_).  % non-positive N --> exit
   65
   66sandbox:safe_meta(clpBNR_toolkit:iterate_until(_N,Test,Goal), [Test, Goal]).
   67
   68%
   69% contractor to split largest interval of a list on midpoint
   70%
   71mid_split_one(Xs) :-
   72	select_split(Xs,X),  % select largest interval with largest width
   73	mid_split(X).        % split it
   74%
   75% contractor to split a single interval on midpoint if sufficiently wide (creates a choicepoint)
   76%
   77mid_split(X) :- small(X), !. % too narrow to split
   78mid_split(X) :-
   79	M is midpoint(X),
   80	({X=<M} ; {M=<X}).
   81%
   82% select interval with largest width
   83%
   84select_split([X],X) :- !.       % select last remaining element
   85select_split([X1,X2|Xs],X) :-   % compare widths and discard one interval
   86	delta(X1,D1),
   87	delta(X2,D2),
   88	(D1 >= D2
   89	 -> select_split([X1|Xs],X)
   90	 ;  select_split([X2|Xs],X)
   91	).
   92%
   93% centred form contractor
   94%
   95% Bind the values of As to the midpoints of Xs. To support repetitive application 
   96% of the contractor (required by the iterator), the contractor should not permanently 
   97% bind anything so findall/3 will be used to achieve this "forward checking" 
   98% (as suggested in [CLIP]). After the call to findall, 'XDs' is a list of narrowed domains
   99% of Xs and are then applied to Xs.
  100%
  101% This contractor can be used with any "centred form", e.g., Newton or Krawczyk, since it
  102% only depends on intervals and their midpoints, hence its name `cf_contractor`. The 
  103% details which distinguish the variety of centred form are built into the variables' 
  104% constraints.
  105%
  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.
  117%
  118% build a cf_contractor for a multivariate expression based on Taylor expansion
  119%
  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}.      % apply constraints
  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),        % copy expression and replace Xs by Zs
  147	term_variables(NExp,AllZs),
  148	partial_derivative(NExp,Z,DZ),  % differentiate wrt. Z and add to generated expression
  149	sum_diffs(Xs, As, Zs, AllZs, Exp, AExp+DZ*(X-A), DEQ).
  150
  151% map expression Exp to an equation equivalent to Exp==0 with numeric RHS
  152make_EQ_(Exp,LHS==RHS) :-    % turn expression into equation equivalent to Exp==0.
  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).        % default (non +/- subexpression)
  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).
  168%
  169% build a cf_contractor by merging a list of cf_contractor's
  170%
  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
  203%
  204%	lin_minimum/3         % find minimum of linear problem using library(simplex)
  205%	lin_maximum/3         % find maximum of linear problem using library(simplex)
  206%
  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).
  212%
  213%	lin_minimize/3,       % lin_minimum/3 and bind vars to generated solution
  214%	lin_maximize/3,       % lin_maximum/3 and bind vars to generated solution
  215%	
  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}  % apply 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}  % apply 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
  260% use an attribute to associate a var with a unique simplex variable name
  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
  267% Name attribute removed before exit.
  268remove_names_([]).
  269remove_names_([V|Vs]) :-
  270	del_attr(V,clpBNR_toolkit),
  271	remove_names_(Vs).
  272		 
  273attr_unify_hook(var(_),_).     % unification always does nothing and succeeds
  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	% Note: skip anything nonvar (already bound due to active constraints)
  298	(simplex_var_(V,SV) -> variable_value(S,SV,V) ; true),
  299	bind_vars_(Vs,S).
  300
  301% clpBNR constraints have already been applied so worst errors have been detected
  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).  % from simplex
  312
  313sim_constraint_(C,SC) :-
  314	C=..[Op,LHS,RHS],              % decompose
  315	constraint_op(Op,COp),         % acceptable to simplex
  316	number(RHS), RHS >= 0,         % requirement of simplex
  317	map_simplex_(LHS,T/T,Sim/[]),  % map to simplex list of terms
  318	!,
  319	SC=..[COp,Sim,RHS].            % recompose
  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
  344%
  345%	local_minima/1,       % apply KT constraints for objective function expression (OFE)
  346%	local_maxima/1,       % semantically equivalent to local_minima/1
  347%
  348local_minima(ObjExp) :-
  349	local_optima_(min,ObjExp,[]).
  350
  351local_maxima(ObjExp) :-
  352	local_optima_(max,ObjExp,[]).
  353%
  354%	local_minima/2,       % apply KT constraints for minima with constraints
  355%	local_maxima/2        % apply KT constraints for maxima with constraints
  356%
  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),  % generate constraints
  366	{Cs}.                                       % then apply
  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)),  % no constraint if PD is a constant
  377	gradient_constraints_(Xs,Cs,Exp).
  378	
  379% for each constraint add to Lagrangian expression with auxiliary KKT constraints
  380lagrangian_(C,MinMax,Exp,LExp, LC) :- nonvar(C),
  381	kt_constraint_(C,CExp, LC), % generate langrange term with multiplier
  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,_).  % positive multiplier and KKT non-negativity condition
  399kt_constraint_(LHS >= RHS, Exp,         LC) :-
  400	kt_constraint_(RHS =< LHS, Exp, LC).  % >= convert to =<