1/*  $Id$
    2
    3    Part of CLP(R) (Constraint Logic Programming over Reals)
    4
    5    Author:        Leslie De Koninck
    6    E-mail:        Leslie.DeKoninck@cs.kuleuven.be
    7    WWW:           http://www.swi-prolog.org
    8		   http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09
    9    Copyright (C): 2004, K.U. Leuven and
   10		   1992-1995, Austrian Research Institute for
   11		              Artificial Intelligence (OFAI),
   12			      Vienna, Austria
   13
   14    This software is part of Leslie De Koninck's master thesis, supervised
   15    by Bart Demoen and daily advisor Tom Schrijvers.  It is based on CLP(Q,R)
   16    by Christian Holzbaur for SICStus Prolog and distributed under the
   17    license details below with permission from all mentioned authors.
   18
   19    This program is free software; you can redistribute it and/or
   20    modify it under the terms of the GNU General Public License
   21    as published by the Free Software Foundation; either version 2
   22    of the License, or (at your option) any later version.
   23
   24    This program is distributed in the hope that it will be useful,
   25    but WITHOUT ANY WARRANTY; without even the implied warranty of
   26    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   27    GNU General Public License for more details.
   28
   29    You should have received a copy of the GNU Lesser General Public
   30    License along with this library; if not, write to the Free Software
   31    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
   32
   33    As a special exception, if you link this library with other files,
   34    compiled with a Free Software compiler, to produce an executable, this
   35    library does not by itself cause the resulting executable to be covered
   36    by the GNU General Public License. This exception does not however
   37    invalidate any other reasons why the executable file might be covered by
   38    the GNU General Public License.
   39*/
   40
   41:- module(store_r,
   42	[
   43	    add_linear_11/3,
   44	    add_linear_f1/4,
   45	    add_linear_ff/5,
   46	    normalize_scalar/2,
   47	    delete_factor/4,
   48	    mult_linear_factor/3,
   49	    nf_rhs_x/4,
   50	    indep/2,
   51	    isolate/3,
   52	    nf_substitute/4,
   53	    mult_hom/3,
   54	    nf2sum/3,
   55	    nf_coeff_of/3,
   56	    renormalize/2	
   57	]).   58
   59% normalize_scalar(S,[N,Z])
   60%
   61% Transforms a scalar S into a linear expression [S,0]
   62
   63normalize_scalar(S,[S,0.0]).
   64
   65% renormalize(List,Lin)
   66%
   67% Renormalizes the not normalized linear expression in List into
   68% a normalized one. It does so to take care of unifications.
   69% (e.g. when a variable X is bound to a constant, the constant is added to
   70% the constant part of the linear expression; when a variable X is bound to
   71% another variable Y, the scalars of both are added)
   72
   73renormalize([I,R|Hom],Lin) :-
   74	length(Hom,Len),
   75	renormalize_log(Len,Hom,[],Lin0),
   76	add_linear_11([I,R],Lin0,Lin).
   77
   78% renormalize_log(Len,Hom,HomTail,Lin)
   79%
   80% Logarithmically renormalizes the homogene part of a not normalized
   81% linear expression. See also renormalize/2.
   82
   83renormalize_log(1,[Term|Xs],Xs,Lin) :-
   84	!,
   85	Term = l(X*_,_),
   86	renormalize_log_one(X,Term,Lin).
   87renormalize_log(2,[A,B|Xs],Xs,Lin) :-
   88	!,
   89	A = l(X*_,_),
   90	B = l(Y*_,_),
   91	renormalize_log_one(X,A,LinA),
   92	renormalize_log_one(Y,B,LinB),
   93	add_linear_11(LinA,LinB,Lin).
   94renormalize_log(N,L0,L2,Lin) :-
   95	P is N>>1,
   96	Q is N-P,
   97	renormalize_log(P,L0,L1,Lp),
   98	renormalize_log(Q,L1,L2,Lq),
   99	add_linear_11(Lp,Lq,Lin).
  100
  101% renormalize_log_one(X,Term,Res)
  102%
  103% Renormalizes a term in X: if X is a nonvar, the term becomes a scalar.
  104
  105renormalize_log_one(X,Term,Res) :-
  106	Term = l(X*K,_),
  107	(   var(X)
  108	->  get_attr(X,clpqr_itf,Att),
  109	    arg(5,Att,order(OrdX)), % Order might have changed
  110	    Res = [0.0,0.0,l(X*K,OrdX)]
  111	;   Xk is X*K,
  112	    normalize_scalar(Xk,Res)
  113	).
  114
  115% ----------------------------- sparse vector stuff ---------------------------- %
  116
  117% add_linear_ff(LinA,Ka,LinB,Kb,LinC)
  118%
  119% Linear expression LinC is the result of the addition of the 2 linear expressions
  120% LinA and LinB, each one multiplied by a scalar (Ka for LinA and Kb for LinB).
  121
  122add_linear_ff(LinA,Ka,LinB,Kb,LinC) :-
  123	LinA = [Ia,Ra|Ha],
  124	LinB = [Ib,Rb|Hb],
  125	LinC = [Ic,Rc|Hc],
  126	Ic is Ia*Ka+Ib*Kb,
  127	Rc is Ra*Ka+Rb*Kb,
  128 	add_linear_ffh(Ha,Ka,Hb,Kb,Hc).
  129
  130% add_linear_ffh(Ha,Ka,Hb,Kb,Hc)
  131%
  132% Homogene part Hc is the result of the addition of the 2 homogene parts Ha and Hb,
  133% each one multiplied by a scalar (Ka for Ha and Kb for Hb)
  134
  135add_linear_ffh([],_,Ys,Kb,Zs) :- mult_hom(Ys,Kb,Zs).
  136add_linear_ffh([l(X*Kx,OrdX)|Xs],Ka,Ys,Kb,Zs) :-
  137	add_linear_ffh(Ys,X,Kx,OrdX,Xs,Zs,Ka,Kb).
  138
  139% add_linear_ffh(Ys,X,Kx,OrdX,Xs,Zs,Ka,Kb)
  140%
  141% Homogene part Zs is the result of the addition of the 2 homogene parts Ys and
  142% [l(X*Kx,OrdX)|Xs], each one multiplied by a scalar (Ka for [l(X*Kx,OrdX)|Xs] and Kb for Ys)
  143
  144add_linear_ffh([],X,Kx,OrdX,Xs,Zs,Ka,_) :- mult_hom([l(X*Kx,OrdX)|Xs],Ka,Zs).
  145add_linear_ffh([l(Y*Ky,OrdY)|Ys],X,Kx,OrdX,Xs,Zs,Ka,Kb) :-
  146	compare(Rel,OrdX,OrdY),
  147	(   Rel = (=)
  148	->  Kz is Kx*Ka+Ky*Kb,
  149	    (   % Kz =:= 0
  150		Kz =< 1.0e-10,
  151		Kz >= -1.0e-10
  152	    ->  add_linear_ffh(Xs,Ka,Ys,Kb,Zs)
  153	    ;   Zs = [l(X*Kz,OrdX)|Ztail],
  154		add_linear_ffh(Xs,Ka,Ys,Kb,Ztail)
  155	    )
  156	;   Rel = (<)
  157	->  Zs = [l(X*Kz,OrdX)|Ztail],
  158	    Kz is Kx*Ka,
  159	    add_linear_ffh(Xs,Y,Ky,OrdY,Ys,Ztail,Kb,Ka)
  160	;   Rel = (>)
  161	->  Zs = [l(Y*Kz,OrdY)|Ztail],
  162	    Kz is Ky*Kb,
  163	    add_linear_ffh(Ys,X,Kx,OrdX,Xs,Ztail,Ka,Kb)
  164     	).
  165
  166% add_linear_f1(LinA,Ka,LinB,LinC)
  167%
  168% special case of add_linear_ff with Kb = 1
  169
  170add_linear_f1(LinA,Ka,LinB,LinC) :-
  171	LinA = [Ia,Ra|Ha],
  172	LinB = [Ib,Rb|Hb],
  173	LinC = [Ic,Rc|Hc],
  174	Ic is Ia*Ka+Ib,
  175	Rc is Ra*Ka+Rb,
  176	add_linear_f1h(Ha,Ka,Hb,Hc).
  177
  178% add_linear_f1h(Ha,Ka,Hb,Hc)
  179%
  180% special case of add_linear_ffh/5 with Kb = 1
  181
  182add_linear_f1h([],_,Ys,Ys).
  183add_linear_f1h([l(X*Kx,OrdX)|Xs],Ka,Ys,Zs) :-
  184	add_linear_f1h(Ys,X,Kx,OrdX,Xs,Zs,Ka).
  185
  186% add_linear_f1h(Ys,X,Kx,OrdX,Xs,Zs,Ka)
  187%
  188% special case of add_linear_ffh/8 with Kb = 1
  189
  190add_linear_f1h([],X,Kx,OrdX,Xs,Zs,Ka) :- mult_hom([l(X*Kx,OrdX)|Xs],Ka,Zs).
  191add_linear_f1h([l(Y*Ky,OrdY)|Ys],X,Kx,OrdX,Xs,Zs,Ka) :-
  192	compare(Rel,OrdX,OrdY),
  193	(   Rel = (=)
  194	->  Kz is Kx*Ka+Ky,
  195	    (   % Kz =:= 0.0
  196		Kz =< 1.0e-10,
  197		Kz >= -1.0e-10
  198	    ->  add_linear_f1h(Xs,Ka,Ys,Zs)
  199	    ;   Zs = [l(X*Kz,OrdX)|Ztail],
  200		add_linear_f1h(Xs,Ka,Ys,Ztail)
  201	    )
  202	;   Rel = (<)
  203	->  Zs = [l(X*Kz,OrdX)|Ztail],
  204	    Kz is Kx*Ka,
  205	    add_linear_f1h(Xs,Ka,[l(Y*Ky,OrdY)|Ys],Ztail)
  206 	;   Rel = (>)
  207	->  Zs = [l(Y*Ky,OrdY)|Ztail],
  208	    add_linear_f1h(Ys,X,Kx,OrdX,Xs,Ztail,Ka)
  209	).
  210
  211% add_linear_11(LinA,LinB,LinC)
  212%
  213% special case of add_linear_ff with Ka = 1 and Kb = 1
  214
  215add_linear_11(LinA,LinB,LinC) :-
  216	LinA = [Ia,Ra|Ha],
  217	LinB = [Ib,Rb|Hb],
  218	LinC = [Ic,Rc|Hc],
  219	Ic is Ia+Ib,
  220	Rc is Ra+Rb,
  221	add_linear_11h(Ha,Hb,Hc).
  222
  223% add_linear_11h(Ha,Hb,Hc)
  224%
  225% special case of add_linear_ffh/5 with Ka = 1 and Kb = 1
  226
  227add_linear_11h([],Ys,Ys).
  228add_linear_11h([l(X*Kx,OrdX)|Xs],Ys,Zs) :-
  229	add_linear_11h(Ys,X,Kx,OrdX,Xs,Zs).
  230
  231% add_linear_11h(Ys,X,Kx,OrdX,Xs,Zs)
  232%
  233% special case of add_linear_ffh/8 with Ka = 1 and Kb = 1
  234
  235add_linear_11h([],X,Kx,OrdX,Xs,[l(X*Kx,OrdX)|Xs]).
  236add_linear_11h([l(Y*Ky,OrdY)|Ys],X,Kx,OrdX,Xs,Zs) :-
  237	compare(Rel,OrdX,OrdY),
  238	(   Rel = (=)
  239	->  Kz is Kx+Ky,
  240	    (   % Kz =:= 0.0
  241		Kz =< 1.0e-10,
  242		Kz >= -1.0e-10
  243	    ->  add_linear_11h(Xs,Ys,Zs)
  244	    ;   Zs = [l(X*Kz,OrdX)|Ztail],
  245		add_linear_11h(Xs,Ys,Ztail)
  246	    )
  247	;   Rel = (<)
  248	->  Zs = [l(X*Kx,OrdX)|Ztail],
  249	    add_linear_11h(Xs,Y,Ky,OrdY,Ys,Ztail)
  250	;   Rel = (>)
  251	->  Zs = [l(Y*Ky,OrdY)|Ztail],
  252	    add_linear_11h(Ys,X,Kx,OrdX,Xs,Ztail)
  253	).
  254
  255% mult_linear_factor(Lin,K,Res)
  256%
  257% Linear expression Res is the result of multiplication of linear
  258% expression Lin by scalar K
  259
  260mult_linear_factor(Lin,K,Mult) :-
  261	TestK is K - 1.0,	% K =:= 1
  262	TestK =< 1.0e-10,
  263	TestK >= -1.0e-10,	% avoid copy
  264	!,
  265	Mult = Lin.
  266mult_linear_factor(Lin,K,Res) :-
  267	Lin = [I,R|Hom],
  268	Res = [Ik,Rk|Mult],
  269	Ik is I*K,
  270	Rk is R*K,
  271	mult_hom(Hom,K,Mult).
  272
  273% mult_hom(Hom,K,Res)
  274%
  275% Homogene part Res is the result of multiplication of homogene part
  276% Hom by scalar K
  277
  278mult_hom([],_,[]).
  279mult_hom([l(A*Fa,OrdA)|As],F,[l(A*Fan,OrdA)|Afs]) :-
  280	Fan is F*Fa,
  281	mult_hom(As,F,Afs).
  282
  283% nf_substitute(Ord,Def,Lin,Res)
  284%
  285% Linear expression Res is the result of substitution of Var in
  286% linear expression Lin, by its definition in the form of linear
  287% expression Def
  288
  289nf_substitute(OrdV,LinV,LinX,LinX1) :-
  290	delete_factor(OrdV,LinX,LinW,K),
  291	add_linear_f1(LinV,K,LinW,LinX1).
  292
  293% delete_factor(Ord,Lin,Res,Coeff)
  294%
  295% Linear expression Res is the result of the deletion of the term
  296% Var*Coeff where Var has ordering Ord from linear expression Lin
  297
  298delete_factor(OrdV,Lin,Res,Coeff) :-
  299	Lin = [I,R|Hom],
  300	Res = [I,R|Hdel],
  301	delete_factor_hom(OrdV,Hom,Hdel,Coeff).
  302
  303% delete_factor_hom(Ord,Hom,Res,Coeff)
  304%
  305% Homogene part Res is the result of the deletion of the term
  306% Var*Coeff from homogene part Hom
  307
  308delete_factor_hom(VOrd,[Car|Cdr],RCdr,RKoeff) :-
  309	Car = l(_*Koeff,Ord),
  310	compare(Rel,VOrd,Ord),
  311	(   Rel= (=)
  312	->  RCdr = Cdr,
  313	    RKoeff=Koeff
  314	;   Rel= (>)
  315	->  RCdr = [Car|RCdr1],
  316	    delete_factor_hom(VOrd,Cdr,RCdr1,RKoeff)
  317	).
  318
  319
  320% nf_coeff_of(Lin,OrdX,Coeff)
  321%
  322% Linear expression Lin contains the term l(X*Coeff,OrdX)
  323
  324nf_coeff_of([_,_|Hom],VOrd,Coeff) :-
  325	nf_coeff_hom(Hom,VOrd,Coeff).
  326
  327% nf_coeff_hom(Lin,OrdX,Coeff)
  328%
  329% Linear expression Lin contains the term l(X*Coeff,OrdX) where the
  330% order attribute of X = OrdX
  331
  332nf_coeff_hom([l(_*K,OVar)|Vs],OVid,Coeff) :-
  333	compare(Rel,OVid,OVar),
  334	(   Rel = (=)
  335	->  Coeff = K
  336	;   Rel = (>)
  337	->  nf_coeff_hom(Vs,OVid,Coeff)
  338	).
  339
  340% nf_rhs_x(Lin,OrdX,Rhs,K)
  341%
  342% Rhs = R + I where Lin = [I,R|Hom] and l(X*K,OrdX) is a term of Hom
  343
  344nf_rhs_x(Lin,OrdX,Rhs,K) :-
  345	Lin = [I,R|Tail],
  346	nf_coeff_hom(Tail,OrdX,K),
  347	Rhs is R+I.	% late because X may not occur in H
  348
  349% isolate(OrdN,Lin,Lin1)
  350%
  351% Linear expression Lin1 is the result of the transformation of linear expression
  352% Lin = 0 which contains the term l(New*K,OrdN) into an equivalent expression Lin1 = New.
  353
  354isolate(OrdN,Lin,Lin1) :-
  355	delete_factor(OrdN,Lin,Lin0,Coeff),
  356	K is -1.0/Coeff,
  357	mult_linear_factor(Lin0,K,Lin1).
  358
  359% indep(Lin,OrdX)
  360%
  361% succeeds if Lin = [0,_|[l(X*1,OrdX)]]
  362
  363indep(Lin,OrdX) :-
  364	Lin = [I,_|[l(_*K,OrdY)]],
  365	OrdX == OrdY,
  366	% K =:= 1.0
  367	TestK is K - 1.0,
  368	TestK =< 1.0e-10,
  369	TestK >= -1.0e-10,
  370	% I =:= 0
  371	I =< 1.0e-10,
  372	I >= -1.0e-10.
  373
  374% nf2sum(Lin,Sofar,Term)
  375%
  376% Transforms a linear expression into a sum
  377% (e.g. the expression [5,_,[l(X*2,OrdX),l(Y*-1,OrdY)]] gets transformed into 5 + 2*X - Y)
  378
  379nf2sum([],I,I).
  380nf2sum([X|Xs],I,Sum) :-
  381	(   % I =:= 0.0
  382	    I =< 1.0e-10,
  383	    I >= -1.0e-10
  384	->  X = l(Var*K,_),
  385 	    (   % K =:= 1.0
  386		TestK is K - 1.0,
  387		TestK =< 1.0e-10,
  388		TestK >= -1.0e-10
  389	    ->  hom2sum(Xs,Var,Sum)
  390	    ;   % K =:= -1.0
  391		TestK is K + 1.0,
  392		TestK =< 1.0e-10,
  393		TestK >= -1.0e-10
  394	    ->  hom2sum(Xs,-Var,Sum)
  395	    ;	hom2sum(Xs,K*Var,Sum)
  396	    )
  397	;   hom2sum([X|Xs],I,Sum)
  398 	).
  399
  400% hom2sum(Hom,Sofar,Term)
  401%
  402% Transforms a linear expression into a sum
  403% this predicate handles all but the first term
  404% (the first term does not need a concatenation symbol + or -)
  405% see also nf2sum/3
  406
  407hom2sum([],Term,Term).
  408hom2sum([l(Var*K,_)|Cs],Sofar,Term) :-
  409	(   % K =:= 1.0
  410	    TestK is K - 1.0,
  411	    TestK =< 1.0e-10,
  412	    TestK >= -1.0e-10
  413	->  Next = Sofar + Var
  414	;   % K =:= -1.0
  415	    TestK is K + 1.0,
  416	    TestK =< 1.0e-10,
  417	    TestK >= -1.0e-10
  418	->  Next = Sofar - Var
  419	;   % K < 0.0
  420	    K < -1.0e-10
  421	->  Ka is -K,
  422	    Next = Sofar - Ka*Var
  423	;   Next = Sofar + K*Var
  424	),
  425	hom2sum(Cs,Next,Term)