%nf,ineq,bv,store %;nf %:- use_module(library(error), [instantiation_error/1, type_error/2]). :- [class]. goal_expansion(geler(X,Y),geler(clpr,X,Y)). % ------------------------------------------------------------------------- % {Constraint} % % Adds the constraint Constraint to the constraint store. % % First rule is to prevent binding with other rules when a variable is input % Constraints are converted to normal form and if necessary, submitted to the linear % equality/inequality solver (bv + ineq) or to the non-linear store (geler) {Rel} :- var(Rel), !, instantiation_error(Rel). {R,Rs} :- !, {R},{Rs}. {R;Rs} :- !, ({R};{Rs}). % for entailment checking {L < R} :- !, nf(L-R,Nf), submit_lt(Nf). {L > R} :- !, nf(R-L,Nf), submit_lt(Nf). {L =< R} :- !, nf(L-R,Nf), submit_le( Nf). {<=(L,R)} :- !, nf(L-R,Nf), submit_le(Nf). {L >= R} :- !, nf(R-L,Nf), submit_le(Nf). {L =\= R} :- !, nf(L-R,Nf), submit_ne(Nf). {L =:= R} :- !, nf(L-R,Nf), submit_eq(Nf). {L = R} :- !, nf(L-R,Nf), submit_eq(Nf). {Rel} :- type_error(clpr_constraint, Rel). % entailed(C) % % s -> c = ~s v c = ~(s /\ ~c) % where s is the store and c is the constraint for which % we want to know whether it is entailed. % C is negated and added to the store. If this fails, then c is entailed by s entailed(C) :- negate(C,Cn), \+ {Cn}. % negate(C,Res). % % Res is the negation of constraint C % first rule is to prevent binding with other rules when a variable is input negate(Rel,_) :- var(Rel), !, throw(instantiation_error(entailed(Rel),1)). negate((A,B),(Na;Nb)) :- !, negate(A,Na), negate(B,Nb). negate((A;B),(Na,Nb)) :- !, negate(A,Na), negate(B,Nb). negate(A=B) :- !. negate(A>B,A=B) :- !. negate(A>=B,A A = 0 % b4) nonlinear -> geler % c) Nf=[A,B|Rest] % c1) A=k % c11) (B=c*X^+1 or B=c*X^-1), Rest=[] -> B=-k/c or B=-c/k % c12) invertible(A,B) % c13) linear(B|Rest) % c14) geler % c2) linear(Nf) % c3) nonlinear -> geler submit_eq([]). % trivial success: case a submit_eq([T|Ts]) :- submit_eq(Ts,T). submit_eq([],A) :- submit_eq_b(A). % case b submit_eq([B|Bs],A) :- submit_eq_c(A,B,Bs). % case c % submit_eq_b(A) % % Handles case b of submit_eq/1 % case b1: A is a constant (non-zero) submit_eq_b(v(_,[])) :- !, fail. % case b2/b3: A is n*X^P => X = 0 submit_eq_b(v(_,[X^P])) :- var(X), P > 0, !, export_binding(X,0.0). % case b2: non-linear is invertible: NL(X) = 0 => X - inv(NL)(0) = 0 submit_eq_b(v(_,[NL^1])) :- nonvar(NL), nl_invertible(NL,X,0.0,Inv), !, nf(-Inv,S), nf_add(X,S,New), submit_eq(New). % case b4: A is non-linear and not invertible => submit equality to geler submit_eq_b(Term) :- term_variables(Term,Vs), geler(Vs,nf_r:resubmit_eq([Term])). % submit_eq_c(A,B,Rest) % % Handles case c of submit_eq/1 % case c1: A is a constant submit_eq_c(v(I,[]),B,Rest) :- !, submit_eq_c1(Rest,B,I). % case c2: A,B and Rest are linear submit_eq_c(A,B,Rest) :- % c2 A = v(_,[X^1]), var(X), B = v(_,[Y^1]), var(Y), linear(Rest), !, Hom = [A,B|Rest], % 'solve_='(Hom). nf_length(Hom,0,Len), log_deref(Len,Hom,[],HomD), solve(HomD). % case c3: A, B or Rest is non-linear => geler submit_eq_c(A,B,Rest) :- Norm = [A,B|Rest], term_variables(Norm,Vs), geler(Vs,nf_r:resubmit_eq(Norm)). % submit_eq_c1(Rest,B,K) % % Handles case c1 of submit_eq/1 % case c11a: % i+kX^p=0 if p is an odd integer % special case: one solution if P is negativ but also for a negative X submit_eq_c1([], v(K,[X^P]), I) :- var(X), P =\= 0, 0 > (-I/K), integer(P) =:= P, 1 =:= integer(P) mod 2, !, Val is -((I/K) ** (1/P)), export_binding(X,Val). % case c11b: % i+kX^p=0 for p =\= 0, integer(P) =:= P % special case: generate 2 solutions if p is a positive even integer submit_eq_c1([], v(K,[X^P]), I) :- var(X), P =\= 0, 0 =< (-I/K), integer(P) =:= P, 0 =:= integer(P) mod 2, !, Val is (-I/K) ** (1/P), ( export_binding(X,Val) ; ValNeg is -Val, export_binding(X, ValNeg) ). % case c11c: % i+kX^p=0 for p =\= 0, 0 =< (-I/K) submit_eq_c1([], v(K,[X^P]), I) :- var(X), P =\= 0, 0 =< (-I/K), !, Val is (-I/K) ** (1/P), export_binding(X,Val). % case c11d: fail if var(X) and none of the above. submit_eq_c1([], v(_K,[X^_P]), _I) :- var(X), !, fail. % case c11e: fail for { -25 = _X^2.5 } and { -25 = _X^(-2.5) } and may be others! % if you uncomment this case { -25 = _X^2.5 } throw an error(evaluation_error(undefined)) % and { -25 = _X^(-2.5) } succeed with an unbound X submit_eq_c1([], v(K,[X^P]), I) :- nonvar(X), X = exp(_,_), % TLS added 03/12 1 =:= abs(P), 0 >= I, 0 >= K, !, fail. % case c12: non-linear, invertible: cNL(X)^1+k=0 => inv(NL)(-k/c) = 0 ; % cNL(X)^-1+k=0 => inv(NL)(-c/k) = 0 submit_eq_c1([],v(K,[NL^P]),I) :- nonvar(NL), ( P =:= 1, Y is -I/K ; P =:= -1, Y is -K/I ), nl_invertible(NL,X,Y,Inv), !, nf(-Inv,S), nf_add(X,S,New), submit_eq(New). % case c13: linear: X + Y + Z + c = 0 => submit_eq_c1(Rest,B,I) :- B = v(_,[Y^1]), var(Y), linear(Rest), !, % 'solve_='( [v(I,[]),B|Rest]). Hom = [B|Rest], nf_length(Hom,0,Len), normalize_scalar(I,Nonvar), log_deref(Len,Hom,[],HomD), add_linear_11(Nonvar,HomD,LinD), solve(LinD). % case c14: other cases => geler submit_eq_c1(Rest,B,I) :- Norm = [v(I,[]),B|Rest], term_variables(Norm,Vs), geler(Vs,nf_r:resubmit_eq(Norm)). % ----------------------------------------------------------------------- % submit_lt(Nf) % % Submits the inequality Nf<0 to the constraint store, where Nf is in normal form. % 0 < 0 => fail submit_lt([]) :- fail. % A + B < 0 submit_lt([A|As]) :- submit_lt(As,A). % submit_lt(As,A) % % Does what submit_lt/1 does where Nf = [A|As] % v(K,P) < 0 submit_lt([],v(K,P)) :- submit_lt_b(P,K). % A + B + Bs < 0 submit_lt([B|Bs],A) :- submit_lt_c(Bs,A,B). % submit_lt_b(P,K) % % Does what submit_lt/2 does where A = [v(K,P)] and As = [] % c < 0 submit_lt_b([],I) :- !, I < -1.0e-10. % cX^1 < 0 : if c < 0 then X > 0, else X < 0 submit_lt_b([X^1],K) :- var(X), !, ( K > 1.0e-10 -> ineq_one_s_p_0(X) % X is strictly negative ; ineq_one_s_n_0(X) % X is strictly positive ). % non-linear => geler submit_lt_b(P,K) :- term_variables(P,Vs), geler(Vs,nf_r:resubmit_lt([v(K,P)])). % submit_lt_c(Bs,A,B) % % Does what submit_lt/2 does where As = [B|Bs]. % c + kX < 0 => kX < c submit_lt_c([],A,B) :- A = v(I,[]), B = v(K,[Y^1]), var(Y), !, ineq_one(strict,Y,K,I). % linear < 0 => solve, non-linear < 0 => geler submit_lt_c(Rest,A,B) :- Norm = [A,B|Rest], ( linear(Norm) -> 'solve_<'(Norm) ; term_variables(Norm,Vs), geler(Vs,nf_r:resubmit_lt(Norm)) ). % submit_le(Nf) % % Submits the inequality Nf =< 0 to the constraint store, where Nf is in normal form. % See also submit_lt/1 % 0 =< 0 => success submit_le([]). % A + B =< 0 submit_le([A|As]) :- submit_le(As,A). % submit_le(As,A) % % See submit_lt/2. This handles less or equal. % v(K,P) =< 0 submit_le([],v(K,P)) :- submit_le_b(P,K). % A + B + Bs =< 0 submit_le([B|Bs],A) :- submit_le_c(Bs,A,B). % submit_le_b(P,K) % % See submit_lt_b/2. This handles less or equal. % c =< 0 submit_le_b([],I) :- !, I < 1.0e-10. % cX^1 =< 0: if c < 0 then X >= 0, else X =< 0 submit_le_b([X^1],K) :- var(X), !, ( K > 1.0e-10 -> ineq_one_n_p_0(X) % X is non-strictly negative ; ineq_one_n_n_0(X) % X is non-strictly positive ). % cX^P =< 0 => geler submit_le_b(P,K) :- term_variables(P,Vs), geler(Vs,nf_r:resubmit_le([v(K,P)])). % submit_le_c(Bs,A,B) % % See submit_lt_c/3. This handles less or equal. % c + kX^1 =< 0 => kX =< 0 submit_le_c([],A,B) :- A = v(I,[]), B = v(K,[Y^1]), var(Y), !, ineq_one(nonstrict,Y,K,I). % A, B & Rest are linear => solve, otherwise => geler submit_le_c(Rest,A,B) :- Norm = [A,B|Rest], ( linear(Norm) -> 'solve_=<'(Norm) ; term_variables(Norm,Vs), geler(Vs,nf_r:resubmit_le(Norm)) ). % submit_ne(Nf) % % Submits the inequality Nf =\= 0 to the constraint store, where Nf is in normal form. % if Nf is a constant => check constant = 0, else if Nf is linear => solve else => geler submit_ne(Norm1) :- ( nf_constant(Norm1,K) -> \+ (K >= -1.0e-10, K =< 1.0e-10) % K =\= 0 ; linear(Norm1) -> 'solve_=\\='(Norm1) ; term_variables(Norm1,Vs), geler(Vs,nf_r:resubmit_ne(Norm1)) ). % linear(A) % % succeeds when A is linear: all elements are of the form v(_,[]) or v(_,[X^1]) linear([]). linear(v(_,Ps)) :- linear_ps(Ps). linear([A|As]) :- linear(A), linear(As). % linear_ps(A) % % Succeeds when A = V^1 with V a variable. % This reflects the linearity of v(_,A). linear_ps([]). linear_ps([V^1]) :- var(V). % excludes sin(_), ... % % Goal delays until Term gets linear. % At this time, Var will be bound to the normalform of Term. % %:- meta_predicate wait_linear( ?, ?, :). % wait_linear(Term,Var,Goal) :- nf(Term,Nf), ( linear(Nf) -> Var = Nf, call(Goal) ; term_variables(Nf,Vars), geler(Vars,nf_r:wait_linear_retry(Nf,Var,Goal)) ). % % geler clients % resubmit_eq(N) :- repair(N,Norm), submit_eq(Norm). resubmit_lt(N) :- repair(N,Norm), submit_lt(Norm). resubmit_le(N) :- repair(N,Norm), submit_le(Norm). resubmit_ne(N) :- repair(N,Norm), submit_ne(Norm). wait_linear_retry(Nf0,Var,Goal) :- repair(Nf0,Nf), ( linear(Nf) -> Var = Nf, call(Goal) ; term_variables(Nf,Vars), geler(Vars,nf_r:wait_linear_retry(Nf,Var,Goal)) ). % ----------------------------------------------------------------------- % nl_invertible(F,X,Y,Res) % % Res is the evaluation of the inverse of nonlinear function F in variable X % where X is Y nl_invertible(sin(X),X,Y,Res) :- Res is asin(Y). nl_invertible(cos(X),X,Y,Res) :- Res is acos(Y). nl_invertible(tan(X),X,Y,Res) :- Res is atan(Y). nl_invertible(exp(B,C),X,A,Res) :- ( nf_constant(B,Kb) -> A > 1.0e-10, Kb > 1.0e-10, TestKb is Kb - 1.0, % Kb =\= 1.0 \+ (TestKb >= -1.0e-10, TestKb =< 1.0e-10), X = C, % note delayed unification Res is log(A)/log(Kb) ; nf_constant(C,Kc), \+ (A >= -1.0e-10, A =< 1.0e-10), % A =\= 0 Kc > 1.0e-10, % Kc > 0 X = B, % note delayed unification Res is A**(1.0/Kc) ). % ----------------------------------------------------------------------- % nf(Exp,Nf) % % Returns in Nf, the normal form of expression Exp % % v(A,[B^C,D^E|...]) means A*B^C*D^E*... where A is a scalar (number) % v(A,[]) means scalar A % variable X => 1*X^1 nf(X,Norm) :- var(X), !, Norm = [v(1.0,[X^1])]. nf(X,Norm) :- number(X), !, nf_number(X,Norm). % nf(#(Const),Norm) :- monash_constant(Const,Value), !, Norm = [v(Value,[])]. % nf(-A,Norm) :- !, nf(A,An), nf_mul_factor(v(-1.0,[]),An,Norm). nf(+A,Norm) :- !, nf(A,Norm). % nf(A+B,Norm) :- !, nf(A,An), nf(B,Bn), nf_add(An,Bn,Norm). nf(A-B,Norm) :- !, nf(A,An), nf(-B,Bn), nf_add(An,Bn,Norm). % nf(A*B,Norm) :- !, nf(A,An), nf(B,Bn), nf_mul(An,Bn,Norm). nf(A/B,Norm) :- !, nf(A,An), nf(B,Bn), nf_div(Bn,An,Norm). % non-linear function, one argument: Term = f(Arg) equals f'(Sa1) = Skel nf(Term,Norm) :- nonlin_1(Term,Arg,Skel,Sa1), !, nf(Arg,An), nf_nonlin_1(Skel,An,Sa1,Norm). % non-linear function, two arguments: Term = f(A1,A2) equals f'(Sa1,Sa2) = Skel nf(Term,Norm) :- nonlin_2(Term,A1,A2,Skel,Sa1,Sa2), !, nf(A1,A1n), nf(A2,A2n), nf_nonlin_2(Skel,A1n,A2n,Sa1,Sa2,Norm). % nf(Term,_) :- type_error(clpr_expression, Term). % nf_number(N,Res) % % If N is a number, N is normalized nf_number(N,Res) :- number(N), ( (N >= -1.0e-10, N =< 1.0e-10) % N =:= 0 -> Res = [] ; Res = [v(N,[])] ). nonlin_1(abs(X),X,abs(Y),Y). nonlin_1(sin(X),X,sin(Y),Y). nonlin_1(cos(X),X,cos(Y),Y). nonlin_1(tan(X),X,tan(Y),Y). nonlin_2(min(A,B),A,B,min(X,Y),X,Y). nonlin_2(max(A,B),A,B,max(X,Y),X,Y). nonlin_2(exp(A,B),A,B,exp(X,Y),X,Y). nonlin_2(pow(A,B),A,B,exp(X,Y),X,Y). % pow->exp nonlin_2(A^B,A,B,exp(X,Y),X,Y). nf_nonlin_1(Skel,An,S1,Norm) :- ( nf_constant(An,S1) -> nl_eval(Skel,Res), nf_number(Res,Norm) ; S1 = An, Norm = [v(1.0,[Skel^1])]). nf_nonlin_2(Skel,A1n,A2n,S1,S2,Norm) :- ( nf_constant(A1n,S1), nf_constant(A2n,S2) -> nl_eval(Skel,Res), nf_number(Res,Norm) ; Skel=exp(_,_), nf_constant(A2n,Exp), integerp(Exp,I) -> nf_power(I,A1n,Norm) ; S1 = A1n, S2 = A2n, Norm = [v(1.0,[Skel^1])] ). % evaluates non-linear functions in one variable where the variable is bound nl_eval(abs(X),R) :- R is abs(X). nl_eval(sin(X),R) :- R is sin(X). nl_eval(cos(X),R) :- R is cos(X). nl_eval(tan(X),R) :- R is tan(X). % evaluates non-linear functions in two variables where both variables are % bound nl_eval(min(X,Y),R) :- R is min(X,Y). nl_eval(max(X,Y),R) :- R is max(X,Y). nl_eval(exp(X,Y),R) :- R is X**Y. monash_constant(X,_) :- var(X), !, fail. monash_constant(p,3.14259265). monash_constant(pi,3.14259265). monash_constant(e,2.71828182). monash_constant(zero,1.0e-10). % % check if a Nf consists of just a constant % nf_constant([],0.0). nf_constant([v(K,[])],K). % split(NF,SNF,C) % % splits a normalform expression NF into two parts: % - a constant term C (which might be 0) % - the homogene part of the expression % % this method depends on the polynf ordering, i.e. [] < [X^1] ... split([],[],0.0). split([First|T],H,I) :- ( First = v(I,[]) -> H = T ; I = 0.0, H = [First|T] ). % nf_add(A,B,C): merges two normalized additions into a new normalized addition % % a normalized addition is one where the terms are ordered, e.g. X^1 < Y^1, X^1 < X^2 etc. % terms in the same variable with the same exponent are added, % e.g. when A contains v(5,[X^1]) and B contains v(4,[X^1]) then C contains v(9,[X^1]). nf_add([],Bs,Bs). nf_add([A|As],Bs,Cs) :- nf_add(Bs,A,As,Cs). nf_add([],A,As,Cs) :- Cs = [A|As]. nf_add([B|Bs],A,As,Cs) :- A = v(Ka,Pa), B = v(Kb,Pb), compare(Rel,Pa,Pb), nf_add_case(Rel,A,As,Cs,B,Bs,Ka,Kb,Pa). % nf_add_case(Rel,A,As,Cs,B,Bs,Ka,Kb,Pa) % % merges sorted lists [A|As] and [B|Bs] into new sorted list Cs % A = v(Ka,Pa) and B = v(Kb,_) % Rel is the ordering relation (<, > or =) between A and B. % when Rel is =, Ka and Kb are added to form a new scalar for Pa nf_add_case(<,A,As,Cs,B,Bs,_,_,_) :- Cs = [A|Rest], nf_add(As,B,Bs,Rest). nf_add_case(>,A,As,Cs,B,Bs,_,_,_) :- Cs = [B|Rest], nf_add(Bs,A,As,Rest). nf_add_case(=,_,As,Cs,_,Bs,Ka,Kb,Pa) :- Kc is Ka + Kb, ( (Kc >= -1.0e-10, Kc =< 1.0e-10) % Kc =:= 0.0 -> nf_add(As,Bs,Cs) ; Cs = [v(Kc,Pa)|Rest], nf_add(As,Bs,Rest) ). nf_mul(A,B,Res) :- nf_length(A,0,LenA), nf_length(B,0,LenB), nf_mul_log(LenA,A,[],LenB,B,Res). nf_mul_log(0,As,As,_,_,[]) :- !. nf_mul_log(1,[A|As],As,Lb,B,R) :- !, nf_mul_factor_log(Lb,B,[],A,R). nf_mul_log(2,[A1,A2|As],As,Lb,B,R) :- !, nf_mul_factor_log(Lb,B,[],A1,A1b), nf_mul_factor_log(Lb,B,[],A2,A2b), nf_add(A1b,A2b,R). nf_mul_log(N,A0,A2,Lb,B,R) :- P is N>>1, Q is N-P, nf_mul_log(P,A0,A1,Lb,B,Rp), nf_mul_log(Q,A1,A2,Lb,B,Rq), nf_add(Rp,Rq,R). % nf_add_2: does the same thing as nf_add, but only has 2 elements to combine. nf_add_2(Af,Bf,Res) :- % unfold: nf_add([Af],[Bf],Res). Af = v(Ka,Pa), Bf = v(Kb,Pb), compare(Rel,Pa,Pb), nf_add_2_case(Rel,Af,Bf,Res,Ka,Kb,Pa). nf_add_2_case(<,Af,Bf,[Af,Bf],_,_,_). nf_add_2_case(>,Af,Bf,[Bf,Af],_,_,_). nf_add_2_case(=,_, _,Res,Ka,Kb,Pa) :- Kc is Ka + Kb, ( (Kc >= -1.0e-10, Kc =< 1.0e-10) % Kc =:= 0 -> Res = [] ; Res = [v(Kc,Pa)] ). % nf_mul_k(A,B,C) % % C is the result of the multiplication of each element of A (of the form v(_,_)) with scalar B (which shouldn't be 0) nf_mul_k([],_,[]). nf_mul_k([v(I,P)|Vs],K,[v(Ki,P)|Vks]) :- Ki is K*I, nf_mul_k(Vs,K,Vks). % nf_mul_factor(A,Sum,Res) % % multiplies each element of the list Sum with factor A which is of the form v(_,_) % and puts the result in the sorted list Res. nf_mul_factor(v(K,[]),Sum,Res) :- !, nf_mul_k(Sum,K,Res). nf_mul_factor(F,Sum,Res) :- nf_length(Sum,0,Len), nf_mul_factor_log(Len,Sum,[],F,Res). % nf_mul_factor_log(Len,[Sum|SumTail],SumTail,F,Res) % % multiplies each element of Sum with F and puts the result in the sorted list Res % Len is the length of Sum % Sum is split logarithmically each step nf_mul_factor_log(0,As,As,_,[]) :- !. nf_mul_factor_log(1,[A|As],As,F,[R]) :- !, mult(A,F,R). nf_mul_factor_log(2,[A,B|As],As,F,Res) :- !, mult(A,F,Af), mult(B,F,Bf), nf_add_2(Af,Bf,Res). nf_mul_factor_log(N,A0,A2,F,R) :- P is N>>1, % P is rounded(N/2) Q is N-P, nf_mul_factor_log(P,A0,A1,F,Rp), nf_mul_factor_log(Q,A1,A2,F,Rq), nf_add(Rp,Rq,R). % mult(A,B,C) % % multiplies A and B into C each of the form v(_,_) mult(v(Ka,La),v(Kb,Lb),v(Kc,Lc)) :- Kc is Ka*Kb, pmerge(La,Lb,Lc). % pmerge(A,B,C) % % multiplies A and B into sorted C, where each is of the form of the second argument of v(_,_) pmerge([],Bs,Bs). pmerge([A|As],Bs,Cs) :- pmerge(Bs,A,As,Cs). pmerge([],A,As,Res) :- Res = [A|As]. pmerge([B|Bs],A,As,Res) :- A = Xa^Ka, B = Xb^Kb, compare(R,Xa,Xb), pmerge_case(R,A,As,Res,B,Bs,Ka,Kb,Xa). % pmerge_case(Rel,A,As,Res,B,Bs,Ka,Kb,Xa) % % multiplies and sorts [A|As] with [B|Bs] into Res where each is of the form of % the second argument of v(_,_) % % A is Xa^Ka and B is Xb^Kb, Rel is ordening relation between Xa and Xb pmerge_case(<,A,As,Res,B,Bs,_,_,_) :- Res = [A|Tail], pmerge(As,B,Bs,Tail). pmerge_case(>,A,As,Res,B,Bs,_,_,_) :- Res = [B|Tail], pmerge(Bs,A,As,Tail). pmerge_case(=,_,As,Res,_,Bs,Ka,Kb,Xa) :- Kc is Ka + Kb, ( Kc =:= 0 -> pmerge(As,Bs,Res) ; Res = [Xa^Kc|Tail], pmerge(As,Bs,Tail) ). % nf_div(Factor,In,Out) % % Out is the result of the division of each element in In (which is of the form v(_,_)) by Factor. % division by zero nf_div([],_,_) :- !, zero_division. % division by v(K,P) => multiplication by v(1/K,P^-1) nf_div([v(K,P)],Sum,Res) :- !, Ki is 1.0/K, mult_exp(P,-1,Pi), nf_mul_factor(v(Ki,Pi),Sum,Res). nf_div(D,A,[v(1.0,[(A/D)^1])]). % zero_division % % called when a division by zero is performed zero_division :- fail. % raise_exception(_) ? % mult_exp(In,Factor,Out) % % Out is the result of the multiplication of the exponents of the elements in In % (which are of the form X^Exp by Factor. mult_exp([],_,[]). mult_exp([X^P|Xs],K,[X^I|Tail]) :- I is K*P, mult_exp(Xs,K,Tail). % % raise to integer powers % % | ?- time({(1+X+Y+Z)^15=0}). (sicstus, try with SWI) % Timing 00:00:02.610 2.610 iterative % Timing 00:00:00.660 0.660 binomial nf_power(N,Sum,Norm) :- integer(N), compare(Rel,N,0), ( Rel = (<) -> Pn is -N, % nf_power_pos(Pn,Sum,Inorm), binom(Sum,Pn,Inorm), nf_div(Inorm,[v(1.0,[])],Norm) ; Rel = (>) -> % nf_power_pos(N,Sum,Norm) binom(Sum,N,Norm) ; Rel = (=) -> % 0^0 is indeterminate but we say 1 Norm = [v(1.0,[])] ). % % N>0 % % iterative method: X^N = X*(X^N-1) nf_power_pos(1,Sum,Norm) :- !, Sum = Norm. nf_power_pos(N,Sum,Norm) :- N1 is N-1, nf_power_pos(N1,Sum,Pn1), nf_mul(Sum,Pn1,Norm). % % N>0 % % binomial method binom(Sum,1,Power) :- !, Power = Sum. binom([],_,[]). binom([A|Bs],N,Power) :- ( Bs = [] -> nf_power_factor(A,N,Ap), Power = [Ap] ; Bs = [_|_] -> factor_powers(N,A,v(1.0,[]),Pas), sum_powers(N,Bs,[v(1.0,[])],Pbs,[]), combine_powers(Pas,Pbs,0,N,1,[],Power) ). combine_powers([],[],_,_,_,Pi,Pi). combine_powers([A|As],[B|Bs],L,R,C,Pi,Po) :- nf_mul(A,B,Ab), nf_mul_k(Ab,C,Abc), nf_add(Abc,Pi,Pii), L1 is L+1, R1 is R-1, C1 is C*R//L1, combine_powers(As,Bs,L1,R1,C1,Pii,Po). nf_power_factor(v(K,P),N,v(Kn,Pn)) :- Kn is K**N, mult_exp(P,N,Pn). factor_powers(0,_,Prev,[[Prev]]) :- !. factor_powers(N,F,Prev,[[Prev]|Ps]) :- N1 is N-1, mult(Prev,F,Next), factor_powers(N1,F,Next,Ps). sum_powers(0,_,Prev,[Prev|Lt],Lt) :- !. sum_powers(N,S,Prev,L0,Lt) :- N1 is N-1, nf_mul(S,Prev,Next), sum_powers(N1,S,Next,L0,[Prev|Lt]). % ------------------------------------------------------------------------------ repair(Sum,Norm) :- nf_length(Sum,0,Len), repair_log(Len,Sum,[],Norm). repair_log(0,As,As,[]) :- !. repair_log(1,[v(Ka,Pa)|As],As,R) :- !, repair_term(Ka,Pa,R). repair_log(2,[v(Ka,Pa),v(Kb,Pb)|As],As,R) :- !, repair_term(Ka,Pa,Ar), repair_term(Kb,Pb,Br), nf_add(Ar,Br,R). repair_log(N,A0,A2,R) :- P is N>>1, Q is N-P, repair_log(P,A0,A1,Rp), repair_log(Q,A1,A2,Rq), nf_add(Rp,Rq,R). repair_term(K,P,Norm) :- length(P,Len), repair_p_log(Len,P,[],Pr,[v(1.0,[])],Sum), nf_mul_factor(v(K,Pr),Sum,Norm). repair_p_log(0,Ps,Ps,[],L0,L0) :- !. repair_p_log(1,[X^P|Ps],Ps,R,L0,L1) :- !, repair_p(X,P,R,L0,L1). repair_p_log(2,[X^Px,Y^Py|Ps],Ps,R,L0,L2) :- !, repair_p(X,Px,Rx,L0,L1), repair_p(Y,Py,Ry,L1,L2), pmerge(Rx,Ry,R). repair_p_log(N,P0,P2,R,L0,L2) :- P is N>>1, Q is N-P, repair_p_log(P,P0,P1,Rp,L0,L1), repair_p_log(Q,P1,P2,Rq,L1,L2), pmerge(Rp,Rq,R). repair_p(Term,P,[Term^P],L0,L0) :- var(Term). repair_p(Term,P,[],L0,L1) :- nonvar(Term), repair_p_one(Term,TermN), nf_power(P,TermN,TermNP), nf_mul(TermNP,L0,L1). % % An undigested term a/b is distinguished from an % digested one by the fact that its arguments are % digested -> cuts after repair of args! % repair_p_one(Term,TermN) :- nf_number(Term,TermN), % freq. shortcut for nf/2 case below !. repair_p_one(A1/A2,TermN) :- repair(A1,A1n), repair(A2,A2n), !, nf_div(A2n,A1n,TermN). repair_p_one(Term,TermN) :- nonlin_1(Term,Arg,Skel,Sa), repair(Arg,An), !, nf_nonlin_1(Skel,An,Sa,TermN). repair_p_one(Term,TermN) :- nonlin_2(Term,A1,A2,Skel,Sa1,Sa2), repair(A1,A1n), repair(A2,A2n), !, nf_nonlin_2(Skel,A1n,A2n,Sa1,Sa2,TermN). repair_p_one(Term,TermN) :- nf(Term,TermN). nf_length([],Li,Li). nf_length([_|R],Li,Lo) :- Lii is Li+1, nf_length(R,Lii,Lo). % ------------------------------------------------------------------------------ % nf2term(NF,Term) % % transforms a normal form into a readable term % empty normal form = 0 nf2term([],0.0). % term is first element (+ next elements) nf2term([F|Fs],T) :- f02t(F,T0), % first element yfx(Fs,T0,T). % next elements yfx([],T0,T0). yfx([F|Fs],T0,TN) :- fn2t(F,Ft,Op), T1 =.. [Op,T0,Ft], yfx(Fs,T1,TN). % f02t(v(K,P),T) % % transforms the first element of the normal form (something of the form v(K,P)) % into a readable term f02t(v(K,P),T) :- ( % just a constant P = [] -> T = K ; TestK is K - 1.0, % K =:= 1 (TestK >= -1.0e-10, TestK =< 1.0e-10) -> p2term(P,T) ; TestK is K + 1.0, % K =:= -1 (TestK >= -1.0e-10, TestK =< 1.0e-10) -> T = -Pt, p2term(P,Pt) ; T = K*Pt, p2term(P,Pt) ). % f02t(v(K,P),T,Op) % % transforms a next element of the normal form (something of the form v(K,P)) % into a readable term fn2t(v(K,P),Term,Op) :- ( TestK is K - 1.0, % K =:= 1 (TestK >= -1.0e-10, TestK =< 1.0e-10) -> Term = Pt, Op = + ; TestK is K + 1.0, % K =:= -1 (TestK >= -1.0e-10, TestK =< 1.0e-10) -> Term = Pt, Op = - ; K < -1.0e-10 % K < 0 -> Kf is -K, Term = Kf*Pt, Op = - ; % K > 0 Term = K*Pt, Op = + ), p2term(P,Pt). % transforms the P part in v(_,P) into a readable term p2term([X^P|Xs],Term) :- ( Xs = [] -> pe2term(X,Xt), exp2term(P,Xt,Term) ; Xs = [_|_] -> Term = Xst*Xtp, pe2term(X,Xt), exp2term(P,Xt,Xtp), p2term(Xs,Xst) ). % exp2term(1,X,X) :- !. exp2term(-1,X,1.0/X) :- !. exp2term(P,X,Term) :- % Term = exp(X,Pn) Term = X^P. pe2term(X,Term) :- var(X), Term = X. pe2term(X,Term) :- nonvar(X), X =.. [F|Args], pe2term_args(Args,Argst), Term =.. [F|Argst]. pe2term_args([],[]). pe2term_args([A|As],[T|Ts]) :- nf2term(A,T), pe2term_args(As,Ts). % transg(Goal,[OutList|OutListTail],OutListTail) % % puts the equalities and inequalities that are implied by the elements in Goal % in the difference list OutList % % called by geler.pl for project.pl transg(resubmit_eq(Nf)) --> { nf2term([],Z), nf2term(Nf,Term) }, [clpr:{Term=Z}]. transg(resubmit_lt(Nf)) --> { nf2term([],Z), nf2term(Nf,Term) }, [clpr:{Term { nf2term([],Z), nf2term(Nf,Term) }, [clpr:{Term= { nf2term([],Z), nf2term(Nf,Term) }, [clpr:{Term=\=Z}]. transg(wait_linear_retry(Nf,Res,Goal)) --> { nf2term(Nf,Term) }, [clpr:{Term=Res},Goal]. integerp(X) :- floor(X)=:=X. integerp(X,I) :- floor(X)=:=X, I is integer(X). /******************************* * SANDBOX * *******************************/ %;ineq % ineq(H,I,Nf,Strictness) % % Solves the inequality Nf < 0 or Nf =< 0 where Nf is in normal form % and H and I are the homogene and inhomogene parts of Nf. ineq([],I,_,Strictness) :- ineq_ground(Strictness,I). ineq([v(K,[X^1])|Tail],I,Lin,Strictness) :- ineq_cases(Tail,I,Lin,Strictness,X,K). ineq_cases([],I,_,Strictness,X,K) :- % K*X + I < 0 or K*X + I =< 0 ineq_one(Strictness,X,K,I). ineq_cases([_|_],_,Lin,Strictness,_,_) :- deref(Lin,Lind), % Id+Hd =< 0 Lind = [Inhom,_|Hom], ineq_more(Hom,Inhom,Lind,Strictness). % ineq_ground(Strictness,I) % % Checks whether a grounded inequality I < 0 or I =< 0 is satisfied. ineq_ground(strict,I) :- I < -1.0e-10. % I < 0 ineq_ground(nonstrict,I) :- I < 1.0e-10. % I =< 0 % ineq_one(Strictness,X,K,I) % % Solves the inequality K*X + I < 0 or K*X + I =< 0 ineq_one(strict,X,K,I) :- ( K > 1.0e-10 % K > 0.0 -> ( I >= -1.0e-10, % I =:= 0.0 I =< 1.0e-10 -> ineq_one_s_p_0(X) % K*X < 0, K > 0 => X < 0 ; Inhom is I/K, ineq_one_s_p_i(X,Inhom) % K*X + I < 0, K > 0 => X + I/K < 0 ) ; ( I >= -1.0e-10, % I =:= 0.0 I =< 1.0e-10 -> ineq_one_s_n_0(X) % K*X < 0, K < 0 => -X < 0 ; Inhom is -I/K, ineq_one_s_n_i(X,Inhom) % K*X + I < 0, K < 0 => -X - I/K < 0 ) ). ineq_one(nonstrict,X,K,I) :- ( K > 1.0e-10 % K > 0.0 -> ( I >= -1.0e-10, % I =:= 0 I =< 1.0e-10 -> ineq_one_n_p_0(X) % K*X =< 0, K > 0 => X =< 0 ; Inhom is I/K, ineq_one_n_p_i(X,Inhom) % K*X + I =< 0, K > 0 => X + I/K =< 0 ) ; ( I >= -1.0e-10, % I =:= 0 I =< 1.0e-10 -> ineq_one_n_n_0(X) % K*X =< 0, K < 0 => -X =< 0 ; Inhom is -I/K, ineq_one_n_n_i(X,Inhom) % K*X + I =< 0, K < 0 => -X - I/K =< 0 ) ). % --------------------------- strict ---------------------------- % ineq_one_s_p_0(X) % % Solves the inequality X < 0 ineq_one_s_p_0(X) :- get_attr(X,itf,Att), arg(4,Att,lin([Ix,_|OrdX])), !, % old variable, this is deref ( \+ arg(1,Att,clpr) -> throw(error(permission_error('mix CLP(Q) variables with', 'CLP(R) variables:',X),context(_))) ; ineq_one_old_s_p_0(OrdX,X,Ix) ). ineq_one_s_p_0(X) :- % new variable, nothing depends on it var_intern(t_u(0.0),X,1). % put a strict inactive upperbound on the variable % ineq_one_s_n_0(X) % % Solves the inequality X > 0 ineq_one_s_n_0(X) :- get_attr(X,itf,Att), arg(4,Att,lin([Ix,_|OrdX])), !, ( \+ arg(1,Att,clpr) -> throw(error(permission_error('mix CLP(Q) variables with', 'CLP(R) variables:',X),context(_))) ; ineq_one_old_s_n_0(OrdX,X,Ix) ). ineq_one_s_n_0(X) :- var_intern(t_l(0.0),X,2). % puts a strict inactive lowerbound on the variable % ineq_one_s_p_i(X,I) % % Solves the inequality X < -I ineq_one_s_p_i(X,I) :- get_attr(X,itf,Att), arg(4,Att,lin([Ix,_|OrdX])), !, ( \+ arg(1,Att,clpr) -> throw(error(permission_error('mix CLP(Q) variables with', 'CLP(R) variables:',X),context(_))) ; ineq_one_old_s_p_i(OrdX,I,X,Ix) ). ineq_one_s_p_i(X,I) :- Bound is -I, var_intern(t_u(Bound),X,1). % puts a strict inactive upperbound on the variable % ineq_one_s_n_i(X,I) % % Solves the inequality X > I ineq_one_s_n_i(X,I) :- get_attr(X,itf,Att), arg(4,Att,lin([Ix,_|OrdX])), !, ( \+ arg(1,Att,clpr) -> throw(error(permission_error('mix CLP(Q) variables with', 'CLP(R) variables:',X),context(_))) ; ineq_one_old_s_n_i(OrdX,I,X,Ix) ). ineq_one_s_n_i(X,I) :- var_intern(t_l(I),X,2). % puts a strict inactive lowerbound on the variable % ineq_one_old_s_p_0(Hom,X,Inhom) % % Solves the inequality X < 0 where X has linear equation Hom + Inhom ineq_one_old_s_p_0([],_,Ix) :- Ix < -1.0e-10. % X = I: Ix < 0 ineq_one_old_s_p_0([l(Y*Ky,_)|Tail],X,Ix) :- ( Tail = [] % X = K*Y + I -> Bound is -Ix/Ky, update_indep(strict,Y,Ky,Bound) % X < 0, X = K*Y + I => Y < -I/K or Y > -I/K (depending on K) ; Tail = [_|_] -> get_attr(X,itf,Att), arg(2,Att,type(Type)), arg(3,Att,strictness(Old)), arg(4,Att,lin(Lin)), udus(Type,X,Lin,0.0,Old) % update strict upperbound ). % ineq_one_old_s_p_0(Hom,X,Inhom) % % Solves the inequality X > 0 where X has linear equation Hom + Inhom ineq_one_old_s_n_0([],_,Ix) :- Ix > 1.0e-10. % X = I: Ix > 0 ineq_one_old_s_n_0([l(Y*Ky,_)|Tail], X, Ix) :- ( Tail = [] % X = K*Y + I -> Coeff is -Ky, Bound is Ix/Coeff, update_indep(strict,Y,Coeff,Bound) ; Tail = [_|_] -> get_attr(X,itf,Att), arg(2,Att,type(Type)), arg(3,Att,strictness(Old)), arg(4,Att,lin(Lin)), udls(Type,X,Lin,0.0,Old) % update strict lowerbound ). % ineq_one_old_s_p_i(Hom,C,X,Inhom) % % Solves the inequality X + C < 0 where X has linear equation Hom + Inhom ineq_one_old_s_p_i([],I,_,Ix) :- Ix + I < -1.0e-10. % X = I ineq_one_old_s_p_i([l(Y*Ky,_)|Tail],I,X,Ix) :- ( Tail = [] % X = K*Y + I -> Bound is -(Ix + I)/Ky, update_indep(strict,Y,Ky,Bound) ; Tail = [_|_] -> Bound is -I, get_attr(X,itf,Att), arg(2,Att,type(Type)), arg(3,Att,strictness(Old)), arg(4,Att,lin(Lin)), udus(Type,X,Lin,Bound,Old) % update strict upperbound ). % ineq_one_old_s_n_i(Hom,C,X,Inhom) % % Solves the inequality X - C > 0 where X has linear equation Hom + Inhom ineq_one_old_s_n_i([],I,_,Ix) :- -Ix + I < -1.0e-10. % X = I ineq_one_old_s_n_i([l(Y*Ky,_)|Tail],I,X,Ix) :- ( Tail = [] % X = K*Y + I -> Coeff is -Ky, Bound is (Ix - I)/Coeff, update_indep(strict,Y,Coeff,Bound) ; Tail = [_|_] -> get_attr(X,itf,Att), arg(2,Att,type(Type)), arg(3,Att,strictness(Old)), arg(4,Att,lin(Lin)), udls(Type,X,Lin,I,Old) % update strict lowerbound ). % -------------------------- nonstrict -------------------------- % ineq_one_n_p_0(X) % % Solves the inequality X =< 0 ineq_one_n_p_0(X) :- get_attr(X,itf,Att), arg(4,Att,lin([Ix,_|OrdX])), !, % old variable, this is deref ( \+ arg(1,Att,clpr) -> throw(error(permission_error('mix CLP(Q) variables with', 'CLP(R) variables:',X),context(_))) ; ineq_one_old_n_p_0(OrdX,X,Ix) ). ineq_one_n_p_0(X) :- % new variable, nothing depends on it var_intern(t_u(0.0),X,0). % nonstrict upperbound % ineq_one_n_n_0(X) % % Solves the inequality X >= 0 ineq_one_n_n_0(X) :- get_attr(X,itf,Att), arg(4,Att,lin([Ix,_|OrdX])), !, ( \+ arg(1,Att,clpr) -> throw(error(permission_error('mix CLP(Q) variables with', 'CLP(R) variables:',X),context(_))) ; ineq_one_old_n_n_0(OrdX,X,Ix) ). ineq_one_n_n_0(X) :- var_intern(t_l(0.0),X,0). % nonstrict lowerbound % ineq_one_n_p_i(X,I) % % Solves the inequality X =< -I ineq_one_n_p_i(X,I) :- get_attr(X,itf,Att), arg(4,Att,lin([Ix,_|OrdX])), !, ( \+ arg(1,Att,clpr) -> throw(error(permission_error('mix CLP(Q) variables with', 'CLP(R) variables:',X),context(_))) ; ineq_one_old_n_p_i(OrdX,I,X,Ix) ). ineq_one_n_p_i(X,I) :- Bound is -I, var_intern(t_u(Bound),X,0). % nonstrict upperbound % ineq_one_n_n_i(X,I) % % Solves the inequality X >= I ineq_one_n_n_i(X,I) :- get_attr(X,itf,Att), arg(4,Att,lin([Ix,_|OrdX])), !, ( \+ arg(1,Att,clpr) -> throw(error(permission_error('mix CLP(Q) variables with', 'CLP(R) variables:',X),context(_))) ; ineq_one_old_n_n_i(OrdX,I,X,Ix) ). ineq_one_n_n_i(X,I) :- var_intern(t_l(I),X,0). % nonstrict lowerbound % ineq_one_old_n_p_0(Hom,X,Inhom) % % Solves the inequality X =< 0 where X has linear equation Hom + Inhom ineq_one_old_n_p_0([],_,Ix) :- Ix < 1.0e-10. % X =I ineq_one_old_n_p_0([l(Y*Ky,_)|Tail],X,Ix) :- ( Tail = [] % X = K*Y + I -> Bound is -Ix/Ky, update_indep(nonstrict,Y,Ky,Bound) ; Tail = [_|_] -> get_attr(X,itf,Att), arg(2,Att,type(Type)), arg(3,Att,strictness(Old)), arg(4,Att,lin(Lin)), udu(Type,X,Lin,0.0,Old) % update nonstrict upperbound ). % ineq_one_old_n_n_0(Hom,X,Inhom) % % Solves the inequality X >= 0 where X has linear equation Hom + Inhom ineq_one_old_n_n_0([],_,Ix) :- Ix > -1.0e-10. % X = I ineq_one_old_n_n_0([l(Y*Ky,_)|Tail], X, Ix) :- ( Tail = [] % X = K*Y + I -> Coeff is -Ky, Bound is Ix/Coeff, update_indep(nonstrict,Y,Coeff,Bound) ; Tail = [_|_] -> get_attr(X,itf,Att), arg(2,Att,type(Type)), arg(3,Att,strictness(Old)), arg(4,Att,lin(Lin)), udl(Type,X,Lin,0.0,Old) % update nonstrict lowerbound ). % ineq_one_old_n_p_i(Hom,C,X,Inhom) % % Solves the inequality X + C =< 0 where X has linear equation Hom + Inhom ineq_one_old_n_p_i([],I,_,Ix) :- Ix + I < 1.0e-10. % X = I ineq_one_old_n_p_i([l(Y*Ky,_)|Tail],I,X,Ix) :- ( Tail = [] % X = K*Y + I -> Bound is -(Ix + I)/Ky, update_indep(nonstrict,Y,Ky,Bound) ; Tail = [_|_] -> Bound is -I, get_attr(X,itf,Att), arg(2,Att,type(Type)), arg(3,Att,strictness(Old)), arg(4,Att,lin(Lin)), udu(Type,X,Lin,Bound,Old) % update nonstrict upperbound ). % ineq_one_old_n_n_i(Hom,C,X,Inhom) % % Solves the inequality X - C >= 0 where X has linear equation Hom + Inhom ineq_one_old_n_n_i([],I,_,Ix) :- -Ix + I < 1.0e-10. % X = I ineq_one_old_n_n_i([l(Y*Ky,_)|Tail],I,X,Ix) :- ( Tail = [] -> Coeff is -Ky, Bound is (Ix - I)/Coeff, update_indep(nonstrict,Y,Coeff,Bound) ; Tail = [_|_] -> get_attr(X,itf,Att), arg(2,Att,type(Type)), arg(3,Att,strictness(Old)), arg(4,Att,lin(Lin)), udl(Type,X,Lin,I,Old) ). % --------------------------------------------------------------- % ineq_more(Hom,Inhom,Lin,Strictness) % % Solves the inequality Lin < 0 or Lin =< 0 with Lin = Hom + Inhom ineq_more([],I,_,Strictness) :- ineq_ground(Strictness,I). % I < 0 or I =< 0 ineq_more([l(X*K,_)|Tail],Id,Lind,Strictness) :- ( Tail = [] -> % X*K < Id or X*K =< Id % one var: update bound instead of slack introduction get_or_add_class(X,_), % makes sure X belongs to a class Bound is -Id/K, update_indep(Strictness,X,K,Bound) % new bound ; Tail = [_|_] -> ineq_more(Strictness,Lind) ). % ineq_more(Strictness,Lin) % % Solves the inequality Lin < 0 or Lin =< 0 ineq_more(strict,Lind) :- ( unconstrained(Lind,U,K,Rest) -> % never fails, no implied value % Lind < 0 => Rest < -K*U where U has no bounds var_intern(t_l(0.0),S,2), % create slack variable S get_attr(S,itf,AttS), arg(5,AttS,order(OrdS)), Ki is -1.0/K, add_linear_ff(Rest,Ki,[0.0,0.0,l(S*1.0,OrdS)],Ki,LinU), % U = (-1/K)*Rest + (-1/K)*S LinU = [_,_|Hu], get_or_add_class(U,Class), same_class(Hu,Class), % put all variables of new lin. eq. of U in the same class get_attr(U,itf,AttU), arg(5,AttU,order(OrdU)), arg(6,AttU,class(ClassU)), backsubst(ClassU,OrdU,LinU) % substitute U by new lin. eq. everywhere in the class ; var_with_def_intern(t_u(0.0),S,Lind,1), % Lind < 0 => Lind = S with S < 0 basis_add(S,_), % adds S to the basis determine_active_dec(Lind), % activate bounds reconsider(S) % reconsider basis ). ineq_more(nonstrict,Lind) :- ( unconstrained(Lind,U,K,Rest) -> % never fails, no implied value % Lind =< 0 => Rest =< -K*U where U has no bounds var_intern(t_l(0.0),S,0), % create slack variable S Ki is -1.0/K, get_attr(S,itf,AttS), arg(5,AttS,order(OrdS)), add_linear_ff(Rest,Ki,[0.0,0.0,l(S*1.0,OrdS)],Ki,LinU), % U = (-1K)*Rest + (-1/K)*S LinU = [_,_|Hu], get_or_add_class(U,Class), same_class(Hu,Class), % put all variables of new lin. eq of U in the same class get_attr(U,itf,AttU), arg(5,AttU,order(OrdU)), arg(6,AttU,class(ClassU)), backsubst(ClassU,OrdU,LinU) % substitute U by new lin. eq. everywhere in the class ; % all variables are constrained var_with_def_intern(t_u(0.0),S,Lind,0), % Lind =< 0 => Lind = S with S =< 0 basis_add(S,_), % adds S to the basis determine_active_dec(Lind), reconsider(S) ). % update_indep(Strictness,X,K,Bound) % % Updates the bound of independent variable X where X < Bound or X =< Bound % or X > Bound or X >= Bound, depending on Strictness and K. update_indep(strict,X,K,Bound) :- get_attr(X,itf,Att), arg(2,Att,type(Type)), arg(3,Att,strictness(Old)), arg(4,Att,lin(Lin)), ( K < -1.0e-10 -> uils(Type,X,Lin,Bound,Old) % update independent lowerbound strict ; uius(Type,X,Lin,Bound,Old) % update independent upperbound strict ). update_indep(nonstrict,X,K,Bound) :- get_attr(X,itf,Att), arg(2,Att,type(Type)), arg(3,Att,strictness(Old)), arg(4,Att,lin(Lin)), ( K < -1.0e-10 -> uil(Type,X,Lin,Bound,Old) % update independent lowerbound nonstrict ; uiu(Type,X,Lin,Bound,Old) % update independent upperbound nonstrict ). % --------------------------------------------------------------------------------------- % % Update a bound on a var xi % % a) independent variable % % a1) update inactive bound: done % % a2) update active bound: % Determine [lu]b including most constraining row R % If we are within: done % else pivot(R,xi) and introduce bound via (b) % % a3) introduce a bound on an unconstrained var: % All vars that depend on xi are unconstrained (invariant) -> % the bound cannot invalidate any Lhs % % b) dependent variable % % repair upper or lower (maybe just swap with an unconstrained var from Rhs) % % % Sign = 1,0,-1 means inside,at,outside % % Read following predicates as update (dependent/independent) (lowerbound/upperbound) (strict) % udl(Type,X,Lin,Bound,Strict) % % Updates lower bound of dependent variable X with linear equation % Lin that had type Type and strictness Strict, to the new non-strict % bound Bound. udl(t_none,X,Lin,Bound,_Sold) :- get_attr(X,itf,AttX), arg(5,AttX,order(Ord)), setarg(2,AttX,type(t_l(Bound))), setarg(3,AttX,strictness(0)), ( unconstrained(Lin,Uc,Kuc,Rest) -> % X = Lin => -1/K*Rest + 1/K*X = U where U has no bounds Ki is -1.0/Kuc, add_linear_ff(Rest,Ki,[0.0,0.0,l(X* -1.0,Ord)],Ki,LinU), get_attr(Uc,itf,AttU), arg(5,AttU,order(OrdU)), arg(6,AttU,class(Class)), backsubst(Class,OrdU,LinU) ; % no unconstrained variables in Lin: make X part of basis and reconsider basis_add(X,_), determine_active_inc(Lin), reconsider(X) ). udl(t_l(L),X,Lin,Bound,Sold) :- TestBL is Bound - L, ( TestBL < -1.0e-10 -> true % new bound is smaller than old one: keep old ; TestBL > 1.0e-10 -> % new bound is larger than old one: use new and reconsider basis Strict is Sold /\ 1, get_attr(X,itf,Att), setarg(2,Att,type(t_l(Bound))), setarg(3,Att,strictness(Strict)), reconsider_lower(X,Lin,Bound) % makes sure that Lin still satisfies lowerbound Bound ; true % new bound is equal to old one, new one is nonstrict: keep old ). udl(t_u(U),X,Lin,Bound,_Sold) :- TestUB is U - Bound, ( TestUB < -1.0e-10 -> fail % new bound is larger than upperbound: fail ; TestUB > 1.0e-10 -> % new bound is smaller than upperbound: add new and reconsider basis get_attr(X,itf,Att), setarg(2,Att,type(t_lu(Bound,U))), reconsider_lower(X,Lin,Bound) % makes sure that Lin still satisfies lowerbound Bound ; solve_bound(Lin,Bound) % new bound is equal to upperbound: solve ). udl(t_lu(L,U),X,Lin,Bound,Sold) :- TestBL is Bound - L, ( TestBL < -1.0e-10 -> true % smaller than lowerbound: keep ; TestBL > 1.0e-10 -> % larger than lowerbound: check upperbound TestUB is U - Bound, ( TestUB < -1.0e-10 -> fail % larger than upperbound: fail ; TestUB > 1.0e-10 -> % smaller than upperbound: use new and reconsider basis Strict is Sold /\ 1, get_attr(X,itf,Att), setarg(2,Att,type(t_lu(Bound,U))), setarg(3,Att,strictness(Strict)), reconsider_lower(X,Lin,Bound) ; % equal to upperbound: if strictness matches => solve Sold /\ 1 =:= 0, solve_bound(Lin,Bound) ) ; true % equal to lowerbound and nonstrict: keep ). % udls(Type,X,Lin,Bound,Strict) % % Updates lower bound of dependent variable X with linear equation % Lin that had type Type and strictness Strict, to the new strict % bound Bound. udls(t_none,X,Lin,Bound,_Sold) :- get_attr(X,itf,AttX), arg(5,AttX,order(Ord)), setarg(2,AttX,type(t_l(Bound))), setarg(3,AttX,strictness(2)), ( unconstrained(Lin,Uc,Kuc,Rest) -> % X = Lin => U = -1/K*Rest + 1/K*X with U an unconstrained variable Ki is -1.0/Kuc, add_linear_ff(Rest,Ki,[0.0,0.0,l(X* -1.0,Ord)],Ki,LinU), get_attr(Uc,itf,AttU), arg(5,AttU,order(OrdU)), arg(6,AttU,class(Class)), backsubst(Class,OrdU,LinU) ; % no unconstrained variables: add X to basis and reconsider basis basis_add(X,_), determine_active_inc(Lin), reconsider(X) ). udls(t_l(L),X,Lin,Bound,Sold) :- TestBL is Bound - L, ( TestBL < -1.0e-10 -> true % smaller than lowerbound: keep ; TestBL > 1.0e-10 -> % larger than lowerbound: use new and reconsider basis Strict is Sold \/ 2, get_attr(X,itf,Att), setarg(2,Att,type(t_l(Bound))), setarg(3,Att,strictness(Strict)), reconsider_lower(X,Lin,Bound) ; % equal to lowerbound: check strictness Strict is Sold \/ 2, get_attr(X,itf,Att), setarg(3,Att,strictness(Strict)) ). udls(t_u(U),X,Lin,Bound,Sold) :- U - Bound > 1.0e-10, % smaller than upperbound: set new bound Strict is Sold \/ 2, get_attr(X,itf,Att), setarg(2,Att,type(t_lu(Bound,U))), setarg(3,Att,strictness(Strict)), reconsider_lower(X,Lin,Bound). udls(t_lu(L,U),X,Lin,Bound,Sold) :- TestBL is Bound - L, ( TestBL < -1.0e-10 -> true % smaller than lowerbound: keep ; TestBL > 1.0e-10 -> % larger than lowerbound: check upperbound and possibly use new and reconsider basis U - Bound > 1.0e-10, Strict is Sold \/ 2, get_attr(X,itf,Att), setarg(2,Att,type(t_lu(Bound,U))), setarg(3,Att,strictness(Strict)), reconsider_lower(X,Lin,Bound) ; % equal to lowerbound: put new strictness Strict is Sold \/ 2, get_attr(X,itf,Att), setarg(3,Att,strictness(Strict)) ). % udu(Type,X,Lin,Bound,Strict) % % Updates upper bound of dependent variable X with linear equation % Lin that had type Type and strictness Strict, to the new non-strict % bound Bound. udu(t_none,X,Lin,Bound,_Sold) :- get_attr(X,itf,AttX), arg(5,AttX,order(Ord)), setarg(2,AttX,type(t_u(Bound))), setarg(3,AttX,strictness(0)), ( unconstrained(Lin,Uc,Kuc,Rest) -> % X = Lin => U = -1/K*Rest + 1/K*X with U an unconstrained variable Ki is -1.0/Kuc, add_linear_ff(Rest,Ki,[0.0,0.0,l(X* -1.0,Ord)],Ki,LinU), get_attr(Uc,itf,AttU), arg(5,AttU,order(OrdU)), arg(6,AttU,class(Class)), backsubst(Class,OrdU,LinU) ; % no unconstrained variables: add X to basis and reconsider basis basis_add(X,_), determine_active_dec(Lin), % try to lower R reconsider(X) ). udu(t_u(U),X,Lin,Bound,Sold) :- TestUB is U - Bound, ( TestUB < -1.0e-10 -> true % larger than upperbound: keep ; TestUB > 1.0e-10 -> % smaller than upperbound: update and reconsider basis Strict is Sold /\ 2, get_attr(X,itf,Att), setarg(2,Att,type(t_u(Bound))), setarg(3,Att,strictness(Strict)), reconsider_upper(X,Lin,Bound) ; true % equal to upperbound and nonstrict: keep ). udu(t_l(L),X,Lin,Bound,_Sold) :- TestBL is Bound - L, ( TestBL < -1.0e-10 -> fail % smaller than lowerbound: fail ; TestBL > 1.0e-10 -> % larger than lowerbound: use new and reconsider basis get_attr(X,itf,Att), setarg(2,Att,type(t_lu(L,Bound))), reconsider_upper(X,Lin,Bound) ; solve_bound(Lin,Bound) % equal to lowerbound: solve ). udu(t_lu(L,U),X,Lin,Bound,Sold) :- TestUB is U - Bound, ( TestUB < -1.0e-10 -> true % larger than upperbound: keep ; TestUB > 1.0e-10 -> % smaller than upperbound: check lowerbound TestBL is Bound - L, ( TestBL < -1.0e-10 -> fail % smaller than lowerbound: fail ; TestBL > 1.0e-10 -> % larger than lowerbound: update and reconsider basis Strict is Sold /\ 2, get_attr(X,itf,Att), setarg(2,Att,type(t_lu(L,Bound))), setarg(3,Att,strictness(Strict)), reconsider_upper(X,Lin,Bound) ; % equal to lowerbound: check strictness and possibly solve Sold /\ 2 =:= 0, solve_bound(Lin,Bound) ) ; true % equal to upperbound and nonstrict: keep ). % udus(Type,X,Lin,Bound,Strict) % % Updates upper bound of dependent variable X with linear equation % Lin that had type Type and strictness Strict, to the new strict % bound Bound. udus(t_none,X,Lin,Bound,_Sold) :- get_attr(X,itf,AttX), arg(5,AttX,order(Ord)), setarg(2,AttX,type(t_u(Bound))), setarg(3,AttX,strictness(1)), ( unconstrained(Lin,Uc,Kuc,Rest) -> % X = Lin => U = -1/K*Rest + 1/K*X with U an unconstrained variable Ki is -1.0/Kuc, add_linear_ff(Rest,Ki,[0.0,0.0,l(X* -1.0,Ord)],Ki,LinU), get_attr(Uc,itf,AttU), arg(5,AttU,order(OrdU)), arg(6,AttU,class(Class)), backsubst(Class,OrdU,LinU) ; % no unconstrained variables: add X to basis and reconsider basis basis_add(X,_), determine_active_dec(Lin), reconsider(X) ). udus(t_u(U),X,Lin,Bound,Sold) :- TestUB is U - Bound, ( TestUB < -1.0e-10 -> true % larger than upperbound: keep ; TestUB > 1.0e-10 -> % smaller than upperbound: update bound and reconsider basis Strict is Sold \/ 1, get_attr(X,itf,Att), setarg(2,Att,type(t_u(Bound))), setarg(3,Att,strictness(Strict)), reconsider_upper(X,Lin,Bound) ; % equal to upperbound: set new strictness Strict is Sold \/ 1, get_attr(X,itf,Att), setarg(3,Att,strictness(Strict)) ). udus(t_l(L),X,Lin,Bound,Sold) :- Bound - L > 1.0e-10, % larger than lowerbound: update and reconsider basis Strict is Sold \/ 1, get_attr(X,itf,Att), setarg(2,Att,type(t_lu(L,Bound))), setarg(3,Att,strictness(Strict)), reconsider_upper(X,Lin,Bound). udus(t_lu(L,U),X,Lin,Bound,Sold) :- TestUB is U - Bound, ( TestUB < -1.0e-10 -> true % larger than upperbound: keep ; TestUB > 1.0e-10 -> % smaller than upperbound: check lowerbound, possibly update and reconsider basis Bound - L > 1.0e-10, Strict is Sold \/ 1, get_attr(X,itf,Att), setarg(2,Att,type(t_lu(L,Bound))), setarg(3,Att,strictness(Strict)), reconsider_upper(X,Lin,Bound) ; % equal to upperbound: update strictness Strict is Sold \/ 1, get_attr(X,itf,Att), setarg(3,Att,strictness(Strict)) ). % uiu(Type,X,Lin,Bound,Strict) % % Updates upper bound of independent variable X with linear equation % Lin that had type Type and strictness Strict, to the new non-strict % bound Bound. uiu(t_none,X,_Lin,Bound,_) :- % X had no bounds get_attr(X,itf,Att), setarg(2,Att,type(t_u(Bound))), setarg(3,Att,strictness(0)). uiu(t_u(U),X,_Lin,Bound,Sold) :- TestUB is U - Bound, ( TestUB < -1.0e-10 -> true % larger than upperbound: keep ; TestUB > 1.0e-10 -> % smaller than upperbound: update. Strict is Sold /\ 2, % update strictness: strictness of lowerbound is kept, % strictness of upperbound is set to non-strict get_attr(X,itf,Att), setarg(2,Att,type(t_u(Bound))), setarg(3,Att,strictness(Strict)) ; true % equal to upperbound and nonstrict: keep ). uiu(t_l(L),X,Lin,Bound,_Sold) :- TestBL is Bound - L, ( TestBL < -1.0e-10 -> fail % Lowerbound was smaller than new upperbound: fail ; TestBL > 1.0e-10 -> % Upperbound is larger than lowerbound: store new bound get_attr(X,itf,Att), setarg(2,Att,type(t_lu(L,Bound))) ; solve_bound(Lin,Bound) % Lowerbound was equal to new upperbound: solve ). uiu(t_L(L),X,Lin,Bound,_Sold) :- TestBL is Bound - L, ( TestBL < -1.0e-10 -> fail % Same as for t_l ; TestBL > 1.0e-10 -> % Same as for t_l (new bound becomes t_Lu) get_attr(X,itf,Att), setarg(2,Att,type(t_Lu(L,Bound))) ; solve_bound(Lin,Bound) % Same as for t_l ). uiu(t_lu(L,U),X,Lin,Bound,Sold) :- TestUB is U - Bound, ( TestUB < -1.0e-10 -> true % Upperbound was smaller than new bound: keep ; TestUB > 1.0e-10 -> TestBL is Bound - L, % Upperbound was larger than new bound: check lowerbound ( TestBL < -1.0e-10 -> fail % Lowerbound was larger than new bound: fail ; TestBL > 1.0e-10 -> % Lowerbound was smaller than new bound: store new bound Strict is Sold /\ 2, get_attr(X,itf,Att), setarg(2,Att,type(t_lu(L,Bound))), setarg(3,Att,strictness(Strict)) ; % Lowerbound was equal to new bound: solve Sold /\ 2 =:= 0, % Only solve when strictness matches solve_bound(Lin,Bound) ) ; true % Upperbound was equal to new bound and new bound non-strict: keep ). uiu(t_Lu(L,U),X,Lin,Bound,Sold) :- % See t_lu case TestUB is U - Bound, ( TestUB < -1.0e-10 -> true ; TestUB > 1.0e-10 -> TestBL is Bound - L, ( TestBL < -1.0e-10 -> fail ; TestBL > 1.0e-10 -> Strict is Sold /\ 2, get_attr(X,itf,Att), setarg(2,Att,type(t_Lu(L,Bound))), setarg(3,Att,strictness(Strict)) ; Sold /\ 2 =:= 0, solve_bound(Lin,Bound) ) ; true ). uiu(t_U(U),X,_Lin,Bound,Sold) :- TestUB is U - Bound, ( TestUB < -1.0e-10 -> true % larger than upperbound: keep ; TestUB > 1.0e-10 -> % smaller than active upperbound: check how much active upperbound can be lowered. % if enough, just lower bound, otherwise update the bound, make X dependent and reconsider basis Strict is Sold /\ 2, ( get_attr(X,itf,Att), arg(5,Att,order(OrdX)), arg(6,Att,class(ClassX)), lb(ClassX,OrdX,Vlb-Vb-Lb), Bound - (Lb + U) < 1.0e-10 -> get_attr(X,itf,Att2), % changed? setarg(2,Att2,type(t_U(Bound))), setarg(3,Att2,strictness(Strict)), pivot_a(Vlb,X,Vb,t_u(Bound)), reconsider(X) ; get_attr(X,itf,Att), arg(5,Att,order(OrdX)), arg(6,Att,class(ClassX)), setarg(2,Att,type(t_U(Bound))), setarg(3,Att,strictness(Strict)), Delta is Bound - U, backsubst_delta(ClassX,OrdX,X,Delta) ) ; true % equal to upperbound and non-strict: keep ). uiu(t_lU(L,U),X,Lin,Bound,Sold) :- TestUB is U - Bound, ( TestUB < -1.0e-10 -> true % larger than upperbound: keep ; TestUB > 1.0e-10 -> TestBL is Bound-L, ( TestBL < -1.0e-10 -> fail % smaller than lowerbound: fail ; TestBL > 1.0e-10 -> % larger than lowerbound: see t_U case for rest Strict is Sold /\ 2, ( get_attr(X,itf,Att), arg(5,Att,order(OrdX)), arg(6,Att,class(ClassX)), lb(ClassX,OrdX,Vlb-Vb-Lb), Bound - (Lb + U) < 1.0e-10 -> get_attr(X,itf,Att2), % changed? setarg(2,Att2,type(t_lU(L,Bound))), setarg(3,Att2,strictness(Strict)), pivot_a(Vlb,X,Vb,t_lu(L,Bound)), reconsider(X) ; get_attr(X,itf,Att), arg(5,Att,order(OrdX)), arg(6,Att,class(ClassX)), setarg(2,Att,type(t_lU(L,Bound))), setarg(3,Att,strictness(Strict)), Delta is Bound - U, backsubst_delta(ClassX,OrdX,X,Delta) ) ; % equal to lowerbound: check strictness and solve Sold /\ 2 =:= 0, solve_bound(Lin,Bound) ) ; true % equal to upperbound and non-strict: keep % smaller than upperbound: check lowerbound ). % uius(Type,X,Lin,Bound,Strict) % % Updates upper bound of independent variable X with linear equation % Lin that had type Type and strictness Strict, to the new strict % bound Bound. (see also uiu/5) uius(t_none,X,_Lin,Bound,_Sold) :- get_attr(X,itf,Att), setarg(2,Att,type(t_u(Bound))), setarg(3,Att,strictness(1)). uius(t_u(U),X,_Lin,Bound,Sold) :- TestUB is U - Bound, ( TestUB < -1.0e-10 -> true ; TestUB > 1.0e-10 -> Strict is Sold \/ 1, get_attr(X,itf,Att), setarg(2,Att,type(t_u(Bound))), setarg(3,Att,strictness(Strict)) ; Strict is Sold \/ 1, get_attr(X,itf,Att), setarg(3,Att,strictness(Strict)) ). uius(t_l(L),X,_Lin,Bound,Sold) :- Bound - L > 1.0e-10, Strict is Sold \/ 1, get_attr(X,itf,Att), setarg(2,Att,type(t_lu(L,Bound))), setarg(3,Att,strictness(Strict)). uius(t_L(L),X,_Lin,Bound,Sold) :- Bound - L > 1.0e-10, Strict is Sold \/ 1, get_attr(X,itf,Att), setarg(2,Att,type(t_Lu(L,Bound))), setarg(3,Att,strictness(Strict)). uius(t_lu(L,U),X,_Lin,Bound,Sold) :- TestUB is U - Bound, ( TestUB < -1.0e-10 -> true ; TestUB > 1.0e-10 -> Bound - L > 1.0e-10, Strict is Sold \/ 1, get_attr(X,itf,Att), setarg(2,Att,type(t_lu(L,Bound))), setarg(3,Att,strictness(Strict)) ; Strict is Sold \/ 1, get_attr(X,itf,Att), setarg(3,Att,strictness(Strict)) ). uius(t_Lu(L,U),X,_Lin,Bound,Sold) :- TestUB is U - Bound, ( TestUB < -1.0e-10 -> true ; TestUB > 1.0e-10 -> Bound - L > 1.0e-10, Strict is Sold \/ 1, get_attr(X,itf,Att), setarg(2,Att,type(t_Lu(L,Bound))), setarg(3,Att,strictness(Strict)) ; Strict is Sold \/ 1, get_attr(X,itf,Att), setarg(3,Att,strictness(Strict)) ). uius(t_U(U),X,_Lin,Bound,Sold) :- TestUB is U - Bound, ( TestUB < -1.0e-10 -> true ; TestUB > 1.0e-10 -> Strict is Sold \/ 1, ( get_attr(X,itf,Att), arg(5,Att,order(OrdX)), arg(6,Att,class(ClassX)), lb(ClassX,OrdX,Vlb-Vb-Lb), Bound - (Lb + U) < 1.0e-10 -> get_attr(X,itf,Att2), % changed? setarg(2,Att2,type(t_U(Bound))), setarg(3,Att2,strictness(Strict)), pivot_a(Vlb,X,Vb,t_u(Bound)), reconsider(X) ; get_attr(X,itf,Att), arg(5,Att,order(OrdX)), arg(6,Att,class(ClassX)), setarg(2,Att,type(t_U(Bound))), setarg(3,Att,strictness(Strict)), Delta is Bound - U, backsubst_delta(ClassX,OrdX,X,Delta) ) ; Strict is Sold \/ 1, get_attr(X,itf,Att), setarg(3,Att,strictness(Strict)) ). uius(t_lU(L,U),X,_Lin,Bound,Sold) :- TestUB is U - Bound, ( TestUB < -1.0e-10 -> true ; TestUB > 1.0e-10 -> Bound - L > 1.0e-10, Strict is Sold \/ 1, ( get_attr(X,itf,Att), arg(5,Att,order(OrdX)), arg(6,Att,class(ClassX)), lb(ClassX,OrdX,Vlb-Vb-Lb), Bound - (Lb + U) < 1.0e-10 -> get_attr(X,itf,Att2), % changed? setarg(2,Att2,type(t_lU(L,Bound))), setarg(3,Att2,strictness(Strict)), pivot_a(Vlb,X,Vb,t_lu(L,Bound)), reconsider(X) ; get_attr(X,itf,Att), arg(5,Att,order(OrdX)), arg(6,Att,class(ClassX)), setarg(2,Att,type(t_lU(L,Bound))), setarg(3,Att,strictness(Strict)), Delta is Bound - U, backsubst_delta(ClassX,OrdX,X,Delta) ) ; Strict is Sold \/ 1, get_attr(X,itf,Att), setarg(3,Att,strictness(Strict)) ). % uil(Type,X,Lin,Bound,Strict) % % Updates lower bound of independent variable X with linear equation % Lin that had type Type and strictness Strict, to the new non-strict % bound Bound. (see also uiu/5) uil(t_none,X,_Lin,Bound,_Sold) :- get_attr(X,itf,Att), setarg(2,Att,type(t_l(Bound))), setarg(3,Att,strictness(0)). uil(t_l(L),X,_Lin,Bound,Sold) :- TestBL is Bound - L, ( TestBL < -1.0e-10 -> true ; TestBL > 1.0e-10 -> Strict is Sold /\ 1, get_attr(X,itf,Att), setarg(2,Att,type(t_l(Bound))), setarg(3,Att,strictness(Strict)) ; true ). uil(t_u(U),X,Lin,Bound,_Sold) :- TestUB is U - Bound, ( TestUB < -1.0e-10 -> fail ; TestUB > 1.0e-10 -> get_attr(X,itf,Att), setarg(2,Att,type(t_lu(Bound,U))) ; solve_bound(Lin,Bound) ). uil(t_U(U),X,Lin,Bound,_Sold) :- TestUB is U - Bound, ( TestUB < -1.0e-10 -> fail ; TestUB > 1.0e-10 -> get_attr(X,itf,Att), setarg(2,Att,type(t_lU(Bound,U))) ; solve_bound(Lin,Bound) ). uil(t_lu(L,U),X,Lin,Bound,Sold) :- TestBL is Bound - L, ( TestBL < -1.0e-10 -> true ; TestBL > 1.0e-10 -> TestUB is U - Bound, ( TestUB < -1.0e-10 -> fail ; TestUB > 1.0e-10 -> Strict is Sold /\ 1, get_attr(X,itf,Att), setarg(2,Att,type(t_lu(Bound,U))), setarg(3,Att,strictness(Strict)) ; Sold /\ 1 =:= 0, solve_bound(Lin,Bound) ) ; true ). uil(t_lU(L,U),X,Lin,Bound,Sold) :- TestBL is Bound - L, ( TestBL < -1.0e-10 -> true ; TestBL > 1.0e-10 -> TestUB is U - Bound, ( TestUB < -1.0e-10 -> fail ; TestUB > 1.0e-10 -> Strict is Sold /\ 1, get_attr(X,itf,Att), setarg(2,Att,type(t_lU(Bound,U))), setarg(3,Att,strictness(Strict)) ; Sold /\ 1 =:= 0, solve_bound(Lin,Bound) ) ; true ). uil(t_L(L),X,_Lin,Bound,Sold) :- TestBL is Bound - L, ( TestBL < -1.0e-10 -> true ; TestBL > 1.0e-10 -> Strict is Sold /\ 1, ( get_attr(X,itf,Att), arg(5,Att,order(OrdX)), arg(6,Att,class(ClassX)), ub(ClassX,OrdX,Vub-Vb-Ub), Bound - (Ub + L) > -1.0e-10 -> get_attr(X,itf,Att2), % changed? setarg(2,Att2,type(t_L(Bound))), setarg(3,Att2,strictness(Strict)), pivot_a(Vub,X,Vb,t_l(Bound)), reconsider(X) ; get_attr(X,itf,Att), arg(5,Att,order(OrdX)), arg(6,Att,class(ClassX)), setarg(2,Att,type(t_L(Bound))), setarg(3,Att,strictness(Strict)), Delta is Bound - L, backsubst_delta(ClassX,OrdX,X,Delta) ) ; true ). uil(t_Lu(L,U),X,Lin,Bound,Sold) :- TestBL is Bound - L, ( TestBL < -1.0e-10 -> true ; TestBL > 1.0e-10 -> TestUB is U - Bound, ( TestUB < -1.0e-10 -> fail ; TestUB > 1.0e-10 -> Strict is Sold /\ 1, ( get_attr(X,itf,Att), arg(5,Att,order(OrdX)), arg(6,Att,class(ClassX)), ub(ClassX,OrdX,Vub-Vb-Ub), Bound - (Ub + L) > -1.0e-10 -> get_attr(X,itf,Att2), % changed? setarg(2,Att2,t_Lu(Bound,U)), setarg(3,Att2,strictness(Strict)), pivot_a(Vub,X,Vb,t_lu(Bound,U)), reconsider(X) ; get_attr(X,itf,Att), arg(5,Att,order(OrdX)), arg(6,Att,class(ClassX)), setarg(2,Att,type(t_Lu(Bound,U))), setarg(3,Att,strictness(Strict)), Delta is Bound - L, backsubst_delta(ClassX,OrdX,X,Delta) ) ; Sold /\ 1 =:= 0, solve_bound(Lin,Bound) ) ; true ). % uils(Type,X,Lin,Bound,Strict) % % Updates lower bound of independent variable X with linear equation % Lin that had type Type and strictness Strict, to the new strict % bound Bound. (see also uiu/5) uils(t_none,X,_Lin,Bound,_Sold) :- get_attr(X,itf,Att), setarg(2,Att,type(t_l(Bound))), setarg(3,Att,strictness(2)). uils(t_l(L),X,_Lin,Bound,Sold) :- TestBL is Bound - L, ( TestBL < -1.0e-10 -> true ; TestBL > 1.0e-10 -> Strict is Sold \/ 2, get_attr(X,itf,Att), setarg(2,Att,type(t_l(Bound))), setarg(3,Att,strictness(Strict)) ; Strict is Sold \/ 2, get_attr(X,itf,Att), setarg(3,Att,strictness(Strict)) ). uils(t_u(U),X,_Lin,Bound,Sold) :- U - Bound > 1.0e-10, Strict is Sold \/ 2, get_attr(X,itf,Att), setarg(2,Att,type(t_lu(Bound,U))), setarg(3,Att,strictness(Strict)). uils(t_U(U),X,_Lin,Bound,Sold) :- U - Bound > 1.0e-10, Strict is Sold \/ 2, get_attr(X,itf,Att), setarg(2,Att,type(t_lU(Bound,U))), setarg(3,Att,strictness(Strict)). uils(t_lu(L,U),X,_Lin,Bound,Sold) :- TestBL is Bound - L, ( TestBL < -1.0e-10 -> true ; TestBL > 1.0e-10 -> U - Bound > 1.0e-10, Strict is Sold \/ 2, get_attr(X,itf,Att), setarg(2,Att,type(t_lu(Bound,U))), setarg(3,Att,strictness(Strict)) ; Strict is Sold \/ 2, get_attr(X,itf,Att), setarg(3,Att,strictness(Strict)) ). uils(t_lU(L,U),X,_Lin,Bound,Sold) :- TestBL is Bound - L, ( TestBL < -1.0e-10 -> true ; TestBL > 1.0e-10 -> U - Bound > 1.0e-10, Strict is Sold \/ 2, get_attr(X,itf,Att), setarg(2,Att,type(t_lU(Bound,U))), setarg(3,Att,strictness(Strict)) ; Strict is Sold \/ 2, get_attr(X,itf,Att), setarg(3,Att,strictness(Strict)) ). uils(t_L(L),X,_Lin,Bound,Sold) :- TestBL is Bound - L, ( TestBL < -1.0e-10 -> true ; TestBL > 1.0e-10 -> Strict is Sold \/ 2, ( get_attr(X,itf,Att), arg(5,Att,order(OrdX)), arg(6,Att,class(ClassX)), ub(ClassX,OrdX,Vub-Vb-Ub), Bound - (Ub + L) > -1.0e-10 -> get_attr(X,itf,Att2), % changed? setarg(2,Att2,type(t_L(Bound))), setarg(3,Att2,strictness(Strict)), pivot_a(Vub,X,Vb,t_l(Bound)), reconsider(X) ; get_attr(X,itf,Att), arg(5,Att,order(OrdX)), arg(6,Att,class(ClassX)), setarg(2,Att,type(t_L(Bound))), setarg(3,Att,strictness(Strict)), Delta is Bound - L, backsubst_delta(ClassX,OrdX,X,Delta) ) ; Strict is Sold \/ 2, get_attr(X,itf,Att), setarg(3,Att,strictness(Strict)) ). uils(t_Lu(L,U),X,_Lin,Bound,Sold) :- TestBL is Bound - L, ( TestBL < -1.0e-10 -> true ; TestBL > 1.0e-10 -> U - Bound > 1.0e-10, Strict is Sold \/ 2, ( get_attr(X,itf,Att), arg(5,Att,order(OrdX)), arg(6,Att,class(ClassX)), ub(ClassX,OrdX,Vub-Vb-Ub), Bound - (Ub + L) > -1.0e-10 -> get_attr(X,itf,Att2), % changed? setarg(2,Att2,type(t_Lu(Bound,U))), setarg(3,Att2,strictness(Strict)), pivot_a(Vub,X,Vb,t_lu(Bound,U)), reconsider(X) ; get_attr(X,itf,Att), arg(5,Att,order(OrdX)), arg(6,Att,class(ClassX)), setarg(2,Att,type(t_Lu(Bound,U))), setarg(3,Att,strictness(Strict)), Delta is Bound - L, backsubst_delta(ClassX,OrdX,X,Delta) ) ; Strict is Sold \/ 2, get_attr(X,itf,Att), setarg(3,Att,strictness(Strict)) ). % reconsider_upper(X,Lin,U) % % Checks if the upperbound of X which is U, satisfies the bounds % of the variables in Lin: let R be the sum of all the bounds on % the variables in Lin, and I be the inhomogene part of Lin, then % upperbound U should be larger or equal to R + I (R may contain % lowerbounds). % See also rcb/3 in bv.pl reconsider_upper(X,[I,R|H],U) :- R + I - U > -1.0e-10, % violation !, dec_step(H,Status), % we want to decrement R rcbl_status(Status,X,[],Binds,[],u(U)), export_binding(Binds). reconsider_upper( _, _, _). % reconsider_lower(X,Lin,L) % % Checks if the lowerbound of X which is L, satisfies the bounds % of the variables in Lin: let R be the sum of all the bounds on % the variables in Lin, and I be the inhomogene part of Lin, then % lowerbound L should be smaller or equal to R + I (R may contain % upperbounds). % See also rcb/3 in bv.pl reconsider_lower(X,[I,R|H],L) :- R + I - L < 1.0e-10, % violation !, inc_step(H,Status), % we want to increment R rcbl_status(Status,X,[],Binds,[],l(L)), export_binding(Binds). reconsider_lower(_,_,_). % % lin is dereferenced % % solve_bound(Lin,Bound) % % Solves the linear equation Lin - Bound = 0 % Lin is the linear equation of X, a variable whose bounds have narrowed to value Bound solve_bound(Lin,Bound) :- Bound >= -1.0e-10, Bound =< 1.0e-10, !, solve(Lin). solve_bound(Lin,Bound) :- Nb is -Bound, normalize_scalar(Nb,Nbs), add_linear_11(Nbs,Lin,Eq), solve(Eq). %;bv_r % For the rhs maint. the following events are important: % % -) introduction of an indep var at active bound B % -) narrowing of active bound % -) swap active bound % -) pivot % % a variables bound (L/U) can have the states: % % -) t_none no bounds % -) t_l inactive lower bound % -) t_u inactive upper bound % -) t_L active lower bound % -) t_U active upper bound % -) t_lu inactive lower and upper bound % -) t_Lu active lower bound and inactive upper bound % -) t_lU inactive lower bound and active upper bound % ----------------------------------- deref ----------------------------------- % % deref(Lin,Lind) % % Makes a linear equation of the form [v(I,[])|H] into a solvable linear % equation. % If the variables are new, they are initialized with the linear equation X=X. deref(Lin,Lind) :- split(Lin,H,I), normalize_scalar(I,Nonvar), length(H,Len), log_deref(Len,H,[],Restd), add_linear_11(Nonvar,Restd,Lind). % log_deref(Len,[Vs|VsTail],VsTail,Res) % % Logarithmically converts a linear equation in normal form ([v(_,_)|_]) into a % linear equation in solver form ([I,R,K*X|_]). Res contains the result, Len is % the length of the part to convert and [Vs|VsTail] is a difference list % containing the equation in normal form. log_deref(0,Vs,Vs,Lin) :- !, Lin = [0.0,0.0]. log_deref(1,[v(K,[X^1])|Vs],Vs,Lin) :- !, deref_var(X,Lx), mult_linear_factor(Lx,K,Lin). log_deref(2,[v(Kx,[X^1]),v(Ky,[Y^1])|Vs],Vs,Lin) :- !, deref_var(X,Lx), deref_var(Y,Ly), add_linear_ff(Lx,Kx,Ly,Ky,Lin). log_deref(N,V0,V2,Lin) :- P is N >> 1, Q is N - P, log_deref(P,V0,V1,Lp), log_deref(Q,V1,V2,Lq), add_linear_11(Lp,Lq,Lin). % deref_var(X,Lin) % % Returns the equation of variable X. If X is a new variable, a new equation % X = X is made. deref_var(X,Lin) :- ( get_attr(X,itf,Att) -> ( \+ arg(1,Att,clpr) -> throw(error(permission_error('mix CLP(Q) variables with', 'CLP(R) variables:',X),context(_))) ; arg(4,Att,lin(Lin)) -> true ; setarg(2,Att,type(t_none)), setarg(3,Att,strictness(0)), Lin = [0.0,0.0,l(X*1.0,Ord)], setarg(4,Att,lin(Lin)), setarg(5,Att,order(Ord)) ) ; Lin = [0.0,0.0,l(X*1.0,Ord)], put_attr(X,itf,t(clpr,type(t_none),strictness(0), lin(Lin),order(Ord),n,n,n,n,n,n)) ). % TODO % % var_with_def_assign(Var,Lin) :- Lin = [I,_|Hom], ( Hom = [] -> % X=k Var = I ; Hom = [l(V*K,_)|Cs] -> ( Cs = [], TestK is K - 1.0, % K =:= 1 TestK =< 1.0e-10, TestK >= -1.0e-10, I >= -1.0e-010, % I =:= 0 I =< 1.0e-010 -> % X=Y Var = V ; % general case var_with_def_intern(t_none,Var,Lin,0) ) ). % var_with_def_intern(Type,Var,Lin,Strictness) % % Makes Lin the linear equation of new variable Var, makes all variables of % Lin, and Var of the same class and bounds Var by type(Type) and % strictness(Strictness) var_with_def_intern(Type,Var,Lin,Strict) :- put_attr(Var,itf,t(clpr,type(Type),strictness(Strict),lin(Lin), order(_),n,n,n,n,n,n)), % check uses Lin = [_,_|Hom], get_or_add_class(Var,Class), same_class(Hom,Class). % TODO % % var_intern(Type,Var,Strict) :- put_attr(Var,itf,t(clpr,type(Type),strictness(Strict), lin([0.0,0.0,l(Var*1.0,Ord)]),order(Ord),n,n,n,n,n,n)), get_or_add_class(Var,_Class). % TODO % % var_intern(Var,Class) :- % for ordered/1 but otherwise free vars get_attr(Var,itf,Att), arg(2,Att,type(_)), arg(4,Att,lin(_)), !, get_or_add_class(Var,Class). var_intern(Var,Class) :- put_attr(Var,itf,t(clpr,type(t_none),strictness(0), lin([0.0,0.0,l(Var*1.0,Ord)]),order(Ord),n,n,n,n,n,n)), get_or_add_class(Var,Class). % ----------------------------------------------------------------------------- % export_binding(Lst) % % Binds variables X to Y where Lst contains elements of the form [X-Y]. export_binding([]). export_binding([X-Y|Gs]) :- export_binding(Y,X), export_binding(Gs). % export_binding(Y,X) % % Binds variable X to Y. If Y is a nonvar and equals 0, then X is set to 0 % (numerically more stable) export_binding(Y,X) :- var(Y), Y = X. export_binding(Y,X) :- nonvar(Y), ( Y >= -1.0e-10, % Y =:= 0 Y =< 1.0e-10 -> X = 0.0 ; Y = X ). % 'solve_='(Nf) % % Solves linear equation Nf = 0 where Nf is in normal form. 'solve_='(Nf) :- deref(Nf,Nfd), % dereferences and turns Nf into solvable form Nfd solve(Nfd). % 'solve_=\\='(Nf) % % Solves linear inequality Nf =\= 0 where Nf is in normal form. 'solve_=\\='(Nf) :- deref(Nf,Lind), % dereferences and turns Nf into solvable form Lind Lind = [Inhom,_|Hom], ( Hom = [] -> % Lind = Inhom => check Inhom =\= 0 \+ (Inhom >= -1.0e-10, Inhom =< 1.0e-10) % Inhom =\= 0 ; % make new variable Nz = Lind var_with_def_intern(t_none,Nz,Lind,0), % make Nz nonzero get_attr(Nz,itf,Att), setarg(8,Att,nonzero) ). % 'solve_<'(Nf) % % Solves linear inequality Nf < 0 where Nf is in normal form. 'solve_<'(Nf) :- split(Nf,H,I), ineq(H,I,Nf,strict). % 'solve_=<'(Nf) % % Solves linear inequality Nf =< 0 where Nf is in normal form. 'solve_=<'(Nf) :- split(Nf,H,I), ineq(H,I,Nf,nonstrict). maximize(Term) :- minimize(-Term). % % This is NOT coded as minimize(Expr) :- inf(Expr,Expr). % % because the new version of inf/2 only visits % the vertex where the infimum is assumed and returns % to the 'current' vertex via backtracking. % The rationale behind this construction is to eliminate % all garbage in the solver data structures produced by % the pivots on the way to the extremal point caused by % {inf,sup}/{2,4}. % % If we are after the infimum/supremum for minimizing/maximizing, % this strategy may have adverse effects on performance because % the simplex algorithm is forced to re-discover the % extremal vertex through the equation {Inf =:= Expr}. % % Thus the extra code for {minimize,maximize}/1. % % In case someone comes up with an example where % % inf(Expr,Expr) % % outperforms the provided formulation for minimize - so be it. % Both forms are available to the user. % minimize(Term) :- wait_linear(Term,Nf,minimize_lin(Nf)). % minimize_lin(Lin) % % Minimizes the linear expression Lin. It does so by making a new % variable Dep and minimizes its value. minimize_lin(Lin) :- deref(Lin,Lind), var_with_def_intern(t_none,Dep,Lind,0), determine_active_dec(Lind), iterate_dec(Dep,Inf), { Dep =:= Inf }. sup(Expression,Sup) :- sup(Expression,Sup,[],[]). sup(Expression,Sup,Vector,Vertex) :- inf(-Expression,-Sup,Vector,Vertex). inf(Expression,Inf) :- inf(Expression,Inf,[],[]). inf(Expression,Inf,Vector,Vertex) :- % wait until Expression becomes linear, Nf contains linear Expression % in normal form wait_linear(Expression,Nf,inf_lin(Nf,Inf,Vector,Vertex)). inf_lin(Lin,_,Vector,_) :- deref(Lin,Lind), var_with_def_intern(t_none,Dep,Lind,0), % make new variable Dep = Lind determine_active_dec(Lind), % minimizes Lind iterate_dec(Dep,Inf), vertex_value(Vector,Values), nb_setval(inf,[Inf|Values]), fail. inf_lin(_,Infimum,_,Vertex) :- catch(nb_getval(inf,L),_,fail), nb_delete(inf), assign([Infimum|Vertex],L). % assign(L1,L2) % % The elements of L1 are pairwise assigned to the elements of L2 % by means of asserting {X =:= Y} where X is an element of L1 and Y % is the corresponding element of L2. assign([],[]). assign([X|Xs],[Y|Ys]) :- {X =:= Y}, % more defensive/expressive than X=Y assign(Xs,Ys). % --------------------------------- optimization ------------------------------ % % The _sn(S) =< 0 row might be temporarily infeasible. % We use reconsider/1 to fix this. % % s(S) e [_,0] = d +xi ... -xj, Rhs > 0 so we want to decrease s(S) % % positive xi would have to be moved towards their lower bound, % negative xj would have to be moved towards their upper bound, % % the row s(S) does not limit the lower bound of xi % the row s(S) does not limit the upper bound of xj % % a) if some other row R is limiting xk, we pivot(R,xk), % s(S) will decrease and get more feasible until (b) % b) if there is no limiting row for some xi: we pivot(s(S),xi) % xj: we pivot(s(S),xj) % which cures the infeasibility in one step % % iterate_dec(OptVar,Opt) % % Decreases the bound on the variables of the linear equation of OptVar as much % as possible and returns the resulting optimal bound in Opt. Fails if for some % variable, a status of unlimited is found. iterate_dec(OptVar,Opt) :- get_attr(OptVar,itf,Att), arg(4,Att,lin([I,R|H])), dec_step(H,Status), ( Status = applied -> iterate_dec(OptVar,Opt) ; Status = optimum, Opt is R + I ). % iterate_inc(OptVar,Opt) % % Increases the bound on the variables of the linear equation of OptVar as much % as possible and returns the resulting optimal bound in Opt. Fails if for some % variable, a status of unlimited is found. iterate_inc(OptVar,Opt) :- get_attr(OptVar,itf,Att), arg(4,Att,lin([I,R|H])), inc_step(H,Status), ( Status = applied -> iterate_inc(OptVar,Opt) ; Status = optimum, Opt is R + I ). % % Status = {optimum,unlimited(Indep,DepT),applied} % If Status = optimum, the tables have not been changed at all. % Searches left to right, does not try to find the 'best' pivot % Therefore we might discover unboundedness only after a few pivots % dec_step_cont([],optimum,Cont,Cont). dec_step_cont([l(V*K,OrdV)|Vs],Status,ContIn,ContOut) :- get_attr(V,itf,Att), arg(2,Att,type(W)), arg(6,Att,class(Class)), ( dec_step_2_cont(W,l(V*K,OrdV),Class,Status,ContIn,ContOut) -> true ; dec_step_cont(Vs,Status,ContIn,ContOut) ). inc_step_cont([],optimum,Cont,Cont). inc_step_cont([l(V*K,OrdV)|Vs],Status,ContIn,ContOut) :- get_attr(V,itf,Att), arg(2,Att,type(W)), arg(6,Att,class(Class)), ( inc_step_2_cont(W,l(V*K,OrdV),Class,Status,ContIn,ContOut) -> true ; inc_step_cont(Vs,Status,ContIn,ContOut) ). dec_step_2_cont(t_U(U),l(V*K,OrdV),Class,Status,ContIn,ContOut) :- K > 1.0e-10, ( lb(Class,OrdV,Vub-Vb-_) -> % found a lower bound Status = applied, pivot_a(Vub,V,Vb,t_u(U)), replace_in_cont(ContIn,Vub,V,ContOut) ; Status = unlimited(V,t_u(U)), ContIn = ContOut ). dec_step_2_cont(t_lU(L,U),l(V*K,OrdV),Class,applied,ContIn,ContOut) :- K > 1.0e-10, Init is L - U, class_basis(Class,Deps), lb(Deps,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_), pivot_b(Vub,V,Vb,t_lu(L,U)), replace_in_cont(ContIn,Vub,V,ContOut). dec_step_2_cont(t_L(L),l(V*K,OrdV),Class,Status,ContIn,ContOut) :- K < -1.0e-10, ( ub(Class,OrdV,Vub-Vb-_) -> Status = applied, pivot_a(Vub,V,Vb,t_l(L)), replace_in_cont(ContIn,Vub,V,ContOut) ; Status = unlimited(V,t_l(L)), ContIn = ContOut ). dec_step_2_cont(t_Lu(L,U),l(V*K,OrdV),Class,applied,ContIn,ContOut) :- K < -1.0e-10, Init is U - L, class_basis(Class,Deps), ub(Deps,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_), pivot_b(Vub,V,Vb,t_lu(L,U)), replace_in_cont(ContIn,Vub,V,ContOut). dec_step_2_cont(t_none,l(V*_,_),_,unlimited(V,t_none),Cont,Cont). inc_step_2_cont(t_U(U),l(V*K,OrdV),Class,Status,ContIn,ContOut) :- K < -1.0e-10, ( lb(Class,OrdV,Vub-Vb-_) -> Status = applied, pivot_a(Vub,V,Vb,t_u(U)), replace_in_cont(ContIn,Vub,V,ContOut) ; Status = unlimited(V,t_u(U)), ContIn = ContOut ). inc_step_2_cont(t_lU(L,U),l(V*K,OrdV),Class,applied,ContIn,ContOut) :- K < -1.0e-10, Init is L - U, class_basis(Class,Deps), lb(Deps,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_), pivot_b(Vub,V,Vb,t_lu(L,U)), replace_in_cont(ContIn,Vub,V,ContOut). inc_step_2_cont(t_L(L),l(V*K,OrdV),Class,Status,ContIn,ContOut) :- K > 1.0e-10, ( ub(Class,OrdV,Vub-Vb-_) -> Status = applied, pivot_a(Vub,V,Vb,t_l(L)), replace_in_cont(ContIn,Vub,V,ContOut) ; Status = unlimited(V,t_l(L)), ContIn = ContOut ). inc_step_2_cont(t_Lu(L,U),l(V*K,OrdV),Class,applied,ContIn,ContOut) :- K > 1.0e-10, Init is U - L, class_basis(Class,Deps), ub(Deps,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_), pivot_b(Vub,V,Vb,t_lu(L,U)), replace_in_cont(ContIn,Vub,V,ContOut). inc_step_2_cont(t_none,l(V*_,_),_,unlimited(V,t_none),Cont,Cont). replace_in_cont([],_,_,[]). replace_in_cont([H1|T1],X,Y,[H2|T2]) :- ( H1 == X -> H2 = Y, T1 = T2 ; H2 = H1, replace_in_cont(T1,X,Y,T2) ). dec_step([],optimum). dec_step([l(V*K,OrdV)|Vs],Status) :- get_attr(V,itf,Att), arg(2,Att,type(W)), arg(6,Att,class(Class)), ( dec_step_2(W,l(V*K,OrdV),Class,Status) -> true ; dec_step(Vs,Status) ). dec_step_2(t_U(U),l(V*K,OrdV),Class,Status) :- K > 1.0e-10, ( lb(Class,OrdV,Vub-Vb-_) -> % found a lower bound Status = applied, pivot_a(Vub,V,Vb,t_u(U)) ; Status = unlimited(V,t_u(U)) ). dec_step_2(t_lU(L,U),l(V*K,OrdV),Class,applied) :- K > 1.0e-10, Init is L - U, class_basis(Class,Deps), lb(Deps,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_), pivot_b(Vub,V,Vb,t_lu(L,U)). dec_step_2(t_L(L),l(V*K,OrdV),Class,Status) :- K < -1.0e-10, ( ub(Class,OrdV,Vub-Vb-_) -> Status = applied, pivot_a(Vub,V,Vb,t_l(L)) ; Status = unlimited(V,t_l(L)) ). dec_step_2(t_Lu(L,U),l(V*K,OrdV),Class,applied) :- K < -1.0e-10, Init is U - L, class_basis(Class,Deps), ub(Deps,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_), pivot_b(Vub,V,Vb,t_lu(L,U)). dec_step_2(t_none,l(V*_,_),_,unlimited(V,t_none)). inc_step([],optimum). % if status has not been set yet: no changes inc_step([l(V*K,OrdV)|Vs],Status) :- get_attr(V,itf,Att), arg(2,Att,type(W)), arg(6,Att,class(Class)), ( inc_step_2(W,l(V*K,OrdV),Class,Status) -> true ; inc_step(Vs,Status) ). inc_step_2(t_U(U),l(V*K,OrdV),Class,Status) :- K < -1.0e-10, ( lb(Class,OrdV,Vub-Vb-_) -> Status = applied, pivot_a(Vub,V,Vb,t_u(U)) ; Status = unlimited(V,t_u(U)) ). inc_step_2(t_lU(L,U),l(V*K,OrdV),Class,applied) :- K < -1.0e-10, Init is L - U, class_basis(Class,Deps), lb(Deps,OrdV,V-t_Lu(L,U)-Init,Vub-Vb-_), pivot_b(Vub,V,Vb,t_lu(L,U)). inc_step_2(t_L(L),l(V*K,OrdV),Class,Status) :- K > 1.0e-10, ( ub(Class,OrdV,Vub-Vb-_) -> Status = applied, pivot_a(Vub,V,Vb,t_l(L)) ; Status = unlimited(V,t_l(L)) ). inc_step_2(t_Lu(L,U),l(V*K,OrdV),Class,applied) :- K > 1.0e-10, Init is U - L, class_basis(Class,Deps), ub(Deps,OrdV,V-t_lU(L,U)-Init,Vub-Vb-_), pivot_b(Vub,V,Vb,t_lu(L,U)). inc_step_2(t_none,l(V*_,_),_,unlimited(V,t_none)). % ------------------------- find the most constraining row -------------------- % % The code for the lower and the upper bound are dual versions of each other. % The only difference is in the orientation of the comparisons. % Indeps are ruled out by their types. % If there is no bound, this fails. % % *** The actual lb and ub on an indep variable X are [lu]b + b(X), where b(X) % is the value of the active bound. % % Nota bene: We must NOT consider infeasible rows as candidates to % leave the basis! % % ub(Class,OrdX,Ub) % % See lb/3: this is similar ub(Class,OrdX,Ub) :- class_basis(Class,Deps), ub_first(Deps,OrdX,Ub). % ub_first(Deps,X,Dep-W-Ub) % % Finds the tightest upperbound for variable X from the linear equations of % basis variables Deps, and puts the resulting bound in Ub. Dep is the basis % variable that generates the bound, and W is bound of that variable that has % to be activated to achieve this. ub_first([Dep|Deps],OrdX,Tightest) :- ( get_attr(Dep,itf,Att), arg(2,Att,type(Type)), arg(4,Att,lin(Lin)), ub_inner(Type,OrdX,Lin,W,Ub), Ub > -1.0e-10 % Ub >= 0 -> ub(Deps,OrdX,Dep-W-Ub,Tightest) ; ub_first(Deps,OrdX,Tightest) ). % ub(Deps,OrdX,TightestIn,TightestOut) % % See lb/4: this is similar ub([],_,T0,T0). ub([Dep|Deps],OrdX,T0,T1) :- ( get_attr(Dep,itf,Att), arg(2,Att,type(Type)), arg(4,Att,lin(Lin)), ub_inner(Type,OrdX,Lin,W,Ub), T0 = _-Ubb, % Ub < Ubb: tighter upper bound is a smaller one Ub - Ubb < -1.0e-10, % Ub >= 0: upperbound should be larger than 0; rare failure Ub > -1.0e-10 -> ub(Deps,OrdX,Dep-W-Ub,T1) % tighter bound, use new bound ; ub(Deps,OrdX,T0,T1) % no tighter bound, keep current one ). % ub_inner(Type,OrdX,Lin,W,Ub) % % See lb_inner/5: this is similar ub_inner(t_l(L),OrdX,Lin,t_L(L),Ub) :- nf_rhs_x(Lin,OrdX,Rhs,K), % Rhs is right hand side of lin. eq. Lin containing term X*K K < -1.0e-10, % K < 0 Ub is (L-Rhs)/K. ub_inner(t_u(U),OrdX,Lin,t_U(U),Ub) :- nf_rhs_x(Lin,OrdX,Rhs,K), K > 1.0e-10, % K > 0 Ub is (U-Rhs)/K. ub_inner(t_lu(L,U),OrdX,Lin,W,Ub) :- nf_rhs_x(Lin,OrdX,Rhs,K), ( K < -1.0e-10 % K < 0, use lowerbound -> W = t_Lu(L,U), Ub = (L-Rhs)/K ; K > 1.0e-10 % K > 0, use upperbound -> W = t_lU(L,U), Ub = (U-Rhs)/K ). % lb(Class,OrdX,Lb) % % Returns in Lb how much we can lower the upperbound of X without violating % a bound of the basisvariables. % Lb has the form Dep-W-Lb with Dep the variable whose bound is violated when % lowering the bound for X more, W the actual bound that has to be activated % and Lb the amount that the upperbound can be lowered. % X has ordering OrdX and class Class. lb(Class,OrdX,Lb) :- class_basis(Class,Deps), lb_first(Deps,OrdX,Lb). % lb_first(Deps,OrdX,Tightest) % % Returns in Tightest how much we can lower the upperbound of X without % violating a bound of Deps. % Tightest has the form Dep-W-Lb with Dep the variable whose bound is violated % when lowering the bound for X more, W the actual bound that has to be % activated and Lb the amount that the upperbound can be lowered. X has % ordering attribute OrdX. lb_first([Dep|Deps],OrdX,Tightest) :- ( get_attr(Dep,itf,Att), arg(2,Att,type(Type)), arg(4,Att,lin(Lin)), lb_inner(Type,OrdX,Lin,W,Lb), Lb < 1.0e-10 % Lb =< 0: Lb > 0 means a violated bound -> lb(Deps,OrdX,Dep-W-Lb,Tightest) ; lb_first(Deps,OrdX,Tightest) ). % lb(Deps,OrdX,TightestIn,TightestOut) % % See lb_first/3: this one does the same thing, but is used for the steps after % the first one and remembers the tightest bound so far. lb([],_,T0,T0). lb([Dep|Deps],OrdX,T0,T1) :- ( get_attr(Dep,itf,Att), arg(2,Att,type(Type)), arg(4,Att,lin(Lin)), lb_inner(Type,OrdX,Lin,W,Lb), T0 = _-Lbb, Lb - Lbb > 1.0e-10, % Lb > Lbb: choose the least lowering, others % might violate bounds Lb < 1.0e-10 % Lb =< 0: violation of a bound (without lowering) -> lb(Deps,OrdX,Dep-W-Lb,T1) ; lb(Deps,OrdX,T0,T1) ). % lb_inner(Type,X,Lin,W,Lb) % % Returns in Lb how much lower we can make X without violating a bound % by using the linear equation Lin of basis variable B which has type % Type and which has to activate a bound (type W) to do so. % % E.g. when B has a lowerbound L, then L should always be smaller than I + R. % So a lowerbound of X (which has scalar K in Lin), could be at most % (L-(I+R))/K lower than its upperbound (if K is positive). % Also note that Lb should always be smaller than 0, otherwise the row is % not feasible. % X has ordering attribute OrdX. lb_inner(t_l(L),OrdX,Lin,t_L(L),Lb) :- nf_rhs_x(Lin,OrdX,Rhs,K), % if linear equation Lin contains the term % X*K, Rhs is the right hand side of that % equation K > 1.0e-10, % K > 0 Lb is (L-Rhs)/K. lb_inner(t_u(U),OrdX,Lin,t_U(U),Lb) :- nf_rhs_x(Lin,OrdX,Rhs,K), K < -1.0e-10, % K < 0 Lb is (U-Rhs)/K. lb_inner(t_lu(L,U),OrdX,Lin,W,Lb) :- nf_rhs_x(Lin,OrdX,Rhs,K), ( K < -1.0e-10 -> W = t_lU(L,U), Lb is (U-Rhs)/K ; K > 1.0e-10 -> W = t_Lu(L,U), Lb is (L-Rhs)/K ). % ---------------------------------- equations -------------------------------- % % backsubstitution will not make the system infeasible, if the bounds on the % indep vars are obeyed, but some implied values might pop up in rows where X % occurs % -) special case X=Y during bs -> get rid of dependend var(s), alias % solve(Lin) :- Lin = [I,_|H], solve(H,Lin,I,Bindings,[]), export_binding(Bindings). % solve(Hom,Lin,I,Bind,BindT) % % Solves a linear equation Lin = [I,_|H] = 0 and exports the generated bindings solve([],_,I,Bind0,Bind0) :- !, I >= -1.0e-10, % I =:= 0: redundant or trivially unsat I =< 1.0e-10. solve(H,Lin,_,Bind0,BindT) :- sd(H,[],ClassesUniq,9-9-0,Category-Selected-_,NV,NVT), get_attr(Selected,itf,Att), arg(5,Att,order(Ord)), isolate(Ord,Lin,Lin1), % Lin = 0 => Selected = Lin1 ( Category = 1 % classless variable, no bounds -> setarg(4,Att,lin(Lin1)), Lin1 = [Inhom,_|Hom], bs_collect_binding(Hom,Selected,Inhom,Bind0,BindT), eq_classes(NV,NVT,ClassesUniq) ; Category = 2 % class variable, no bounds -> arg(6,Att,class(NewC)), class_allvars(NewC,Deps), ( ClassesUniq = [_] % rank increasing -> bs_collect_bindings(Deps,Ord,Lin1,Bind0,BindT) ; Bind0 = BindT, bs(Deps,Ord,Lin1) ), eq_classes(NV,NVT,ClassesUniq) ; Category = 3 % classless variable, all variables in Lin and % Selected are bounded -> arg(2,Att,type(Type)), setarg(4,Att,lin(Lin1)), deactivate_bound(Type,Selected), eq_classes(NV,NVT,ClassesUniq), basis_add(Selected,Basis), undet_active(Lin1), % we can't tell which bound will likely be a % problem at this point Lin1 = [Inhom,_|Hom], bs_collect_binding(Hom,Selected,Inhom,Bind0,Bind1), % only if % Hom = [] rcbl(Basis,Bind1,BindT) % reconsider entire basis ; Category = 4 % class variable, all variables in Lin and Selected % are bounded -> arg(2,Att,type(Type)), arg(6,Att,class(NewC)), class_allvars(NewC,Deps), ( ClassesUniq = [_] % rank increasing -> bs_collect_bindings(Deps,Ord,Lin1,Bind0,Bind1) ; Bind0 = Bind1, bs(Deps,Ord,Lin1) ), deactivate_bound(Type,Selected), basis_add(Selected,Basis), % eq_classes( NV, NVT, ClassesUniq), % 4 -> var(NV) equate(ClassesUniq,_), undet_active(Lin1), rcbl(Basis,Bind1,BindT) ). % % Much like solve, but we solve for a particular variable of type t_none % % solve_x(H,Lin,I,X,[Bind|BindT],BindT) % % solve_x(Lin,X) :- Lin = [I,_|H], solve_x(H,Lin,I,X,Bindings,[]), export_binding(Bindings). solve_x([],_,I,_,Bind0,Bind0) :- !, I >= -1.0e-10, % I =:= 0: redundant or trivially unsat I =< 1.0e-10. solve_x(H,Lin,_,X,Bind0,BindT) :- sd(H,[],ClassesUniq,9-9-0,_,NV,NVT), get_attr(X,itf,Att), arg(5,Att,order(OrdX)), isolate(OrdX,Lin,Lin1), ( arg(6,Att,class(NewC)) -> class_allvars(NewC,Deps), ( ClassesUniq = [_] % rank increasing -> bs_collect_bindings(Deps,OrdX,Lin1,Bind0,BindT) ; Bind0 = BindT, bs(Deps,OrdX,Lin1) ), eq_classes(NV,NVT,ClassesUniq) ; setarg(4,Att,lin(Lin1)), Lin1 = [Inhom,_|Hom], bs_collect_binding(Hom,X,Inhom,Bind0,BindT), eq_classes(NV,NVT,ClassesUniq) ). % solve_ord_x(Lin,OrdX,ClassX) % % Does the same thing as solve_x/2, but has the ordering of X and its class as % input. This also means that X has a class which is not sure in solve_x/2. solve_ord_x(Lin,OrdX,ClassX) :- Lin = [I,_|H], solve_ord_x(H,Lin,I,OrdX,ClassX,Bindings,[]), export_binding(Bindings). solve_ord_x([],_,I,_,_,Bind0,Bind0) :- I >= -1.0e-10, % I =:= 0 I =< 1.0e-10. solve_ord_x([_|_],Lin,_,OrdX,ClassX,Bind0,BindT) :- isolate(OrdX,Lin,Lin1), Lin1 = [_,_|H1], sd(H1,[],ClassesUniq1,9-9-0,_,NV,NVT), % do sd on Lin without X, then % add class of X ord_add_element(ClassesUniq1,ClassX,ClassesUniq), class_allvars(ClassX,Deps), ( ClassesUniq = [_] % rank increasing -> bs_collect_bindings(Deps,OrdX,Lin1,Bind0,BindT) ; Bind0 = BindT, bs(Deps,OrdX,Lin1) ), eq_classes(NV,NVT,ClassesUniq). % sd(H,[],ClassesUniq,9-9-0,Category-Selected-_,NV,NVT) % sd(Hom,ClassesIn,ClassesOut,PreferenceIn,PreferenceOut,[NV|NVTail],NVTail) % % ClassesOut is a sorted list of the different classes that are either in % ClassesIn or that are the classes of the variables in Hom. Variables that do % not belong to a class yet, are put in the difference list NV. sd([],Class0,Class0,Preference0,Preference0,NV0,NV0). sd([l(X*K,_)|Xs],Class0,ClassN,Preference0,PreferenceN,NV0,NVt) :- get_attr(X,itf,Att), ( arg(6,Att,class(Xc)) % old: has class -> NV0 = NV1, ord_add_element(Class0,Xc,Class1), ( arg(2,Att,type(t_none)) -> preference(Preference0,2-X-K,Preference1) % has class, no bounds => category 2 ; preference(Preference0,4-X-K,Preference1) % has class, is bounded => category 4 ) ; % new: has no class Class1 = Class0, NV0 = [X|NV1], % X has no class yet, add to list of new variables ( arg(2,Att,type(t_none)) -> preference(Preference0,1-X-K,Preference1) % no class, no bounds => category 1 ; preference(Preference0,3-X-K,Preference1) % no class, is bounded => category 3 ) ), sd(Xs,Class1,ClassN,Preference1,PreferenceN,NV1,NVt). % % A is best sofar, B is current % smallest prefered preference(A,B,Pref) :- A = Px-_-_, B = Py-_-_, ( Px < Py -> Pref = A ; Pref = B ). % eq_classes(NV,NVTail,Cs) % % Attaches all classless variables NV to a new class and equates all other % classes with this class. The equate operation only happens after attach_class % because the unification of classes can bind the tail of the AllVars attribute % to a nonvar and then the attach_class operation wouldn't work. eq_classes(NV,_,Cs) :- var(NV), !, equate(Cs,_). eq_classes(NV,NVT,Cs) :- class_new(Su,clpr,NV,NVT,[]), % make a new class Su with NV as the variables attach_class(NV,Su), % attach the variables NV to Su equate(Cs,Su). equate([],_). equate([X|Xs],X) :- equate(Xs,X). % % assert: none of the Vars has a class attribute yet % attach_class(Xs,_) :- var(Xs), % Tail !. attach_class([X|Xs],Class) :- get_attr(X,itf,Att), setarg(6,Att,class(Class)), attach_class(Xs,Class). % unconstrained(Lin,Uc,Kuc,Rest) % % Finds an unconstrained variable Uc (type(t_none)) in Lin with scalar Kuc and % removes it from Lin to return Rest. unconstrained(Lin,Uc,Kuc,Rest) :- Lin = [_,_|H], sd(H,[],_,9-9-0,Category-Uc-_,_,_), Category =< 2, get_attr(Uc,itf,Att), arg(5,Att,order(OrdUc)), delete_factor(OrdUc,Lin,Rest,Kuc). % % point the vars in Lin into the same equivalence class % maybe join some global data % same_class([],_). same_class([l(X*_,_)|Xs],Class) :- get_or_add_class(X,Class), same_class(Xs,Class). % get_or_add_class(X,Class) % % Returns in Class the class of X if X has one, or a new class where X now % belongs to if X didn't have one. get_or_add_class(X,Class) :- get_attr(X,itf,Att), arg(1,Att,CLP), ( arg(6,Att,class(ClassX)) -> ClassX = Class ; setarg(6,Att,class(Class)), class_new(Class,CLP,[X|Tail],Tail,[]) ). % allvars(X,Allvars) % % Allvars is a list of all variables in the class to which X belongs. allvars(X,Allvars) :- get_attr(X,itf,Att), arg(6,Att,class(C)), class_allvars(C,Allvars). % deactivate_bound(Type,Variable) % % The Type of the variable is changed to reflect the deactivation of its % bounds. % t_L(_) becomes t_l(_), t_lU(_,_) becomes t_lu(_,_) and so on. deactivate_bound(t_l(_),_). deactivate_bound(t_u(_),_). deactivate_bound(t_lu(_,_),_). deactivate_bound(t_L(L),X) :- get_attr(X,itf,Att), setarg(2,Att,type(t_l(L))). deactivate_bound(t_Lu(L,U),X) :- get_attr(X,itf,Att), setarg(2,Att,type(t_lu(L,U))). deactivate_bound(t_U(U),X) :- get_attr(X,itf,Att), setarg(2,Att,type(t_u(U))). deactivate_bound(t_lU(L,U),X) :- get_attr(X,itf,Att), setarg(2,Att,type(t_lu(L,U))). % intro_at(X,Value,Type) % % Variable X gets new type Type which reflects the activation of a bound with % value Value. In the linear equations of all the variables belonging to the % same class as X, X is substituted by [0,Value,X] to reflect the new active % bound. intro_at(X,Value,Type) :- get_attr(X,itf,Att), arg(5,Att,order(Ord)), arg(6,Att,class(Class)), setarg(2,Att,type(Type)), ( Value >= -1.0e-10, % Value =:= 0 Value =< 1.0e-010 -> true ; backsubst_delta(Class,Ord,X,Value) ). % undet_active(Lin) % % For each variable in the homogene part of Lin, a bound is activated % if an inactive bound exists. (t_l(L) becomes t_L(L) and so on) undet_active([_,_|H]) :- undet_active_h(H). % undet_active_h(Hom) % % For each variable in homogene part Hom, a bound is activated if an % inactive bound exists (t_l(L) becomes t_L(L) and so on) undet_active_h([]). undet_active_h([l(X*_,_)|Xs]) :- get_attr(X,itf,Att), arg(2,Att,type(Type)), undet_active(Type,X), undet_active_h(Xs). % undet_active(Type,Var) % % An inactive bound of Var is activated if such exists % t_lu(L,U) is arbitrarily chosen to become t_Lu(L,U) undet_active(t_none,_). % type_activity undet_active(t_L(_),_). undet_active(t_Lu(_,_),_). undet_active(t_U(_),_). undet_active(t_lU(_,_),_). undet_active(t_l(L),X) :- intro_at(X,L,t_L(L)). undet_active(t_u(U),X) :- intro_at(X,U,t_U(U)). undet_active(t_lu(L,U),X) :- intro_at(X,L,t_Lu(L,U)). % determine_active_dec(Lin) % % Activates inactive bounds on the variables of Lin if such bounds exist. % If the type of a variable is t_none, this fails. This version is aimed % to make the R component of Lin as small as possible in order not to violate % an upperbound (see reconsider/1) determine_active_dec([_,_|H]) :- determine_active(H,-1). % determine_active_inc(Lin) % % Activates inactive bounds on the variables of Lin if such bounds exist. % If the type of a variable is t_none, this fails. This version is aimed % to make the R component of Lin as large as possible in order not to violate % a lowerbound (see reconsider/1) determine_active_inc([_,_|H]) :- determine_active(H,1). % determine_active(Hom,S) % % For each variable in Hom, activates its bound if it is not yet activated. % For the case of t_lu(_,_) the lower or upper bound is activated depending on % K and S: % If sign of K*S is negative, then lowerbound, otherwise upperbound. determine_active([],_). determine_active([l(X*K,_)|Xs],S) :- get_attr(X,itf,Att), arg(2,Att,type(Type)), determine_active(Type,X,K,S), determine_active(Xs,S). determine_active(t_L(_),_,_,_). determine_active(t_Lu(_,_),_,_,_). determine_active(t_U(_),_,_,_). determine_active(t_lU(_,_),_,_,_). determine_active(t_l(L),X,_,_) :- intro_at(X,L,t_L(L)). determine_active(t_u(U),X,_,_) :- intro_at(X,U,t_U(U)). determine_active(t_lu(L,U),X,K,S) :- TestKs is K*S, ( TestKs < -1.0e-10 % K*S < 0 -> intro_at(X,L,t_Lu(L,U)) ; TestKs > 1.0e-10 -> intro_at(X,U,t_lU(L,U)) ). % % Careful when an indep turns into t_none !!! % detach_bounds(V) :- get_attr(V,itf,Att), arg(2,Att,type(Type)), arg(4,Att,lin(Lin)), arg(5,Att,order(OrdV)), arg(6,Att,class(Class)), setarg(2,Att,type(t_none)), setarg(3,Att,strictness(0)), ( indep(Lin,OrdV) -> ( ub(Class,OrdV,Vub-Vb-_) -> % exchange against thightest class_basis_drop(Class,Vub), pivot(Vub,Class,OrdV,Vb,Type) ; lb(Class,OrdV,Vlb-Vb-_) -> class_basis_drop(Class,Vlb), pivot(Vlb,Class,OrdV,Vb,Type) ; true ) ; class_basis_drop(Class,V) ). detach_bounds_vlv(OrdV,Lin,Class,Var,NewLin) :- ( indep(Lin,OrdV) -> Lin = [_,R|_], ( ub(Class,OrdV,Vub-Vb-_) -> % in verify_lin, class might contain two occurrences of Var, % but it doesn't matter which one we delete class_basis_drop(Class,Var), pivot_vlv(Vub,Class,OrdV,Vb,R,NewLin) ; lb(Class,OrdV,Vlb-Vb-_) -> class_basis_drop(Class,Var), pivot_vlv(Vlb,Class,OrdV,Vb,R,NewLin) ; NewLin = Lin ) ; NewLin = Lin, class_basis_drop(Class,Var) ). % ----------------------------- manipulate the basis -------------------------- % basis_drop(X) % % Removes X from the basis of the class to which X belongs. basis_drop(X) :- get_attr(X,itf,Att), arg(6,Att,class(Cv)), class_basis_drop(Cv,X). % basis(X,Basis) % % Basis is the basis of the class to which X belongs. basis(X,Basis) :- get_attr(X,itf,Att), arg(6,Att,class(Cv)), class_basis(Cv,Basis). % basis_add(X,NewBasis) % % NewBasis is the result of adding X to the basis of the class to which X % belongs. basis_add(X,NewBasis) :- get_attr(X,itf,Att), arg(6,Att,class(Cv)), class_basis_add(Cv,X,NewBasis). % basis_pivot(Leave,Enter) % % Removes Leave from the basis of the class to which it belongs, and adds % Enter to that basis. basis_pivot(Leave,Enter) :- get_attr(Leave,itf,Att), arg(6,Att,class(Cv)), class_basis_pivot(Cv,Enter,Leave). % ----------------------------------- pivot ----------------------------------- % pivot(Dep,Indep) % % The linear equation of variable Dep, is transformed into one of variable % Indep, containing Dep. Then, all occurrences of Indep in linear equations are % substituted by this new definition % % Pivot ignoring rhs and active states % pivot(Dep,Indep) :- get_attr(Dep,itf,AttD), arg(4,AttD,lin(H)), arg(5,AttD,order(OrdDep)), get_attr(Indep,itf,AttI), arg(5,AttI,order(Ord)), arg(5,AttI,class(Class)), delete_factor(Ord,H,H0,Coeff), K is -1.0/Coeff, add_linear_ff(H0,K,[0.0,0.0,l(Dep* -1.0,OrdDep)],K,Lin), backsubst(Class,Ord,Lin). % pivot_a(Dep,Indep,IndepT,DepT) % % Removes Dep from the basis, puts Indep in, and pivots the equation of % Dep to become one of Indep. The type of Dep becomes DepT (which means % it gets deactivated), the type of Indep becomes IndepT (which means it % gets activated) pivot_a(Dep,Indep,Vb,Wd) :- basis_pivot(Dep,Indep), get_attr(Indep,itf,Att), arg(2,Att,type(Type)), arg(5,Att,order(Ord)), arg(6,Att,class(Class)), pivot(Dep,Class,Ord,Vb,Type), get_attr(Indep,itf,Att2), %changed? setarg(2,Att2,type(Wd)). pivot_b(Vub,V,Vb,Wd) :- ( Vub == V -> get_attr(V,itf,Att), arg(5,Att,order(Ord)), arg(6,Att,class(Class)), setarg(2,Att,type(Vb)), pivot_b_delta(Vb,Delta), % nonzero(Delta) backsubst_delta(Class,Ord,V,Delta) ; pivot_a(Vub,V,Vb,Wd) ). pivot_b_delta(t_Lu(L,U),Delta) :- Delta is L-U. pivot_b_delta(t_lU(L,U),Delta) :- Delta is U-L. % select_active_bound(Type,Bound) % % Returns the bound that is active in Type (if such exists, 0 otherwise) select_active_bound(t_L(L),L). select_active_bound(t_Lu(L,_),L). select_active_bound(t_U(U),U). select_active_bound(t_lU(_,U),U). select_active_bound(t_none,0.0). % % for project.pl % select_active_bound(t_l(_),0.0). select_active_bound(t_u(_),0.0). select_active_bound(t_lu(_,_),0.0). % pivot(Dep,Class,IndepOrd,DepAct,IndAct) % % See pivot/2. % In addition, variable Indep with ordering IndepOrd has an active bound IndAct % % % Pivot taking care of rhs and active states % pivot(Dep,Class,IndepOrd,DepAct,IndAct) :- get_attr(Dep,itf,Att), arg(4,Att,lin(H)), arg(5,Att,order(DepOrd)), setarg(2,Att,type(DepAct)), select_active_bound(DepAct,AbvD), % New current value for Dep select_active_bound(IndAct,AbvI), % New current value of Indep delete_factor(IndepOrd,H,H0,Coeff), % Dep = ... + Coeff*Indep + ... AbvDm is -AbvD, AbvIm is -AbvI, add_linear_f1([0.0,AbvIm],Coeff,H0,H1), K is -1.0/Coeff, add_linear_ff(H1,K,[0.0,AbvDm,l(Dep* -1.0,DepOrd)],K,H2), % Indep = -1/Coeff*... + 1/Coeff*Dep add_linear_11(H2,[0.0,AbvIm],Lin), backsubst(Class,IndepOrd,Lin). pivot_vlv(Dep,Class,IndepOrd,DepAct,AbvI,Lin) :- get_attr(Dep,itf,Att), arg(4,Att,lin(H)), arg(5,Att,order(DepOrd)), setarg(2,Att,type(DepAct)), select_active_bound(DepAct,AbvD), % New current value for Dep delete_factor(IndepOrd,H,H0,Coeff), % Dep = ... + Coeff*Indep + ... AbvDm is -AbvD, AbvIm is -AbvI, add_linear_f1([0.0,AbvIm],Coeff,H0,H1), K is -1.0/Coeff, add_linear_ff(H1,K,[0.0,AbvDm,l(Dep* -1.0,DepOrd)],K,Lin), % Indep = -1/Coeff*... + 1/Coeff*Dep add_linear_11(Lin,[0.0,AbvIm],SubstLin), backsubst(Class,IndepOrd,SubstLin). % backsubst_delta(Class,OrdX,X,Delta) % % X with ordering attribute OrdX, is substituted in all linear equations of % variables in the class Class, by linear equation [0,Delta,l(X*1,OrdX)]. This % reflects the activation of a bound. backsubst_delta(Class,OrdX,X,Delta) :- backsubst(Class,OrdX,[0.0,Delta,l(X*1.0,OrdX)]). % backsubst(Class,OrdX,Lin) % % X with ordering OrdX is substituted in all linear equations of variables in % the class Class, by linear equation Lin backsubst(Class,OrdX,Lin) :- class_allvars(Class,Allvars), bs(Allvars,OrdX,Lin). % bs(Vars,OrdV,Lin) % % In all linear equations of the variables Vars, variable V with ordering % attribute OrdV is substituted by linear equation Lin. % % valid if nothing will go ground % bs(Xs,_,_) :- var(Xs), !. bs([X|Xs],OrdV,Lin) :- ( get_attr(X,itf,Att), arg(4,Att,lin(LinX)), nf_substitute(OrdV,Lin,LinX,LinX1) % does not change attributes -> setarg(4,Att,lin(LinX1)), bs(Xs,OrdV,Lin) ; bs(Xs,OrdV,Lin) ). % % rank increasing backsubstitution % % bs_collect_bindings(Deps,SelectedOrd,Lin,Bind,BindT) % % Collects bindings (of the form [X-I] where X = I is the binding) by % substituting Selected in all linear equations of the variables Deps (which % are of the same class), by Lin. Selected has ordering attribute SelectedOrd. % % E.g. when V = 2X + 3Y + 4, X = 3V + 2Z and Y = 4X + 3 % we can substitute V in the linear equation of X: X = 6X + 9Y + 2Z + 12 % we can't substitute V in the linear equation of Y of course. bs_collect_bindings(Xs,_,_,Bind0,BindT) :- var(Xs), !, Bind0 = BindT. bs_collect_bindings([X|Xs],OrdV,Lin,Bind0,BindT) :- ( get_attr(X,itf,Att), arg(4,Att,lin(LinX)), nf_substitute(OrdV,Lin,LinX,LinX1) % does not change attributes -> setarg(4,Att,lin(LinX1)), LinX1 = [Inhom,_|Hom], bs_collect_binding(Hom,X,Inhom,Bind0,Bind1), bs_collect_bindings(Xs,OrdV,Lin,Bind1,BindT) ; bs_collect_bindings(Xs,OrdV,Lin,Bind0,BindT) ). % bs_collect_binding(Hom,Selected,Inhom,Bind,BindT) % % Collects binding following from Selected = Hom + Inhom. % If Hom = [], returns the binding Selected-Inhom (=0) % bs_collect_binding([],X,Inhom) --> [X-Inhom]. bs_collect_binding([_|_],_,_) --> []. % % reconsider the basis % % rcbl(Basis,Bind,BindT) % % rcbl([],Bind0,Bind0). rcbl([X|Continuation],Bind0,BindT) :- ( rcb_cont(X,Status,Violated,Continuation,NewContinuation) % have a culprit -> rcbl_status(Status,X,NewContinuation,Bind0,BindT,Violated) ; rcbl(Continuation,Bind0,BindT) ). rcb_cont(X,Status,Violated,ContIn,ContOut) :- get_attr(X,itf,Att), arg(2,Att,type(Type)), arg(4,Att,lin([I,R|H])), ( Type = t_l(L) % case 1: lowerbound: R + I should always be larger % than the lowerbound -> R + I - L < 1.0e-10, Violated = l(L), inc_step_cont(H,Status,ContIn,ContOut) ; Type = t_u(U) % case 2: upperbound: R + I should always be smaller % than the upperbound -> R + I - U > -1.0e-10, Violated = u(U), dec_step_cont(H,Status,ContIn,ContOut) ; Type = t_lu(L,U) % case 3: check both -> At is R + I, ( At - L < 1.0e-10 -> Violated = l(L), inc_step_cont(H,Status,ContIn,ContOut) ; At - U > -1.0e-10 -> Violated = u(U), dec_step_cont(H,Status,ContIn,ContOut) ) ). % other types imply nonbasic variable or unbounded variable % % reconsider one element of the basis % later: lift the binds % reconsider(X) :- rcb(X,Status,Violated), !, rcbl_status(Status,X,[],Binds,[],Violated), export_binding(Binds). reconsider(_). % % Find a basis variable out of its bound or at its bound % Try to move it into whithin its bound % a) impossible -> fail % b) optimum at the bound -> implied value % c) else look at the remaining basis variables % % % Idea: consider a variable V with linear equation Lin. % When a bound on a variable X of Lin gets activated, its value, multiplied % with the scalar of X, is added to the R component of Lin. % When we consider the lowerbound of V, it must be smaller than R + I, since R % contains at best the lowerbounds of the variables in Lin (but could contain % upperbounds, which are of course larger). So checking this can show the % violation of a bound of V. A similar case works for the upperbound. rcb(X,Status,Violated) :- get_attr(X,itf,Att), arg(2,Att,type(Type)), arg(4,Att,lin([I,R|H])), ( Type = t_l(L) % case 1: lowerbound: R + I should always be larger % than the lowerbound -> R + I - L < 1.0e-10, % R + I =< L Violated = l(L), inc_step(H,Status) ; Type = t_u(U) % case 2: upperbound: R + I should always be smaller % than the upperbound -> R + I - U > -1.0e-10, % R + I >= U Violated = u(U), dec_step(H,Status) ; Type = t_lu(L,U) % case 3: check both -> At is R + I, ( At - L < 1.0e-10 % At =< L -> Violated = l(L), inc_step(H,Status) ; At - U > -1.0e-10 % At >= U -> Violated = u(U), dec_step(H,Status) ) ). % other types imply nonbasic variable or unbounded variable % rcbl_status(Status,X,Continuation,[Bind|BindT],BindT,Violated) % % rcbl_status(optimum,X,Cont,B0,Bt,Violated) :- rcbl_opt(Violated,X,Cont,B0,Bt). rcbl_status(applied,X,Cont,B0,Bt,Violated) :- rcbl_app(Violated,X,Cont,B0,Bt). rcbl_status(unlimited(Indep,DepT),X,Cont,B0,Bt,Violated) :- rcbl_unl(Violated,X,Cont,B0,Bt,Indep,DepT). % % Might reach optimum immediately without changing the basis, % but in general we must assume that there were pivots. % If the optimum meets the bound, we backsubstitute the implied % value, solve will call us again to check for further implied % values or unsatisfiability in the rank increased system. % rcbl_opt(l(L),X,Continuation,B0,B1) :- get_attr(X,itf,Att), arg(2,Att,type(Type)), arg(3,Att,strictness(Strict)), arg(4,Att,lin(Lin)), Lin = [I,R|_], Opt is R + I, TestLO is L - Opt, ( TestLO < -1.0e-10 % L < Opt -> narrow_u(Type,X,Opt), % { X =< Opt } rcbl(Continuation,B0,B1) ; TestLO =< 1.0e-10, % L = Opt Strict /\ 2 =:= 0, % meets lower Mop is -Opt, normalize_scalar(Mop,MopN), add_linear_11(MopN,Lin,Lin1), Lin1 = [Inhom,_|Hom], ( Hom = [] -> rcbl(Continuation,B0,B1) % would not callback ; solve(Hom,Lin1,Inhom,B0,B1) ) ). rcbl_opt(u(U),X,Continuation,B0,B1) :- get_attr(X,itf,Att), arg(2,Att,type(Type)), arg(3,Att,strictness(Strict)), arg(4,Att,lin(Lin)), Lin = [I,R|_], Opt is R + I, TestUO is U - Opt, ( TestUO > 1.0e-10 % U > Opt -> narrow_l(Type,X,Opt), % { X >= Opt } rcbl(Continuation,B0,B1) ; TestUO >= -1.0e-10, % U = Opt Strict /\ 1 =:= 0, % meets upper Mop is -Opt, normalize_scalar(Mop,MopN), add_linear_11(MopN,Lin,Lin1), Lin1 = [Inhom,_|Hom], ( Hom = [] -> rcbl(Continuation,B0,B1) % would not callback ; solve(Hom,Lin1,Inhom,B0,B1) ) ). % % Basis has already changed when this is called % rcbl_app(l(L),X,Continuation,B0,B1) :- get_attr(X,itf,Att), arg(4,Att,lin([I,R|H])), ( R + I - L > 1.0e-10 % R+I > L: within bound now -> rcbl(Continuation,B0,B1) ; inc_step(H,Status), rcbl_status(Status,X,Continuation,B0,B1,l(L)) ). rcbl_app(u(U),X,Continuation,B0,B1) :- get_attr(X,itf,Att), arg(4,Att,lin([I,R|H])), ( R + I - U < -1.0e-10 % R+I < U: within bound now -> rcbl(Continuation,B0,B1) ; dec_step(H,Status), rcbl_status(Status,X,Continuation,B0,B1,u(U)) ). % % This is never called for a t_lu culprit % rcbl_unl(l(L),X,Continuation,B0,B1,Indep,DepT) :- pivot_a(X,Indep,t_L(L),DepT), % changes the basis rcbl(Continuation,B0,B1). rcbl_unl(u(U),X,Continuation,B0,B1,Indep,DepT) :- pivot_a(X,Indep,t_U(U),DepT), % changes the basis rcbl(Continuation,B0,B1). % narrow_u(Type,X,U) % % Narrows down the upperbound of X (type Type) to U. % Fails if Type is not t_u(_) or t_lu(_) narrow_u(t_u(_),X,U) :- get_attr(X,itf,Att), setarg(2,Att,type(t_u(U))). narrow_u(t_lu(L,_),X,U) :- get_attr(X,itf,Att), setarg(2,Att,type(t_lu(L,U))). % narrow_l(Type,X,L) % % Narrows down the lowerbound of X (type Type) to L. % Fails if Type is not t_l(_) or t_lu(_) narrow_l( t_l(_), X, L) :- get_attr(X,itf,Att), setarg(2,Att,type(t_l(L))). narrow_l( t_lu(_,U), X, L) :- get_attr(X,itf,Att), setarg(2,Att,type(t_lu(L,U))). % ----------------------------------- dump ------------------------------------ % dump_var(Type,Var,I,H,Dump,DumpTail) % % Returns in Dump a representation of the linear constraint on variable % Var which has linear equation H + I and has type Type. dump_var(t_none,V,I,H) --> !, ( { H = [l(W*K,_)], V == W, I >= -1.0e-10, % I=:=0 I =< 1.0e-010, TestK is K - 1.0, % K=:=1 TestK >= -1.0e-10, TestK =< 1.0e-10 } -> % indep var [] ; {nf2sum(H,I,Sum)}, [V = Sum] ). dump_var(t_L(L),V,I,H) --> !, dump_var(t_l(L),V,I,H). % case lowerbound: V >= L or V > L % say V >= L, and V = K*V1 + ... + I, then K*V1 + ... + I >= L % and K*V1 + ... >= L-I and V1 + .../K = (L-I)/K dump_var(t_l(L),V,I,H) --> !, { H = [l(_*K,_)|_], % avoid 1 >= 0 get_attr(V,itf,Att), arg(3,Att,strictness(Strict)), Sm is Strict /\ 2, Kr is 1.0/K, Li is Kr*(L - I), mult_hom(H,Kr,H1), nf2sum(H1,0.0,Sum), ( K > 1.0e-10 % K > 0 -> dump_strict(Sm,Sum >= Li,Sum > Li,Result) ; dump_strict(Sm,Sum =< Li,Sum < Li,Result) ) }, [Result]. dump_var(t_U(U),V,I,H) --> !, dump_var(t_u(U),V,I,H). dump_var(t_u(U),V,I,H) --> !, { H = [l(_*K,_)|_], % avoid 0 =< 1 get_attr(V,itf,Att), arg(3,Att,strictness(Strict)), Sm is Strict /\ 1, Kr is 1.0/K, Ui is Kr*(U-I), mult_hom(H,Kr,H1), nf2sum(H1,0.0,Sum), ( K > 1.0e-10 % K > 0 -> dump_strict(Sm,Sum =< Ui,Sum < Ui,Result) ; dump_strict(Sm,Sum >= Ui,Sum > Ui,Result) ) }, [Result]. dump_var(t_Lu(L,U),V,I,H) --> !, dump_var(t_l(L),V,I,H), dump_var(t_u(U),V,I,H). dump_var(t_lU(L,U),V,I,H) --> !, dump_var(t_l(L),V,I,H), dump_var(t_u(U),V,I,H). dump_var(t_lu(L,U),V,I,H) --> !, dump_var(t_l(L),V,I,H), dump_var(t_U(U),V,I,H). dump_var(T,V,I,H) --> % should not happen [V:T:I+H]. % dump_strict(FilteredStrictness,Nonstrict,Strict,Res) % % Unifies Res with either Nonstrict or Strict depending on FilteredStrictness. % FilteredStrictness is the component of strictness related to the bound: 0 % means nonstrict, 1 means strict upperbound, 2 means strict lowerbound, % 3 is filtered out to either 1 or 2. dump_strict(0,Result,_,Result). dump_strict(1,_,Result,Result). dump_strict(2,_,Result,Result). % dump_nz(V,H,I,Dump,DumpTail) % % Returns in Dump a representation of the nonzero constraint of variable V % which has linear % equation H + I. dump_nz(_,H,I) --> { H = [l(_*K,_)|_], Kr is 1.0/K, I1 is -Kr*I, mult_hom(H,Kr,H1), nf2sum(H1,0.0,Sum) }, [Sum =\= I1]. /******************************* * SANDBOX * *******************************/ %;store_r % normalize_scalar(S,[N,Z]) % % Transforms a scalar S into a linear expression [S,0] normalize_scalar(S,[S,0.0]). % renormalize(List,Lin) % % Renormalizes the not normalized linear expression in List into % a normalized one. It does so to take care of unifications. % (e.g. when a variable X is bound to a constant, the constant is added to % the constant part of the linear expression; when a variable X is bound to % another variable Y, the scalars of both are added) renormalize([I,R|Hom],Lin) :- length(Hom,Len), renormalize_log(Len,Hom,[],Lin0), add_linear_11([I,R],Lin0,Lin). % renormalize_log(Len,Hom,HomTail,Lin) % % Logarithmically renormalizes the homogene part of a not normalized % linear expression. See also renormalize/2. renormalize_log(1,[Term|Xs],Xs,Lin) :- !, Term = l(X*_,_), renormalize_log_one(X,Term,Lin). renormalize_log(2,[A,B|Xs],Xs,Lin) :- !, A = l(X*_,_), B = l(Y*_,_), renormalize_log_one(X,A,LinA), renormalize_log_one(Y,B,LinB), add_linear_11(LinA,LinB,Lin). renormalize_log(N,L0,L2,Lin) :- P is N>>1, Q is N-P, renormalize_log(P,L0,L1,Lp), renormalize_log(Q,L1,L2,Lq), add_linear_11(Lp,Lq,Lin). % renormalize_log_one(X,Term,Res) % % Renormalizes a term in X: if X is a nonvar, the term becomes a scalar. renormalize_log_one(X,Term,Res) :- var(X), Term = l(X*K,_), get_attr(X,itf,Att), arg(5,Att,order(OrdX)), % Order might have changed Res = [0.0,0.0,l(X*K,OrdX)]. renormalize_log_one(X,Term,Res) :- nonvar(X), Term = l(X*K,_), Xk is X*K, normalize_scalar(Xk,Res). % ----------------------------- sparse vector stuff ---------------------------- % % add_linear_ff(LinA,Ka,LinB,Kb,LinC) % % Linear expression LinC is the result of the addition of the 2 linear expressions % LinA and LinB, each one multiplied by a scalar (Ka for LinA and Kb for LinB). add_linear_ff(LinA,Ka,LinB,Kb,LinC) :- LinA = [Ia,Ra|Ha], LinB = [Ib,Rb|Hb], LinC = [Ic,Rc|Hc], Ic is Ia*Ka+Ib*Kb, Rc is Ra*Ka+Rb*Kb, add_linear_ffh(Ha,Ka,Hb,Kb,Hc). % add_linear_ffh(Ha,Ka,Hb,Kb,Hc) % % Homogene part Hc is the result of the addition of the 2 homogene parts Ha and Hb, % each one multiplied by a scalar (Ka for Ha and Kb for Hb) add_linear_ffh([],_,Ys,Kb,Zs) :- mult_hom(Ys,Kb,Zs). add_linear_ffh([l(X*Kx,OrdX)|Xs],Ka,Ys,Kb,Zs) :- add_linear_ffh(Ys,X,Kx,OrdX,Xs,Zs,Ka,Kb). % add_linear_ffh(Ys,X,Kx,OrdX,Xs,Zs,Ka,Kb) % % Homogene part Zs is the result of the addition of the 2 homogene parts Ys and % [l(X*Kx,OrdX)|Xs], each one multiplied by a scalar (Ka for [l(X*Kx,OrdX)|Xs] and Kb for Ys) add_linear_ffh([],X,Kx,OrdX,Xs,Zs,Ka,_) :- mult_hom([l(X*Kx,OrdX)|Xs],Ka,Zs). add_linear_ffh([l(Y*Ky,OrdY)|Ys],X,Kx,OrdX,Xs,Zs,Ka,Kb) :- compare(Rel,OrdX,OrdY), ( Rel = (=) -> Kz is Kx*Ka+Ky*Kb, ( % Kz =:= 0 Kz =< 1.0e-10, Kz >= -1.0e-10 -> add_linear_ffh(Xs,Ka,Ys,Kb,Zs) ; Zs = [l(X*Kz,OrdX)|Ztail], add_linear_ffh(Xs,Ka,Ys,Kb,Ztail) ) ; Rel = (<) -> Zs = [l(X*Kz,OrdX)|Ztail], Kz is Kx*Ka, add_linear_ffh(Xs,Y,Ky,OrdY,Ys,Ztail,Kb,Ka) ; Rel = (>) -> Zs = [l(Y*Kz,OrdY)|Ztail], Kz is Ky*Kb, add_linear_ffh(Ys,X,Kx,OrdX,Xs,Ztail,Ka,Kb) ). % add_linear_f1(LinA,Ka,LinB,LinC) % % special case of add_linear_ff with Kb = 1 add_linear_f1(LinA,Ka,LinB,LinC) :- LinA = [Ia,Ra|Ha], LinB = [Ib,Rb|Hb], LinC = [Ic,Rc|Hc], Ic is Ia*Ka+Ib, Rc is Ra*Ka+Rb, add_linear_f1h(Ha,Ka,Hb,Hc). % add_linear_f1h(Ha,Ka,Hb,Hc) % % special case of add_linear_ffh/5 with Kb = 1 add_linear_f1h([],_,Ys,Ys). add_linear_f1h([l(X*Kx,OrdX)|Xs],Ka,Ys,Zs) :- add_linear_f1h(Ys,X,Kx,OrdX,Xs,Zs,Ka). % add_linear_f1h(Ys,X,Kx,OrdX,Xs,Zs,Ka) % % special case of add_linear_ffh/8 with Kb = 1 add_linear_f1h([],X,Kx,OrdX,Xs,Zs,Ka) :- mult_hom([l(X*Kx,OrdX)|Xs],Ka,Zs). add_linear_f1h([l(Y*Ky,OrdY)|Ys],X,Kx,OrdX,Xs,Zs,Ka) :- compare(Rel,OrdX,OrdY), ( Rel = (=) -> Kz is Kx*Ka+Ky, ( % Kz =:= 0.0 Kz =< 1.0e-10, Kz >= -1.0e-10 -> add_linear_f1h(Xs,Ka,Ys,Zs) ; Zs = [l(X*Kz,OrdX)|Ztail], add_linear_f1h(Xs,Ka,Ys,Ztail) ) ; Rel = (<) -> Zs = [l(X*Kz,OrdX)|Ztail], Kz is Kx*Ka, add_linear_f1h(Xs,Ka,[l(Y*Ky,OrdY)|Ys],Ztail) ; Rel = (>) -> Zs = [l(Y*Ky,OrdY)|Ztail], add_linear_f1h(Ys,X,Kx,OrdX,Xs,Ztail,Ka) ). % add_linear_11(LinA,LinB,LinC) % % special case of add_linear_ff with Ka = 1 and Kb = 1 add_linear_11(LinA,LinB,LinC) :- LinA = [Ia,Ra|Ha], LinB = [Ib,Rb|Hb], LinC = [Ic,Rc|Hc], Ic is Ia+Ib, Rc is Ra+Rb, add_linear_11h(Ha,Hb,Hc). % add_linear_11h(Ha,Hb,Hc) % % special case of add_linear_ffh/5 with Ka = 1 and Kb = 1 add_linear_11h([],Ys,Ys). add_linear_11h([l(X*Kx,OrdX)|Xs],Ys,Zs) :- add_linear_11h(Ys,X,Kx,OrdX,Xs,Zs). % add_linear_11h(Ys,X,Kx,OrdX,Xs,Zs) % % special case of add_linear_ffh/8 with Ka = 1 and Kb = 1 add_linear_11h([],X,Kx,OrdX,Xs,[l(X*Kx,OrdX)|Xs]). add_linear_11h([l(Y*Ky,OrdY)|Ys],X,Kx,OrdX,Xs,Zs) :- compare(Rel,OrdX,OrdY), ( Rel = (=) -> Kz is Kx+Ky, ( % Kz =:= 0.0 Kz =< 1.0e-10, Kz >= -1.0e-10 -> add_linear_11h(Xs,Ys,Zs) ; Zs = [l(X*Kz,OrdX)|Ztail], add_linear_11h(Xs,Ys,Ztail) ) ; Rel = (<) -> Zs = [l(X*Kx,OrdX)|Ztail], add_linear_11h(Xs,Y,Ky,OrdY,Ys,Ztail) ; Rel = (>) -> Zs = [l(Y*Ky,OrdY)|Ztail], add_linear_11h(Ys,X,Kx,OrdX,Xs,Ztail) ). % mult_linear_factor(Lin,K,Res) % % Linear expression Res is the result of multiplication of linear % expression Lin by scalar K mult_linear_factor(Lin,K,Mult) :- TestK is K - 1.0, % K =:= 1 TestK =< 1.0e-10, TestK >= -1.0e-10, % avoid copy !, Mult = Lin. mult_linear_factor(Lin,K,Res) :- Lin = [I,R|Hom], Res = [Ik,Rk|Mult], Ik is I*K, Rk is R*K, mult_hom(Hom,K,Mult). % mult_hom(Hom,K,Res) % % Homogene part Res is the result of multiplication of homogene part % Hom by scalar K mult_hom([],_,[]). mult_hom([l(A*Fa,OrdA)|As],F,[l(A*Fan,OrdA)|Afs]) :- Fan is F*Fa, mult_hom(As,F,Afs). % nf_substitute(Ord,Def,Lin,Res) % % Linear expression Res is the result of substitution of Var in % linear expression Lin, by its definition in the form of linear % expression Def nf_substitute(OrdV,LinV,LinX,LinX1) :- delete_factor(OrdV,LinX,LinW,K), add_linear_f1(LinV,K,LinW,LinX1). % delete_factor(Ord,Lin,Res,Coeff) % % Linear expression Res is the result of the deletion of the term % Var*Coeff where Var has ordering Ord from linear expression Lin delete_factor(OrdV,Lin,Res,Coeff) :- Lin = [I,R|Hom], Res = [I,R|Hdel], delete_factor_hom(OrdV,Hom,Hdel,Coeff). % delete_factor_hom(Ord,Hom,Res,Coeff) % % Homogene part Res is the result of the deletion of the term % Var*Coeff from homogene part Hom delete_factor_hom(VOrd,[Car|Cdr],RCdr,RKoeff) :- Car = l(_*Koeff,Ord), compare(Rel,VOrd,Ord), ( Rel= (=) -> RCdr = Cdr, RKoeff=Koeff ; Rel= (>) -> RCdr = [Car|RCdr1], delete_factor_hom(VOrd,Cdr,RCdr1,RKoeff) ). % nf_coeff_of(Lin,OrdX,Coeff) % % Linear expression Lin contains the term l(X*Coeff,OrdX) nf_coeff_of([_,_|Hom],VOrd,Coeff) :- nf_coeff_hom(Hom,VOrd,Coeff). % nf_coeff_hom(Lin,OrdX,Coeff) % % Linear expression Lin contains the term l(X*Coeff,OrdX) where the % order attribute of X = OrdX nf_coeff_hom([l(_*K,OVar)|Vs],OVid,Coeff) :- compare(Rel,OVid,OVar), ( Rel = (=) -> Coeff = K ; Rel = (>) -> nf_coeff_hom(Vs,OVid,Coeff) ). % nf_rhs_x(Lin,OrdX,Rhs,K) % % Rhs = R + I where Lin = [I,R|Hom] and l(X*K,OrdX) is a term of Hom nf_rhs_x(Lin,OrdX,Rhs,K) :- Lin = [I,R|Tail], nf_coeff_hom(Tail,OrdX,K), Rhs is R+I. % late because X may not occur in H % isolate(OrdN,Lin,Lin1) % % Linear expression Lin1 is the result of the transformation of linear expression % Lin = 0 which contains the term l(New*K,OrdN) into an equivalent expression Lin1 = New. isolate(OrdN,Lin,Lin1) :- delete_factor(OrdN,Lin,Lin0,Coeff), K is -1.0/Coeff, mult_linear_factor(Lin0,K,Lin1). % indep(Lin,OrdX) % % succeeds if Lin = [0,_|[l(X*1,OrdX)]] indep(Lin,OrdX) :- Lin = [I,_|[l(_*K,OrdY)]], OrdX == OrdY, % K =:= 1.0 TestK is K - 1.0, TestK =< 1.0e-10, TestK >= -1.0e-10, % I =:= 0 I =< 1.0e-10, I >= -1.0e-10. % nf2sum(Lin,Sofar,Term) % % Transforms a linear expression into a sum % (e.g. the expression [5,_,[l(X*2,OrdX),l(Y*-1,OrdY)]] gets transformed into 5 + 2*X - Y) nf2sum([],I,I). nf2sum([X|Xs],I,Sum) :- ( % I =:= 0.0 I =< 1.0e-10, I >= -1.0e-10 -> X = l(Var*K,_), ( % K =:= 1.0 TestK is K - 1.0, TestK =< 1.0e-10, TestK >= -1.0e-10 -> hom2sum(Xs,Var,Sum) ; % K =:= -1.0 TestK is K + 1.0, TestK =< 1.0e-10, TestK >= -1.0e-10 -> hom2sum(Xs,-Var,Sum) ; hom2sum(Xs,K*Var,Sum) ) ; hom2sum([X|Xs],I,Sum) ). % hom2sum(Hom,Sofar,Term) % % Transforms a linear expression into a sum % this predicate handles all but the first term % (the first term does not need a concatenation symbol + or -) % see also nf2sum/3 hom2sum([],Term,Term). hom2sum([l(Var*K,_)|Cs],Sofar,Term) :- ( % K =:= 1.0 TestK is K - 1.0, TestK =< 1.0e-10, TestK >= -1.0e-10 -> Next = Sofar + Var ; % K =:= -1.0 TestK is K + 1.0, TestK =< 1.0e-10, TestK >= -1.0e-10 -> Next = Sofar - Var ; % K < 0.0 K < -1.0e-10 -> Ka is -K, Next = Sofar - Ka*Var ; Next = Sofar + K*Var ), hom2sum(Cs,Next,Term). %;itf /* Part of CLP(R) (Constraint Logic Programming over Reals) Author: Leslie De Koninck E-mail: Leslie.DeKoninck@cs.kuleuven.be WWW: http://www.swi-prolog.org http://www.ai.univie.ac.at/cgi-bin/tr-online?number+95-09 Copyright (C): 2004, K.U. Leuven and 1992-1995, Austrian Research Institute for Artificial Intelligence (OFAI), Vienna, Austria This software is part of Leslie De Koninck's master thesis, supervised by Bart Demoen and daily advisor Tom Schrijvers. It is based on CLP(Q,R) by Christian Holzbaur for SICStus Prolog and distributed under the license details below with permission from all mentioned authors. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA As a special exception, if you link this library with other files, compiled with a Free Software compiler, to produce an executable, this library does not by itself cause the resulting executable to be covered by the GNU General Public License. This exception does not however invalidate any other reasons why the executable file might be covered by the GNU General Public License. */ do_checks(Y,Ty,St,Li,Or,Cl,No,Later) :- numbers_only(Y), verify_nonzero(No,Y), verify_type(Ty,St,Y,Later,[]), verify_lin(Or,Cl,Li,Y), maplist(call,Later). numbers_only(Y) :- ( var(Y) -> true ; integer(Y) -> true ; float(Y) -> true ; type_error(real, Y) ). % verify_nonzero(Nonzero,Y) % % if Nonzero = nonzero, then verify that Y is not zero % (if possible, otherwise set Y to be nonzero) verify_nonzero(nonzero,Y) :- ( var(Y) -> ( get_attr(Y,itf,Att) -> setarg(8,Att,nonzero) ; put_attr(Y,itf,t(clpr,n,n,n,n,n,n,nonzero,n,n,n)) ) ; ( Y < -1.0e-10 -> true ; Y > 1.0e-10 ) ). verify_nonzero(n,_). % X is not nonzero % verify_type(type(Type),strictness(Strict),Y,[OL|OLT],OLT) % % if possible verifies whether Y satisfies the type and strictness of X % if not possible to verify, then returns the constraints that follow from % the type and strictness verify_type(type(Type),strictness(Strict),Y) --> verify_type2(Y,Type,Strict). verify_type(n,n,_) --> []. verify_type2(Y,TypeX,StrictX) --> {var(Y)}, !, verify_type_var(TypeX,Y,StrictX). verify_type2(Y,TypeX,StrictX) --> {verify_type_nonvar(TypeX,Y,StrictX)}. % verify_type_nonvar(Type,Nonvar,Strictness) % % verifies whether the type and strictness are satisfied with the Nonvar verify_type_nonvar(t_none,_,_). verify_type_nonvar(t_l(L),Value,S) :- ilb(S,L,Value). verify_type_nonvar(t_u(U),Value,S) :- iub(S,U,Value). verify_type_nonvar(t_lu(L,U),Value,S) :- ilb(S,L,Value), iub(S,U,Value). verify_type_nonvar(t_L(L),Value,S) :- ilb(S,L,Value). verify_type_nonvar(t_U(U),Value,S) :- iub(S,U,Value). verify_type_nonvar(t_Lu(L,U),Value,S) :- ilb(S,L,Value), iub(S,U,Value). verify_type_nonvar(t_lU(L,U),Value,S) :- ilb(S,L,Value), iub(S,U,Value). % ilb(Strict,Lower,Value) & iub(Strict,Upper,Value) % % check whether Value is satisfiable with the given lower/upper bound and % strictness. % strictness is encoded as follows: % 2 = strict lower bound % 1 = strict upper bound % 3 = strict lower and upper bound % 0 = no strict bounds ilb(S,L,V) :- S /\ 2 =:= 0, !, L - V < 1.0e-10. % non-strict ilb(_,L,V) :- L - V < -1.0e-10. % strict iub(S,U,V) :- S /\ 1 =:= 0, !, V - U < 1.0e-10. % non-strict iub(_,U,V) :- V - U < -1.0e-10. % strict % % Running some goals after X=Y simplifies the coding. It should be possible % to run the goals here and taking care not to put_atts/2 on X ... % % verify_type_var(Type,Var,Strictness,[OutList|OutListTail],OutListTail) % % returns the inequalities following from a type and strictness satisfaction % test with Var verify_type_var(t_none,_,_) --> []. verify_type_var(t_l(L),Y,S) --> llb(S,L,Y). verify_type_var(t_u(U),Y,S) --> lub(S,U,Y). verify_type_var(t_lu(L,U),Y,S) --> llb(S,L,Y), lub(S,U,Y). verify_type_var(t_L(L),Y,S) --> llb(S,L,Y). verify_type_var(t_U(U),Y,S) --> lub(S,U,Y). verify_type_var(t_Lu(L,U),Y,S) --> llb(S,L,Y), lub(S,U,Y). verify_type_var(t_lU(L,U),Y,S) --> llb(S,L,Y), lub(S,U,Y). % llb(Strict,Lower,Value,[OL|OLT],OLT) and lub(Strict,Upper,Value,[OL|OLT],OLT) % % returns the inequalities following from the lower and upper bounds and the % strictness see also lb and ub llb(S,L,V) --> {S /\ 2 =:= 0}, !, [clpr:{L =< V}]. llb(_,L,V) --> [clpr:{L < V}]. lub(S,U,V) --> {S /\ 1 =:= 0}, !, [clpr:{V =< U}]. lub(_,U,V) --> [clpr:{V < U}]. % % We used to drop X from the class/basis to avoid trouble with subsequent % put_atts/2 on X. Now we could let these dead but harmless updates happen. % In R however, exported bindings might conflict, e.g. 0 \== 0.0 % % If X is indep and we do _not_ solve for it, we are in deep shit % because the ordering is violated. % verify_lin(order(OrdX),class(Class),lin(LinX),Y) :- !, ( indep(LinX,OrdX) -> detach_bounds_vlv(OrdX,LinX,Class,Y,NewLinX), % if there were bounds, they are requeued already class_drop(Class,Y), nf(-Y,NfY), deref(NfY,LinY), add_linear_11(NewLinX,LinY,Lind), ( nf_coeff_of(Lind,OrdX,_) -> % X is element of Lind solve_ord_x(Lind,OrdX,Class) ; solve(Lind) % X is gone, can safely solve Lind ) ; class_drop(Class,Y), nf(-Y,NfY), deref(NfY,LinY), add_linear_11(LinX,LinY,Lind), solve(Lind) ). verify_lin(_,_,_,_).