home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / CLIPPER / MISC / AIE8908.ZIP / BAYES next >
Encoding:
Text File  |  1988-05-04  |  12.1 KB  |  470 lines

  1. program BayesNet(NodeFile,LinkFile,Output);
  2.  
  3.   This code is a basic implementation of Judea Pearl's belief
  4.   propagation algorithm for tree-structured Bayesian belief
  5.   networks.  The procedures and functions can be divided into
  6.   three basic groups:
  7.   
  8.   Math Support:
  9.      Normalize
  10.      MakeIdentityVector
  11.      TermProduct
  12.      TermQuotient
  13.      MatMult
  14.   Core:
  15.      ReviseBelief
  16.      UpdateNode
  17.      SubmitEvidence
  18.   General Support:
  19.      ReadString
  20.      FindNode
  21.      DumpNetwork
  22.         DumpNode
  23.      ReadNet
  24.         ReadNodes
  25.         ReadLinks
  26.      
  27.   The Core routines are described in the August AI Expert article.
  28.   The main program is set up to run the example from the May AI
  29.   Expert article.  It reads the net from two data files which are
  30.   described in ReadNodes and ReadLinks.  Be sure to figure out how
  31.   to RESET these files so that they get picked up correctly by those
  32.   procedures.   
  33.  
  34. const
  35.   MaxString      = 15;
  36.   MaxValues      =  5;
  37.         
  38. ype
  39.   StringRange = 1..MaxString;
  40.   ValueRange  = 1..MaxValues;
  41.   StringType  = packed array[StringRange] of char;
  42.   NetVector   = record
  43.                    Data:  array [ValueRange] of real;
  44.                    NVals: ValueRange
  45.                 end;
  46.   CPType      = record
  47.                    Data:        array[ValueRange,ValueRange] of real;
  48.                    NRows,NCols: ValueRange
  49.                 end;
  50.   NetNodePtr  = ^NetNode;
  51.   NetNode     = record
  52.                    Name:                       StringType;
  53.                    NumValues:                  ValueRange;
  54.                    Values:                     array[ValueRange] of 
  55. tringType;
  56.                    Belief,Pi,IncomingPi,
  57.                    ExternalLambda,
  58.                    Lambda,OutgoingLambda:      NetVector;
  59.                    Parent,NextNode,
  60.                    NextSibling,FirstChild:     NetNodePtr;
  61.                    CPMatrix,TransCPMatrix:     CPType
  62.                 end;
  63.                        
  64. ar      NodeFile,LinkFile:            Text;
  65.         NetRoot,NodeList:             NetNodePtr;
  66.         EvidenceVector:               NetVector;
  67. { ******************** Math Support ******************** }
  68.        
  69. rocedure Normalize(var Vector: NetVector);
  70.  Scales incoming Vector so that it sums to unity }
  71. ar i:   ValueRange;
  72.    Sum: real;
  73. begin
  74. um := 0;
  75. ith Vector do
  76.   begin
  77.   for i := 1 to NVals do
  78.      Sum := Sum + Data[i];
  79.   for i := 1 to NVals do
  80.      Data[i] := Data[i] / Sum
  81.   end
  82. nd;
  83.         
  84. rocedure MakeIdentityVector(var Vector: NetVector;Length: ValueRange);
  85.  Makes incoming Vector into an identity vector of specified length}
  86. ar i: ValueRange;
  87. begin
  88. ith Vector do
  89.   begin   
  90.   NVals := Length;
  91.   for i := 1 to Length do
  92.      Data[i] := 1.0   
  93.   end   
  94. nd;
  95. procedure TermProduct(var V1,V2,Result: NetVector);
  96.  Returns term product of V1 and V2 in Result }                      
  97. ar i:  ValueRange;
  98. begin
  99. f v1.NVals <> v2.Nvals then
  100.   writeln('*** Dimension error in TermProduct ***');
  101. ith Result do
  102.   begin
  103.   Nvals := V1.Nvals;
  104.   for i := 1 to NVals do
  105.      Data[i] := V1.Data[i] * V2.Data[i]
  106.   end
  107. end;
  108. procedure TermQuotient(var V1,V2,Result: NetVector);
  109.  Returns term quotient of V1 and V2 in Result }                               
  110.             
  111. ar i:  ValueRange;
  112. begin
  113. f v1.NVals <> v2.Nvals then
  114.   writeln('*** Dimension error in TermQuotient ***');
  115. ith Result do
  116.   begin
  117.   Nvals := V1.Nvals;
  118.   for i := 1 to NVals do
  119.      Data[i] := V1.Data[i] / V2.Data[i]
  120.   end
  121. nd;
  122. procedure MatMult(var InMat:  CPType;var InVec:  NetVector;var OutVec: 
  123. etVector);
  124.  Simplified matrix multiplication matrix routine.  Multiplies InMat * InVec
  125.  to produce OutVec.  Interprets InVec to be a NVals X 1 matrix. }
  126. ar Row,Col: ValueRange;
  127. begin
  128. f InMat.NCols <> InVec.NVals then
  129.   writeln('*** Dimension error in MatMult ***');
  130. ith InMat do
  131.   begin
  132.   OutVec.NVals := NRows;
  133.   for Row := 1 to NRows do
  134.      begin
  135.      OutVec.Data[Row] := 0.0;
  136.      for Col := 1 to NCols do
  137.         OutVec.Data[Row] := OutVec.Data[Row] + Data[Row,Col] * InVec.Data[Col]
  138.       end
  139.   end
  140. nd;
  141. { ******************** Core ******************** }
  142. procedure ReviseBelief(Node: NetNodePtr);
  143. ar Child:  NetNodePtr;
  144. egin
  145. ith Node^ do
  146.   begin
  147.   { Part (a) of Figure 4 }
  148.   if Parent <> nil then
  149.      MatMult(TransCPMatrix,IncomingPi,Pi);
  150.   { Part (b) of Figure 4 }
  151.   Lambda := ExternalLambda;
  152.   Child := FirstChild;
  153.   while Child <> nil do
  154.      begin
  155.      TermProduct(Child^.OutgoingLambda,Lambda,Lambda);
  156.      Child := Child^.NextSibling
  157.      end;
  158.   { Shaded part of Figure 4 }
  159.   TermProduct(Lambda,Pi,Belief);   
  160.   Normalize(Belief)
  161.   end
  162. nd;
  163. procedure UpdateNode(Node,Sender: NetNodePtr);
  164. ar Child:  NetNodePtr;
  165. egin
  166. ith Node^ do
  167.   begin
  168.   ReviseBelief(Node);
  169.   { Update OutgoingLambda & send update message to parent
  170.     (part (c) of Figure 4) }
  171.   if (Parent <> Sender) and (Parent <> nil) then
  172.      begin
  173.      MatMult(CPMatrix,Lambda,OutgoingLambda);
  174.      UpdateNode(Parent,Node)
  175.      end;
  176.   { Update IncomingPi and send update message to children
  177.     (part (d) of Figure 4) }
  178.   Child := FirstChild;
  179.   while Child <> nil do
  180.      begin
  181.      if Child <> Sender then
  182.         begin
  183.         TermQuotient(Belief,Child^.OutgoingLambda,Child^.IncomingPi);
  184.         UpdateNode(Child,Node)
  185.         end;
  186.      Child := Child^.NextSibling
  187.      end
  188.   end
  189. nd;
  190. procedure SubmitEvidence(Node: NetNodePtr;var Evidence: NetVector);
  191. ar i: ValueRange;
  192. egin
  193. ith node^ do
  194.   begin
  195.   writeln('Submitting evidence to ',Node^.Name,', evidence is:');
  196.   for i := 1 to Evidence.NVals do
  197.      writeln('[',Values[i],'] = ',Evidence.Data[i]);
  198.   TermProduct(Evidence,ExternalLambda,ExternalLambda);
  199.   UpdateNode(Node,nil)
  200.   end
  201. nd;
  202. { ******************** General Support ******************** }
  203. function ReadString(var InFile: Text;var InString: StringType): boolean;
  204.  Reads InFile, returning next string in InString.  Returns FALSE upon
  205.  encountering end of file, otherwise returns TRUE. }
  206. ar i,j:  StringRange;
  207. begin
  208. f eof(InFile) then
  209.   ReadString := false
  210. lse
  211.   begin
  212.   i := 1;
  213.   while not eoln(InFile) do
  214.      begin
  215.      read(InFile,InString[i]);
  216.      i := i + 1
  217.      end;
  218.   readln(InFile);
  219.   for j := i to MaxString do
  220.      InString[j] := ' ';
  221.   ReadString := true   
  222.   end;   
  223. nd;
  224.         
  225. unction FindNode(NodeName: StringType):NetNodePtr;
  226.  Searches network for node having specified NodeName. }                  
  227. ar CurrentNode:   NetNodePtr;
  228. begin
  229. urrentNode := NodeList;
  230. hile (CurrentNode^.Name <> NodeName) and (CurrentNode <> nil) do
  231.   CurrentNode := CurrentNode^.NextNode;
  232. f CurrentNode = nil then
  233.   begin
  234.   writeln('*** Error in FindNode -- cannot find ',NodeName);
  235.   FindNode := nil
  236.   end
  237. lse
  238.   FindNode := CurrentNode
  239. nd;
  240.        
  241. rocedure DumpNetwork(Node: NetNodePtr);
  242.  Recursive procedure to dump network, given pointer to root }
  243. procedure DumpNode(Node: NetNodePtr);
  244.  Simple procedure to dump a single node }
  245. onst Stars = '*************************************************';
  246. var CurrentValue,NumRows,NumCols,Row,Col:  ValueRange;
  247. begin
  248. riteln(Stars);
  249. ith Node^ do
  250.   begin
  251.   writeln('Dumping ',Name);
  252.   for CurrentValue := 1 to NumValues do
  253.      writeln('Pi[',Values[CurrentValue],'] = ',Pi.Data[CurrentValue]);
  254.   for CurrentValue := 1 to NumValues do
  255.      writeln('Lambda[',Values[CurrentValue],'] = ',Lambda.Data[CurrentValue]);
  256.   for CurrentValue := 1 to NumValues do
  257.      writeln('Belief[',Values[CurrentValue],'] = ',Belief.Data[CurrentValue]);
  258.   if Parent <> nil then
  259.      begin
  260.      writeln;
  261.      writeln('CP Matrix:');
  262.      for Row := 1 to CPMatrix.NRows do
  263.         begin
  264.         for Col := 1 to CPMatrix.NCols do
  265.            write(CPMatrix.Data[Row,Col]);
  266.         writeln
  267.         end
  268.      end
  269.   end;
  270. riteln(Stars);
  271. riteln('Type <cr> to continue...');
  272. eadln
  273. nd;   { of DumpNode }
  274. var CurrentNode: NetNodePtr;
  275. begin
  276. f Node <> nil then
  277.   begin
  278.   DumpNode(Node);
  279.   CurrentNode := Node^.FirstChild;
  280.   while CurrentNode <> nil do
  281.      begin
  282.      DumpNetwork(CurrentNode);
  283.      CurrentNode := CurrentNode^.NextSibling
  284.      end 
  285.   end  
  286. nd;
  287. procedure ReadNet(var NodeFile,LinkFile: Text);
  288.       
  289. rocedure ReadNodes(Var NodeFile: Text);
  290.  This procedure reads the NodeFile.  Format of file is as follows:
  291.    Node 1 name
  292.   Node 1 number of values
  293.   Node 1 value 1 name
  294.   Node 1 value 1 prior probability (ignored except for root node)
  295.   Node 1 value 2 name
  296.   Node 1 value 2 prior probability (ignored except for root node)
  297.           .....
  298.   Node 1 value n name
  299.   Node 1 value n prior probability (ignored except for root node)
  300.   Node 2 name
  301.           .....
  302.           etc.
  303.    
  304. ar NodeName:      StringType;
  305.    CurrentValue:  ValueRange;
  306.    eofStatus:     boolean;
  307.    CurrentNode:   NetNodePtr;
  308.                    
  309. egin
  310. eset(NodeFile);
  311. odeList := nil;
  312. hile ReadString(NodeFile,NodeName) do
  313.   begin
  314.   new(CurrentNode);
  315.   with CurrentNode^ do
  316.      begin
  317.      Name := NodeName;
  318.      readln(NodeFile,NumValues);
  319.      for CurrentValue := 1 to NumValues do
  320.         begin
  321.         eofStatus := ReadString(NodeFile,Values[CurrentValue]);
  322.         readln(NodeFile,Pi.Data[CurrentValue])
  323.         end;
  324.      Pi.NVals := NumValues;
  325.      Parent := nil;
  326.      NextSibling := nil;
  327.      FirstChild := nil;
  328.      NextNode := NodeList;
  329.      NodeList := CurrentNode;
  330.      MakeIdentityVector(ExternalLambda,NumValues);
  331.      MakeIdentityVector(Lambda,NumValues)
  332.      end
  333.   end;
  334. lose(NodeFile)
  335. nd;   { or ReadNodes }
  336. procedure ReadLinks(var LinkFile: Text);
  337.  This procedure reads the NodeFile.  Be careful here, upper/lower case
  338.  must match identically the node names in NodeFile.  Format of file is
  339.  as follows:
  340.    Top Node name for first link
  341.   Bottom Node name for first link
  342.   1st row of CP matrix
  343.   2nd row of CP matrix
  344.          ....
  345.   nth row of CP matrix
  346.   Top Node name for second link
  347.   Bottom Node name for second link
  348.   1st row of CP matrix
  349.   2nd row of CP matrix
  350.          ....
  351.   nth row of CP matrix
  352.       etc.
  353.  
  354. ar TopNodeName,BottomNodeName:      StringType;
  355.    TopNode,BottomNode:              NetNodePtr;
  356.    Row,Col:                         ValueRange;
  357.    eofStatus:                       boolean;
  358.                    
  359. egin
  360. eset(LinkFile);
  361. hile ReadString(LinkFile,TopNodeName) do
  362.   begin
  363.   TopNode := FindNode(TopNodeName);
  364.   eofStatus := ReadString(LinkFile,BottomNodeName);
  365.   BottomNode := FindNode(BottomNodeName);
  366.   with BottomNode^ do
  367.      begin
  368.      CPMatrix.NRows := TopNode^.NumValues;
  369.      CPMatrix.NCols := NumValues;
  370.      TransCPMatrix.NRows := CPMatrix.Ncols;
  371.      TransCPMatrix.NCols := CPMatrix.NRows;
  372.      for Row := 1 to CPMatrix.NRows do
  373.         begin
  374.         for Col := 1 to CPMatrix.Ncols do
  375.            begin
  376.            read(LinkFile,CPMatrix.Data[Row,Col]);
  377.            TransCPMatrix.Data[Col,Row] := CPMatrix.Data[Row,Col]
  378.            end;
  379.         readln(LinkFile)
  380.         end;
  381.      NextSibling := TopNode^.FirstChild;
  382.      Parent := TopNode;
  383.      MakeIdentityVector(OutgoingLambda,TopNode^.NumValues)
  384.      end;
  385.   TopNode^.FirstChild := BottomNode
  386.   end
  387. nd;   { of ReadLinks }
  388. begin
  389. eadNodes(NodeFile);
  390. eadLinks(LinkFile);
  391.  Find root of network. }
  392. etRoot := NodeList;
  393. hile NetRoot^.Parent <> nil do
  394.   NetRoot := NetRoot^.NextNode;
  395.  Initialize network }   
  396. pdateNode(NetRoot,nil)
  397. nd;
  398.     
  399. egin
  400. { Read network in }
  401. eadNet(NodeFile,LinkFile);
  402. { Take a look }
  403. umpNetwork(NetRoot);
  404. { Store evidence from rain alarm in EvidenceVector }
  405. ith EvidenceVector do
  406.   begin
  407.   Data[1] := 0.8;
  408.   Data[2] := 0.04;
  409.   NVals := 2
  410.   end;
  411. { Submit EvidenceVector to Rain node }   
  412. ubmitEvidence(FindNode('Rain           '),EvidenceVector);
  413. { Take a look }
  414. umpNetwork(NetRoot);
  415. { Store evidence from telephone call in EvidenceVector }
  416. ith EvidenceVector do
  417.   begin
  418.   Data[1] := 1.0;
  419.   Data[2] := 0.02;
  420.   NVals := 2
  421.   end;
  422. { Submit EvidenceVector to Sunburn node }   
  423. ubmitEvidence(FindNode('Sunburn        '),EvidenceVector);
  424. { Take a look }
  425. umpNetwork(NetRoot)
  426. end.
  427. Clouds
  428. ain
  429. .6 0.4
  430. .0 1.0
  431. ain
  432. lay Game
  433. .05 0.95
  434. .00 0.00
  435. louds
  436. unburn
  437. .1 0.9
  438. .7 0.3
  439. louds
  440.  
  441. resent
  442. 1
  443. bsent
  444. 9
  445. ain
  446.  
  447. resent
  448.  
  449. bsent
  450.  
  451. lay Game
  452.  
  453. es
  454.  
  455. o
  456.  
  457. unburn
  458.  
  459. es
  460.  
  461. o
  462.  
  463. s
  464. unburn
  465. .1 0.9
  466. .7 0.3
  467. louds
  468.  
  469. resent
  470.