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	var(X),
  107	Term = l(X*K,_),
  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)].
  111renormalize_log_one(X,Term,Res) :-
  112	nonvar(X),
  113	Term = l(X*K,_),
  114	Xk is X*K,
  115	normalize_scalar(Xk,Res).
  116
  117% ----------------------------- sparse vector stuff ---------------------------- %
  118
  119% add_linear_ff(LinA,Ka,LinB,Kb,LinC)
  120%
  121% Linear expression LinC is the result of the addition of the 2 linear expressions
  122% LinA and LinB, each one multiplied by a scalar (Ka for LinA and Kb for LinB).
  123
  124add_linear_ff(LinA,Ka,LinB,Kb,LinC) :-
  125	LinA = [Ia,Ra|Ha],
  126	LinB = [Ib,Rb|Hb],
  127	LinC = [Ic,Rc|Hc],
  128	Ic is Ia*Ka+Ib*Kb,
  129	Rc is Ra*Ka+Rb*Kb,
  130 	add_linear_ffh(Ha,Ka,Hb,Kb,Hc).
  131
  132% add_linear_ffh(Ha,Ka,Hb,Kb,Hc)
  133%
  134% Homogene part Hc is the result of the addition of the 2 homogene parts Ha and Hb,
  135% each one multiplied by a scalar (Ka for Ha and Kb for Hb)
  136
  137add_linear_ffh([],_,Ys,Kb,Zs) :- mult_hom(Ys,Kb,Zs).
  138add_linear_ffh([l(X*Kx,OrdX)|Xs],Ka,Ys,Kb,Zs) :-
  139	add_linear_ffh(Ys,X,Kx,OrdX,Xs,Zs,Ka,Kb).
  140
  141% add_linear_ffh(Ys,X,Kx,OrdX,Xs,Zs,Ka,Kb)
  142%
  143% Homogene part Zs is the result of the addition of the 2 homogene parts Ys and
  144% [l(X*Kx,OrdX)|Xs], each one multiplied by a scalar (Ka for [l(X*Kx,OrdX)|Xs] and Kb for Ys)
  145
  146add_linear_ffh([],X,Kx,OrdX,Xs,Zs,Ka,_) :- mult_hom([l(X*Kx,OrdX)|Xs],Ka,Zs).
  147add_linear_ffh([l(Y*Ky,OrdY)|Ys],X,Kx,OrdX,Xs,Zs,Ka,Kb) :-
  148	compare(Rel,OrdX,OrdY),
  149	(   Rel = (=)
  150	->  Kz is Kx*Ka+Ky*Kb,
  151	    (   % Kz =:= 0
  152		Kz =< 1.0e-10,
  153		Kz >= -1.0e-10
  154	    ->  add_linear_ffh(Xs,Ka,Ys,Kb,Zs)
  155	    ;   Zs = [l(X*Kz,OrdX)|Ztail],
  156		add_linear_ffh(Xs,Ka,Ys,Kb,Ztail)
  157	    )
  158	;   Rel = (<)
  159	->  Zs = [l(X*Kz,OrdX)|Ztail],
  160	    Kz is Kx*Ka,
  161	    add_linear_ffh(Xs,Y,Ky,OrdY,Ys,Ztail,Kb,Ka)
  162	;   Rel = (>)
  163	->  Zs = [l(Y*Kz,OrdY)|Ztail],
  164	    Kz is Ky*Kb,
  165	    add_linear_ffh(Ys,X,Kx,OrdX,Xs,Ztail,Ka,Kb)
  166     	).
  167
  168% add_linear_f1(LinA,Ka,LinB,LinC)
  169%
  170% special case of add_linear_ff with Kb = 1
  171
  172add_linear_f1(LinA,Ka,LinB,LinC) :-
  173	LinA = [Ia,Ra|Ha],
  174	LinB = [Ib,Rb|Hb],
  175	LinC = [Ic,Rc|Hc],
  176	Ic is Ia*Ka+Ib,
  177	Rc is Ra*Ka+Rb,
  178	add_linear_f1h(Ha,Ka,Hb,Hc).
  179
  180% add_linear_f1h(Ha,Ka,Hb,Hc)
  181%
  182% special case of add_linear_ffh/5 with Kb = 1
  183
  184add_linear_f1h([],_,Ys,Ys).
  185add_linear_f1h([l(X*Kx,OrdX)|Xs],Ka,Ys,Zs) :-
  186	add_linear_f1h(Ys,X,Kx,OrdX,Xs,Zs,Ka).
  187
  188% add_linear_f1h(Ys,X,Kx,OrdX,Xs,Zs,Ka)
  189%
  190% special case of add_linear_ffh/8 with Kb = 1
  191
  192add_linear_f1h([],X,Kx,OrdX,Xs,Zs,Ka) :- mult_hom([l(X*Kx,OrdX)|Xs],Ka,Zs).
  193add_linear_f1h([l(Y*Ky,OrdY)|Ys],X,Kx,OrdX,Xs,Zs,Ka) :-
  194	compare(Rel,OrdX,OrdY),
  195	(   Rel = (=)
  196	->  Kz is Kx*Ka+Ky,
  197	    (   % Kz =:= 0.0
  198		Kz =< 1.0e-10,
  199		Kz >= -1.0e-10
  200	    ->  add_linear_f1h(Xs,Ka,Ys,Zs)
  201	    ;   Zs = [l(X*Kz,OrdX)|Ztail],
  202		add_linear_f1h(Xs,Ka,Ys,Ztail)
  203	    )
  204	;   Rel = (<)
  205	->  Zs = [l(X*Kz,OrdX)|Ztail],
  206	    Kz is Kx*Ka,
  207	    add_linear_f1h(Xs,Ka,[l(Y*Ky,OrdY)|Ys],Ztail)
  208 	;   Rel = (>)
  209	->  Zs = [l(Y*Ky,OrdY)|Ztail],
  210	    add_linear_f1h(Ys,X,Kx,OrdX,Xs,Ztail,Ka)
  211	).
  212
  213% add_linear_11(LinA,LinB,LinC)
  214%
  215% special case of add_linear_ff with Ka = 1 and Kb = 1
  216
  217add_linear_11(LinA,LinB,LinC) :-
  218	LinA = [Ia,Ra|Ha],
  219	LinB = [Ib,Rb|Hb],
  220	LinC = [Ic,Rc|Hc],
  221	Ic is Ia+Ib,
  222	Rc is Ra+Rb,
  223	add_linear_11h(Ha,Hb,Hc).
  224
  225% add_linear_11h(Ha,Hb,Hc)
  226%
  227% special case of add_linear_ffh/5 with Ka = 1 and Kb = 1
  228
  229add_linear_11h([],Ys,Ys).
  230add_linear_11h([l(X*Kx,OrdX)|Xs],Ys,Zs) :-
  231	add_linear_11h(Ys,X,Kx,OrdX,Xs,Zs).
  232
  233% add_linear_11h(Ys,X,Kx,OrdX,Xs,Zs)
  234%
  235% special case of add_linear_ffh/8 with Ka = 1 and Kb = 1
  236
  237add_linear_11h([],X,Kx,OrdX,Xs,[l(X*Kx,OrdX)|Xs]).
  238add_linear_11h([l(Y*Ky,OrdY)|Ys],X,Kx,OrdX,Xs,Zs) :-
  239	compare(Rel,OrdX,OrdY),
  240	(   Rel = (=)
  241	->  Kz is Kx+Ky,
  242	    (   % Kz =:= 0.0
  243		Kz =< 1.0e-10,
  244		Kz >= -1.0e-10
  245	    ->  add_linear_11h(Xs,Ys,Zs)
  246	    ;   Zs = [l(X*Kz,OrdX)|Ztail],
  247		add_linear_11h(Xs,Ys,Ztail)
  248	    )
  249	;   Rel = (<)
  250	->  Zs = [l(X*Kx,OrdX)|Ztail],
  251	    add_linear_11h(Xs,Y,Ky,OrdY,Ys,Ztail)
  252	;   Rel = (>)
  253	->  Zs = [l(Y*Ky,OrdY)|Ztail],
  254	    add_linear_11h(Ys,X,Kx,OrdX,Xs,Ztail)
  255	).
  256
  257% mult_linear_factor(Lin,K,Res)
  258%
  259% Linear expression Res is the result of multiplication of linear
  260% expression Lin by scalar K
  261
  262mult_linear_factor(Lin,K,Mult) :-
  263	TestK is K - 1.0,	% K =:= 1
  264	TestK =< 1.0e-10,
  265	TestK >= -1.0e-10,	% avoid copy
  266	!,
  267	Mult = Lin.
  268mult_linear_factor(Lin,K,Res) :-
  269	Lin = [I,R|Hom],
  270	Res = [Ik,Rk|Mult],
  271	Ik is I*K,
  272	Rk is R*K,
  273	mult_hom(Hom,K,Mult).
  274
  275% mult_hom(Hom,K,Res)
  276%
  277% Homogene part Res is the result of multiplication of homogene part
  278% Hom by scalar K
  279
  280mult_hom([],_,[]).
  281mult_hom([l(A*Fa,OrdA)|As],F,[l(A*Fan,OrdA)|Afs]) :-
  282	Fan is F*Fa,
  283	mult_hom(As,F,Afs).
  284
  285% nf_substitute(Ord,Def,Lin,Res)
  286%
  287% Linear expression Res is the result of substitution of Var in
  288% linear expression Lin, by its definition in the form of linear
  289% expression Def
  290
  291nf_substitute(OrdV,LinV,LinX,LinX1) :-
  292	delete_factor(OrdV,LinX,LinW,K),
  293	add_linear_f1(LinV,K,LinW,LinX1).
  294
  295% delete_factor(Ord,Lin,Res,Coeff)
  296%
  297% Linear expression Res is the result of the deletion of the term
  298% Var*Coeff where Var has ordering Ord from linear expression Lin
  299
  300delete_factor(OrdV,Lin,Res,Coeff) :-
  301	Lin = [I,R|Hom],
  302	Res = [I,R|Hdel],
  303	delete_factor_hom(OrdV,Hom,Hdel,Coeff).
  304
  305% delete_factor_hom(Ord,Hom,Res,Coeff)
  306%
  307% Homogene part Res is the result of the deletion of the term
  308% Var*Coeff from homogene part Hom
  309
  310delete_factor_hom(VOrd,[Car|Cdr],RCdr,RKoeff) :-
  311	Car = l(_*Koeff,Ord),
  312	compare(Rel,VOrd,Ord),
  313	(   Rel= (=)
  314	->  RCdr = Cdr,
  315	    RKoeff=Koeff
  316	;   Rel= (>)
  317	->  RCdr = [Car|RCdr1],
  318	    delete_factor_hom(VOrd,Cdr,RCdr1,RKoeff)
  319	).
  320
  321
  322% nf_coeff_of(Lin,OrdX,Coeff)
  323%
  324% Linear expression Lin contains the term l(X*Coeff,OrdX)
  325
  326nf_coeff_of([_,_|Hom],VOrd,Coeff) :-
  327	nf_coeff_hom(Hom,VOrd,Coeff).
  328
  329% nf_coeff_hom(Lin,OrdX,Coeff)
  330%
  331% Linear expression Lin contains the term l(X*Coeff,OrdX) where the
  332% order attribute of X = OrdX
  333
  334nf_coeff_hom([l(_*K,OVar)|Vs],OVid,Coeff) :-
  335	compare(Rel,OVid,OVar),
  336	(   Rel = (=)
  337	->  Coeff = K
  338	;   Rel = (>)
  339	->  nf_coeff_hom(Vs,OVid,Coeff)
  340	).
  341
  342% nf_rhs_x(Lin,OrdX,Rhs,K)
  343%
  344% Rhs = R + I where Lin = [I,R|Hom] and l(X*K,OrdX) is a term of Hom
  345
  346nf_rhs_x(Lin,OrdX,Rhs,K) :-
  347	Lin = [I,R|Tail],
  348	nf_coeff_hom(Tail,OrdX,K),
  349	Rhs is R+I.	% late because X may not occur in H
  350
  351% isolate(OrdN,Lin,Lin1)
  352%
  353% Linear expression Lin1 is the result of the transformation of linear expression
  354% Lin = 0 which contains the term l(New*K,OrdN) into an equivalent expression Lin1 = New.
  355
  356isolate(OrdN,Lin,Lin1) :-
  357	delete_factor(OrdN,Lin,Lin0,Coeff),
  358	K is -1.0/Coeff,
  359	mult_linear_factor(Lin0,K,Lin1).
  360
  361% indep(Lin,OrdX)
  362%
  363% succeeds if Lin = [0,_|[l(X*1,OrdX)]]
  364
  365indep(Lin,OrdX) :-
  366	Lin = [I,_|[l(_*K,OrdY)]],
  367	OrdX == OrdY,
  368	% K =:= 1.0
  369	TestK is K - 1.0,
  370	TestK =< 1.0e-10,
  371	TestK >= -1.0e-10,
  372	% I =:= 0
  373	I =< 1.0e-10,
  374	I >= -1.0e-10.
  375
  376% nf2sum(Lin,Sofar,Term)
  377%
  378% Transforms a linear expression into a sum
  379% (e.g. the expression [5,_,[l(X*2,OrdX),l(Y*-1,OrdY)]] gets transformed into 5 + 2*X - Y)
  380
  381nf2sum([],I,I).
  382nf2sum([X|Xs],I,Sum) :-
  383	(   % I =:= 0.0
  384	    I =< 1.0e-10,
  385	    I >= -1.0e-10
  386	->  X = l(Var*K,_),
  387 	    (   % K =:= 1.0
  388		TestK is K - 1.0,
  389		TestK =< 1.0e-10,
  390		TestK >= -1.0e-10
  391	    ->  hom2sum(Xs,Var,Sum)
  392	    ;   % K =:= -1.0
  393		TestK is K + 1.0,
  394		TestK =< 1.0e-10,
  395		TestK >= -1.0e-10
  396	    ->  hom2sum(Xs,-Var,Sum)
  397	    ;	hom2sum(Xs,K*Var,Sum)
  398	    )
  399	;   hom2sum([X|Xs],I,Sum)
  400 	).
  401
  402% hom2sum(Hom,Sofar,Term)
  403%
  404% Transforms a linear expression into a sum
  405% this predicate handles all but the first term
  406% (the first term does not need a concatenation symbol + or -)
  407% see also nf2sum/3
  408
  409hom2sum([],Term,Term).
  410hom2sum([l(Var*K,_)|Cs],Sofar,Term) :-
  411	(   % K =:= 1.0
  412	    TestK is K - 1.0,
  413	    TestK =< 1.0e-10,
  414	    TestK >= -1.0e-10
  415	->  Next = Sofar + Var
  416	;   % K =:= -1.0
  417	    TestK is K + 1.0,
  418	    TestK =< 1.0e-10,
  419	    TestK >= -1.0e-10
  420	->  Next = Sofar - Var
  421	;   % K < 0.0
  422	    K < -1.0e-10
  423	->  Ka is -K,
  424	    Next = Sofar - Ka*Var
  425	;   Next = Sofar + K*Var
  426	),
  427	hom2sum(Cs,Next,Term)