loopybp.java

这是一个matlab的java实现。里面有许多内容。请大家慢慢捉摸。

/* Copyright (C) 2003 Univ. of Massachusetts Amherst, Computer Science Dept.
   This file is part of "MALLET" (MAchine Learning for LanguagE Toolkit).
   http://www.cs.umass.edu/~mccallum/mallet
   This software is provided under the terms of the Common Public License,
   version 1.0, as published by http://www.opensource.org.  For further
   information, see the file `LICENSE' included with this distribution. */

package edu.umass.cs.mallet.grmm;

import salvo.jesus.graph.Edge;
import java.util.Iterator;



/**
 * The loopy belief propagation algorithm for approximate inference in
 *  general graphical models.
 *
 * Created: Wed Nov  5 19:30:15 2003
 *
 * @author <a href="mailto:[e-mail]">Charles Sutton</a>
 * @version $Id: LoopyBP.java,v 1.1 2004/07/15 17:53:31 casutton Exp $
 */
public class LoopyBP extends BeliefPropagation {

	public static final int DEFAULT_MAX_ITER = 1000;

	private int iterUsed;
	private int maxIter;
	private double threshold = 0.0001;

	DiscretePotential[][] oldMessages;

	public int iterationsUsed () { return iterUsed; }

	// Note this does not have the sophisticated Terminator interface
	// that we've got in TRP.
	public LoopyBP() { this (DEFAULT_MAX_ITER); }
	public LoopyBP(int maxIter) { this.maxIter = maxIter; }
	
	private void initOldMessages (UndirectedModel mdl)
	{
		int n = mdl.getVerticesCount ();
		oldMessages = new DiscretePotential [n][n];
		for (Iterator it = mdl.getEdgeSet().iterator(); it.hasNext();) {
			Edge edge = (Edge) it.next();
			Variable v1 = (Variable) edge.getVertexA ();
			Variable v2 = (Variable) edge.getVertexB ();
			int i = mdl.getIndex (v1);
			int j = mdl.getIndex (v2);
			oldMessages [i][j] = new MultinomialPotential (mdl.get (i));
		}
	}

	private boolean hasConverged () 
	{
		for (int i = 0; i < oldMessages.length; i++) {
			for (int j = 0; j < oldMessages[i].length; j++) {
				DiscretePotential ptl1 = oldMessages [i][j];
				DiscretePotential ptl2  = messages [i][j];
				if (oldMessages [i][j] != null) {
					assert messages [i][j] != null 
						: "Message went from nonnull to null "+i+" --> "+j; 
					for (Iterator it = ptl1.assignmentIterator(); it.hasNext();) {
						Assignment assn = (Assignment) it.next();
						double val1 = ptl1.phi (assn);
						double val2 = ptl2.phi (assn);
						if (Math.abs (val1 - val2) > threshold) {
							return false;
						}
					}
				}
			}
		}
		return true;
	}

	public void computeMarginals (UndirectedModel mdl) 
	{
		initOldMessages (mdl);

		int iter;
		for (iter = 0; iter < maxIter; iter++) {
			logger.finer ("***LoopyBP iteration "+iter);
			super.initForGraph (mdl);
			propagate (mdl);
			if (hasConverged ()) break;
			oldMessages = messages;
		}
		iterUsed = iter;
		if (iter >= maxIter) {
			logger.info ("***Loopy BP quitting: not converged after "+maxIter+" iterations.");
		} else {
			iterUsed++;  // there's an off-by-one b/c of location of above break
			logger.info ("***LoopyBP converged: "+iterUsed+" iterations");
		}
	}
		
	private void propagate (UndirectedModel mdl)
	{
		// Send all messages. In this implementation, we send messages
		//  synchronously, so it doesn't matter what order we send
		//  messages in.
		for (Iterator it = mdl.getEdgeSet().iterator(); it.hasNext();) {
			Edge edge = (Edge) it.next();
			Variable v1 = (Variable) edge.getVertexA ();
			Variable v2 = (Variable) edge.getVertexB ();
			sendMessage (mdl, v1, v2);
			sendMessage (mdl, v2, v1);
		}
	}

	protected DiscretePotential msgProduct (DiscretePotential product, int idx, int excludeMsgFrom)
	{
		logger.finest ("Multipling messages to "+idx+" (except for "
									  + excludeMsgFrom + ")");
		if (product == null) {
			product = new MultinomialPotential (mdlCurrent.get (idx));
		}

		// Assumes that oldMessages [idx][idx] == null;
		for (int j = 0; j < oldMessages[idx].length; j++) {
			if ((oldMessages[idx][j] != null) && (j != excludeMsgFrom)) {
				logger.finest ("Multiplying in oldMessages["+idx+"]["+j+"] ="+
											 oldMessages[idx][j]);
				product.multiplyBy (oldMessages [idx][j]);
			}
		}
		return product;
	}

	public DiscretePotential lookupMarginal (Variable v1, Variable v2)
	{
		int idx1 = mdlCurrent.getIndex (v1);
		int idx2 = mdlCurrent.getIndex (v2);

		DiscretePotential product1 = msgProduct (null, idx1, idx2);
		DiscretePotential product2 = msgProduct (null, idx2, idx1);
		product1.multiplyBy (product2);

		DiscretePotential edgePtl = mdlCurrent.potentialOfEdge (v1, v2);
		product1.multiplyBy (edgePtl);

		DiscretePotential marg = product1.marginalize (new Variable[] { v1, v2 });
		marg.normalize ();

		return marg;
	}

	public DiscretePotential lookupMarginal (Clique c)
	{
		switch (c.size ()) {
			case 1: return lookupMarginal (c.get (0));
			case 2: return lookupMarginal (c.get (0), c.get (1));
			default: throw new IllegalArgumentException 
				("LoopyBP currently only supports node and edge cliques.");
		}
	}

	// xxx Assumes UndirectedModel
	public double lookupLogJoint (Assignment assn)
	{
		double accum = 0.0;

		// Compute using BP-factorization 
		// prod_s (p(x_s))^-(deg(s)-1) * ...
		for (Iterator it = mdlCurrent.getVerticesIterator(); it.hasNext();) {
			Variable var = (Variable) it.next();
			DiscretePotential ptl = lookupMarginal (var);
			int deg = mdlCurrent.getDegree(var);
			if (deg > 1)
				accum -= (deg - 1) * Math.log (ptl.phi (assn));
		}
		 
		// ... * prod_{st} p(x_s, x_t)
		for (Iterator it = mdlCurrent.getEdgeSet().iterator(); it.hasNext();) {
			Edge edge = (Edge) it.next();
			Variable v1 = (Variable) edge.getVertexA ();
			Variable v2 = (Variable) edge.getVertexB ();
			DiscretePotential p12 = lookupMarginal (v1, v2);
			DiscretePotential p1 = lookupMarginal (v1);
			DiscretePotential p2 = lookupMarginal (v2);
			accum += Math.log (p12.phi (assn));
		}
		
		return accum;
	}

}