You selected belief02.pl

headline:-
  wn('% -----------------------------------------------------------  %'),
  wn('%   ambiguous beliefs and decision making: Choquet rationality %'),
  wn('% -----------------------------------------------------------  %'),
  h0.
h0:-
   WL=[
'%  main predicates in this program: ',
'%  set_model/0,1: select a belief system from the model base. ',
'%  ------------ game theory under uncertain beliefs ------------ ',
'%  support/1: support of current belief model bel/3.  ',
'%  best_response/2: nash strategy under ambiguous beliefs.  ',
'%  bleifs_equilibrium/2: equilibrium of the game in beliefs.  ',
'%  -------------- Choquet rational decision making ------------ ',
'%  test_ceu/3: Choquet expected utility (CEU) with belief function.',
'%  payoff(choquet(A),B,C):  another CEU construction ',
'%  payoff(expected_conservative(A),B,C): additive decomposition.',
'%  --------------- belief functions ---------------------------- ',
'%  sbel/3,spos/3,bpa/3:  belief/possibility functions and b.p.a. ',
'%  modularity/1,3,4, covexity/1:  check for modularity and convexity. ',
'%  update_mass/3,update_sbel/3,update_spos/3:  Dempster-Shafer rule. ',
'%  update_bel_ul/3,update_pos_ul/3: Upper-lower rule. ',
'%  --------------- references ---------------------------------- ',
'%  references1,references2: 1:belief and decision, 2:game theory.',
'%  me: edit history and about the author.',
'%  h0/0:  show this command help again.'
   ],
   display_all_lines(WL),nl.
me:-
   M=[
'%  file: belief02.pl',
'%  privious version: belief01.pl(29 Jun 2003)',
'%  edited: 14--23 Jan 2004 (belief02.pl)',
'%  author: Kenryo INDO (Kanto Gakuen University)',
'%  modified: not yet'
   ],
   display_all_lines(M).
references1:-
   R1=[
'% references:',
'% [1] Shafer,G.(1976). Mathematical Theory of Evidence. Princeton ',
'%  University Press.',
'% [2] Dempster, A.P.(1967). Upper and lower probabilities induced by a',
'%  multivalued mapping. Annals of Mathematical Statistics 38: 325-339.',
'% [3] Moral,S.and De Campos, L.M.(1991). Updating uncertain information.',
'%  In the Proceedings of IPMU90, LNCS 521, Springer, pp.58-67. ',
'% [4] Dow, J. and S. R. da Costa Werlang (1992). Excess volatility of ',
'%  stock prices and Knightian uncertainty. European Economic Review ',
'%   36: 631-638.',
'% [5] Jaffray J-.Y. and P. Wakker (1994). Decision making with belief ',
'%  functions: compatibility and incompatibility with the sure-thing ',
'%  principle. Journal of Risk and Uncertainty 8: 255-271.',
'% [6] Gilboa,I. and D. Schmeidler (1993). Updating ambiguous beliefs.',
'%  Journal of Economic Theory 59: 33-49.',
'% [7] Mukerji, S. (1997). Understanding the nonadditive probability ',
'%  decision model. Economic Theory 9: 23-46.'
   ],
   display_all_lines(R1).
references2:-
   R2=[
'% References: ',
'% [8] Gilboa, I. and D. Schmeidler (1989). JME 18: 141-153.',
'% [9] Dow, J. and S.R.C. Werlang (1994). JET 64: 305-324.',
'% [10] Lo, K.C. (1996). JET 71: 443-484.',
'% [11] Groes, E. et al. (1998). Theory and Decision 44: 37-66.',
'% [12] Eichberger, J. and D. Kelsey (2000). GEB 30: 183-215.',
'% [13] Marinacci, M. (2000). GEB 31: 191-219.'
   ],
   display_all_lines(R2).

wn(Z):-
   write(Z),nl.
display_all_lines(A):-
   forall(member(B,A),(nl,write(B))).

:- headline.


%---------------------------------------------------
%   User Interface : Model Setting
%---------------------------------------------------

:- dynamic current_model/1.
:- dynamic states/1.
:- dynamic bpa0/2.
:- dynamic bel0/2.
:- dynamic pos0/2.
:- dynamic payoff0/3.
:- dynamic act/1.

% default
current_model(ipd2).

models([ipd2]).

set_model:-
  set_model(_EX).

set_model(EX):-
   models(M),
   (
    var(EX)
    ->
     (
      nl,write((models:M)),
      nl,write('please specify the model >'),
      read(EX)
     )
    ;
   true),
   (member(EX,M)->true;(write('no such model. '),fail)),
   nl,write((EX,'Would you like to set this model ? (y) >')),
   read(y),
   update_current_model(EX).

update_current_model(EX):-
   models(M),
   member(EX,M),
   initialize_model,
   forall(
     (model_predicate(F,Arity),F1=..[F,EX|Arity],clause(F1,_),F1),
     (G=..[F|Arity],assert(G))
   ).
   %display_current_model.
      
initialize_model:-
   forall(
     (model_predicate(F,Arity),length(Arity,N)),
     abolish(F/N)
   ).

display_current_model:-
   current_model(EX),
   (nl,write((model:EX))),
   forall(
     (model_predicate(F,Arity),F1=..[F|Arity],clause(F1,_),F1),
     (nl,write(F1))
   ).

model_predicate(states,[_]).
model_predicate(bpa0,[_,_]).
model_predicate(bel0,[_,_]).
model_predicate(pos0,[_,_]).
model_predicate(payoff0,[_,_,_]).
model_predicate(act,[_]).



:- initialize_model.



%------------------
%   Examples
%------------------
%
% example 1. Ellsberg 3 color problem
% (Gilboa and Schmeidler, 1993, p.36)
%---------------------------------------------------
/*
states([r,b,y]).

% Basic probability assignment of the example
bpa0([],0).
bpa0([r],1/3).
bpa0([b,y],2/3).
%bpa0([r,b,y],1/3).  % violate to normalization.
*/

/*
bel0([],0).
bel0(E,1/3):- member(E,[[r],[b,r],[r,y]]).
bel0([b,y],2/3).
bel0([b,r,y],1).

% note: 
% see also the results of update_sbel(A/[b,y],_,B) for example 1.

*/


% example 2. Dow and Werlang(1992)'s trader
%---------------------------------------------------
/*
states([1,2,3]).

asset(stock).
asset(cash).
time(T):-member(T,[0,1,2]).
return(cash,0).
return(stock,state(1),1).
return(stock,state(2),1/2).
return(stock,state(3),0).
partition(h1,[1]).
partition(h2,[2,3]).
know(S,H,E):-
   partition(_,H),
   member(S,H),
   sevent(E),
   subset_of(H,E,_).
%
belief(time(0),[],0).
belief(time(0),X,1/4):- member(X,[[1],[2],[3]]).
belief(time(0),X,1/2):- member(X,[[1,2],[2,3],[1,3]]).
belief(time(0),[1,2,3],1).

% belief function for example 3.
bel0(F,Yq,Y):-
  event(F),
  sort(F,F1),
  belief(time(0),F1,Yq),
  Y is Yq.

*/

% note: the core of v(.)=bel(.) in example 2 is 
% the convex hull of additive measures.
%  Core(v(.))=[p | p(.)>=v(.)].
%
%  p*([1])=1/2, p*([2])=1/4, p*([3])=1/4,
%  p*([1])=1/4, p*([2])=1/2, p*([3])=1/4,
%  p*([1])=1/4, p*([2])=1/4, p*([3])=1/2.
%

% example 3. TV set sales
% in Jaffray and Wakker(1994), example 2.1.
%---------------------------------------------------
/*

states([qL,qM,qH]).

bpa0([],0).
bpa0([qL],0).
bpa0([qM],0).
bpa0([qH],0.1).
bpa0([qL,qM],0.6).
bpa0([qL,qH],0).
bpa0([qM,qH],0.3).
bpa0([qL,qM,qH],0).

bel0([],0).
bel0([qL],0).
bel0([qM],0).
bel0([qH],0.1).
bel0([qL,qM],0.6).
bel0([qL,qH],0.1).
bel0([qM,qH],0.4).
bel0([qL,qM,qH],1).
*/

/*
pos0([],0).
pos0([qL],0.6).
pos0([qM],0.9).
pos0([qH],0.4).
pos0([qL,qM],0.9).
pos0([qL,qH],1).
pos0([qM,qH],1).
pos0([qL,qM,qH],1).
*/


/*
Table 1. Beliefs in the TV sales example.

           []    [L]   [M]   [H]  [L,M]  [L,H]  [M,H]   S
------------------------------------------------------------
Bel(=v)     0    0.0   0.0   0.1   0.6    0.1    0.4    1
Pl (=v~)    0	 0.6   0.9   0.4   0.9    1.0    1.0    1
m           0	 0.0   0.0   0.1   0.6    0.0    0.3    0
------------------------------------------------------------
 (Cited from Jaffray and Wakker(1994), Table 3.1).

% mass and belief of example 2, 
% the TV set sales problem (Jaffray and Wakker, 1994)

(a) S=[L,M,H]
 m([])=0.
 m([L])=v([L])=0. 
 m([M])=v([M])=0.
 m([H])=v([H])=0.1.
 m([L,M])=-v([L])-v([M])+v([L,M])=0.6.
 m([L,H])=-v([L])-v([H])+v([L,H]) =-0.1+0.1=0.
 m([M,H])=-v([M])-v([H])+v([M,H])=-0.1+0.4=0.3.
 m(S)=v([L])+v([M])+v([H])-v([L,M])-v([L,H])-v([M,H])+v([S]) 
     =0.1-1.1+1=1-(0.1+0.6+0.3)=0.
(b) bpa->bel
 v([])= m([])=0.
 v([L])= m([])+m([L])=0. 
 v([M])= m([])+m([M])=0. 
 v([H])= m([])+m([H])=0.1.
 v([L,M])
    = m([])+m([L])+m([M])+m([L,M])=0.6.
 v([L,H])
    = m([])+m([L])+m([H])+m([L,H])=0.1.
 v([M,H])
    = m([])+m([M])+m([H])+m([M,H])=0.4.
 v(S)
    = m([])+m([L])+m([M])+m([H])
     +m([L,M])+m([L,H])+m([M,H])
     +m([L,M,H])=1.
(c) bpa->pos
 v~([])=0.
 v~([L])
   = m([L])+m([L,M])+m([L,H])+m(S)=0.6. 
 v~([M])
   = m([M])+m([L,M])+m([M,H])+m(S)=0.9. 
 v~([H])
   = m([H])+m([L,H])+m([M,H])+m(S)=0.4. 
 v~([L,M])
   = m([L])+m([M])+m([L,M])
     +m([L,H])+m([M,H])+m(S)=0.9.
 v~([L,H])
   = m([L])+m([H])+m([L,H])
    +m([L,M])+m([M,H])+m(S)=1.0.
 v~([M,H])
   = m([M])+m([H])+m([M,H])
    +m([L,M])+m([L,H])+m(S)=1.0.
 v~(S)
   = m([])+m([L])+m([M])+m([H])
    +m([L,M])+m([L,H])+m([M,H])
    +m([L,M,H])=1=v(S).
*/

% a brief explanation for TV set example.
% ---------------------------------------------------
/*
A story: Suppose you are the sales person of a TV set.
This product has the possible quality classes, qH, qM, or qL,
each of which brings about earning size in accordance with this 
order of quality. 

You may use the following uncertain evidences of eight types 
listed below in order to predict the purchase behavior of 
customer you now supposed to faced.

% Although it will not be utilized hereafter in this version, ... 
% available evidences and law of movement.

% (E1) share w.r.t. quality: 
current_share_by_quality(qH, 0.1).
current_share_by_quality(qM, 0.3).
current_share_by_quality(qL, 0.6).

% (E2) possible behavior conditionalized on current state.
possible_purchase_of_customer(qH / [qH ]).
possible_purchase_of_customer(qM / [qL, qM]).
possible_purchase_of_customer(qL / [qL]).

% (E3)--(E8) boundary conditions of possible sales.
possible_sales_share([qH], lower(0.1), upper(0.4)).
possible_sales_share([qM], lower(0), upper(0.9)).
possible_sales_share([qL], lower(0),upper(0.6)).
possible_sales_share([qM,qH], lower(0.4),upper(1)).
possible_sales_share([qL,qM], lower(0.6),upper(0.9)).
possible_sales_share([qL,qH], lower(0.1),upper(1)).

possible_sales_share(E, upper(Ub), lower(Lb)):-
   sevent(E),
   findall((Uq,Lq),
     (
      member(Q,E),
      maximum_sales_share_of(Q,Uq),
      minimum_sales_share_of(Q,Lq)
     ),
   D),
   findall(Uq,member((Uq,_),D),Us),
   sum(Us,Ub),
   findall(Lq,member((_,Lq),D),Ls),
   sum(Ls,Lb).

maximum_sales_share_of(Q,Y):-
   state(Q),
   findall(B,
     (
      possible_purchase_of_customer(Q/C),
      member(X,C),
      current_share_by_quality(X,B)
     ),
   A),
   sum(A,Y).

minimum_sales_share_of(qL,0).
minimum_sales_share_of(qM,0).
minimum_sales_share_of(qH,B):-current_share_by_quality(qH,B).
*/


% example 4. contingency trade (Mukerji,1997)
%---------------------------------------------------

% the default mode: states/1, bpa0/2, payoff0/3, act/1.

states([s1,s2,s3,s4]).
primitive_states([w1,w2,w3,w4,w5]).
prior(w1,0.1).
prior(w2,0.2).
prior(w3,0.3).
prior(w4,0.1).
prior(w5,0.3).
know(at(W),map(H),E):-implicate(W,H),sevent(E),subset_of(H,E,_).
%
bpa0(S,P):-sevent(S),implicate(W,S),prior(W,P).
bpa0(S,0):-sevent(S),\+ implicate(_W,S).

/*
% table 2 of Mukerji(1997, p.29).
% to verify bel/3 and sbel/3. 
bel0([],0).
bel0([s1],0.1).
bel0([s2],0).
bel0([s3],0.1).
bel0([s4],0.3).
bel0([s1,s2],0.1).
bel0([s1,s3],0.2).
bel0([s1,s4],0.4).
bel0([s2,s3],0.3).
bel0([s2,s4],0.3).
bel0([s3,s4],0.4).
bel0([s1,s2,s3],0.4).
bel0([s1,s2,s4],0.4).
bel0([s2,s3,s4],0.6).
bel0([s1,s3,s4],0.5).
bel0([s1,s2,s3,s4],1).
*/

% the implication mapping.
implicate(w1,[s1]).
implicate(w2,[s2,s3]).
implicate(w3,[s4]).
implicate(w4,[s3]).
implicate(w5,[s1,s2,s3,s4]).

% inverse implication.
inverse_implication(X,Y):-
   event(X),
   findall(W,
     (
      implicate(W,S),wn(implicate(W,S)),
      subset(S,X)
     ),
   Ws),
   flatten(Ws,Ws1),
   sort(Ws1,Y).

/*
% sample execution.

?- inverse_implication([s1,s2,s4],B).
implicate(w1, [s1])
implicate(w2, [s2, s3])
implicate(w3, [s4])
implicate(w4, [s3])
implicate(w5, [s1, s2, s3, s4])

B = [w1, w3] ;

No
?- 
*/

% decision part of Mukerji's example
%------------------------------------

act(f).

payoff0(f,s1,10).
payoff0(f,s2,7).
payoff0(f,s3,4).
payoff0(f,s4,15).

act(mukerji,f).

payoff0(mukerji,f,s1,10).
payoff0(mukerji,f,s2,7).
payoff0(mukerji,f,s3,4).
payoff0(mukerji,f,s4,15).


%-----------------------------------------
% sample executions for Mukerji's example.
%-----------------------------------------

/*

?- payoff(expected_conservative(F),P,Q).
w1, 0.1, 10
w2, 0.2, 4
w3, 0.3, 15
w4, 0.1, 4
w5, 0.3, 4

F = f
P = [0.1*10, 0.2*4, 0.3*15, 0.1*4, 0.3*4]
Q = 7.9 

Yes
?- payoff(choquet(A),B,C).

A = [s1, s2, s3, s4]
B1 = []
B2 = [0]
B3 = [0]
C = 0 ;

A = [s1, s2, s4]
B1 = [s3]
B2 = [4, 0]
B3 = [ (4-0)*1, 0]
C = 4 ;

A = [s1, s4]
B1 = [s2, s3]
B2 = [7, 4, 0]
B3 = [ (7-4)*0.4, (4-0)*1, 0]
C = 5.2 ;

A = [s4]
B1 = [s1, s2, s3]
B2 = [10, 7, 4, 0]
B3 = [ (10-7)*0.4, (7-4)*0.4, (4-0)*1, 0]
C = 6.4 ;

A = []
B1 = [s4, s1, s2, s3]
B2 = [15, 10, 7, 4, 0]
B3 = [ (15-10)*0.3, (10-7)*0.4, (7-4)*0.4, (4-0)*1, 0]
C = 7.9 ;

No
?- payoff(choquet1(f),[A,B1,B2,B3],C).

A = [s1, s2, s3, s4]
B1 = []
B2 = [0]
B3 = [0]
C = 0 ;

A = [s1, s2, s3]
B1 = [s4]
B2 = [15, 0]
B3 = [15* (0.3-0), 0]
C = 4.5 ;

A = [s2, s3]
B1 = [s4, s1]
B2 = [10, 15, 0]
B3 = [10* (0.4-0.3), 15* (0.3-0), 0]
C = 5.5 ;

A = [s3]
B1 = [s4, s1, s2]
B2 = [7, 10, 15, 0]
B3 = [7* (0.4-0.4), 10* (0.4-0.3), 15* (0.3-0), 0]
C = 5.5 ;

A = []
B1 = [s4, s1, s2, s3]
B2 = [4, 7, 10, 15, 0]
B3 = [4* (1-0.4), 7* (0.4-0.4), 10* (0.4-0.3), 15* (0.3-0), 0]
C = 7.9 ;


No
?- 
*/



%---------------------------------------------------
% normal form game theory under uncertainty and 
% equilibrium in beliefs (Ref. 9-13).
%---------------------------------------------------
% added: Jan 2004.


% example 5. A symmetric two-stage prisoners dilemma
% (Dow and Werlang,1994) Ref. 9.
%---------------------------------------------------

act(ipd2,F):-
   strategy(ipd2,F,_).

payoff0(ipd2,F,S,P):-
   act(ipd2,F),
   game(ipd2,payoff, [F, S], [P, _]).


% the meaning of the states:
% strategies in two step prisoners dilemma game.

% four unconditional strategies.
strategy(ipd2,ff,[(1,f),(2,f->f),(2,c->f)]). 
strategy(ipd2,fc,[(1,f),(2,f->c),(2,c->c)]).
strategy(ipd2,cf,[(1,c),(2,f->f),(2,c->f)]).
strategy(ipd2,cc,[(1,c),(2,f->c),(2,c->c)]).
% four conditional strategies.
strategy(ipd2,w,[(1,f),(2,f->f),(2,c->c)]).
strategy(ipd2,x,[(1,f),(2,f->c),(2,c->f)]).
strategy(ipd2,y,[(1,c),(2,f->f),(2,c->c)]). % Tit For Tat.
strategy(ipd2,z,[(1,c),(2,f->c),(2,c->f)]).


%
% display figures (game trees)
% -----------------------------------------------------------  %

figure(K/G):-game(G,figure,K,Figure),
   forall(member(L,Figure),(nl,write(L))).

% payoff
game(pd,paramater,payoff(a,A)):- A is 1.25.  % A >1.
game(pd,paramater,payoff(b,B)):- B is -0.5. % B <0.

game(pd,payoff, [f, f], [0, 0]).
game(pd,payoff, [c, c], [1, 1]).
game(pd,payoff, [f, c], [A, B]):-
   game(pd,paramater,payoff(a,A)),
   game(pd,paramater,payoff(b,B)).
game(pd,payoff, [c, f], [B, A]):-
   game(pd,paramater,payoff(a,A)),
   game(pd,paramater,payoff(b,B)).

game(pd,figure,1,Figure):-
   Figure=[
     '%        f       c       ',
     '%    +-------+-------+   ',
     '%    |   0   |   a   |   ',
     '%  f |   0   |   b   l   ',
     '%    +-------+-------+   ',
     '%    |   b   |   1   |   ',
     '%  c |   a   |   1   l   ',
     '%    +-------+-------+   ',
     '% Fig. one-shot game of prisoners dilemma.'
   ].


% payoffs in two step prisoners dilemma game without discount.

game(ipd2,payoff, [A, B], [P1, P2]):-
   strategy(ipd2,A,ASP), 
   strategy(ipd2,B,BSP),
   ASP=[(1,A1),(2,f->_A2F),(2,c->_A2C)],
   BSP=[(1,B1),(2,f->_B2F),(2,c->_B2C)], 
   game(pd,payoff, [A1, B1], [X1, Y1]),
   member((2,B1->A2),ASP),
   member((2,A1->B2),BSP),
   game(pd,payoff, [A2, B2], [X2, Y2]),
   P1 is X1 + X2,
   P2 is Y1 + Y2.


/*
?- game(ipd2,payoff,A,P).

A = [ff, ff]
P = [0, 0] ;

A = [ff, fc]
P = [1.25, -0.5] ;

A = [ff, cf]
P = [1.25, -0.5] ;

A = [ff, cc]
P = [2.5, -1] ;

A = [ff, w]
P = [0, 0] ;

A = [ff, x]
P = [1.25, -0.5] ;

A = [ff, y]
P = [1.25, -0.5] ;

A = [ff, z]
P = [2.5, -1] ;

A = [fc, ff]
P = [-0.5, 1.25] 

(...)

*/

% the belief system for ipd2.
%----------------------------------------------------

states(ipd2,[ff,fc,cf,cc,w,x,y,z]).

event(ipd2,E):-
  states(ipd2,S),
  subset_of(E1,_,S),
  sort_by_list(E1,S,E).

bel_ipd2([],0).
bel_ipd2([cf],0.4).
bel_ipd2(E,0):-
   member(E,[[w],[cc],[ff],[ff,cc],[cc,w],[ff,cc,w]]).
bel_ipd2(E,0.4):-
   member(E,[[ff,cf],[cf,w],[cf,cc]]).
bel_ipd2(E,0.4):-
   member(E,[[ff,cf,w],[ff,cf,cc]]).
bel_ipd2([cc,cf,w],0.8).
bel_ipd2([ff,cf,cc,w],1).
bel_ipd2(S,1):-states(ipd2,S).

bel0(ipd2,B,P):-
   event(ipd2,B),
   intersection(B,[ff,cf,cc,w],C),
   bel_ipd2_pattern(C,P).

bel_ipd2_pattern(C,P):-
   bel_ipd2(C1,P),
   seteq(C,C1),
   !.
bel_ipd2_pattern(C,0):-
   bel_ipd2(D,0),
   subset(C,D),
   !.



% sample executions as for Dow and Werlang's IPD game.
%---------------------------------------------------------

/*

?- bel0(ipd2,[A],P),P>0.

A = cf
P = 0.4 ;

No
?- length(A,8),bel0(ipd2,A,P),P>0.

A = [ff, fc, cf, cc, w, x, y, z]
P = 1 ;

No
*/



% -----------------------------------------------------------  %
%  Choquet integral (CEU) under belief functions  
% -----------------------------------------------------------  %
% privious version: belief01.pl(29 Jun 2003)
% modified: 20 Jan 2004.

% a specified model: the ranking of the outcomes of ipd2.

ranked_states_in_ipd2(2,strategy(A),replies(BL),payoffs(P1L)):-
   strategy(ipd2,A,_), 
   findall((P1,B),game(ipd2,payoff, [A, B], [P1, _]),R),
   sort(R,R1),
   reverse(R1,Q),
   findall(B,member((P1,B),Q),BL),
   findall(P1,member((P1,B),Q),P1L).


% sample execution
%---------------------------------------------------------

/*

?- ranked_states_in_ipd2(2,strategy(A),replies(BL),payoffs(P1L)).

A = ff
BL = [z, cc, y, x, fc, cf, w, ff]
P1L = [2.5, 2.5, 1.25, 1.25, 1.25, 1.25, 0, 0] ;

A = fc
BL = [z, cc, x, fc, y, cf, w, ff]
P1L = [2.25, 2.25, 1, 1, 0.75, 0.75, -0.5, -0.5] ;

A = cf
BL = [y, cc, z, cf, w, fc, x, ff]
P1L = [2.25, 2.25, 1, 1, 0.75, 0.75, -0.5, -0.5] ;

A = cc
BL = [y, cc, z, w, fc, cf, x, ff]
P1L = [2, 2, 0.5, 0.5, 0.5, 0.5, -1, -1] ;

A = w
BL = [z, cc, x, fc, y, cf, w, ff]
P1L = [2.25, 2.25, 1.25, 1.25, 0.75, 0.75, 0, 0] ;

A = x
BL = [z, cc, y, cf, x, fc, w, ff]
P1L = [2.5, 2.5, 1.25, 1.25, 1, 1, -0.5, -0.5] ;

A = y
BL = [y, cc, w, fc, z, cf, x, ff]
P1L = [2, 2, 0.75, 0.75, 0.5, 0.5, -0.5, -0.5] 

A = z
BL = [y, cc, z, cf, w, fc, x, ff]
P1L = [2.25, 2.25, 1, 1, 0.5, 0.5, -1, -1] ;

No

*/



/*********************************/
/*    model 1 :                  */
/*      ranking and cumulating   */
/*      by states                */
/*********************************/


%  outcomes ranked by states 
%-----------------------------------------------------
%  The following programs cause multiple cumulative path.

payoff(ranked(F),SS,[]):-
   act(F),
   states(SS).

payoff(ranked(F),Remains,[(X,S)|Z]):-
   act(F),
   payoff(ranked(F),R1,Z),
   (R1=[]->!,fail;true),
   state(S),
   payoff0(F,S,X),
   member(S,R1),
   \+ (
     member(S1,R1),
     payoff0(F,S1,X1),
     X1 < X
   ),
   subtract(R1,[S],Remains).


% CEU : Choquet integral representation of EU under uncertainty.
%-------------------------------------------------------------

payoff(choquet(F),[Y,Z,W],E1):-
   act(F),
   payoff(choquet_1(F),[[],Y,Z,W],E1).

payoff(choquet_1(F),[SS,[],[0],[0]],0):-
   act(F),
   states(SS).
payoff(choquet_1(F),[Remains,Y1,[X|Z],[V|W]],E):-
   act(F),
   payoff(choquet_1(F),[R1,Y,Z,W],E1),
   (R1=[]->!,fail;true),
   findall(S,
     max(X,
       (
        state(S),
        member(S,R1),
        payoff0(F,S,X)
       )
     ),
   A),
   subtract(R1,A,Remains),
   A=[T|_],
   payoff0(F,T,X),
   bel(Y,_,B0),
   append(Y,A,Y1),
   bel(Y1,_,B1),
   V = X * (B1-B0),
   %nl,write(X*(bel(Y1)-bel(Y))),write('='),write(V),
   E is E1 + V.

% CEU : an equivalent formalization.
%---------------------------------

payoff(choquet0(F),[Y,Z,W],E1):-
   act(F),
   payoff(choquet_2(F),[[],Y,Z,W],E1).

payoff(choquet_2(F),[SS,[],[0],[0]],0):-
   act(F),
   states(SS).
payoff(choquet_2(F),[Remains,[S|Y],[X|Z],[V|W]],E):-
   act(F),
   payoff(choquet_2(F),[R1,Y,Z,W],E1),
   (R1=[]->!,fail;true),
   state(S),
   payoff0(F,S,X),
   member(S,R1),
   \+ (
     member(S1,R1),
     payoff0(F,S1,X1),
     X1 < X
   ),
   subtract(R1,[S],Remains),
   bel(R1,_,B),
   Z = [Z1|_],
   V = (X-Z1) * B,
   %nl,write('dX*_bel'(R1)=(X-Z1)*B),write('='(V)),
   E is E1 + V.


% conservative (and cautious DM's) extension of act 
% and its EU representation (Mukerji, p.41).
%-----------------------------------------------------
payoff(conservative(F),W,X):- 
   act(F),
   implicate(W,H),
   findall(Y,
     (
      member(S,H),
      payoff0(F,S,Y)
     ),
   Z),
   min_of(X,Z).


payoff(expected_conservative(F),Es,X):- 
   act(F),
   findall((W,Q,Z),
     (
      prior(W,Q),
      payoff(conservative(F),W,Z),wn((W,Q,Z))
     ),
   D),
   findall(E,
     (
      member((W,Q,Z),D),
      E = Q * Z
     ),
   Es),
   sum(Es,X).


% sample execution for ipd2
%---------------------------------------------------------

/*

?- payoff(choquet(A),B,C).

A = ff
B = [[cc, z, cf, fc, x, y, ff, w], [0, 1.25, 2.5, 0], [0* (1-0.4), 1.25* (0.4-0), 2.5* (0-0), 0]]
C = 0.5 ;

A = fc
B = [[cc, z, fc, x, cf, y, ff, w], [-0.5, 0.75, 1, 2.25, 0], [-0.5* (1-0.4), 0.75* (0.4-0), 1* (0-0), 2.25* (0-0), 0]]
C = 5.55112e-017 ;

A = cf
B = [[cc, y, cf, z, fc, w, ff, x], [-0.5, 0.75, 1, 2.25, 0], [-0.5* (1-0.8), 0.75* (0.8-0.4), 1* (0.4-0), 2.25* (0-0), 0]]
C = 0.6 ;

A = cc
B = [[cc, y, cf, fc, w, z, ff, x], [-1, 0.5, 2, 0], [-1* (1-0.8), 0.5* (0.8-0), 2* (0-0), 0]]
C = 0.2 ;

A = w
B = [[cc, z, fc, x, cf, y, ff, w], [0, 0.75, 1.25, 2.25, 0], [0* (1-0.4), 0.75* (0.4-0), 1.25* (0-0), 2.25* (0-0), 0]]
C = 0.3 ;

A = x
B = [[cc, z, cf, y, fc, x, ff, w], [-0.5, 1, 1.25, 2.5, 0], [-0.5* (1-0.4), 1* (0.4-0.4), 1.25* (0.4-0), 2.5* (0-0), 0]]
C = 0.2 ;

A = y
B = [[cc, y, fc, w, cf, z, ff, x], [-0.5, 0.5, 0.75, 2, 0], [-0.5* (1-0.8), 0.5* (0.8-0), 0.75* (0-0), 2* (0-0), 0]]
C = 0.3 ;

A = z
B = [[cc, y, cf, z, fc, w, ff, x], [-1, 0.5, 1, 2.25, 0], [-1* (1-0.8), 0.5* (0.8-0.4), 1* (0.4-0), 2.25* (0-0), 0]]
C = 0.4 ;

No
?- 

*/

/*********************************/
/*    model 2 :                  */
/*      ranking and cumulating   */
/*      by events                */
/*********************************/


ranked_states(act(A),states(SL),payoffs(PL)):-
   act(A), 
   setof((P,S),P^S^payoff0(A,S,P),R),
   reverse(R,R1),
   findall(S,member((P,S),R1),SL),
   findall(P,member((P,S),R1),PL).



leveled_event(act(A),event(E),payoff(P),rank(K/N)):-
   act(A), 
   setof(P,S^payoff0(A,S,P),R),
   length(R,N),
   reverse(R,R1),
   nth1(K,R1,P),
   findall(S,payoff0(A,S,P),E).


% sample execution for ipd2 
%---------------------------------------------------------

/*

?- set_model(ipd2).

ipd2, Would you like to set this model ? (y) >y.

Yes
?- ranked_states(A,B,C).

A = act(ff)
B = states([y, cc, z, w, fc, cf, x, ff])
C = payoffs([2.5, 2.5, 1.25, 1.25, 1.25, 1.25, 0, 0]) ;

A = act(fc)
B = states([y, cc, w, fc, z, cf, x, ff])
C = payoffs([2.25, 2.25, 1, 1, 0.75, 0.75, -0.5, -0.5]) 

Yes
?- leveled_event(A,B,C,D).

A = act(ff)
B = event([cc, z])
C = payoff(2.5)
D = rank(1/3) ;

A = act(ff)
B = event([fc, cf, x, y])
C = payoff(1.25)
D = rank(2/3) ;

A = act(ff)
B = event([ff, w])
C = payoff(0)
D = rank(3/3) ;

A = act(fc)
B = event([cc, z])
C = payoff(2.25)
D = rank(1/4) 

Yes
?- 
*/


% pessimistic agent.

ceu_payoff(act(F),rank(1/N),[A,[A],[X],[B],[X*B]],V):-
   act(F),
   % the best outcome.
   leveled_event(act(F),event(A),payoff(X),rank(1/N)),
   bel(A,_,B),
   V is X * B.

ceu_payoff(act(F),rank(K/N),[Y,[A|O],[X|Z],[B|C],[V|U]],E):-
   act(F),
   ceu_payoff(act(F),rank(K0/N),[W,O,Z,C,U],E0),
   (K0=N->!,fail;K is K0 + 1),
   leveled_event(act(F),event(A),payoff(X),rank(K/N)),
   C=[B0|_],
   append(A,W,Y),
   bel(Y,_,B),
   V = X * (B-B0),
   E is E0 + V.

test_ceu(Act,[Levels,Events,Payoffs,Beliefs,Terms_of_CEU],CEU_Value):-
   act(Act),
   ceu_payoff(
     act(Act),
     rank(Levels/Levels),
     [_,Events,Payoffs,Beliefs,Terms_of_CEU],
     CEU_Value
   ).


% optimistic agent.

dceu_payoff(act(F),rank(N/N),[A,[A],[X],[B],[X*B]],V):-
   act(F),
   % the best outcome.
   leveled_event(act(F),event(A),payoff(X),rank(N/N)),
   bel(A,_,B),
   V is X * B.

dceu_payoff(act(F),rank(K/N),[Y,[A|O],[X|Z],[B|C],[V|U]],E):-
   act(F),
   dceu_payoff(act(F),rank(K0/N),[W,O,Z,C,U],E0),
   (K0=1->!,fail;K is K0 - 1),
   leveled_event(act(F),event(A),payoff(X),rank(K/N)),
   C=[B0|_],
   append(A,W,Y),
   bel(Y,_,B),
   V = X * (B-B0),
   E is E0 + V.

test_dceu(Act,[Levels,Events,Payoffs,Beliefs,Terms_of_CEU],CEU_Value):-
   act(Act),
   dceu_payoff(
     act(Act),
     rank(1/Levels),
     [_,Events,Payoffs,Beliefs,Terms_of_CEU],
     CEU_Value
   ).



%---------------------------------------------------
%  test run: computing CEU for examples 
%---------------------------------------------------

/*

% CEU for the default model (Mukerji's trade example):
%---------------------------------------------------

?- test_ceu(Act,[Levels,Events,Payoffs,Beliefs,Terms_of_CEU],CEU_Value).

Act = f
Levels = 4
Events = [[s3], [s2], [s1], [s4]]
Payoffs = [4, 7, 10, 15]
Beliefs = [1, 0.4, 0.4, 0.3]
Terms_of_CEU = [4* (1-0.4), 7* (0.4-0.4), 10* (0.4-0.3), 15*0.3]
CEU_Value = 7.9 

Yes

% CEU for the ipd2 model(Dow and Werlang's game example):
%---------------------------------------------------
% Cf., see remarks for the proposition 5.1 in Groes et al.(1998, p.56).


?- set_model(ipd2).

ipd2, Would you like to set this model ? (y) >y.

Yes
?- test_ceu(Act,[Levels,Events,Payoffs,Beliefs,Terms_of_CEU],CEU_Value).

Act = ff
Levels = 3
Events = [[ff, w], [fc, cf, x, y], [cc, z]]
Payoffs = [0, 1.25, 2.5]
Beliefs = [1, 0.4, 0]
Terms_of_CEU = [0* (1-0.4), 1.25* (0.4-0), 2.5*0]
CEU_Value = 0.5 ;

Act = fc
Levels = 4
Events = [[ff, w], [cf, y], [fc, x], [cc, z]]
Payoffs = [-0.5, 0.75, 1, 2.25]
Beliefs = [1, 0.4, 0, 0]
Terms_of_CEU = [-0.5* (1-0.4), 0.75* (0.4-0), 1* (0-0), 2.25*0]
CEU_Value = 5.55112e-017 ;

Act = cf
Levels = 4
Events = [[ff, x], [fc, w], [cf, z], [cc, y]]
Payoffs = [-0.5, 0.75, 1, 2.25]
Beliefs = [1, 0.8, 0.4, 0]
Terms_of_CEU = [-0.5* (1-0.8), 0.75* (0.8-0.4), 1* (0.4-0), 2.25*0]
CEU_Value = 0.6 ;

Act = cc
Levels = 3
Events = [[ff, x], [fc, cf, w, z], [cc, y]]
Payoffs = [-1, 0.5, 2]
Beliefs = [1, 0.8, 0]
Terms_of_CEU = [-1* (1-0.8), 0.5* (0.8-0), 2*0]
CEU_Value = 0.2 ;

Act = w
Levels = 4
Events = [[ff, w], [cf, y], [fc, x], [cc, z]]
Payoffs = [0, 0.75, 1.25, 2.25]
Beliefs = [1, 0.4, 0, 0]
Terms_of_CEU = [0* (1-0.4), 0.75* (0.4-0), 1.25* (0-0), 2.25*0]
CEU_Value = 0.3 ;

Act = x
Levels = 4
Events = [[ff, w], [fc, x], [cf, y], [cc, z]]
Payoffs = [-0.5, 1, 1.25, 2.5]
Beliefs = [1, 0.4, 0.4, 0]
Terms_of_CEU = [-0.5* (1-0.4), 1* (0.4-0.4), 1.25* (0.4-0), 2.5*0]
CEU_Value = 0.2 ;

Act = y
Levels = 4
Events = [[ff, x], [cf, z], [fc, w], [cc, y]]
Payoffs = [-0.5, 0.5, 0.75, 2]
Beliefs = [1, 0.8, 0, 0]
Terms_of_CEU = [-0.5* (1-0.8), 0.5* (0.8-0), 0.75* (0-0), 2*0]
CEU_Value = 0.3 ;

Act = z
Levels = 4
Events = [[ff, x], [fc, w], [cf, z], [cc, y]]
Payoffs = [-1, 0.5, 1, 2.25]
Beliefs = [1, 0.8, 0.4, 0]
Terms_of_CEU = [-1* (1-0.8), 0.5* (0.8-0.4), 1* (0.4-0), 2.25*0]
CEU_Value = 0.4 ;

No
?- 
*/


%---------------------------------------------------
% as if EU : assuming the additivity erroneously.
%---------------------------------------------------

as_if_eu_payoff(act(F),rank(1/N),[A,[A],[X],[B],[X*B]],V):-
   act(F),
   % the best outcome.
   leveled_event(act(F),event(A),payoff(X),rank(1/N)),
   bel(A,_,B),
   V is X * B.

as_if_eu_payoff(act(F),rank(K/N),[Y,[A|O],[X|Z],[B|C],[V|U]],E):-
   act(F),
   as_if_eu_payoff(act(F),rank(K0/N),[W,O,Z,C,U],E0),
   (K0=N->!,fail;K is K0 + 1),
   leveled_event(act(F),event(A),payoff(X),rank(K/N)),
   append(A,W,Y),
   bel(A,_,B),
   V = X * B,
   E is E0 + V.

test_as_if_eu(Act,[Levels,Events,Payoffs,Beliefs,Terms_of_CEU],CEU_Value):-
   act(Act),
   as_if_eu_payoff(
     act(Act),
     rank(Levels/Levels),
     [_,Events,Payoffs,Beliefs,Terms_of_CEU],
     CEU_Value
   ).

/*

% test run for ipd2.
% -----------------------------------------------------------  %

?- test_as_if_eu(Act,[Levels,Events,Payoffs,Beliefs,Terms_of_CEU],CEU_Value).

Act = ff
Levels = 3
Events = [[ff, w], [fc, cf, x, y], [cc, z]]
Payoffs = [0, 1.25, 2.5]
Beliefs = [0, 0.4, 0]
Terms_of_CEU = [0*0, 1.25*0.4, 2.5*0]
CEU_Value = 0.5 ;

Act = fc
Levels = 4
Events = [[ff, w], [cf, y], [fc, x], [cc, z]]
Payoffs = [-0.5, 0.75, 1, 2.25]
Beliefs = [0, 0.4, 0, 0]
Terms_of_CEU = [-0.5*0, 0.75*0.4, 1*0, 2.25*0]
CEU_Value = 0.3 ;

Act = cf
Levels = 4
Events = [[ff, x], [fc, w], [cf, z], [cc, y]]
Payoffs = [-0.5, 0.75, 1, 2.25]
Beliefs = [0, 0, 0.4, 0]
Terms_of_CEU = [-0.5*0, 0.75*0, 1*0.4, 2.25*0]
CEU_Value = 0.4 ;

Act = cc
Levels = 3
Events = [[ff, x], [fc, cf, w, z], [cc, y]]
Payoffs = [-1, 0.5, 2]
Beliefs = [0, 0.4, 0]
Terms_of_CEU = [-1*0, 0.5*0.4, 2*0]
CEU_Value = 0.2 ;

Act = w
Levels = 4
Events = [[ff, w], [cf, y], [fc, x], [cc, z]]
Payoffs = [0, 0.75, 1.25, 2.25]
Beliefs = [0, 0.4, 0, 0]
Terms_of_CEU = [0*0, 0.75*0.4, 1.25*0, 2.25*0]
CEU_Value = 0.3 ;

Act = x
Levels = 4
Events = [[ff, w], [fc, x], [cf, y], [cc, z]]
Payoffs = [-0.5, 1, 1.25, 2.5]
Beliefs = [0, 0, 0.4, 0]
Terms_of_CEU = [-0.5*0, 1*0, 1.25*0.4, 2.5*0]
CEU_Value = 0.5 ;

Act = y
Levels = 4
Events = [[ff, x], [cf, z], [fc, w], [cc, y]]
Payoffs = [-0.5, 0.5, 0.75, 2]
Beliefs = [0, 0.4, 0, 0]
Terms_of_CEU = [-0.5*0, 0.5*0.4, 0.75*0, 2*0]
CEU_Value = 0.2 ;

Act = z
Levels = 4
Events = [[ff, x], [fc, w], [cf, z], [cc, y]]
Payoffs = [-1, 0.5, 1, 2.25]
Beliefs = [0, 0, 0.4, 0]
Terms_of_CEU = [-1*0, 0.5*0, 1*0.4, 2.25*0]
CEU_Value = 0.4 ;

No
?- 
*/


%--------------------------------------------------
% support of (nonadditive-)probability measure.
%--------------------------------------------------
% added: 20-22 Jan 2004.
%  definition of support by Dow and Werlang(1994).
%  Ref. 9. Also in Ref. 12.
%  supp(P):=A s.t. P(A*)=0 and for all subset B of A, P(B*)>0.
%  where C* denotes complement of event C. i.e.,
%  the minimal event event whose complement is measure zero. 

proper_subset_of(S,Y):-
   var(Y),
   sevent(Y),
   subset_of(S,_,Y),
   Y \= S.

proper_subset_of(S,Y):-
   \+ var(Y),
   subset_of(S,_,Y),
   Y \= S.

support0(S):-
   findall(X,
     (
      sbel([X],_,B),
      B>0
     ),
   S).

event_whose_complements_is_measure_zero(S):-
   states(W),
   sbel(Y,_,0),
   subtract(W,Y,S).

support1(S):-
   states(W),
   event_whose_complements_is_measure_zero(S),
   forall(
     proper_subset_of(Y,S),
     (
      subtract(W,Y,X),
      bel(X,_,B),
      B>0
     )
   ).

support(S):-
   states(W),
   event_whose_complements_is_measure_zero(S),
   \+ (
     proper_subset_of(Y,S),
     subtract(W,Y,X),  % complementary event of Y
     bel(X,_,0)
   ).

/*
?- support(S).

S = [s1, s3, s4] ;

No
?- set_model(ipd2).

ipd2, Would you like to set this model ? (y) >y.

Yes
?- support(S).

S = [cf] ;

No
?-

*/


%---------------------------------------------------
%   Beliefs Equilibrium 
%---------------------------------------------------
% added: 23 Jan 2004.
% Ref. 9 and 12.


%  best_response/2: nash strategy under ambiguous beliefs. 

best_response(BR):-
   findall(A,
     max(U,test_ceu(A,_,U)),
   BR).


%  bliefs_equilibrium/2: equilibrium of the game in beliefs. 

% definition (Equilibrium in Beliefs) Ref.s 9 and 12.
% support(B) is a subset of intersection of best responses. 

equilibrium_in_beliefs(Sup,BR,Out):-
   support(Sup),
   best_response(BR),
   subtract(Sup,BR,Out).


/*

?- set_model(ipd2).

ipd2, Would you like to set this model ? (y) >y.

Yes
?- best_response(A).

A = [cf] ;

No
?- equilibrium_in_beliefs(A,B,C).

A = [cf]
B = [cf]
C = [] ;

No
?- 

*/


%---------------------------------------------------
%   Part I.   Belief Function
%---------------------------------------------------
%  privious version: belief01.pl(29 Jun 2003)
%  modified: 19--23 Jan 2004.


% probabilistic state space: state space, event
%---------------------------------------------------

state(S):- states(A),member(S,A).

event(E):-
   states(A),
  %subset_of(E,_,A).   % to take the easy way out for computation.
   bag1(E,A,_N).

% modified: 15 Jan 2004. to use sort_by_list/3 instead of sort/2.

sevent(X):-
   event(X),
   states(S),
   sort_by_list(X,S,X).

all_event(S):-setof(E,event(E),S).

all_sevent(S):-setof(E,sevent(E),S).

cevent(E,X):-
   complementary_event(E,X).

complementary_event(E,X):-
   states(A),
   event(X),
   subtract(A,X,C),
   sevent(E),
   seteq(C,E).



% Basic probability assignment (bpa) 
% computed from database bpa0.
%---------------------------------------------------
bpa(E,B):-
   event(E),
   bpa0(E1,B),
   seteq(E,E1).
bpa(E,0):-
   event(E),
   \+ (
     bpa0(E1,_),
     seteq(E,E1)
   ).


% Belief functions
%---------------------------------------------------
% modified: 14 Jan 2004. bug fix. (bel0/3-->bel0/2 and to use event/1.)

bel(E,X,X):-
   (clause(bel0(_,_),_)->true;fail),
   event(E),
   bel0(E1,X),
   seteq(E1,E).

% modified: 29 Jun 2003. bug fix. (setof was used.)

bel(E,Xq,X):-
   (clause(bel0(_,_),_)->fail;true),
   (clause(bpa0(_,_),_)->true;fail),
   event(E),
   bagof(B,
     F^(
       sevent(F),
       subset(F,E),
       bpa(F,B)
       %,nl,write(bpa(F,B))
     ),
   G),
   sum_eq(G,Xq,X).

% belief function for sorted event.
sbel(E,Xq,X):-
   sevent(E),
   bel(E,Xq,X).

% Possibility function (conjugate belief function)
%---------------------------------------------------
pos(E,Xq,X):-
   (clause(pos0(_,_),_)->fail;true),
   (clause(bpa0(_,_),_)->true;fail),
   event(E),
   findall(B,
     (
      sevent(F),
      intersection(F,E,M),
      M \= [],
      bpa(F,B)
     ),
   G),
   sum_eq(G,Xq,X).

spos(E,Xq,X):-
   sevent(E),
   pos(E,Xq,X).

% a comparison.
pos_1(E,1-Xq,X):-
   event(E),
   cevent(C,E),
   bel(C,Xq,_),
   X is 1 - Xq.

% a verification program as for v(A)=1-v~(-A).
%
% ?- sbel(E,_,B),cevent(F,E),pos_1(F,1-X,_),Y is X.


% b.p.a. via Mevious inversion
%---------------------------------------------------
%  mass(B) = sum(for(subset(A,B)), (-1)^|B-A| * v(A)).

mass([],0,0).
mass(E,Yq,Y):-
   event(E),
   E \= [],
   findall(A,
     (
      sevent(F),
      subset_of(F,_,E),
      bel(F,Bq,_B),
      %nl,write(bel(F,Bq,_B)),
      mevious(F,E,K),
      %nl,write(mevious(F,E,K)),
      A = K * Bq
     ),
   G),
   sum_eq(G,Yq,Y1),
   Y is Y1.

% Coefficients in the Mevious inversion formula.
%---------------------------------------------------
mevious(X,Y,Z):-
   event(X),
   event(Y),
   subtract(Y,X,W),
   length(W,M),
   Z = (-1)^M.



%---------------------------------------------------
% Checking the modularity (i.e., 2-monotonicity).
%---------------------------------------------------
% sevent by using sort_by_list of the all states. 

pairwise_element_of([K,K1],B,[A1,A2]):-
   nth1(K,B,A1),
   nth1(K1,B,A2),
   K1 >= K.
   
pairwise_sevent(E,F):-
   setof(A,sevent(A),B),
   pairwise_element_of(_K,B,[E,F]).
   %F\=E,

:- dynamic temp_pairwise_sevent/2.

generate_all_pairwise_sevents:-
   confirm_user_1,
   abolish(temp_pairwise_sevent/2),
   pairwise_sevent(E,F),
   assert(temp_pairwise_sevent(E,F)),
   fail.

generate_all_pairwise_sevents:-
   nl,write('Now you hold data. Continue run ? It may cost minutes.(y)>'),
   read(y).

confirm_user_1:-
   nl,write('Abolish old data and create new ? (y)>'),
   read(y).

union_and_intesection_of_pairwise_sevent(E,F,G,H):-
   union_of_pairwise_sevent(E,F,G),
   intersection_of_pairwise_sevent(E,F,H).

union_of_pairwise_sevent(E,F,G):-
   union(E,F,G1),
   seteq(G1,G),
   sevent(G).

intersection_of_pairwise_sevent(E,F,H):-
   intersection(E,F,H1),
   seteq(H1,H),
   sevent(H).


%----------------------------------------
%  test of modularity for pairwise event 
%----------------------------------------

modularity([E,F],[G,H],REL,Z):-
   (var(E) ->  pairwise_sevent(E,F); true),
   sbel(E,_Yq1,Y1),
   sbel(F,_Yq2,Y2),
   union_and_intesection_of_pairwise_sevent(E,F,G,H),
   sbel(G,_Yq3,Y3),
   sbel(H,_Yq4,Y4),
   X1 is Y1 + Y2,
   X2 is Y3 + Y4,
   case_of_modularity1(X1,X2,REL0,Z),
   REL0 =.. [OP,X1,X2],
   REL =.. [OP, Y1 + Y2, Y3 + Y4].

case_of_modularity1(X1,X2, REL, modular):-
   num_eq(X1,X2),!, REL=(X1 = X2).
case_of_modularity1(X1,X2, REL, supermodular):-
   \+ num_eq(X1,X2), X1X2,!, REL=(X1 > X2).

num_eq(X,Y):-
   Z is (X - Y)^2,
   Z < 10^(-10).

:- dynamic temp_modularity/4.

update_temp_modularity([E,F],UI,R,Y):-
   modularity([E,F],UI,R,Y),
   retractall(temp_modularity([E,F],_,_,_)),
   assert(temp_modularity([E,F],UI,R,Y)).


%---------------------------------
%  total test of modularity
%---------------------------------

modularity(Z):-
   generate_all_pairwise_sevents,
   findall(Y,
     (
      temp_pairwise_sevent(E,F),
      update_temp_modularity([E,F],_,_,Y)
     ),
   G),
   %sort(G,G1),write(G1),
   case_of_modularity2(G,Z).

case_of_modularity2(G,modular):- sort(G,[modular]),!.
case_of_modularity2(G,supermodular):- \+ member(submodular,G),!.
case_of_modularity2(G,submodular):- \+ member(supermodular,G),!.
case_of_modularity2(_G,nonlinear).


% test run for default example.
%---------------------------------

/*

?- modularity(A).

Abolish old data and create new ? (y)>y.

Now you hold data. Continue run ? It may cost minutes.(y)>y.

A = supermodular 

Yes
?- 

*/


% rephrase into convexity.
convexity(E,F,G,H,I):-
   modularity(E,F,G,H,J),
   A = (concave,submodular),
   B = (convex,supermodular),
   C = (linear,modular),
   C1 = nonlinear,
   D = (I,J),
   member(D,[A,B,C,C1]).

convexity(concave):- modularity(submodular).
convexity(convex):- modularity(supermodular).
convexity(linear):- modularity(modular).


%---------------------------------
%  partial test of modularity
%---------------------------------

modularity(part,W,Z):-
   generate_pairwise_sevents_restricted_to_states(W),
   findall(Y,
     (
      temp_pairwise_sevent(E,F),
      update_temp_modularity([E,F],_,_,Y)
     ),
   G),
   %sort(G,G1),write(G1),
   case_of_modularity2(G,Z).

generate_pairwise_sevents_restricted_to_states(W):-
   event(W),
   confirm_user_2(W),
   abolish(temp_pairwise_sevent/2),
   pairwise_sevent_restricted_to_states(E,F,W),
   assert(temp_pairwise_sevent(E,F)),
   fail.

generate_pairwise_sevents_restricted_to_states(_W):-
   nl,write('Now you hold data. Continue run ? It may cost minutes.(y)>'),
   read(y).

confirm_user_2(E):-
   nl,
   write(
     (
      'selected states':E,
      'ok ? (y) >'
     )
   ),
   read(y),
   confirm_user_1.

pairwise_sevent_restricted_to_states(E,F,W):-
   event(W),
   setof(A,
     (
      sevent(A),
      forall(member(X,A),member(X,W))
     ),
   B),
   pairwise_element_of(_K,B,[E,F]).
   %F\=E,


% test run for ipd2 example.
%---------------------------------

/*

?- set_model(ipd2).

ipd2, Would you like to set this model ? (y) >y.

Yes
?- modularity(part,[ff,cf,cc,w],A).

selected states:[ff, cf, cc, w], ok ? (y) >y.

Abolish old data and create new ? (y)>y.

Now you hold data. Continue run ? It may cost minutes.(y)>y.

A = supermodular ;

No
?- 
*/

%---------------------------------------------------
% Checking additivity.
%---------------------------------------------------
additivity(E,F,Yq1 + Yq2,Yq3,Z):-
   sevent(E),
   sevent(F),
   F \= E,
   union(E,F,G1),
   sevent(G),
   seteq(G1,G),
   intersection(E,F,[]),
   bel(E,Yq1,_Y1),
   bel(F,Yq2,_Y2),
   bel(G,Yq3,_Y3),
   X1 is Yq1 + Yq2,
   X2 is Yq3,
   (num_eq(X1,X2) -> Z = additive; true),
   (X1 > X2 -> Z = subadditive; true),
   (X1 < X2 -> Z = superadditive; true).

additivity(Z):-
   findall(Y,
     (
      sevent(E),
      sevent(F),
      additivity(E,F,_,_,Y)
     ),
   G),%nl,write(G),
   (sort(G,[additive])->Z = additive; true),
   (sort(G,[additive,superadditive])->Z = superadditive; true),
   (sort(G,[additive,subadditive])->Z = subadditive; true),
   (Z = nonadditive; true).
   

% Checking monotonicity in event.
%---------------------------------------------------
monotonicity(E,F,Yq1,Yq2,Z):-
   sevent(E),
   sevent(F),
   F \= E,
   subset(E,F),
   bel(E,Yq1,Y1),
   bel(F,Yq2,Y2),
   (Y1 > Y2 -> Z = nonmonotonic; true),
   (Y1 =< Y2 -> Z = monotone; true).

monotonicity(Z):-
   findall(Y,
     (
      sevent(E),
      sevent(F),
      monotonicity(E,F,_,_,Y)
     ),
   G),%nl,write(G),
   (
    sort(G,[monotone])
     -> Z = monotone;
        Z = nonmonotonic
   ).

% definition of capacity.
is_capcity:-  monotonicity(monotone).


/*
% sample executions.

?- sevent(A),mass(A,B,C),bpa(A,D).

A = [r]
B = 1/3* -1^0
C = 0.333333
D = 1/3 ;

A = [b]
B = 0* -1^0
C = 0
D = 0 ;

A = [y]
B = 0* -1^0
C = 0
D = 0 ;

A = [b, r]
B = 1/3* -1^1+ (1/3+0)* -1^0+0* -1^1
C = 0
D = 0 ;

A = [r, y]
B = 1/3* -1^1+ (1/3+0)* -1^0+0* -1^1
C = 0
D = 0 ;

A = [b, y]
B = (2/3+0)* -1^0+0* -1^1
C = 0.666667
D = 2/3 


A = [b, r, y]
B = 1/3* -1^2+ (2/3+0)* -1^1+ (1/3+0)* -1^1+ (2/3+1/3+0)* -1^0+0* -1^2
C = 0.333333
D = 1/3 ;

No
?- sevent(A),cevent(D,A),pos(A,B,C).

A = [r]
D = [b, y]
B = 1/3+0
C = 0.333333 

Yes
?- sevent(A),cevent(D,A),bel(D,E,F),pos(A,B,C).

A = [r]
D = [b, y]
E = 2/3+0
F = 0.666667
B = 1/3+0
C = 0.333333 

Yes
?- 
*/

%---------------------------------------------------
%   Part II.   Updating Belief Function
%---------------------------------------------------


% Updating bpa by Dempster-Shafer rule
%---------------------------------------------------
mass_of_intersection_with(A,B,C,Yq):-
   % mass of A which has intersection C with B.
   event(A),
   event(B),
   intersection(A,B,D),
   sort(D,C),
   mass(A,Yq,_).

update_mass([]/D,0,0):-event(D).
update_mass(E/D,Yq,Y):-
   event(E),
   \+ E = [],
   event(D),
   \+ D = [],
   findall((C,M),
     (
      sevent(B),
      mass_of_intersection_with(B,D,C,M),
      C \= []
     ),
   X),
   %length(X,NX),
   %nl,write(find(NX,X)),
   findall(Q,
     member((_,Q),X),
   BX0),
   sum(BX0,M0),
   sort(E,E1),
   findall(Q,
     member((E1,Q),X),
   BX1),
   sum(BX1,M1),
   Yq = M1 /M0,
   Y is Yq.

% Updating bel by Dempster-Shafer rule
%---------------------------------------------------
update_bel(E/D,Xq,X):-
   event(E),
   event(D),
%   \+ D = [],
   findall(Yq,
     (
      sevent(F),
      subset(F,E),
      update_mass(F/D,Yq,_Y)
     ),
   G),
   sum_of_fractions(G,Xq,X,_).

% sum of fractions switched by case of zero divisor
sum_of_fractions(G,Xq,X,no):-
   \+ (member(_/D,G),0 is D),
   sum_eq(G,Xq,X1),
   X is X1.

sum_of_fractions(G,Xq,X,yes):-
   member(_/D,G),
   0 is D,
   Xq = sum(G),
   X = indefinite.

% updating bel only when an event is of sorted states.

update_sbel(E/D,Xq,X):-
   sevent(E),
   sevent(D),
%   \+ D = [],
   update_bel(E/D,Xq,X).


% a comparison.
%               v(union(A,-B)) - v(-B)
% v_ds(A|B) =  ------------------------
%               1 - v(-B)

update_bel_1(E/D,Xq,X):-
   event(E),
   event(D),
%   \+ D = [],
   cevent(F,D),
   union(E,F,G),
   bel(F,XFq,_),
   bel(G,XGq,_),
   Yq = (XGq - XFq),
   Zq = (1 - XFq),
   sum_of_fractions([Yq/Zq],Xq,X,_).

%  a verification program code like as below may be convenient to you 
%  in order to verify the equality of three versions of update_bel.
%
%  ?- update_sbel(E/D,_,C),(update_bel_2(E/D,_,F)->true;F=non).

% another comparison.
%               v~(B) - v~(intersect(-A,B))
% v_ds(A|B) =  ----------------------------
%               v~(B)

update_bel_2(E/D,Xq,X):-
   event(E),
   event(D),
%   \+ D = [],
   cevent(F,E),
   intersection(D,F,G),
   pos(D,YDq,_),
   pos(G,YGq,_),
   Yq = (YDq - YGq),
   sum_of_fractions([Yq/YDq],Xq,X,_).


% Updating pos by Dempster-Shafer rule
%---------------------------------------------------
update_pos(E/D,Xq,X):-
   event(E),
   event(D),
   \+ D = [],
   findall(Yq,
     (
       sevent(F),
       intersection(F,E,M),
       M \= [],
       update_mass(F/D,Yq,_Y)
     ),
   G),
   sum_of_fractions(G,Xq,X,_).

update_spos(E/D,Xq,X):-
   sevent(E),
   sevent(D),
%   \+ D = [],
   update_pos(E/D,Xq,X).

% a code for test:
%
% ?- update_sbel(E/D,_,C),_,F),
%      cevent(G,E),(update_pos(G/D,_,F)->true;F=non).


% comparison
%               v~(intersection(A,B))
% v~ds(A|B) =  ------------------------
%               v~(B)

update_pos_1(E/D,Xq,X):-
   event(E),
   event(D),
   \+ D = [],
   intersection(D,E,G),
   pos(G,Yed,_),
   pos(D,Yd,_),
   sum_of_fractions([Yed/Yd],Xq,X,_).

% a code for test:
%
% ?- update_spos(A,_,C),_,F),
%      (update_pos_1(A,_,F)->true;F=non).


% Updating bel by upper-lower probabilities conditioning
%----------------------------------------------------------------
/*
                 v(intersect(A,B))
  v_ul(A|B)  =  -----------------------------------------
                 v(intersect(A,B)) + v~(intersect(-A,B))
*/

update_bel_ul(E/D,Xq,X):-
   event(E),
   event(D),
   cevent(F,E),
   intersection(E,D,G),
   intersection(F,D,H),
   bel(G,YGq,_),
   pos(H,YHq,_),
   Yq = YGq,
   Zq = (YGq + YHq),
   sum_of_fractions([Yq/Zq],Xq,X,_).


% Updating pos by upper-lower probabilities conditioning
%----------------------------------------------------------------
/*
                 v~(intersect(A,B))
  v~ul(A|B)  =  -----------------------------------------
                 v~(intersect(A,B)) + v(intersect(-A,B))
*/

update_pos_ul(E/D,Xq,X):-
   event(E),
   event(D),
   cevent(F,E),
   intersection(E,D,G),
   intersection(F,D,H),
   pos(G,YGq,_),
   bel(H,YHq,_),
   Yq = YGq,
   Zq = (YGq + YHq),
   sum_of_fractions([Yq/Zq],Xq,X,_).

/*

% sample executions of updating for the example 1:

?- update_sbel(A/[b,y],_,B).

A = []
B = 0 ;

A = [r]
B = 0 ;

A = [b]
B = 0 ;

A = [y]
B = 0 ;

A = [b, r]
B = 0 ;

A = [r, y]
B = 0 ;

A = [b, y]
B = 1 ;

A = [b, r, y]
B = 1 ;

No
?- sevent(A),update_bel_ul(A/[b,y],_,B).

A = []
B = 0 ;

A = [r]
B = 0 ;

A = [b]
B = 0 ;

A = [y]
B = 0 ;

A = [b, r]
B = 0 ;

A = [r, y]
B = 0 ;

A = [b, y]
B = 1 ;

A = [b, r, y]
B = 1 ;

No
?- 

% Dempster-Shafer rule for the example 2 (Dow and Werlang,1992):

case 1)  H=[1]

v([1]|[1])
 = [v([1] U [1]c)-v([1]c)]/[1-v([1]c)] 
 = [v([1] U [2,3])-v([2,3])]/[1-v([2,3])]
 = (1-1/2)/(1-1/2)
 = 1
v([2]|[1])  = [v([2,3])-v([2,3])]/[1-v([2,3])]
 = 0.  
v([3]|[1])  = [v([2,3])-v([2,3])]/[1-v([2,3])]
 = 0.  
v([1,2]|[1]) = [v([1,2,3])-v([2,3])]/[1-v([2,3])]
 = 1.  
v([1,3]|[1]) = [v([1,2,3])-v([2,3])]/[1-v([2,3])]
 = 1.  
v([2,3]|[1]) = [v([2,3])-v([2,3])]/[1-v([2,3])]
 = 0.  
v([1,2,3]|[1]) = [v([1,2,3])-v([2,3])]/[1-v([2,3])]
 = 1.  
case 2)   H=[2,3]
v([1]|[2,3])  =[v([1] U [2,3]c)-v([2,3]c)] / [1-v([2,3]c)]
  = [ v([1] U [1])-v([1])]/[1-v([1])]
  =0.  
v([2]|[2,3]) = [v([1,2])-v([1])]/[1-v([1])]
  = 1/3.
v([3]|[2,3]) = [v([1,3])-v([1])]/[1-v([1])] 
  = 1/3.
v([1,2]|[2,3]) = [v([1,2])-v([1])]/[1-v([1])]
  = 1/3.
v([1,3]|[2,3]) = [v([1,3])-v([1])]/[1-v([1])]
  = 1/3.
v([2,3]|[2,3]) = [v([1,2,3])-v([1])]/[1-v([1])]
  = 1.
v([1,2,3]|[2,3]) = [v([1,2,3])-v([1])]/[1-v([1])]
  = 1.


*/

%
% -----------------------------------------------------------  %
%   Utilities for list operations
% -----------------------------------------------------------  %
%  edited: 14 Feb 2003. (imported from: set1.pl)
%  edited: 19 Jan 2004. case of generating equivalent set.
%
% equality for pair of set
% -----------------------------------------------------------  %
seteq(X,Y):-
   \+ var(X),
   length(X,N),
   \+ var(Y),
   length(Y,N),
   sort(X,Sort),
   sort(Y,Sort).

seteq(X,Y):-
   \+ var(X),
   length(X,N),
   var(Y),
   bag1(Y,X,N).


%
%
% 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,
   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([],_A,0).
bag0([C|B],A,N):-
   length([C|B],N),
   bag0(B,A,_N1),
   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],[_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).
%
% subset_of/3 : subset-enumeration 
% -----------------------------------------------------------  %
subset_of(A,N,As):-
   length(As,L),
   length(D,L),
   list_projection(D,As,B),
   length(B,N),
   sort(B,A).
%
% complementary list projection
%--------------------------------------------------
% added: 10 Jan 2003.
c_list_projection(X,Y,Z):-
   complement(X,XC,_N),
   list_projection(XC,Y,Z).

complement(X,XC,N):-
   \+ (var(X),var(N)),
   bag0(X,[1,0],N),
   zeros(Zero,N),
   ones(One,N),
   replace(X,Zero,One,XC).


%
% sort without removal of duplicates
%--------------------------------------------------
asort(A,B):-
   sort(A,C),
   bagof(CK,
     J^K^(
      nth1(J,C,CK),
      nth1(K,A,CK)
     ),
   B).

% sort by a list of ordering
%--------------------------------------------------
% added: 14 Jan 2004. 
% modified: 16 Jan 2004. 

sort_by_list(X,OL,Y):-
   (var(X)->bag1(X,OL,_);true),
   list_projection(_,OL,Y),
   seteq(X,Y).


% -----------------------------------------------------------  %
% Arithmetic 
% -----------------------------------------------------------  %
%
% evaluation of a nummerical value
% -----------------------------------------------------------  %
eval_number(X,X1):-
   X1 is X,
   number(X1).

% maximal solution for given goal clause : a naive solver 
%---------------------------------------------------------

max(X,Goal):-
  % X: the objective variable,
  % Goal: the objective function and constraints,
  setof((X,Goal),Goal,Z),
  member((X,Goal),Z),
  \+ (
    member((Y,_),Z),
    Y > 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):-
    length(L,_),
    findall(M,member(M,L),Mx),
    length(Mx,N).

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

% added: 27 feb 2003.
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.
%
% product
% -----------------------------------------------------------  %
product([],1).
product([X|Members],Z):-
   product(Members,Z1),
  %number(X),
   Z is Z1 * X.



% output to file.
% -----------------------------------------------------------  %
tell_test(Goal):-
  open('tell.txt',write,S),
  tell('tell.txt'),
  Goal,
  current_stream('tell.txt',write,S),
  tell(user),wn(end),
  close(S).
%


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_goal(File,forall,G):-
   G0 = (nl,write(G),write('.')),
   G1 = forall(G,G0),
   tell_goal(File,G1).


% 実行時刻の取得
%--------------------------------

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.



/*********************************
%   Appendix. 
%     Brief Explanations for 
%     Belief Function
%     and Updating Rules
*********************************/



% Section 1. Basic Notions.
%-----------------------------
% (0) Assume a frame of discernment (i.e, a state space with some 
% epistemic interpretation of it).  Then beliefs of agent about uncertain 
% events are represented by belief function stated as follows.
% (1) A belief function is 0-1 normalized totally monotone capacity.
% (2) A belief function is also defined by  the bpa (mass) function. 
% Belief function of an event E summarizes its evidence in that
% it is the sum of all bpas of its focal set F included in E .
% (3) Conversely, a bpa function can be calculated via inversion 
% from belief functions.
% (4) A belief function has its conjugate, a possibility function.
% (5) Updating rule of belief function can be computed by using 
% updated bpa function which is the restriction of prior bpa to 
% the observed event.

% Section 2. Formula (Ref.1-3).
%-----------------------------------
/*
(1. n-monotone capacity)
  v(union(Ak))
   >=  sum(subset(J,[1,..,n]),
         (-1)^(|J|+1) * v(intersect(for(member(k,J)),Ak))
       ).

(2. definition via basic probability assignment; bpa)
  v(A) =  bel(A)
       =  sum(for(subset(B,A),B\=[]),mass(A)).

(3. Mobius inversion formula)
  mass(B) = sum(for(subset(A,B)), (-1)^|B-A| * v(A).
   where 
  mass([]) = 0,
  sum(mass(A))=1.

(4. plausibility function (possibility measure))
  v~(A) = v(All) - v(All - A) = 1 - v(-A).
        = sum(for(intersect(A,B)), mass(B)).

(5. belief updating)
 > Dempster's compostion rule:

  m_12(A) = sum(for(intersection(X,Y,A)), m_1(X) * m_2(Y)).

 > Dempster-Shafer conditioning:

  m_ds(A|B)
   = sum(for(intersection(E,B,A)),  m(E) / (1-Sum) ),
   where 
   Sum = sum(for(intersection(E,B,F),F\=[]), m(E)).

 >The objective part of the interpretation:

                 v(union(A,-B)) - v(-B)
  v_ds(A|B)  =  ------------------------
                 1 - v(-B)

               v~(B) - v~(intersect(-A,B))
           =  -----------------------------
               v~(B)

  > conditioning by Upper-Lower envelope of additive measures:

                 v(intersect(A,B))
  v_ul(A|B)  =  -----------------------------------------
                 v(intersect(A,B)) + v~(intersect(-A,B))

                 v~(intersect(A,B))
  v~ul(A|B)  =  -----------------------------------------
                 v~(intersect(A,B)) + v(difference(B,A))

*/





% end of the program.

return to front page.