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(clpcd_store,
   42	[
   43	    add_linear_11/4,
   44	    add_linear_f1/5,
   45	    add_linear_ff/6,
   46	    normalize_scalar/2,
   47	    delete_factor/4,
   48	    mult_linear_factor/4,
   49	    nf_rhs_x/5,
   50	    isolate/4,
   51	    nf_substitute/5,
   52	    mult_hom/4,
   53	    nf_coeff_of/3,
   54	    renormalize/3	
   55	]).   56
   57:- use_module(library(clpcd/domain_ops)).   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]).
   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(CLP, [I,R|Hom], Lin) :-
   74	length(Hom,Len),
   75	renormalize_log(Len, CLP, Hom, [], Lin0),
   76	add_linear_11(CLP, [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, CLP, [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(CLP, LinA, LinB, Lin).
   94renormalize_log(N, CLP, L0, L2, Lin) :-
   95	P is N>>1,
   96	Q is N-P,
   97	renormalize_log(P, CLP, L0, L1, Lp),
   98	renormalize_log(Q, CLP, L1, L2, Lq),
   99	add_linear_11(CLP, 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,clpcd_itf,Att),
  109	arg(5,Att,order(OrdX)), % Order might have changed
  110	Res = [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(CLP, LinA, Ka, LinB, Kb, LinC) :-
  125	LinA = [Ia,Ra|Ha],
  126	LinB = [Ib,Rb|Hb],
  127	LinC = [Ic,Rc|Hc],
  128	eval_d(CLP, Ia*Ka+Ib*Kb, Ic),
  129	eval_d(CLP, Ra*Ka+Rb*Kb, Rc),
  130 	add_linear_ffh(Ha, CLP, 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([], CLP, _, Ys, Kb, Zs) :- mult_hom(Ys,CLP,Kb,Zs).
  138add_linear_ffh([l(X*Kx,OrdX)|Xs], CLP, Ka, Ys, Kb, Zs) :-
  139	add_linear_ffh(Ys, CLP, 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([], CLP, X, Kx, OrdX, Xs, Zs, Ka, _) :-
  147            mult_hom([l(X*Kx,OrdX)|Xs], CLP, Ka, Zs).
  148add_linear_ffh([l(Y*Ky,OrdY)|Ys], CLP, X, Kx, OrdX, Xs, Zs, Ka, Kb) :-
  149	compare(Rel,OrdX,OrdY),
  150	(   Rel = (=)
  151	->  eval_d(CLP, Kx*Ka+Ky*Kb, Kz),
  152	    (   % Kz =:= 0
  153                compare_d(CLP, =, Kx*Ka, -Ky*Kb)
  154	    ->  add_linear_ffh(Xs, CLP, Ka, Ys, Kb, Zs)
  155	    ;   Zs = [l(X*Kz,OrdX)|Ztail],
  156		add_linear_ffh(Xs, CLP, Ka, Ys, Kb, Ztail)
  157	    )
  158	;   Rel = (<)
  159	->  Zs = [l(X*Kz,OrdX)|Ztail],
  160	    eval_d(CLP, Kx*Ka, Kz),
  161	    add_linear_ffh(Xs, CLP, Y, Ky, OrdY, Ys, Ztail, Kb, Ka)
  162	;   Rel = (>)
  163	->  Zs = [l(Y*Kz,OrdY)|Ztail],
  164	    eval_d(CLP, Ky*Kb, Kz),
  165	    add_linear_ffh(Ys, CLP, 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(CLP, LinA, Ka, LinB, LinC) :-
  173	LinA = [Ia,Ra|Ha],
  174	LinB = [Ib,Rb|Hb],
  175	LinC = [Ic,Rc|Hc],
  176	eval_d(CLP, Ia*Ka+Ib, Ic),
  177	eval_d(CLP, Ra*Ka+Rb, Rc),
  178	add_linear_f1h(Ha, CLP, 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], CLP, Ka, Ys, Zs) :-
  186	add_linear_f1h(Ys, CLP, 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([], CLP, X, Kx, OrdX, Xs, Zs, Ka) :-
  193            mult_hom([l(X*Kx,OrdX)|Xs], CLP, Ka, Zs).
  194add_linear_f1h([l(Y*Ky,OrdY)|Ys], CLP, X, Kx, OrdX, Xs, Zs, Ka) :-
  195	compare(Rel,OrdX,OrdY),
  196	(   Rel = (=)
  197	->  eval_d(CLP, Kx*Ka+Ky, Kz),
  198	    (   % Kz =:= 0.0
  199                compare_d(CLP, =, Kx*Ka, -Ky)
  200	    ->  add_linear_f1h(Xs, CLP, Ka, Ys, Zs)
  201	    ;   Zs = [l(X*Kz,OrdX)|Ztail],
  202		add_linear_f1h(Xs, CLP, Ka, Ys, Ztail)
  203	    )
  204	;   Rel = (<)
  205	->  Zs = [l(X*Kz,OrdX)|Ztail],
  206	    eval_d(CLP, Kx*Ka, Kz),
  207	    add_linear_f1h(Xs, CLP, Ka, [l(Y*Ky,OrdY)|Ys], Ztail)
  208 	;   Rel = (>)
  209	->  Zs = [l(Y*Ky,OrdY)|Ztail],
  210	    add_linear_f1h(Ys, CLP, 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(CLP, LinA, LinB, LinC) :-
  218	LinA = [Ia,Ra|Ha],
  219	LinB = [Ib,Rb|Hb],
  220	LinC = [Ic,Rc|Hc],
  221	eval_d(CLP, Ia+Ib, Ic),
  222	eval_d(CLP, Ra+Rb, Rc),
  223	add_linear_11h(Ha, CLP, 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], CLP, Ys, Zs) :-
  231	add_linear_11h(Ys, CLP, 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], CLP, X, Kx, OrdX, Xs, Zs) :-
  239	compare(Rel,OrdX,OrdY),
  240	(   Rel = (=)
  241	->  eval_d(CLP, Kx+Ky, Kz),
  242	    (   % Kz =:= 0
  243		compare_d(CLP, =, Kx, -Ky)
  244	    ->  add_linear_11h(Xs, CLP, Ys, Zs)
  245	    ;   Zs = [l(X*Kz,OrdX)|Ztail],
  246		add_linear_11h(Xs, CLP, Ys, Ztail)
  247	    )
  248	;   Rel = (<)
  249	->  Zs = [l(X*Kx,OrdX)|Ztail],
  250	    add_linear_11h(Xs, CLP, Y, Ky, OrdY, Ys, Ztail)
  251	;   Rel = (>)
  252	->  Zs = [l(Y*Ky,OrdY)|Ztail],
  253	    add_linear_11h(Ys, CLP, X, Kx, OrdX, Xs, Ztail)
  254	).
  255
  256% mult_linear_factor(Lin,K,Res)
  257%
  258% Linear expression Res is the result of multiplication of linear
  259% expression Lin by scalar K
  260
  261mult_linear_factor(CLP, Lin, K, Mult) :-
  262        compare_d(CLP, =, K, 1),
  263	!,
  264	Mult = Lin.
  265mult_linear_factor(CLP, Lin, K, Res) :-
  266	Lin = [I,R|Hom],
  267	Res = [Ik, Rk|Mult],
  268	eval_d(CLP, I*K, Ik),
  269	eval_d(CLP, R*K, Rk),
  270	mult_hom(Hom, CLP, K, Mult).
  271
  272div_linear_factor(CLP, Lin, K, Mult) :-
  273	compare_d(CLP, =, K, 1),
  274	!,
  275	Mult = Lin.
  276div_linear_factor(CLP, Lin, K, Res) :-
  277	Lin = [I,R|Hom],
  278	Res = [Ik,Rk|Mult],
  279	eval_d(CLP, I/K, Ik),
  280	eval_d(CLP, R/K, Rk),
  281	div_hom(Hom, CLP, K, Mult).
  282
  283% mult_hom(Hom,CLP,K,Res)
  284%
  285% Homogene part Res is the result of multiplication of homogene part
  286% Hom by scalar K
  287
  288mult_hom([], _, _, []).
  289mult_hom([l(A*Fa,OrdA)|As], CLP, F, [l(A*Fan,OrdA)|Afs]) :-
  290	eval_d(CLP, F*Fa, Fan),
  291	mult_hom(As, CLP, F, Afs).
  292
  293div_hom([], _, _, []).
  294div_hom([l(A*Fa,OrdA)|As], CLP, F, [l(A*Fan,OrdA)|Afs]) :-
  295	eval_d(CLP, Fa/F, Fan),
  296	div_hom(As, CLP, F, Afs).
  297
  298% nf_substitute(Ord,Def,Lin,Res)
  299%
  300% Linear expression Res is the result of substitution of Var in
  301% linear expression Lin, by its definition in the form of linear
  302% expression Def
  303
  304nf_substitute(CLP, OrdV, LinV, LinX, LinX1) :-
  305	delete_factor(OrdV,LinX,LinW,K),
  306	add_linear_f1(CLP, LinV, K, LinW, LinX1).
  307
  308% delete_factor(Ord,Lin,Res,Coeff)
  309%
  310% Linear expression Res is the result of the deletion of the term
  311% Var*Coeff where Var has ordering Ord from linear expression Lin
  312
  313delete_factor(OrdV,Lin,Res,Coeff) :-
  314	Lin = [I,R|Hom],
  315	Res = [I,R|Hdel],
  316	delete_factor_hom(OrdV,Hom,Hdel,Coeff).
  317
  318% delete_factor_hom(Ord,Hom,Res,Coeff)
  319%
  320% Homogene part Res is the result of the deletion of the term
  321% Var*Coeff from homogene part Hom
  322
  323delete_factor_hom(VOrd,[Car|Cdr],RCdr,RKoeff) :-
  324	Car = l(_*Koeff,Ord),
  325	compare(Rel,VOrd,Ord),
  326	(   Rel= (=)
  327	->  RCdr = Cdr,
  328	    RKoeff=Koeff
  329	;   Rel= (>)
  330	->  RCdr = [Car|RCdr1],
  331	    delete_factor_hom(VOrd,Cdr,RCdr1,RKoeff)
  332	).
  333
  334
  335% nf_coeff_of(Lin,OrdX,Coeff)
  336%
  337% Linear expression Lin contains the term l(X*Coeff,OrdX)
  338
  339nf_coeff_of([_,_|Hom],VOrd,Coeff) :-
  340	nf_coeff_hom(Hom,VOrd,Coeff).
  341
  342% nf_coeff_hom(Lin,OrdX,Coeff)
  343%
  344% Linear expression Lin contains the term l(X*Coeff,OrdX) where the
  345% order attribute of X = OrdX
  346
  347nf_coeff_hom([l(_*K,OVar)|Vs],OVid,Coeff) :-
  348	compare(Rel,OVid,OVar),
  349	(   Rel = (=)
  350	->  Coeff = K
  351	;   Rel = (>)
  352	->  nf_coeff_hom(Vs,OVid,Coeff)
  353	).
  354
  355% nf_rhs_x(CLP, Lin,OrdX,Rhs,K)
  356%
  357% Rhs = R + I where Lin = [I,R|Hom] and l(X*K,OrdX) is a term of Hom
  358
  359nf_rhs_x(CLP, Lin,OrdX,Rhs,K) :-
  360	Lin = [I,R|Tail],
  361	nf_coeff_hom(Tail,OrdX,K),
  362	eval_d(CLP, R+I, Rhs).	% late because X may not occur in H
  363
  364% isolate(OrdN,Lin,Lin1)
  365%
  366% Linear expression Lin1 is the result of the transformation of linear expression
  367% Lin = 0 which contains the term l(New*K,OrdN) into an equivalent expression Lin1 = New.
  368
  369isolate(CLP, OrdN, Lin, Lin1) :-
  370	delete_factor(OrdN,Lin,Lin0,Coeff),
  371	div_linear_factor(CLP, Lin0, -Coeff, Lin1)