data_of_this_program([
  line('%-------------------------------------------------------%'),
  title(' elimination by aspects (EBA) model by Tversky(1972).'),
  line('%-------------------------------------------------------%'),
  program_name('eba01.pl'),
  created(by('Kenryo INDO <Kanto Gakuen University>')),
  created(date('21-23 May 2003')),
  modified(date('25-26,28 May 2003')),
  privious_program('eba0.pl(14 May 2003)'),
  reference(1,author('Tversky, A.(1972)')),
  reference(1,title('Elimination by aspects')),
  reference(1,journal('Psychological Review 79(4): 281-299.')),
  reference(2,author('Tversky, A.(1972)')),
  reference(2,title('Choice by elimination')),
  reference(2,journal('J Mathematical Psychology 9: 341-367.')),
  line('%-------------------------------------------------------%'),
  main_predicates(aspects/3),
  main_predicates(utility/3),
  main_predicates(eliminate_by_aspects/4),
  main_predicates(iterative_elimination/4)
]).

:- data_of_this_program(X),forall(member(T,X),(nl,write(T))).

:- dynamic mc/3.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    basic data of the alternatives
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

alternatives(model(1),[x,y,z]).
alternatives(model(2),[x,y,z]).
alternatives(model(3),[x,y,z]).

aspects(model(1),base,x,[a,w]).
aspects(model(1),base,y,[b,w]).
aspects(model(1),base,z,[c]).

aspects(model(2),base,x,[a,t,r,w]).
aspects(model(2),base,y,[b,t,s,w]).
aspects(model(2),base,z,[c,r,s,w]).

aspects(model(3),base,x,[b,c,w]).
aspects(model(3),base,y,[b,c,w]).
aspects(model(3),base,z,[a,w]).

aspects(M,subset,S,Ax):-
   alternatives(M,Alts),
   subset_of(S,_,Alts),
   S \= [],
   findall(A,
     (
      member(E,S),
      aspects(M,base,E,Asp),
      member(A,Asp)
     ),
   Ax0),
   sort(Ax0,Ax).
aspects(M,uniqueness,X,Ix):-
   aspects(M,base,X,Ax),
   alternatives(M,Alts),
   subtract(Alts,[X],OTX),
   aspects(M,subset,OTX,Bx),
   subtract(Ax,Bx,Ix).
aspects(M,shared(base),[X,Y],Axy):-
   aspects(M,base,X,Ax),
   aspects(M,base,Y,Ay),
   %Y \= X,
   intersection(Ax,Ay,Axy).
aspects(M,shared(subset),[X,Y],Axy):-
   aspects(M,subset,X,Ax),
   aspects(M,subset,Y,Ay),
   %Y \= X,
   intersection(Ax,Ay,Axy).
aspects(M,difference,S-T,Dst):-
   aspects(M,subset,S,As),
   aspects(M,subset,T,At),
   subtract(As,At,Dst).
aspects(M,vector,S,P):-
   alternatives(M,Alts),
   aspects(M,subset,Alts,Asp0),
   aspects(M,subset,S,Ax),
   list_projection(P,Asp0,Ax).


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    utilities and probabilistic choice
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%----utility for each aspect.

utility(model(1),aspect(a),1).
utility(model(1),aspect(b),1).
utility(model(1),aspect(c),2).
utility(model(1),aspect(w),1).

utility(model(2),aspect(a),1).
utility(model(2),aspect(b),3).
utility(model(2),aspect(c),5).
utility(model(2),aspect(r),2).
utility(model(2),aspect(s),3).
utility(model(2),aspect(t),5).
utility(model(2),aspect(w),1).

utility(model(3),aspect(a),1).
utility(model(3),aspect(b),1).
utility(model(3),aspect(c),2).
utility(model(3),aspect(w),3).

%----utility for a set of alternatives.

utility(M,S=Asp,U=Uq):-
   aspects(M,subset,S,Asp),
   findall(B,
     (
      member(A,Asp),
      utility(M,aspect(A),B)
     ),
   U0),
   sum_eq(U0,Uq,U).

%----utility for a set of aspects.

utility(M,aspects(S),U=Uq):-
   alternatives(M,ALL),
   aspects(M,subset,ALL,Asp),
   subset_of(S,_,Asp),
   findall(UA,
     (
      member(A,S),
      utility(M,aspect(A),UA)
     ),
   US),
   sum_eq(US,Uq,U).


%----probabilistic choice from a set of alternatives.
% note: this not always eliminate the complement.

probabilistic_choice(M,S->T,P=(Pq1=(Pq=u(C)/u(A)))):-
   utility(M,S=A, U0=Uq0),
   U0 \= 0,
   subset_of(T,_,S),
   T \= [],
   aspects(M,shared(subset),[S,T],C),
   utility(M,aspects(C), U1=Uq1),
   Pq1 = U1/U0,
   Pq = Uq1/Uq0,
   P is Pq.


%----probabilistic choice for the set of aspects.
probabilistic_choice(M,S->T,P=(Pq=u(T)/u(S))):-
   utility(M,aspects(S),_),
   findall(US,
     (
      subset_of(T0,_,S),
      utility(M,aspects(T0),_=US)
     ),
   USS),
   sum_eq(USS,Uq0,_),
   subset_of(T,_,S),
   T\=[],
   T\=S,
   utility(M,aspects(T),_=Uq1),
   Pq =Uq1/Uq0,
   P is Pq.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    elimination procedure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

eliminate_by_aspects(Model,S->(T:Deleted),F):-
   aspects(Model,subset,S,AS),
   subset_of(F,_,AS),  F \=[],
   subset_of(Deleted,_,S), Deleted \=[],
   aspects(Model,subset,Deleted,AD),
   intersection(AD,F,[]),
   subtract(S,Deleted,T),
   \+ (
     member(X,T),
     aspects(Model,base,X,AX),
     \+ subset(F,AX)
   ).

eliminate_by_aspects(M,S->(T:D),F,P=Formulas):-
   eliminate_by_aspects(M,S->(T:D),F),
   utility(M,aspects(F), U1=U1eq),
   find_all_alternative_ebas(M,S,EBAs),
   %display_each_eba(EBAs),
   sum_utility_of_aspects(EBAs,_,U0sEq,U0),
   sum_equations_of_utilities(EBAs,_,UF0),
   P is U1eq / U0sEq,
   Input =[U1,U0,U1eq,U0sEq,F,UF0],
   equation_of_choice_probability(Input,Formulas).


/*
%------stock of the modeling error--------
eliminate_by_aspects(Model,S->(T:Deleted),F):-
   aspects(Model,difference,S-Deleted,FocalSet),
   subset_of(Deleted,_,S),
   subtract(S,Deleted,T),
   subset_of(F,_,FocalSet),
   F \=[],
   aspects(Model,shared(subset),[S,T],Common),
   Common \= [],
   subset_of(F,_,Common).
eliminate_by_aspects(Model,S->(T:Deleted),F):-
   aspects(Model,subset,S,AS),
   subset_of(T,_,S), T\=[],
   aspects(Model,subset,T,AT),
   subtract(AS,AT,FocalSet), FocalSet \=[],
   subset_of(F,_,FocalSet), F \=[],
   subtract(S,T,Deleted).
eliminate_by_aspects(Model,S->(T:Deleted),F):-
   aspects(Model,subset,S,AS),
   subset_of(T,_,S), T\=[],
   aspects(Model,subset,T,AT),
   subtract(AS,AT,FocalSet), FocalSet \=[],
   subset_of(F,_,FocalSet), F \=[],
   subtract(S,T,Deleted),
   aspects(Model,shared(subset),[T,Deleted],[]).
*/

find_all_alternative_ebas(M,S,EBAs):-
   EBA = eba(M,S->(T0:D0),F0,U0eq=u(F0)),
   findall(EBA,
     (
      eliminate_by_aspects(M,S->(T0:D0),F0),
      utility(M,aspects(F0),_=U0eq)
     ),
   EBAs).

display_each_eba(EBAs):-
   forall(member(X,EBAs),(nl,write(X))).

sum_utility_of_aspects(EBAs,U0s,U0sEq,U0):-
   X = eba(_,_,_,U0eq=_),
   findall(U0eq,member(X,EBAs),U0s),
   sum_eq(U0s,U0sEq,U0).

sum_equations_of_utilities(EBAs,UF0s,UF0):-
   X = eba(_,_,F0,_=u(F0)),
   findall(u(F0),member(X,EBAs),UF0s),
   eqsum(UF0s,UF0).

equation_of_choice_probability(Input,Formulas):-
   Input =[U1,U0,U1eq,U0sEq,F,UF0],
   Formula1= U1/U0,
   Formula2= U1eq/U0sEq,
   Formula3= u(F)/UF0,
   Formulas= (Formula1=(Formula2=Formula3)).


/*
%-------sample executions----------

?- eliminate_by_aspects(model(1),[x,y,z]->A,B,C=(D=(E=F))).

A = [x, y]:[z]
B = [w]
C = 0.2
D = 1/5
E = 1/ (1+1+2+1)
F = u([w])/ (u([a])+u([b])+u([c])+u([w])) ;

A = [z]:[x, y]
B = [c]
C = 0.4
D = 2/5
E = 2/ (1+1+2+1)
F = u([c])/ (u([a])+u([b])+u([c])+u([w])) ;

A = [y]:[x, z]
B = [b]
C = 0.2
D = 1/5
E = 1/ (1+1+2+1)
F = u([b])/ (u([a])+u([b])+u([c])+u([w])) ;

A = [x]:[y, z]
B = [a]
C = 0.2
D = 1/5
E = 1/ (1+1+2+1)
F = u([a])/ (u([a])+u([b])+u([c])+u([w])) ;

No
?- 
?- eliminate_by_aspects(model(2),A,B,C=(C1=(C2=C3))).

eba(model(1), ([y, z]->[y]:[z]), [c], 2=u([c]))
eba(model(1), ([y, z]->[y]:[z]), [b], 1=u([b]))
eba(model(1), ([y, z]->[y]:[z]), [b, c], 2+1=u([b, c]))
eba(model(1), ([y, z]->[z]:[y]), [a], 1=u([a]))

A = [y, z]->[y]:[z]
B = [c]
C = 0.285714
C1 = 2/7
C2 = 2/ (1+ (2+1)+1+2)
C3 = u([c])/ (u([a])+u([b, c])+u([b])+u([c])) ;

<--omit-->

?- eliminate_by_aspects(model(1),[x,y,z]->[z]:[x,y],B,C=(D=(E=F))).

eba(model(1), ([x, y, z]->[x, y]:[z]), [c], 2=u([c]))
eba(model(1), ([x, y, z]->[x, y]:[z]), [b], 1=u([b]))
eba(model(1), ([x, y, z]->[x, y]:[z]), [b, c], 2+1=u([b, c]))
eba(model(1), ([x, y, z]->[z]:[x, y]), [a], 1=u([a]))

B = [a]
C = 0.142857
D = 1/7
E = 1/ (1+ (2+1)+1+2)
F = u([a])/ (u([a])+u([b, c])+u([b])+u([c])) ;

No
?- eliminate_by_aspects(model(2),[x,y,z]->A,B,C=(D=(E=F))).

A = [x, y]:[z]
B = [c]
C = 0.285714
D = 2/7
E = 2/ (1+ (2+1)+1+2)
F = u([c])/ (u([a])+u([b, c])+u([b])+u([c])) ;

A = [x, y]:[z]
B = [b]
C = 0.142857
D = 1/7
E = 1/ (1+ (2+1)+1+2)
F = u([b])/ (u([a])+u([b, c])+u([b])+u([c])) ;

A = [x, y]:[z]
B = [b, c]
C = 0.428571
D = 3/7
E = (2+1)/ (1+ (2+1)+1+2)
F = u([b, c])/ (u([a])+u([b, c])+u([b])+u([c])) ;

A = [z]:[x, y]
B = [a]
C = 0.142857
D = 1/7
E = 1/ (1+ (2+1)+1+2)
F = u([a])/ (u([a])+u([b, c])+u([b])+u([c])) ;

No
?-
%-------end of sample executions----------
*/

:- dynamic ieba /5.

iterative_elimination_0(M,(S->S),([],[Asp],[1]),1=(1=1)):-
   aspects(M,subset,S,Asp).

iterative_elimination_0(M,(S->T),Paths,P1=(P1eq=P1eq2)):-
   aspects(M,subset,S,_AS),
   subset_of(T,_,S),T\=[],T\=S,
   aspects(M,subset,T,_AT),
   eliminate_by_aspects(M,T1->(T:D),F,P=Peq),
   iterative_elimination_0(M,(S->T1),Paths0,_P0=(P0eq=P0eq2)),
   Paths0=(DelPath,AspPath,PrPath),
   Paths=([D|DelPath],[F|AspPath],[P|PrPath]),
   P1eq2 = P0eq2 * Peq,
   P1eq = P0eq * P,
   P1 is P1eq.


/*
?- iterative_elimination_0(model(2),[x,y,z]->A,B,C=_).

A = [x, y, z]
B = [], [[a, b, c, w]], [1]
C = 1 ;

A = [z]
B = [[x, y]], [[a], [a, b, c, w]], [0.142857, 1]
C = 0.142857 ;

A = [x, y]
B = [[z]], [[c], [a, b, c, w]], [0.285714, 1]
C = 0.285714 ;

A = [x, y]
B = [[z]], [[b], [a, b, c, w]], [0.142857, 1]
C = 0.142857 ;

A = [x, y]
B = [[z]], [[b, c], [a, b, c, w]], [0.428571, 1]
C = 0.428571 ;

No
?- 

*/

iterative_elimination(M,(S->S),([],[Asp],[1]),1=(1)):-
   aspects(M,subset,S,Asp).

iterative_elimination(M,(S->T),Paths,P1=(P1eq)):-
   aspects(M,subset,S,_AS),
   subset_of(T,_,S),T\=[],T\=S,
   aspects(M,subset,T,_AT),
   R = eba(M,T1->(T:D),F,P=Peq),
   findall(R,
     eliminate_by_aspects(M,T1->(T:D),F,P=Peq),
   EBA0),
   W = ieba(M,(S->T),NewPaths,P=Peq,P0=(P0eq)),
   findall(W,
     (
      member(R,EBA0),
      iterative_elimination_0(M,(S->T1),Paths0,P0=(P0eq)),
      Paths0=(DelPath,AspPath,PrPath),
      NewPaths=([D|DelPath],[F|AspPath],[P|PrPath])
     ),
   EBA),
   EBA \=[],
   forall(member(W,EBA),(nl,write(W),write('.')%,assert(W)
)),
   findall(P*P0,member(W,EBA),Ps),
   findall(NewPaths,member(W,EBA),Paths),
   sum_eq(Ps,P1eq,P1).


/*

?- iterative_elimination(model(3),[x,y,z]->A,B,C=_).

A = [x, y, z]
B = [], [[a, b, c, w]], [1]
C = 1 ;

ieba(model(1), ([x, y, z]->[z]), ([[x, y]], [[a], [a, b, c, w]], [0.142857, 1]), 0.142857= (1/7= (1/ (1+ (2+1)+1+2)=u([a])/ (u([a])+u([b, c])+u([b])+u([c])))), 1= (1=1)).

A = [z]
B = [ ([[x, y]], [[a], [a, b, c, w]], [0.142857, 1])]
C = 0.142857*1 ;

ieba(model(1), ([x, y, z]->[x, y]), ([[z]], [[c], [a, b, c, w]], [0.285714, 1]), 0.285714= (2/7= (2/ (1+ (2+1)+1+2)=u([c])/ (u([a])+u([b, c])+u([b])+u([c])))), 1= (1=1)).
ieba(model(1), ([x, y, z]->[x, y]), ([[z]], [[b], [a, b, c, w]], [0.142857, 1]), 0.142857= (1/7= (1/ (1+ (2+1)+1+2)=u([b])/ (u([a])+u([b, c])+u([b])+u([c])))), 1= (1=1)).
ieba(model(1), ([x, y, z]->[x, y]), ([[z]], [[b, c], [a, b, c, w]], [0.428571, 1]), 0.428571= (3/7= ((2+1)/ (1+ (2+1)+1+2)=u([b, c])/ (u([a])+u([b, c])+u([b])+u([c])))), 1= (1=1)).

A = [x, y]
B = [ ([[z]], [[c], [a, b, c, w]], [0.285714, 1]), ([[z]], [[b], [a, b, c|...]], [0.142857, 1]), ([[z]], [[b, c], [a, b|...]], [0.428571, 1])]
C = 0.857143 ;

No
?-
?- tell_goal( 'out_eba.txt',
iterative_elimination(model(2),[x,y,z]->A,[B|B1],(C=(D=E)))
).

A = _G169
B = _G171
B1 = _G172
C = _G174
C1 = _G175 

Yes
?- ['out_eba.txt'].

Yes.
?- iterative_elimination(model(2),[x,y,z]->[A],[B|B1],C=C1).

A = z
B = [[y], [x]], [[r], [s], [a, b, c, r|...]], [0.0666667, 0.0247934, 1]
B1 = [ ([[y], [x]], [[r], [c], [a, b, c|...]], [0.0666667, 0.0413223, 1]), ([[y], [x]], [[r], [c, s], [a, b|...]], [0.0666667, 0.0661157, 1]), ([[y], [x]], [[r], [b], [a|...]], [0.0666667, 0.0247934, 1]), ([[y], [x]], [[r], [b|...], [...|...]], [0.0666667, 0.0495868, 1]), ([[y], [x]], [[r], [...|...]|...], [0.0666667, 0.0661157|...]), ([[y], [...]], [[...]|...], [0.0666667|...]), ([[...]|...], [...|...], [...|...]), ([...|...], ..., ...), (..., ...)|...]
C = 0.362141
C1 = ... +... +... *... +0.166667*0.0247934+0.0666667*0.0909091+0.0666667*0.0661157+0.0666667*0.0495868+0.0666667*0.0247934+0.0666667*0.0661157+0.0666667*0.0413223+0.0666667*0.0247934 ;

% and so on.


*/

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    simulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ----- iterative elimination with mc trials

init_mc(M):-
   retractall(mc(M,_,_)).

eba_process(M,E_step,prob([1])):-
   init_mc(M),
   alternatives(M,S),
   E_step = [(S,S),[],[Asp]],
   iterative_elimination(M,(S,S),[],[Asp]),
   MC = mc(M,[last(0)|E_step],prob(1/1,[])),
   assert(MC).

eba_process(M,E_step,prob([P|Q])):-
   E_step = [(S,T),[_D|Path],[_F|AspPath]],
   E_step0 = [(S,T0),Path,AspPath],
   MC0 = mc(M,[last(N0)|E_step0],prob(_,Q)),
   MC0,
   N is N0 + 1,
   make_mc_data(T0->T,N,MC0),
   eba_selection(M,N,E_step,P).

% ----- monte carlo simulation

make_mc_data(T0->T,N,MC0):-
   MC0 = mc(M,[last(_)|E_step0],prob(P0/_,Q)),
   E_step0 = [(S,T0),Path,AspPath],
   E_step = [(S,T),[D|Path],[F|AspPath]],
   MC = mc(M,[last(N)|E_step],prob(P/_,[P0|Q])),
   EBA = eliminate_by_aspects(M,T0->(T:D),F),
   PROB = probabilistic_choice(M,T0->T,P=_),
   forall(
     (
      EBA,nl,write(EBA),
      T \= T0,
      PROB,
      nl,write(PROB)
     ),
     (
      update_mc(MC)
     )
   ).

update_mc(MC):-
   MC = mc(M,[last(N)|_],prob(P/CP,Q)),
   MC0 = mc(M,[last(N)|E_step0],prob(P0/CP0,Q)),
   clause(MC0,true),
   retract(MC0),
   % nl,write(retract:MC0),
   MC1 = mc(M,[N|E_step0],prob(P0/CP0,Q)),
   assert(MC1),
   % nl,write(assert:MC1),
   CP is CP0 + P,
   assert(MC),
   nl,write(new:MC).

update_mc(MC):-
   MC = mc(M,[last(N)|_],prob(P/P,Q)),
   MC0 = mc(M,[last(N)|_],prob(_,Q)),
   \+ clause(MC0,true),
   assert(MC),
   nl,write(new:MC).

eba_selection(M,N,E_step,P):-
   member(X,[best(N),last(N)]),
   mc(M,[X|E_step],prob(P/_CP,_Q)).




%------------------------------------------------
%  mathematical utiliteis
%------------------------------------------------
max_of(X,[X]).
max_of(Z,[X|Y]):-
   max_of(Z1,Y),
   (X > Z1 -> Z=X; Z=Z1).
min_of(X,[X]).
min_of(Z,[X|Y]):-
   min_of(Z1,Y),
   (X < Z1 -> Z=X; Z=Z1).
dnum_seq([],N):-N<0,!.
dnum_seq([0],1).
dnum_seq([A|Q],N):-
   A is N - 1,
   length(Q,A),
   dnum_seq(Q,A).
anum_seq(Aseq,N):-dnum_seq(Dseq,N),sort(Dseq,Aseq).

dnum_seq1(Q,N):-
   M is N + 1,
   dnum_seq(Q0,M),
   subtract(Q0,[0],Q).
anum_seq1(Q,N):-
   M is N + 1,
   anum_seq(Q0,M),
   subtract(Q0,[0],Q).

% bag0/3 : allow multiplicity
% -----------------------------------------------------------  %
bag0([],_A,0).
bag0([C|B],A,N):-length([C|B],N),bag0(B,A,_N1),%N is N1 + 1,
  member(C,A).

zeros(Zero,N):-bag0(Zero,[0],N).
ones(One,N):-bag0(One,[1],N).

%
% bag1/3 : do not allow multiplicity
% -----------------------------------------------------------  %
bag1([],_A,0).
bag1([C|B],A,N1):-
  \+var(A),
  length(A,L),
  anum_seq(Q,L),
  member(N,Q),
  length(B,N),bag1(B,A,N),N1 is N + 1,
  member(C,A),\+member(C,B).

% ordering/3
% -----------------------------------------------------------  %
ordering(A,B,C):-bag1(A,B,C).

% a sequence of binary choice for a list:
%--------------------------------------------------
list_projection([],[],[]).
list_projection([X|Y],[_|B],C):-
   list_projection(Y,B,C),
   X = 0.
list_projection([X|Y],[A|B],[A|C]):-
   list_projection(Y,B,C),
   X = 1.

% complementary list projection
%--------------------------------------------------
c_list_projection([],[],[]).
c_list_projection([X|Y],[_|B],C):-
   c_list_projection(Y,B,C),
   X = 1.
c_list_projection([X|Y],[A|B],[A|C]):-
   c_list_projection(Y,B,C),
   X = 0.

subset_of(A,N,As):-
   length(As,L),
   length(D,L),
   list_projection(D,As,B),
   length(B,N),
   sort(B,A).

sum([],0).
sum([X|Members],Sum):-
   sum(Members,Sum1),
   Sum is Sum1 + X.

sum_eq([],0,0).
sum_eq([X],X,X).
sum_eq([X|Members],Eq,Sum):-
   Members \= [],
   sum_eq(Members,Eq1,Sum1),
   Eq = Eq1 + X,
   Sum is Sum1 + X.

reqsum(A,B):-
   reverse(A,C),
   eqsum(C,B).
eqsum([],0).
eqsum([X|Members],Sum):-
   eqsum(Members,Sum0),
   (
    X=0 -> Sum = Sum0
     ;
    (
     Sum0=0 -> Sum = X
      ;
     Sum = Sum0 + X
    )
   ).

%
% product
% -----------------------------------------------------------  %
product([],1).
product([X|Members],Z):-
   product(Members,Z1),
  %number(X),
   Z is Z1 * X.
%
% weighted sum
% -----------------------------------------------------------  %
product_sum([],[],[],0).
product_sum([P|Q],[A|B],[E|F],V):-
    length(Q,N),
    length(B,N),
    product_sum(Q,B,F,V1),
    E is P * A,
    V is V1 + E.
%
% product sum value with equational
% -----------------------------------------------------------  %
product_sum_eq([],[],[],0,0).
product_sum_eq([P|Q],[A|B],[E|F],V,Vq):-
    length(Q,N),
    length(B,N),
    product_sum_eq(Q,B,F,V1,Vq1),
    Eq = (P) * A,
    E is Eq,
    (Vq1=0 -> Vq = Eq; Vq = Vq1 + Eq),
    V is V1 + E.


% tell a successful goal
%--------------------------------

tell_goal(File,G):-
   (current_stream(File,write,S0)->close(S0);true),
   open(File,write,S),
   tell(File),
   nl,
   tstamp('% file output start time ',_),
   nl,
   write('%----------  start from here ------------%'),
   nl,
   G,
   nl,
   write('%----------  end of data ------------%'),
   nl,
   tstamp('% file output end time ',_),
   tell(user),
   close(S),
   % The following is to cope with the duplicated stream problem.
   (current_stream(File,write,S1)->close(S1);true).


% tell all successful goal
%--------------------------------

tell_goal(File,forall,G):-
   G0 = (nl,write(G),write('.')),
   G1 = forall(G,G0),
   tell_goal(File,G1).


% time stamp
%--------------------------------

tstamp(no,T):-
   get_time(U),
   convert_time(U,A,B,C,D,E,F,_G),
   T = [date(A/B/C), time(D:E:F)],
   nl.

tstamp(Word,T):-
   \+ var(Word),
   Word \= no,
   get_time(U),
   convert_time(U,A,B,C,D,E,F,_G),
   T = [date(A/B/C), time(D:E:F)],
%   format('~`.t~t~a~30|~`.t~t~70|~n',[Word]),
   write((Word,T)),
   nl.




%----------end of pogram ----------%


/*

%-----------tombeau of bugs -----------------

init_mc(M):-
   retractall(mc(M,_,_)).

eba_process(M,E_step,prob([1])):-
   init_mc(M),
   alternatives(M,S),
   E_step = [(S,S),[],[Asp]],
   iterative_elimination(M,(S,S),[],[Asp]),
   MC = mc(M,[last(0)|E_step],prob(1/1,[])),
   assert(MC).

eba_process(M,E_step,prob([P|Q])):-
   E_step = [(S,T),[_D|Path],[_F|AspPath]],
   E_step0 = [(S,T0),Path,AspPath],
   MC0 = mc(M,[last(N0)|E_step0],prob(_,Q)),
   MC0,
   N is N0 + 1,
   make_mc_data(T0->T,N,MC0),
   eba_selection(M,N,E_step,P).

% ----- monte carlo simulation

make_mc_data(T0->T,N,MC0):-
   MC0 = mc(M,[last(_)|E_step0],prob(P0/_,Q)),
   E_step0 = [(S,T0),Path,AspPath],
   E_step = [(S,T),[D|Path],[F|AspPath]],
   MC = mc(M,[last(N)|E_step],prob(P/_,[P0|Q])),
   EBA = eliminate_by_aspects(M,T0->(T:D),F),
   PROB = probabilistic_choice(M,T0->T,P=_),
   forall(
     (
      EBA,nl,write(EBA),
      T \= T0,
      PROB,
      nl,write(PROB)
     ),
     (
      update_mc(MC)
     )
   ).

update_mc(MC):-
   MC = mc(M,[last(N)|_],prob(P/CP,Q)),
   MC0 = mc(M,[last(N)|E_step0],prob(P0/CP0,Q)),
   clause(MC0,true),
   retract(MC0),
   % nl,write(retract:MC0),
   MC1 = mc(M,[N|E_step0],prob(P0/CP0,Q)),
   assert(MC1),
   % nl,write(assert:MC1),
   CP is CP0 + P,
   assert(MC),
   nl,write(new:MC).

update_mc(MC):-
   MC = mc(M,[last(N)|_],prob(P/P,Q)),
   MC0 = mc(M,[last(N)|_],prob(_,Q)),
   \+ clause(MC0,true),
   assert(MC),
   nl,write(new:MC).

eba_selection(M,N,E_step,P):-
   member(X,[best(N),last(N)]),
   mc(M,[X|E_step],prob(P/_CP,_Q)).

iterative_elimination_0(M,(S,S),([],[Asp],[1]),1=(1=1)):-
   aspects(M,subset,S,Asp).

iterative_elimination_0(M,(S,T),Paths,Q=((Q0*P)=Qeq)):-
   subset_of(T,_,S),T\=[],T\=S,
   aspects(M,subset,T,_AT),
   iterative_elimination_0(M,(S,T1),(DelPath,AspPath,PrPath),Q0=(_=Q0eq)),
   Paths=([D|DelPath],[F|AspPath],[P|PrPath]),
   eliminate_by_aspects(M,T1->(T:D),F,P=Peq),
   Qeq= Q0eq * Peq,
   Q is Qeq.

*/