You selected matrix01.pl

% -----------------------------------------------------------  %
%  Basics of Matrix and Linear Algebra by Prolog
% -----------------------------------------------------------  %
%  file: matrix1b.pl
%  imported:  (math1.pl version: 21 Feb 2003)
%  edited: 25-26 Feb 2003.(matrix1.pl)
%  modified: 18 Apr 2003.  case_of_diagonal(4,Pxx,_,_,_,_).
%  modified: 8-9 Jun 2003.  symbolic equations.
%  modified: 20 Aug 2003.  matrix_is with modulo. 
%  modified: 20 Aug 2003.  the program has rewrited with 
%   some new predicates len/2, kth/3, and rev/2 each of which 
%   is the alternative of SWI system predicates length/2, 
%   nth1/3, and reverse/2 respectively. 


:- dynamic matrix /2.
:- dynamic cell /3.


% -----------------------------------------------------------  %
%  definitions of matrix representation
% -----------------------------------------------------------  %
%
cell(1,1,1).
cell(1,2,2).
cell(2,1,3).
cell(2,2,4).
matrix(i(2),[[1,0],[0,1]]).
matrix(i(3),[[1,0,0],[0,1,0],[0,0,1]]).
matrix(ma,[[0,-1],[1,0]]).
matrix(mb,[[1,-1],[2,0]]).
matrix(mc,[[3,-6,-2],[1,-2,-1],[-2,6,3]]).
matrix(p,[[1,3,2],[0,1,1],[1,0,-2]]).
matrix(q,[[2,-6,-1],[-1,4,1],[1,-3,-1]]).
matrix(r,[[1,-1,-1],[-1,1,1],[1,-1,-1]]).
matrix(s,[[1,1,1],[2,1,-1],[3,-1,1]]).
matrix(t,[[1,1,1],[1,1,1],[3,-1,1]]).
matrix(n2sq,[[a,b],[c,d]]).
matrix(x,[[x,z],[y,w]]).

% identity matrix
matrix(i(N),I):-
   integer(N),
   N > 3,
   matrix_of_size(I,N,N),
   findall(R,
     (
      kth(K,I,R),
      findall(C,
        (
         kth(L,R,C),
         (L = K -> C = 1; C = 0)
        ),
      R)
     ),
   I),!.

% scalar multiplication
% modified: 8 June 2003. symbolic case.
matrix(C * X,B):-
   %eval_number(C,_),
   matrix(X,A),
   matrix_of_size(A,N,M),
   matrix_of_size(B,N,M),
   findall(Br,
     (
      kth(K,A,Ar),
      kth(K,B,Br),
      findall(Brc,
        (
         kth(L,Ar,Arc),
         kth(L,Br,Brc),
         evaluate_if_numerical(Brc is C * Arc)
        ),
      Br)
     ),
   B),
   !.

% zero matrix
matrix(o(N,M),O):-
   matrix(z(N,M,0),O).

% uniform bias
matrix(z(N,M,C),O):-
   \+ C == cell,
   matrix_of_size(O,N,M),
   findall(Ok,
     (
      kth(_K,O,Ok),
      findall(C,
        (
         kth(_L,Ok,C)
        ),
      Ok)
     ),
   O),!.

% indirect definition via cells.
matrix(z(N,M,cell),O):-
   X =cell(_,_,_),
   clause(X,_),
   matrix_of_size(O,N,M),
   findall(Ok,
     (
      kth(K,O,Ok),
      findall(C,
        (
         kth(L,Ok,C),
         cell(K,L,C)
        ),
      Ok)
     ),
   O),!.

% numerical evaluation : a matrix version of `is /2'
% modified: 8 June 2003. for cases of symbolic
matrix_is(C, A):-
   matrix(C is A).

matrix(C is X):-
   \+ var(X),
   (matrix(X,A);X=A),
   matrix_of_size(A,N,M),
   matrix_of_size(C,N,M),
   findall(Ck,
     (
      kth(K,A,Ak),
      kth(K,C,Ck),
      findall(Ckl,
        (
         kth(L,Ak,Akl),
         evaluate_if_numerical(Ckl is Akl), 
         kth(L,Ck,Ckl)
        ),
      Ck)
     ),
   C),!.

% mod
% added: 20 Aug 2003. in order for design0.pl

matrix_mod(C is (A mod 2^2)):-
   !,
   C = A.

matrix_mod(C is (A mod P)):-
   matrix_mod(C, A, mod(P)).

matrix_mod(C, A, mod(P)):-
   \+ var(A),
   integer(P),
   (matrix(X,A);X=A),
   matrix_of_size(A,N,M),
   matrix_of_size(C,N,M),
   findall(Ck,
     (
      kth(K,A,Ak),
      kth(K,C,Ck),
      findall(Ckl,
        (
         kth(L,Ak,Akl),
         Ckl is Akl mod P, 
         kth(L,Ck,Ckl)
        ),
      Ck)
     ),
   C),!.


% -----------------------------------------------------------  %
% basic calculations of matrix and linear algebra
% -----------------------------------------------------------  %
% modified: 8 June 2003.
%
matrix_add(X + Y = Z):-
   (matrix(X,A);X=A),
   (matrix(Y,B);Y=B),
   matrix_of_size(A,N,M),
   matrix_of_size(B,N,M),
   findall(ZK,
     (
      kth(K,A,AK),
      kth(K,B,BK),
      findall(ZKL,
        (
         kth(L,AK,AKL),
         kth(L,BK,BKL),
         evaluate_if_numerical(AKL1 is AKL),
         evaluate_if_numerical(BKL1 is BKL),
         evaluate_if_numerical(ZKL is AKL1 + BKL1)
        ),
      ZK)
     ),
   Z).

matrix_negative(- Y = Z):-
   matrix(Y,B),
   matrix_of_size(B,_N,_M),
   findall(ZK,
     (
      kth(_K,B,BK),
      findall(ZKL,
        (
         kth(_L,BK,BKL),
         evaluate_if_numerical(ZKL is (-1)*BKL)
        ),
      ZK)
     ),
   Z).

matrix_subtract(X - Y = Z):-
   matrix_negative(- Y = Z1),
   retractall(matrix(temp,_)),
   assert(matrix(temp,Z1)),
   matrix_add(X + temp = Z).


% matrix generator
% -----------------------------------------------------------  %

matrix_of_size(A,N,M):-  
   len(A,N),
   kth(1,A,A1),
   len(A1,M),
   findall(B,
     (
      member(B,A),
      len(B,M)
     ),
   A).

% restriction by rows and columns.
% -----------------------------------------------------------  %
restricted_matrix(X,[B],[C],[[D]]):-
   (matrix(X,A);X=A),
   matrix_of_size(A,_N0,_M0),
   kth(B,A,Ab),
   kth(C,Ab,D).
   
restricted_matrix(X,B,C,D):-
   % X: original matrix
   % B: restriction for rows
   % C: restriction for columns
   % D: restricted matrix
   B \= [_],
   C \= [_],
   (matrix(X,A);X=A),
   matrix_of_size(A,N0,M0),
   len(B,N),
   len(C,M),
   N =< N0,
   M =< M0,
   matrix_of_size(D,N,M),
   findall(DK,
     (
      member(K,B),
      kth(K,A,AK),
      findall(AKL,
        (
         member(L,C),
         kth(L,AK,AKL)
        ),
      DK)
     ),
   D).

% vertical / horizontal projections
% -----------------------------------------------------------  %
matrix_rows(A,[R],[B]):-
   kth(R,A,B).
matrix_rows(A1,R,B):-
   \+ R = [_],
   (matrix(A1,A);A=A1),
   matrix_of_size(A,_N,M),
   anum_seq1(Q,M),
   restricted_matrix(A,R,Q,B).

matrix_cols(A1,C,B):-
   (matrix(A1,A);A=A1),
   matrix_of_size(A,N,_M),
   anum_seq1(Q,N),
   restricted_matrix(A,Q,C,B).

% index for matrix
matrix_index(X,A,B,[C,D]):-
   (matrix(X,A);X=A),
   matrix_of_size(A,_N0,_M0),
   restricted_matrix(A,[C],[D],[[B]]).

% cf., collect the L(i)-th element from each row
%  acording to a choice list L.
matrix_collect(X,A,B,L):-
   matrix(X,A),
   index_of_tuple(A,B,L).

% transposition
% -----------------------------------------------------------  %
matrix_transpose(Y -> Z):-
   (matrix(Y,A);Y=A),
   matrix_of_size(A,N,M),
   matrix_of_size(Z,M,N),
   findall(ZK,
     (
      kth(K,Z,ZK),
      findall(ZKL,
        (
         kth(L,ZK,ZKL),
         matrix_index(Y,A,ZKL,[L,K])
        ),
      ZK)
     ),
   Z).


% multiplication
% -----------------------------------------------------------  %
% modified: 8 June 2003.
matrix_multiply(X * Y = Z):-
   \+ number(X),
   (matrix(X,A);X=A),
   (matrix(Y,B);Y=B),
   matrix_of_size(A,N,L0),
   matrix_of_size(B,L0,M),
   matrix_of_size(Z,N,M),
   matrix_transpose(Y -> C),!,
   findall(ZK,
     (
      kth(K,Z,ZK),
      findall(ZKL,
        (
         kth(L,ZK,ZKL),
         matrix_rows(A,[K],[AK]),
         matrix_rows(C,[L],[CL]),
         eq_sum_product(AK,CL,_,ZKL,_)
        ),
      ZK)
     ),
   Z).

matrix_multiply(X * Y = Z):-
   eval_number(X,_),
   (matrix(Y,B);Y=B),
   matrix_of_size(B,_N,M),
   matrix(X * i(M),Ix),
   matrix_multiply(Ix * Y = Z1),
   matrix(Z is Z1).

matrix_multiply_series([X],X,1):-
   matrix(X,A);matrix(A,X).

matrix_multiply_series([X|Y],Z,K):-
   matrix_multiply_series(Y,Z1,K1),
   K is K1 + 1,
   matrix_multiply(X * Z1 = Z).

matrix_powered(X ^ Y = Z):-
   integer(Y),
   matrix(X,_A),
   bag0(B,[X],Y),
   matrix_multiply_series(B,Z,Y).


% expanded matrix 
% -----------------------------------------------------------  %
matrix_expand(X | Y = C):-
   (matrix(X,A);X=A),
   (matrix(Y,B);Y=B),
   matrix_of_size(A,N,M),
   matrix_of_size(B,N,M1),
   M2 is M + M1,
   matrix_of_size(C,N,M2),
   findall(CK,
     (
      kth(K,A,AK),
      kth(K,B,BK),
      append(AK,BK,CK)
     ),
   C),
   !.

% modified: 8 Jun 2003.
% -----------------------------------------------------------  %
row_vector_multiply(C*B=D):-
   len(C,M),
   len(B,M),
   eq_sum_product(B,C,D,_,_).

row_vector_multiply(C0*B=D):-
   evaluate_if_numerical(C is C0),
   len(B,M),
   bag0(CL,[C],M),
   eq_sum_product(B,CL,D,_,_).


%
% the inverse of regular matrix.
% -----------------------------------------------------------  %

matrix_inverse(P ^ -1 = Q):-
   (matrix(P,B);P=B),
   matrix_of_size(B,N,M),
   len(Xs,N),
   eliminate_3(P,Xs,_,V),
   %matrix_of_size(Q,N,M),
   %matrix_expand(i(N)|Q = V). % good sense. but... 
   anum_seq1(Ys1,N),
   matrix(z(1,M,M),O),
   matrix_add([Ys1]+O=R),
   matrix([Ys] is R),
   matrix_cols(V,Ys,Q),
   !.


/*
% a sample execution.

?- matrix_inverse(p^(-1)=Q),matrix_multiply(p*Q=E),matrix(F is E).

Q = [[2, -6, -1], [-1, 4, 1], [1, -3, -1]]
E = [[2*1+3* -1+1*2, 2* -3+3*4+1* -6, 2* -1+3*1+1* -1], [1*1+1* -1+0*2, 1* -3+1*4+0* -6, 1* -1+1*1+0* -1], [-2*1+0* -1+1*2, -2* -3+0*4+1* -6, -2* -1+0*1+1* -1]]
F = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] 

*/

% -----------------------------------------------------------  %
% Gauss elimination
% -----------------------------------------------------------  %
%
%  row operations
% -----------------------------------------------------------  %
% modified: 8 Jun 2003  

row_operation(P,0,[_,_], _, P).

row_operation(P,C0,[K,L], A + C * B, FF):-
   % +P: a matrix
   % K,A: target row (K-th row)
   % L,B: pivot row (L-th row)
   % +C: constant
   % E: outcome
   evaluate_if_numerical(C is C0),
   \+ C = 0, 
   (matrix(P,Q);Q=P),
   matrix_of_size(Q,N,M),
   matrix_of_size(E,N,M),
   kth(K,Q,A),
   kth(L,Q,B),
   %matrix_rows(Q,[K],[A]),
   %matrix_rows(Q,[L],[B]),
   %matrix(A1 is A),
   %matrix(B1 is B),
   row_vector_multiply(C*B=CB),
   matrix_add([A]+[CB]=XXX),
   matrix([ACB] is XXX),
   replace(K/N,Q,ACB,E),
   nl,write(row_operation(row(K)=row(K)+const(C)*row(L))),
   matrix(FF is E),
   nl,write((P)),
   nl,write((->)),
   nl,write((FF)).

row_operation_series(1,P,[[C,K,L]],[A + C * B],[E]):-
   row_operation(P,C,[K,L], A + C * B, E).

row_operation_series(K,P,[[C,K,L]|X], [A + C * B|Y], [E|[F|G]]):-
   row_operation_series(K1,P,X,Y,[F|G]),
   K is K1 + 1,
   row_operation(F,C,[K,L], A + C * B, E).

expand_to_eliminate(A,B):-
   (matrix(A,P);P=A),
   matrix_of_size(P,N,N),
   matrix(i(N),I),
   matrix_expand(P|I=B).

%
% to make the diagonal element c(X,X)=1.
% -----------------------------------------------------------  %

% level 2
diagonalize_2(P,Q,[],[],[]):-
   (matrix(P,Q);P=Q),
   matrix_of_size(Q,_N,_M).
   %expand_to_eliminate(P,B).

diagonalize_2(P,Q,[E|F],[G|H],[X|Y]):-
   diagonalize_2(P,E,F,H,Y),
   len(E,N),
   (len(Y,N)->!,fail;true),
   matrix_index(E,_,Pxx,[X,X]),
   \+ member(X,Y),
   diagonalize_1((E,X,Pxx),G,Q,_).

% base level
diagonalize_1((B,X,Pxx),(L,C,A),P1,E):-
   % X,Pxx: target of diagnolization
   matrix_index(B,B,Pxx,[X,X]),
   case_of_diagonal(Q,Pxx,X,B,C,L),
   row_operation(B,C,[X,L], A, P1),
  %~~~~~~~~~~~~~~~~~~~~~~~!
   matrix(E is P1),
   nl,write(case_dialog_1((Q,X,X))).

% modified: 8 Jun 2003. 
case_of_diagonal(0,Pxx,X,_,C,X):-
   (
    (Pxx=..[_|Sxx],Sxx\=[],member(X,Sxx),\+ number(X));
    (Pxx=..[Sxx],\+ number(Sxx))
   ),
   C1 = 1/Pxx-1,
   evaluate_if_numerical(C is C1),
   !.

case_of_diagonal(1,Pxx,X,_,0,X):-
   Vxx is Pxx,
   Vxx is 1,
   !.

case_of_diagonal(2,Pxx,X,_,-2,X):-
   Vxx is Pxx,
   Vxx is -1,
   !.

case_of_diagonal(3,Pxx,X,_,C,X):-
   Vxx is Pxx,
   (Vxx > 0; Vxx < 0),
   C is (1-Pxx)/Pxx,
   !.

case_of_diagonal(4,Pxx,X,B,1/Plx,L):-
   Vxx is Pxx,
   Vxx is 0, 
   matrix_index(B,B,Plx,[L,X]),
   L > X,    % added: 18 Apr 2003.
   Vlx is Plx,
   (Vlx > 0; Vlx < 0),
   !.

%
% elimination of non diagonal elements.
% -----------------------------------------------------------  %

eliminate_3(P,[X],[Q],V):-
   expand_to_eliminate(P,PE),
   len(PE,N),
   matrix_index(PE,PE,_Pxx,[X,X]),
   len(Xs,N),
   eliminate_2(PE,PE,X,Xs,_,Q),
   matrix(V is Q).

eliminate_3(P,[X|Y],[Q|R],V):-
   eliminate_3(P,Y,[Q1|R],_V1),
   len(Q1,N),
   (len(Y,N)->!,fail;true),
   matrix_index(Q1,Q1,_Pxx,[X,X]),
   \+ member(X,Y),
   len(Xs,N),
   eliminate_2(Q1,Q1,X,Xs,_,Q),
   matrix(V is Q).


%
% eliminate non-zero in the X-th column.
%
eliminate_2(P,B,X,[K|G],[Q|H],V):-
   eliminate_2(P,B,X,G,[Q1|H],_V1),
   len(Q1,N),
   (len(G,N)->!,fail;true),
   matrix_index(Q1,Q1,Pkx,[K,X]),
   \+ member(K,G),
   case_of_elimination(_EC,Pkx,C),
   row_operation(Q1,C,[K,X], _A, Q),
   matrix(V is Q),
   nl,write(elim2(case_of_elim(_EC),[row:K,col:X])).

%
% start with diagonals of 1 (if numerical). 
% modified: 8 Jun 2003.
eliminate_2(P,B,X,[X],[Q],V):-
   eliminate_1(P,B,Q,V),
   matrix_index(Q,Q,Pxx,[X,X]),
   case_of_diagonal(Case,Pxx,X,_,C0,X),
   member(Case,[0,1]),
   nl,write(elim2(case_of_diag(1),[row:X,col:X],C0)).

eliminate_1(P,B,Q,V):-
   (matrix(P,B);P=B),
   matrix_of_size(B,N,_M),
   anum_seq1(Y,N),
   diagonalize_2(P,Q,_,_,Y),
   matrix(V is Q).

% added: 8 Jun 2003.
case_of_elimination(0,Pkx,Qkx):-
   (
    (Pkx=..[_|Skx],Skx\=[],member(X,Skx),\+ number(X));
    (Pkx=..[Skx],\+ number(Skx))
   ),
   evaluate_if_numerical(Qkx is (-Pkx)),
   !.

case_of_elimination(1,Pkx,-1):-
   Vkx is Pkx,
   Vkx is 1,
   !.

case_of_elimination(2,Pkx,1):-
   Vkx is Pkx,
   Vkx is -1,
   !.

case_of_elimination(3,Pkx,Qkx):-
   Vkx is Pkx,
   (Vkx > 0; Vkx < 0),
   Qkx is - Vkx,
   !.

case_of_elimination(4,Pkx,0):-
   Vkx is Pkx,
   Vkx is 0,
   !.

% -----------------------------------------------------------  %
% Arithmetic and so on including probabilistic operators
% -----------------------------------------------------------  %
%  edited: 26 Feb 2003. (imported from: math1.pl)

%
% evaluate the numeric value 'if' it is numeric equation
% -----------------------------------------------------------  %
% added: 8-9 June 2003 
evaluate_if_numerical(X1 is X):-
   functor(X,Functor,Arity),
   (
    (member(Functor,[+,-,*,/]),Arity=2);
    (member(Functor,[-]),Arity=1)
   ),
   X=..[Functor|Args],
   (
    forall( member(Y,Args), number(Y))
     -> X1 is X
     ; 
      (
       eval_equation(step1,X,XX0),
       eval_equation(step1b,XX0,XX1),
       eval_equation(step2,XX1,XX2),
       eval_equation(step3,XX2,X1)
      )
   ),
   !.
evaluate_if_numerical(X is X):-
   functor(X,Functor,Arity),
   \+ member(Functor,[+,-,*,/]),
   Arity=0,
   !.


%
% evaluation of a nummerical value
% -----------------------------------------------------------  %

%eval_equation(step1,X,X1)
eval_equation(step1,_/0, infinite).
eval_equation(step1,0/_, 0).
eval_equation(step1,A/A, 1).
eval_equation(step1,(A/B)*(B/A), 1).
eval_equation(step1,A/(C/A), 1/C).
     %-------
eval_equation(step1,0*_Arg2, 0).
eval_equation(step1,_Arg1*0, 0).
eval_equation(step1,0+Arg2, Arg2).
eval_equation(step1,0-Arg2, Arg2).
eval_equation(step1,Arg1+0, Arg1).
eval_equation(step1,Arg1-0, Arg1).
eval_equation(step1,A+(C-A), C).
eval_equation(step1,A-(C+A), C).
eval_equation(step1,(-A)+(A+C), C).
eval_equation(step1,(-A)+(C+A), C).
eval_equation(step1,(B+C)-B, C).
eval_equation(step1,(C+B)-B, C).
eval_equation(step1,(C-B)+B, C).
eval_equation(step1,((-B)+C)+B, C).
     %-------
eval_equation(step1,1*B, B).
eval_equation(step1,B*1, B).
eval_equation(step1,(-1)*B, -B).
eval_equation(step1,B*(-1), -B).
eval_equation(step1,(-A)+A, 0).
eval_equation(step1,A-A, 0).
eval_equation(step1,(-A)/A, -1).
eval_equation(step1,A/(-A), -1).
     %-------
eval_equation(step1,(-1/A)*A, -1).
eval_equation(step1,A*(-1/A), -1).
eval_equation(step1,-(1/A)*A, -1).
eval_equation(step1,A*(-(1/A)), -1).
eval_equation(step1,((-1)/A)*A, -1).
eval_equation(step1,A*((-1)/A), -1).
eval_equation(step1,(1/A)*(-A), -1).
eval_equation(step1,(-A)*(1/A), -1).
eval_equation(step1,(1/A)*A, 1).
eval_equation(step1,A*(1/A), 1).
eval_equation(step1,(-C/A)*A, -C).
eval_equation(step1,A*(-C/A), -C).
eval_equation(step1,(-(C/A))*A, -C).
eval_equation(step1,A*(-(C/A)), -C).
eval_equation(step1,((-C)/A)*A, -C).
eval_equation(step1,A*((-C)/A), -C).
eval_equation(step1,(C/A)*(-A), -C).
eval_equation(step1,(-A)*(C/A), -C).
eval_equation(step1,(C/A)*A, C).
eval_equation(step1,A*(C/A), C).
     %-------
eval_equation(step1,(-C/A)*(B/D), -(B*C)/(A*D)).
eval_equation(step1,(B/D)*(-C/A), -(B*C)/(A*D)).
eval_equation(step1,(-(C/A)*(B/D)), -(B*C)/(A*D)).
eval_equation(step1,(-(B/D)*(C/A)), -(B*C)/(A*D)).
eval_equation(step1,(-(C/A))*(B/D), -(B*C)/(A*D)).
eval_equation(step1,(B/D)*(-(C/A)), -(B*C)/(A*D)).
eval_equation(step1,((-C)/A)*(B/D), -(B*C)/(A*D)).
eval_equation(step1,(B/D)*((-C)/A), -(B*C)/(A*D)).
eval_equation(step1,(C/A)*(-B/D), -(C*B)/(A*D)).
eval_equation(step1,(-B/D)*(C/A), -(C*B)/(A*D)).
eval_equation(step1,(C/A)*((-B)/D), -(C*B)/(A*D)).
eval_equation(step1,((-B)/D)*(C/A), -(C*B)/(A*D)).
eval_equation(step1,(C/A)*(A/D), C/D).
eval_equation(step1,(A/D)*(C/A), C/D).
eval_equation(step1,(C/A)*(B/D), (C*B)/(A*D)).
     %-------
eval_equation(step1,(-C/A)*(A/D), -C/D).
eval_equation(step1,(A/D)*(-C/A), -C/D).
eval_equation(step1,(-(C/A)*(A/D)), -C/D).
eval_equation(step1,(-(A/D)*(C/A)), -C/D).
eval_equation(step1,(-(C/A))*(A/D), -C/D).
eval_equation(step1,(A/D)*(-(C/A)), -C/D).
eval_equation(step1,((-C)/A)*(A/D), -C/D).
eval_equation(step1,(A/D)*((-C)/A), -C/D).
eval_equation(step1,(C/A)*(-A/D), -C/D).
eval_equation(step1,(-A/D)*(C/A), -C/D).
eval_equation(step1,(C/A)*((-A)/D), -C/D).
eval_equation(step1,((-A)/D)*(C/A), -C/D).
     %-------
eval_equation(step1,(-1/A)*B, -B/A).
eval_equation(step1,B*(-1/A), -B/A).
eval_equation(step1,(-(1/A))*B, -B/A).
eval_equation(step1,B*(-(1/A)), -B/A).
eval_equation(step1,((-1)/A)*B, -B/A).
eval_equation(step1,B*((-1)/A), -B/A).
eval_equation(step1,(1/A)*(-B), -B/A).
eval_equation(step1,(-B)*(1/A), -B/A).
eval_equation(step1,(1/A)*B, B/A).
eval_equation(step1,B*(1/A), B/A).
eval_equation(step1,(-C/A)*B, -C*B/A).
eval_equation(step1,B*(-C/A), -C*B/A).
eval_equation(step1,(-(C/A))*B, -C*B/A).
eval_equation(step1,B*(-(C/A)), -C*B/A).
eval_equation(step1,((-C)/A)*B, -C*B/A).
eval_equation(step1,B*((-C)/A), -C*B/A).
eval_equation(step1,(C/A)*(-B), -C*B/A).
eval_equation(step1,(-B)*(C/A), -C*B/A).
eval_equation(step1,(C/A)*B, C*B/A).
eval_equation(step1,B*(C/A), C*B/A).
     %-------
eval_equation(step1,Arg1+(-B), Arg1-B).
eval_equation(step1,Arg1-(-B), Arg1+B).
eval_equation(step1,Arg1+B, Arg1-C):- number(B),C is abs(B),C>B.
eval_equation(step1,Arg1-B, Arg1+C):- number(B),C is abs(B),C>B.
     %-------
eval_equation(step1,A/(C*A), 1/C).
eval_equation(step1,A/(A*C), 1/C).
eval_equation(step1,(A*C)/A, C).
eval_equation(step1,(A*C)/A, C).
      %-------
eval_equation(step1,(-A)/(C*A), -1/C).
eval_equation(step1,(-A)/(A*C), -1/C).
eval_equation(step1,A/(C*(-A)), -1/C).
eval_equation(step1,A/((-A)*C), -1/C).
eval_equation(step1,(A*C)/(-A), -C).
eval_equation(step1,(A*C)/(-A), -C).
eval_equation(step1,((-A)*C)/A, -C).
eval_equation(step1,((-A)*C)/A, -C).
     %-------
eval_equation(step1,(A*B)/(C/A), B/C).
eval_equation(step1,(B*A)/(C/A), B/C).
eval_equation(step1,A/(C/(A*B)), B/C).
eval_equation(step1,A/(C/(B*A)), B/C).
eval_equation(step1,A/(C/B), A*B/C).
eval_equation(step1,(A/B)/C, A/(B*C)).
eval_equation(step1,A/B/C, A/(B*C)).
eval_equation(step1,(-A/B/C), (-A/(B*C))).
     %-------
eval_equation(step1,((-A)*B)/(C/A), -B/C).
eval_equation(step1,(B*A)/(C/(-A)), -B/C).
eval_equation(step1,(-A)/(C/(A*B)), -B/C).
eval_equation(step1,(-A)/(C/(B*A)), -B/C).
eval_equation(step1,(A*B)/(C/(-A)), -B/C).
eval_equation(step1,(B*(-A))/(C/A), -B/C).
eval_equation(step1,A/(C/((-A)*B)), -B/C).
eval_equation(step1,A/(C/(B*(-A))), -B/C).
     %-------
eval_equation(step1,A+(C*A), (C+1)*A).
eval_equation(step1,A+(A*C), (C+1)*A).
eval_equation(step1,(A*C)+A, (C+1)*A).
eval_equation(step1,(A*C)+A, (C+1)*A).
eval_equation(step1,A-(C*A), -(C-1)*A).
eval_equation(step1,A-(A*C), -(C-1)*A).
eval_equation(step1,(C*a)-A, (C-1)*A).
eval_equation(step1,(A*C)-A, (C-1)*A).
      %-------
eval_equation(step1,(-A)+(C*A), (C-1)*A).
eval_equation(step1,(-A)+(A*C), (C-1)*A).
eval_equation(step1,(A*C)+(-A), (C-1)*A).
eval_equation(step1,(A*C)+(-A), (C-1)*A).
eval_equation(step1,(-A)-(C*A), -(C+1)*A).
eval_equation(step1,(-A)-(A*C), -(C+1)*A).
eval_equation(step1,(A*C)-(-A), (C+1)*A).
eval_equation(step1,(A*C)-(-A), (C+1)*A).
      %-------
eval_equation(step1,(-A)*B, -(A*B)).
eval_equation(step1,A*(-B), -(A*B)).
eval_equation(step1,(-A)/B, -A/B).
eval_equation(step1,A/(-B), -A/B).
eval_equation(step1,(-A)/(-B), A/B).
eval_equation(step1,(-(-A/B)), A/B).
      %-------
eval_equation(step1,X,X).

eval_equation(step1b,A*B, C):-
   eval_equation(step1,A, AX),
   eval_equation(step1,B, BX),
   eval_equation(step1,AX*BX, C).
eval_equation(step1b,A-B, C):-
   eval_equation(step1,A, AX),
   eval_equation(step1,B, BX),
   eval_equation(step1,AX-BX, C).
eval_equation(step1b,A/B, C):-
   eval_equation(step1,A, AX),
   eval_equation(step1,B, BX),
   eval_equation(step1,AX/BX, C).
eval_equation(step1b,A+B, C):-
   eval_equation(step1,A, AX),
   eval_equation(step1,B, BX),
   eval_equation(step1,AX+BX, C).
eval_equation(step1b,X,X).

eval_equation(step2,X,X1):-
       (
        member(X,
          [
           (C*A)+(D*A),
           (C*A)+(A*D),
           (A*C)+(D*A),
           (A*C)+(A*D)
          ]
        )
         ->
        (eval_equation(step1,C+D,Eq1),X1 = Eq1*A)
       );
      %-------
       (
        member(X,
          [
           (C*(-A))+(D*A),
           (C*(-A))+(A*D),
           ((-A)*C)+(D*A),
           ((-A)*C)+(A*D)
          ]
        )
         ->
        (eval_equation(step1,C-D,Eq1),X1 = Eq1*A)
       );
      %-------
       (
        member(X,
          [
           (C*A)+(D*(-A)),
           (C*A)+((-A)*D),
           (A*C)+(D*(-A)),
           (A*C)+((-A)*D)
          ]
        )
         ->
        (eval_equation(step1,C-D,Eq1),X1 = Eq1*A)
       );
     %-------
       (
        member(X,
          [
           (C*A)-(D*A),
           (C*A)-(A*D),
           (A*C)-(D*A),
           (A*C)-(A*D)
          ]
        )
         ->
        (eval_equation(step1,C-D,Eq1),X1 = Eq1*A)
       );
      %-------
       (
        member(X,
          [
           (C*(-A))-(D*A),
           (C*(-A))-(A*D),
           ((-A)*C)-(D*A),
           ((-A)*C)-(A*D)
          ]
        )
         ->
        (eval_equation(step1,(-(C+D)),Eq1),X1 = Eq1*A)
       );

      %-------
       (
        member(X,
          [
           (C*A)-(D*(-A)),
           (C*A)-((-A)*D),
           (A*C)-(D*(-A)),
           (A*C)-((-A)*D)
          ]
        )
         ->
        (eval_equation(step1,(-(C+D)),Eq1),X1 = Eq1*A)
       );
     %-------
       (
        member(X,
          [
           (C*A)/(D*A),
           (C*A)/(A*D),
           (A*C)/(D*A),
           (A*C)/(A*D)
          ]
        )
         ->
        (eval_equation(step1,(C/D),Eq1),X1 = Eq1*A)
       );
      %-------
       (
        member(X,
          [
           (C*(-A))/(D*A),
           (C*(-A))/(A*D),
           ((-A)*C)/(D*A),
           ((-A)*C)/(A*D)
          ]
        )
         ->
        (eval_equation(step1,(-(C/D)),Eq1),X1 = Eq1*A)
       );
      %-------
       (
        member(X,
          [
           (C*A)/(D*(-A)),
           (C*A)/((-A)*D),
           (A*C)/(D*(-A)),
           (A*C)/((-A)*D)
          ]
        )
         ->
        (eval_equation(step1,(-(C/D)),Eq1),X1 = Eq1*A)
       );
     %-------
       X1=X.
eval_equation(step4,Num/Den,X1):-
     eval_equation(step1,Num,NumEV),
     eval_equation(step1,Den,DenEV),
     eval_equation(step1,NumEV/DenEV,X1).
eval_equation(step5,(Num1+Num2)/Den,X1):-
     eval_equation(step1,Num1,NumEV1),
     eval_equation(step1,Num2,NumEV2),
     eval_equation(step1,NumEV1+NumEV2,NumEV),
     eval_equation(step1,Den,DenEV),
     eval_equation(step1,NumEV/DenEV,X1).
eval_equation(step3,X,X1):-
/*
       (X=C+(D/A) -> (eval_equation(step1,(C*A),E),eval_equation(step4,[(E+D)/A],X1)));
       (X=C+(-D/A) -> (eval_equation(step1,(C*A),E),eval_equation(step4,[(E-D)/A],X1)));
       (X=C+(-(D/A)) -> (eval_equation(step1,(C*A),E),eval_equation(step4,[(E-D)/A],X1)));
       (X=C+((-D)/A) -> (eval_equation(step1,(C*A),E),eval_equation(step4,[(E-D)/A],X1)));
       (X=C+(D/(-A)) -> (eval_equation(step1,(C*A),E),eval_equation(step4,[(E-D)/A],X1)));
*/
     %-------
       (X=(C/A)+(D/A) -> eval_equation(step4,[(C+D)/A],X1));
       (X=(-C/A)+(-D/A) -> eval_equation(step4,[(C+D)/A],X1));
       (X=(C/A)+(-D/A) -> eval_equation(step4,[(C-D)/A],X1));
       (X=(C/A)+(-(D/A)) -> eval_equation(step4,[(C-D)/A],X1));
       (X=(C/A)+((-D)/A) -> eval_equation(step4,[(C-D)/A],X1));
       (X=(C/A)+(D/(-A)) -> eval_equation(step4,[(C-D)/A],X1));
       (X=(-C/A)+(D/A) -> eval_equation(step4,[(D-C)/A],X1));
       (X=(-(C/A))+(D/A) -> eval_equation(step4,[(D-C)/A],X1));
       (X=(-(C)/A)+(D/A) -> eval_equation(step4,[(D-C)/A],X1));
       (X=(C/(-A))+(D/A) -> eval_equation(step4,[(D-C)/A],X1));
     %-------
       (X=(C/A)-(D/A)  -> eval_equation(step4,[(C-D)/A],X1));
       (X=(-C/A)-(-D/A)  -> eval_equation(step4,[(D-C)/A],X1));
       (X=(-C/A)-(D/A)  -> eval_equation(step4,[(-C-D)/A],X1));
       (X=(C/A)-(-D/A)  -> eval_equation(step4,[(C-D)/A],X1)); 
       (X=(-(C/A))-(D/A)  -> eval_equation(step4,[(-C-D)/A],X1));
       (X=((-C)/A)-(D/A)  -> eval_equation(step4,[(-C-D)/A],X1));
       (X=(C/(-A))-(D/A)  -> eval_equation(step4,[(-C-D)/A],X1));
       (X=(C/A)-(-(D/A))  -> eval_equation(step4,[(C-D)/A],X1));
       (X=(C/A)-((-D)/A)  -> eval_equation(step4,[(C-D)/A],X1));
       (X=(C/A)-(D/(-A))  -> eval_equation(step4,[(C-D)/A],X1));
      %-------
       (X=(C/A)+(D/B)  -> eval_equation(step5,[(C*D+D*A)/AB],X1));
       (X=(-C/A)+(D/B) -> eval_equation(step5,[((-C*B)+D*A)/AB],X1));
       (X=(-(C/A))+(D/B) -> eval_equation(step5,[((-C*B)+D*A)/AB],X1));
       (X=((-C)/A)+(D/B) -> eval_equation(step5,[((-C*B)+D*A)/AB],X1));
       (X=(C/(-A))+(D/B) -> eval_equation(step5,[((-C*B)+D*A)/AB],X1));
       (X=(C/A)+(-D/B) -> eval_equation(step5,[(C*B-D*A)/AB],X1));
       (X=(C/A)+(-(D/B)) -> eval_equation(step5,[(C*B-D*A)/AB],X1));
       (X=(C/A)+((-D)/B) -> eval_equation(step5,[(C*B-D*A)/AB],X1));
       (X=(C/A)+(D/(-B)) -> eval_equation(step5,[(C*B-D*A)/AB],X1));
     %-------
       (X=(C/A)-(D/B)  -> eval_equation(step5,[(C*D-D*A)/AB],X1));
       (X=(-C/A)-(D/B) -> eval_equation(step5,[((-C*B)-D*A)/AB],X1));
       (X=(-(C/A))-(D/B) -> eval_equation(step5,[((-C*B)-D*A)/AB],X1));
       (X=((-C)/A)-(D/B) -> eval_equation(step5,[((-C*B)-D*A)/AB],X1));
       (X=(C/(-A))-(D/B) -> eval_equation(step5,[((-C*B)-D*A)/AB],X1));
       (X=(C/A)-(-D/B) -> eval_equation(step5,[(C*B+D*A)/AB],X1));
       (X=(C/A)-(-(D/B)) -> eval_equation(step5,[(C*B+D*A)/AB],X1));
       (X=(C/A)-((-D)/B) -> eval_equation(step5,[(C*B+D*A)/AB],X1));
       (X=(C/A)-(D/(-B)) -> eval_equation(step5,[(C*B+D*A)/AB],X1));
     %-------
       X1=X.


%
% evaluation of a nummerical value
% -----------------------------------------------------------  %
eval_number(X,X):-
   number(X),!.


%
% max,min
% -----------------------------------------------------------  %
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).


% count frequency of occurence of the specified value of variable, M.
% -----------------------------------------------------------  %
% note: Both of M and L have to be specified.

counter(N,M,L):-
    len(L,_),
    findall(M,member(M,L),Mx),
    len(Mx,N).

% sum
% -----------------------------------------------------------  %
sum([],0).
sum([X|Members],Sum):-
   sum(Members,Sum1),
  %number(X),
   Sum is Sum1 + 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):-
    len(Q,N),
    len(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):-
    len(Q,N),
    len(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.


%
% symbolic representation of eq_sum_product (product_sum). 
% added: 8-9 June 2003 
% -----------------------------------------------------------  %
eq_sum_product([],[],[],0,0).
eq_sum_product([P|Q],[A|B],[E|F],V,Vq):-
    len(Q,N),
    len(B,N),
    eq_sum_product(Q,B,F,V1,Vq1),
    evaluate_if_numerical(EP is P),
    evaluate_if_numerical(EA is A),
    evaluate_if_numerical(E is EP * EA),
    (Vq1=0 -> Vq = EP*EA; Vq = Vq1 + (EP * EA)),
    evaluate_if_numerical(V is V1+E).
rev_eq_sum_product(P,A,F1,V,Vq):-
    rev(P,Q),
    rev(A,B),
    rev_eq_sum_product(Q,B,F,V,Vq),
    rev(F,F1).


%
% symbolic representation of sum. eqsum/2 and reqsum/2
% cited from: eba01.pl (28 May 2003) 
% -----------------------------------------------------------  %
eqsum([],0).
eqsum([X|Members],Sum):-
   eqsum(Members,Sum0),
   (
    X=0 -> Sum = Sum0
     ;
    (
     Sum0=0 -> Sum = X
      ;
     Sum = Sum0 + X
    )
   ).
reqsum(A,B):-
   rev(A,C),
   eqsum(C,B).


%
% allocation
% -----------------------------------------------------------  %
allocation(N,A,[X|Y]):-
    allocation(N,A,A,[X|Y]).
allocation(0,_,0,[]).
allocation(N,A,B,[X|Y]):-
    integer(A),
    len([X|Y],N),
    allocation(_N1,A,B1,Y),
    % N1 is N - 1,
    % sum(Y,B1),
    K is A - B1 + 1,
    len(L,K),
    nth0(X,L,X),
    B is B1 + X.
%
% probability (percentile) by using allocation
% -----------------------------------------------------------  %
probabilities(0,[]).
probabilities(N,[X|Y]):-
    integer(N),
    len([X|Y],N),
    allocation(N,100,[X|Y]).
% 
% any ratio (weight) can be interpreted into a prob.
% -----------------------------------------------------------  %
scale(W,1/Z,P):-
    findall(Y,(kth(_K,W,X),Y is X/Z),P).
probabilities(W,N,P):-
    len(W,N),
    sum(W,Z),
    scale(W,1/Z,P).
%
% degenerate probability
%---------------------------------------------
degenerate(Ps):-
   kth(K,Ps,P),
   characteristic_vector(K,P,Ps,Ps).
% 
% probability over base set with steps of levels.
% -----------------------------------------------------------  %
make_a_prob(P,base(M),steps(L)):-
    var(P),
    len(P,M),
    allocation(M,L,W),
    probabilities(W,M,P).
make_a_prob(P,base(M),_):-
    \+ var(P),
    len(P,M),
    \+ (
     member(P1,P),
     (
      var(P1);
      P1 > 1;
      P1 < 0
     )
    ),
    sum(P,1).
%
% expected value
% -----------------------------------------------------------  %
expected_value(W,A,E/100):-
    len(A,N),
    probabilities(W,N,P),
    product_sum(P,A,_,E).
%
% expected value with equational
% -----------------------------------------------------------  %
expected_value_eq(W,A,E/100,Eq):-
    len(A,N),
    probabilities(W,N,P),
    product_sum_eq(P,A,_,E,Eq).

%
%---------------------------------------------------
%  net present value (NPV) 
%---------------------------------------------------
%
%  time preference: discount factor
%---------------------------------------------

interest_rate(1.1).

discount_factor(R,Y,DF,DFV):-
   DF = R ^ (-Y),
   DFV is DF.

npv(A,Y,Eq,V):-
   interest_rate(R),
   discount_factor(R,Y,DF,_),
   Eq = DF * A,
   V is Eq.

%
% conditional probabilities
% -----------------------------------------------------------  %
probability_of_event(W,E,P):-
    % conditionalization by event specified directly
    event(E),
    (var(E)->E = E1; sort(E,E1)),
    G = member(S,E1),
    findall(A,(probability(W,S,A),G),Ps),
    sum(Ps,P).
probability_of_event(W,E,P,G):-
    \+ var(G), % conditionalization via constraints indirectly
    G=(Goal,M,[W,S,A]),  % constraints with params
    findall([S1,A1],
      (
       (M=do->(W=W1,S=S1,A=A1);true),
       probability(W1,S1,A1),
       Goal
      ),
    Xs),
    findall(S,member([S,A],Xs),E0),
    findall(A,member([S,A],Xs),Ps),
    sort(E0,E),
    sum(Ps,P).
%
% -----------------------------------------------------------  %
%   Utilities for list operations
% -----------------------------------------------------------  %
%  edited: 14 Feb 2003. (imported from: set1.pl)
%
% index for tuples.
% -----------------------------------------------------------  %
% 1) only mention for a direct product of sets.
index_of_tuple(B,A,Index):-
   \+ var(B),
   \+ var(A),
   len(B,LN),  % base sets
   len(A,LN),  
   len(Index,LN),
   findall(L,
     (
      kth(K,B,BJ), %write(a(K,B,BJ)),
      kth(L,BJ,AJ),%write(b(L,BJ,AJ)),
      kth(K,A,AJ)  %,write(c(K,A,AJ)),nl
     ),
   Index).
index_of_tuple(B,A,Index):-
   \+ var(B),
   \+ var(Index),
   var(A),
   len(B,LN),  % base sets
   len(Index,LN),
   len(A,LN),  
   findall(AJ,
     (
      kth(K,B,BJ),
      kth(K,Index,L),
      kth(L,BJ,AJ)
     ),
   A).
%
% descending/ascending natural number sequence less than N.
% -----------------------------------------------------------  %
dnum_seq([],N):-N<0,!.
dnum_seq([0],1).
dnum_seq([A|Q],N):-
   A is N - 1,
   len(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).
%
% inquire the goal multiplicity
% -----------------------------------------------------------  %
sea_multiple(Goal,Cond,N,M):-
  Clause=..Goal,
  findall(Cond,Clause,Z),len(Z,N),sort(Z,Q),len(Q,M).
%
bag0([],_A,0).
bag0([C|B],A,N):-
   len([C|B],N),
   bag0(B,A,_N0),
   member(C,A).
zeros(Zero,N):-bag0(Zero,[0],N).
ones(One,N):-bag0(One,[1],N).
%
% subset_of/3 : subset-enumeration 
% -----------------------------------------------------------  %
subset_of(A,N,As):-
   len(As,L),
   len(D,L),
   list_projection(D,As,B),
   len(B,N),
   sort(B,A).
% a sequence of binary choice for a list:
%--------------------------------------------------
list_projection([],[],[]).
list_projection([X|Y],[_A|B],C):-
   X = 0,
   list_projection(Y,B,C).
list_projection([X|Y],[A|B],[A|C]):-
   X = 1,
   list_projection(Y,B,C).
%
% characteristic_vector/3
% -----------------------------------------------------------  %
% modified: 8 Feb 2003.  without using nth1.
% modified: 13 Feb 2003.  without using member.
characteristic_vector(X,B,Index):-
   \+ var(B),
   %member(X,B),
   list_projection(Index,B,[X]).
characteristic_vector(1,X,[X|B],[1|DX]):-
   characteristic_vector(X,[X|B],[1|DX]).
characteristic_vector(K,X,[_|B],[0|DX]):-
   characteristic_vector(K1,X,B,DX),
   K is K1 + 1.
%
% replace(Project,Goal,Base,Goal1):-
% -----------------------------------------------------------  %
% added: 15 Oct 2002.
  % a sequence of replacement of a subset of elements in Goal 
  % which specified by a list, Project, 0-1^n, over Base 
  % a list of length n, which brings about Goal1.
  % summary:
  %  X=1 --> preserve the value of Base.
  %  X=0 --> do replace with Goal1.
replace([],[],[],[]).
replace([X|A],[_|B],[Z|C],[Z|D]):-
   X = 0,
   replace(A,B,C,D).
replace([X|A],[Y|B],[_|C],[Y|D]):-
   X = 1,
   replace(A,B,C,D).
%
% replace/4
% -----------------------------------------------------------  %
% modified: 14 Feb 2003. bug fix.
replace(K/N,L,S,L1):-
   \+ var(S),
   \+ var(L),
   len(L,N),
   len(L1,N),
   kth(K,L1,S),
   characteristic_vector(K,_S0,L,V),
   c_replace(V,L,L1,L1).

%
c_replace([],[],[],[]).
c_replace([X|A],[_|B],[Z|C],[Z|D]):-
   X = 1,
   c_replace(A,B,C,D).
c_replace([X|A],[Y|B],[_|C],[Y|D]):-
   X = 0,
   c_replace(A,B,C,D).

%
% sort without removal of duplicates
%--------------------------------------------------
asort(A,B):-
   sort(A,C),
   bagof(CK,
     J^K^(
      kth(J,C,CK),
      kth(K,A,CK)
     ),
   B).
%
% permutation.
% -----------------------------------------------------------  %
% modified: 1 Sep 2002. to be used in is_neutral/2. 
% modified: 15 Oct 2002. add a non-variable constraint for the base set A. 

permutation_of(A,P,APs):-
   \+var(A),
   len(A,M),
   ordering(P,A,M),
   asc_nnseq(Qm,M),
   my_maplist(nth_of_permutation(A,P),Qm,APs).

%caution:
%notations 'A->B's bellow have cause precedence errors in IF-prolog.

nth_poa(P,K,Ak->Pk):-
   alternatives(A),
   nth_of_permutation(A,P,K,Ak->Pk).
nth_of_permutation(A,P,K,Ak->Pk):-
   len(A,M),
   ordering(P,A,M),
   nth_0(K,A,Ak),
   nth_0(K,P,Pk).

permutate_of_order(PoA,[Q->R]):-
   poa(_P,PoA),
   alternatives(A),
   len(A,M),
   ordering(Q,A,M),%wn(Q),
   ordering(R,A,M),%wn(R),
   permutation(Q,PoA,R).

permutation([],[],[]).
permutation(Q,[A->P|PoA1],R):-
   subtract(Q,[A],Q1),nth_1(K,Q,A),
   subtract(R,[P],R1),nth_1(K,R,P),
   permutation(Q1,PoA1,R1).

%
% projection operator via index set. (an exchange economy?)
% -----------------------------------------------------------  %
% choice1 from base1 :: choice2 from base2.
% modified : 15 Oct 2002. to be order-neutral (pending:-)
% 
pcm([Choice1,Base1],[Choice2,Base2]):-
    pairwise_contract_map([Choice1,Base1],[Choice2,Base2]).


pairwise_contract_map([Choice1,Base1],[Choice2,Base2]):-
  %  len(Base1,N2),
  %  len(Base2,N2),
    subset_of(Choice1,N1,Base1),
    subset_of(Choice2,N1,Base2),
    list_projection(Project,Base1,Choice1),
    list_projection(Project,Base2,Choice2).
%
% ppcm /2 using plist_projection 
% added: 15 Oct 2002.
ppcm([Choice1,Base1],[Choice2,Base2]):-
    subset_of(C1,N1,Base1),
    subset_of(C2,N1,Base2),
    plist_projection(Project,Base1,Choice1,C1),
    plist_projection(Project,Base2,Choice2,C2).

%
% -----------------------------------------------------------  %
%   Utilities for outputs
% -----------------------------------------------------------  %

% redirect to file
%--------------------------------

tell_goal(File,G):-
   (current_stream(File,write,S0)->close(S0);true),
   open(File,write,S),
   tell(File),
   nl,
   time_stamp('% file output start time ',_),
   nl,
   write('%----------  start from here ------------%'),
   nl,
   G,
   nl,
   write('%----------  end of data ------------%'),
   nl,
   time_stamp('% 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_goal(File,forall,G):-
   G0 = (nl,write(G),write('.')),
   G1 = forall(G,G0),
   tell_goal(File,G1).


tell_goal(File,forall_such_that,G,Condition):-
   % G should be occurred in the Condition.
   WRITE = (nl,write(G),write('.')),
   G1 = forall(Condition,WRITE),
   tell_goal(File,G1).


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

time_stamp(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.

time_stamp(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.


% len/2 : alternative to the length/2

len(A,B):-
   len_0(A,B).

% length/2 of SWI-prolog fails at the case of both unbound variables.
len_0(A,B):-
   var(A),
   var(B),
   !,
   fail.
len_0([],0).
len_0([_|A],N):-
   ((integer(N),N>0)->N0 is N -1; true),
   len_0(A,N0),
   ((integer(N0))->true; !,fail),
   N is N0 + 1.

% kth/3 : alternative to the nth1/3
kth(K,Y,X):-
   kth_member(K,Y,X).

kth_member(1,[X|Y],X):-
   \+ var(Y).
kth_member(K,[_|Y],X):-
   \+ var(Y),
   kth_member(K1,Y,X),
   K is K1 + 1.

% rev/2 : alternative to the reverse/2
rev(A,B):-
   len(A,L),
   len(B,L),
   rev(A,[],B),
   !.

rev(A,A,[]):-!.
rev(A,B,[C|D]):-
   rev(A,[C|B],D).




return to front page.