Merged SVN head ...
authorJoris Mooij <jorism@marvin.jorismooij.nl>
Mon, 29 Sep 2008 09:50:40 +0000 (11:50 +0200)
committerJoris Mooij <jorism@marvin.jorismooij.nl>
Mon, 29 Sep 2008 09:50:40 +0000 (11:50 +0200)
- Everything from SVN head is now merged (apart from files which
  are only in SVN head, and some stuff described in STATUS)
- Added old history to ChangeLog
- Added Iterations() to each method and to output of tests/test
- Added damping to various algorithms
- Cleaned up headers of algorithms

24 files changed:
ChangeLog
README
STATUS
TODO
include/dai/bp.h
include/dai/daialg.h
include/dai/exactinf.h
include/dai/factorgraph.h
include/dai/hak.h
include/dai/jtree.h
include/dai/lc.h
include/dai/mf.h
include/dai/mr.h
include/dai/regiongraph.h
include/dai/treeep.h
src/bp.cpp
src/hak.cpp
src/jtree.cpp
src/lc.cpp
src/mf.cpp
src/mr.cpp
src/treeep.cpp
tests/test.cpp
tests/testfast.out

index 3cb88ac..d5c9cbc 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -125,3 +125,310 @@ libDAI-0.2.0 (2006-11-30)
 -------------------------
 
 First public release.
+
+
+0.1.5 (2006-11-30)
+------------------
+
+Regressions
+
+- tests/testlcbp and tests/testlcbp are broken.
+- EXACT method does not work anymore.
+- The Properties framework gives a speed penalty because of the lookup
+  costs involved; inner loops should be optimized.
+
+General framework
+
+- DAIAlg is now a template class; typedefs for DAIAlg<FactorGraph> and for 
+  DAIAlg<RegionGraph> are provided. In this way, we do not have to write "wrapper"
+  functions to forward functionality from either FactorGraph or RegionGraph
+  to DAIAlg. Functionality like clamping can be implemented in FactorGraph
+  and in RegionGraph and no explicit interface is needed in descendants.
+- New abstract base class InfAlg added, representing an inference algorithm,
+  from which DAIAlg<T> inherits. This solves the incompatibility problems of
+  DAIAlg<T> for different T (e.g. DAIAlg<FactorGraph> was incompatible with
+  DAIAlg<RegionGraph>). More work is required to reduce code duplication
+  (make FactorGraph part of InfAlg).
+- Added generic interface (nrVars(), Vars(), nrFactors(), factor(size_t),
+  beliefs(), belief(VarSet &), ...) to InfAlg and descendants.
+- Added a saveProbs/undoProbs interface to InfAlg and descendants that enables
+  one to save a few factors, modify them (e.g. clamp them), and then restore them
+  to their old values. Undo should also init the corresponding messages / beliefs.
+  This can be used if a given factor graph repeatedly needs to be clamped in
+  different ways and an approximation method is run for each clamping; using the
+  saveProbs/undoProbs can give a significant speed increase.
+- Switched to a general Properties framework that handles the parameters of
+  all inference methods in a uniform manner. The Properties class is a map of 
+  several properties in boost::any objects, indexed by their names (strings). 
+  It can read from a stream and write to a stream. It is recursive, in the sense
+  that a Properties object can hold a variable of type Properties as well.
+- Added a generic way of constructing inference algorithms from a factor graph,
+  name and properties object. Added the newInfAlg function which constructs 
+  the requested object. This is used by LCBP, the Matlab interface and the
+  command line (test) interface.
+- Added a generic enum framework for enum parameters. Although implemented as a
+  hack, it is quite useful in that it drastically simplifies and reduces the
+  amount of code for handling enum parameters.
+- Provided generic functions for calculating marginals in different ways that
+  work for all approximate inference methods.
+
+Bugfixes
+
+- Fixed GBP free energy.
+- Fixed bug in junctiontree (it didn't init the _vars variable).
+- Corrected two bugs in operator&& and operator|| in VarSet (they returned
+  the logical NOT of what they should return).
+- Fixed bug in RegionGraph::RecomputeOR(s).
+- Fixed bug in utils/create_dreg_fg: 
+  graph structure was not random for given parameters (forgot to call srand()).
+- TreeEP bug workaround: use the complete junction tree instead of a subtree.
+- Fixed bug in JTree::HUGIN() and JTree:ShaferShenoy() in case of junction tree 
+  that consists of one outer region only.
+- Fixed INIT bug in LCBP2::UpdatePancake().
+- Fixed MaxDiffs flow (except for MR).
+
+New functionality
+
+- HAK supports several default cluster choices:
+  minimal (only factors)
+  delta   (Markov blankets)
+  loop    (all loops consisting of loops consisting of <loopdepth> or less variables)
+  Only the maximal clusters are used as outer clusters.
+- Implemented TreeEP. It generalizes the heuristic method described in the
+  Minka & Qi paper for obtaining a tree with the most relevant interactions to
+  higher order interactions. Almost all optimizations described in the Minka & Qi
+  paper are used, except that evidence is passed over the whole tree instead of
+  relevant subsets (the latter is almost implemented but buggy). Also added
+  alternative (worst-case) algorithm that uses a maximum spanning tree on the
+  weighted graph where the weight between neighbours i and j is given by
+  N(psi,i,j), where psi is the product of all the factors involving both i and j
+  (which is an upper bound on the effective interaction between i and j).
+- Implemented MR (MontanariRizzo) based on Bastian's code, but extended it 
+  to be able to handle connectivities larger than 3 (in principle, up to 32).
+  It supports different initialization methods (the original RESPPROP,
+  the CLAMPING method and EXACT which uses JTree) and different update methods
+  (FULL and LINEAR).
+- Implemented LCBP2, an analogon of LCBP which represents pancakes as little 
+  networks and uses some approximate inference method on them for calculating 
+  marginals.
+- Now there are several LCBP variants (LCBP, LCBPI, LCBPJ, LCBPK, LCBPL);
+  LCBPJ works only for pairwise, LCBPK is LCBP improved for higher order
+  interactions and LCBPL is LCBPI improved for higher-order interactions.
+- Wrote one single program utils/createfg for creating various types of
+  random factor graphs.
+- Wrote utility to visualize factor graphs using graphviz.
+  (it uses the BOOST Program Options library)
+- Added fginfo utility that displays some info about a .fg file.
+- Implemented Factor::strength function that calculates the potential strength 
+  N(psi,i,j) between variables i and j as described in cs.IT:0504030
+- Wrote a general MatLab interface matlab/dai (similar to tests/test);
+  this unified the matlab functions dai, dai_bp, dai_mf, dai_jt, dai_tep, dai_cvm.
+- Added MATLAB routine that returns contraction matrix A for BP convergence analysis.
+- Implemented a MATLAB interface ai_potstrength for Factor::strength
+- Added Martijn's x2x
+
+Improvements of existing code
+
+- Reimplemented RegionGraph and descendants: a RegionGraph ISA FactorGraph 
+  and also a BipartiteGraph<FRegion,Region>. It now also keeps a map that 
+  associates outer region indices to factor indices (no powers yet, this 
+  is deemed superfluous) and provides functions to recompute (some of) the 
+  outer regions from the factors.
+- InfAlg descendants run() methods now stop immediately and return NAN in case
+  they detect NANs. Only BP does not do NAN checking for performance reasons.
+- LCBP now works with factors containing zeroes (by defining x/0 = 0).
+- HAK, GBP and DoubleLoop now conform to the standards for verbose reporting,
+  timing and convergence criteria.
+- Implemented logZ() for JTree. It does the calculation during message-passing.
+- Marginal2ndO now optionally divides by the single node beliefs (to the power n-2); 
+  hopefully this will give more accurate approximations.
+- Marginal and Marginal2ndO (optionally) use the new saveProbs/undoProbs functionality 
+  for a faster way of calculating marginals, which does not require a call to init() 
+  nor cloning the whole object for each clamping. This leads to a significant speedup.
+- LCBP (and LCBP2) now have complete flexibility in the specification of the
+  inner method, i.e. the method used to generate the initial cavity approximations. 
+  One can pass two strings, a name and a properties string, and LCBP simply invokes 
+  newInfAlg to construct the corresponding inference algorithm and uses the generic
+  marginal functions to approximate cavity marginals.
+- Replaced the global "method" variable by local properties and removed ai.h
+- Added some methods to Factor (operators *, *=, /, /= with doubles as
+  second argument, operators -, +=, -= with other Factors as second
+  arguments, randomize(), RemoveFirstOrderInteractions) and similar
+  operations to Prob
+- Moving towards boost::program_options for handling command line arguments
+  (tests/test is done).
+- Renamed some FactorGraph methods:
+  nr_vars       -> nrVars
+  nr_factors    -> nrFactors
+  varind        -> findVar
+  factorind     -> findFactor
+  makeFacCavity -> makeFactorCavity
+- LCBP_SEQMAXRES has been removed because it did strange things.
+- Implemented RandomDRegularGraph
+- Implemented JTree::calcMarginal for marginals that are not confined 
+  within one cluster (using cut-set conditioning).
+- Added isConnected() method to FactorGraph (some methods do not work with
+  disconnected factor graphs).
+- Pair beliefs are now calculated in a symmetrical way by calcPairBeliefs
+- Removed single node interaction "correction" code from clamping methods
+- Removed calcCavityDist and calcCavityDist2ndO
+- No longer depends on GSL.
+- Increased portability by combining platform dependant utility functions
+  in util.{h,cpp}.
+- Wrote *.m files providing help
+
+Testing framework
+
+- Made a new and significantly improved testing framework that provides most
+  functionality from the command line.
+- The basis is provided by tests/test, which uses the newInfAlg functionality
+  and enables the user to easily compare from the command line different 
+  inference methods on a given factorgraph. All parameters can be specified. 
+  Output consists of CPU time, average and maximum single variable marginal 
+  errors, relative logZ error and MaxDiff(). 
+- tests/aliases.conf contains aliases for standard combinations of methods
+  and their options (which can be used in tests/test).
+- tests/large contains several bash/python scripts that create random factor 
+  graphs, compare several approximate inference algorithms (using tests/test) and
+  allow for easy visualization of the results using PyX.
+- Added several .fg files for test purposes to /tests (e.g. two ALARM versions
+  alarm.fg and alarm_bnt.fg; testfast.fg, a 4x4 periodic Ising grid for 
+  regression testing).
+- Added a regression test to the Makefile which is included in the standard
+  target.  It compares all inference methods on tests/testfast.fg with the
+  results stored in tests/testfast.out
+
+Miscellaneous
+
+- Expanded all tabs to spaces (":set tabstop 4\n:set expandtab\n:retab" in vim)
+- Experimental MATLAB code added for StarEP approximation on cavity
+- Renamed project to libDAI and changed directory name accordingly.
+- Renamed JunctionTree to JTree.
+- Fixed licensing (now it's officially GPL).
+- Improved README
+
+
+revision 252
+------------
+
+Functionality
+- Added RegionGraph, GBP, CVM and HAK (double-loop).
+- Added JunctionTree (with two update algorithms, HUGIN and Shafer-Shenoy), which is a 
+  RegionGraph.
+- NormType is now chosen automatically (in case of negative factors, Prob::NORMLINF is used,
+  otherwise the default Prob::NORMPROB is used). Also, in case of negative factors, the 
+  RegionGraph constructors assign each Factor to a unique outer region instead of dividing 
+  it over all subsuming outer regions. See README for when negative factors are known to work
+  and when not.
+- FactorGraph::FactorGraph(const vector<Factor>) only gives a warning in case of short loops,
+  it does not automatically merge factors anymore.
+- Removed BP_SEQMAXRESNOCLEAR (all cavity initialization methods now are implicitly NOCLEAR)
+- Added MATLAB interface functions ai_readfg, ai_removeshortloops and ai_bp
+- Added LCBP-III type that should be equivalent to LCBP-II, but can handle zeroes
+  in potentials. Note that it is significantly slower than LCBP-II (and has to be reimplemented
+  such that it does not store the complete pancakes, but represents them as little factor graphs).
+
+Implementation / code
+- Made a seperate type WeightedGraph, which until now only implements Prim's
+  maximal spanning tree algorithm and is only used by the junction tree code. It might
+  be good to make it a class somewhere in the future and extend it's interface.
+- Made a seperate class ClusterGraph, which is only used by the junction tree
+  code. It's main purpose is a graph-theoretical variable elimination algorithm.
+- Implemented the heuristic "minimum-new-edges-in-clique-graph" for variable elimination.
+- Massive code cleanup, moving towards "generic" programming style, using 
+  multiple inheritance and polymorphism.
+  o BP, LCBP, MF, HAK and JunctionTree now inherit from a common DAIAlg class
+  o Made generic functions Marginal, Marginal2ndO, calcCavityDist, calcCavityDist2ndO, clamp
+    that can be used by FactorGraph-based DAIAlgs.
+  o Created TProb<T> class, which stores a probability vector (without the accompanying indexing
+    and VarSet) and provides functionality for it (which is used by TFactor<T>).
+  o Rewrote the VarSet class. It now caches its statespace(). It now privately inherits from set<Var>.
+    I had to overload the insert methods of set<Var> so that they calculate the new statespace.
+  o Rewrote the TFactor class. The TFactor class now HAS a TProb and HAS a VarSet.
+- Rewrote BP to use the new TProb<T> interface. Performance of BP improved up to a factor 6 by:
+  o using Prob's instead of Factor's;
+  o splitting the multiplication of the messages into two parts (thanks to Vicenc!);
+  o optimizing the calculation of the beliefs (only the message calculations were optimized till now).
+  o replacing FactorGraph::_nb1 and _nb2 (which were set<size_t>) by vector<size_t>
+- LCBP now seperately stores cavitydists and pancakes. Added InitPancakes() method
+  that takes the cavitydists and multiplies them with the relevant factors. This
+  resulted in an API change in AI which now accepts and returns initial cavitydists
+  instead of initial pancakes.
+
+Minor changes
+- Started writing DoxyGen documentation
+- Renamed lcptab2fg.m matlab/ai_writefg.m
+- Moved all matlab stuff to matlab/
+- More detailed reporting (also reports clocks used).
+- Marginal and Marginal2ndO now use *differences* in logZ to avoid NaNs.
+- Improved testing suite.
+- Removed logreal support.
+- FactorGraph now also supports input streams and ignores comment lines in .fg files.
+- Added tests/create_full_fg.cpp and tests/create_ising_fg.cpp which create
+  full and periodic 2D Ising networks according to some command line parameters.
+- Now logZ really returns logZ instead of -logZ.
+- Added FactorGraph::WriteToDotFile
+
+
+0.1.4 (2006-04-13)
+------------------
+- Added file IO routines to read and write factorgraphs.
+- Added L-infinity normalization for use with (partially) negative factors.
+- Renamed BetheF, MFF to logZ, which now use complex numbers to be able to
+handle (partially) negative factors.
+- Added test suite.
+- All probabilities are now represented using double instead of LogReal<double>.
+- Combined Alg and Alg3 into LCBP. Several update schemes possible.
+- Combined several variants of BP into doBP. Several update schemes possible.
+Now uses precalculated indices for optimization.
+- Renamed Node -> Var and Potential -> Factor.
+- Extensive code cleanup. More use of OO techniques. Changed API.
+- MaxIter now means maximum number of passes (corresponding to number of
+_parallel_ updates).
+- MaxDiff now means the maximum L-infinity distance between the updated and
+original single variable beliefs, for all AI methods.
+- Implemented RBP which converges faster than BP for difficult problems.
+- Now uses method parameter which is a bitmask combining outmethod and inmethod
+(see ai.h).
+
+
+0.1.3   (2006-03-22)
+--------------------
+- All AI methods now return maxdiff
+- ai.cpp:
+    o Now returns maxdiffinner and maxdiffouter
+    o New BP2ndO innermethod (estimate only 2nd order cavity interactions)
+    o New InitCav outermethod (only create initial cavity distributions)
+- bp.cpp:
+    o New CavityDist2ndO which estimates 2nd order cavity interactions
+- Makefile:
+    o Bugfix: removed dependencies on algwim.*
+
+
+0.1.2   (2006-02-28)
+--------------------
+- Cleaned up alg.cpp (removed Alg2 and its corresponding data structures).
+- Added the possibility to provide initial cavity distributions as an input
+argument to ai (not much error checking is done, so be careful).
+- Potentials2mx now correctly sets the dimensions of the P field (i.e. for
+the output arguments Q, Q0 of ai).
+- Removed algwim.* since it does not work.
+
+
+0.1.1   (2006-02-28)
+--------------------
+- The constructors of (Log)FactorGraph and LogFactorGraph from a
+vector<(Log)Potential> now merge potentials to prevent short loops (of length 
+4) in the factor graph. These are used in ai to construct the factor graphs 
+from the psi argument. If compiled with DEBUG defined, the method calc_nb()
+of BipGraph checks for the existence of short loops.
+- Changed calling syntax of ai (now the actual syntax *does* correspond to its
+description in the help).
+- ai does not hook cout anymore (which caused weird segfaults).
+- Fixed a bug in an assert statement in the matlab interface code in ai.cpp.
+- Removed network.* since it is not useful.
+
+
+0.1.0   (2006-02-28)
+--------------------
+First version worthy a version number.
diff --git a/README b/README
index 1fe39c9..f8e36a7 100644 (file)
--- a/README
+++ b/README
@@ -27,6 +27,18 @@ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 ----------------------------------------------------------------------------------
 
 
+SCIENTISTS: please be aware that the fact that this program is released as Free
+Software does not excuse you from scientific propriety, which obligates you to
+give appropriate credit! If you write a scientific paper describing research
+that made substantive use of this program, it is your moral obligation as a
+scientist to (a) mention the fashion in which this software was used, including
+the version number, with a citation to the literature, to allow replication;
+(b) mention this software in the Acknowledgements section.  The appropriate
+citation is: J. M. Mooij (2008) libDAI: A free/open source C++ library for
+Discrete Approximate Inference methods, http://mloss.org/software/view/77/.
+Moreover, as a personal note, I would appreciate it if you would email me with
+citations of papers referencing this work so I can mention them to my funding
+agent and tenure committee.
 
 What is libDAI?
 ---------------
diff --git a/STATUS b/STATUS
index 66d4ced..b63f5cf 100644 (file)
--- a/STATUS
+++ b/STATUS
@@ -1,34 +1,20 @@
-- Idea: a FactorGraph and a RegionGraph are often equipped with
-extra properties for nodes and edges. The code to initialize those
-is often quite similar; maybe this can be abstracted to classes
-like ExtFactorGraph and ExtRegionGraph (Ext stands for Extended), e.g.
-template <typename NodeProperties, typename EdgeProperties>
-class ExtFactorGraph : public FactorGraph {
-       public:
-               std::vector<NodeProperties>               nodeProps;
-               std::vector<std::vector<EdgeProperties> > edgeProps;
-       // blabla
-}
-A disadvantage of this approach may be worse cachability.
-A problem is if there are nog properties for nodes (type 1), nodes (type 2)
-or edges. Maybe this can be solved using specializations, or using variadac
-template arguments or something similar?
-- BipartiteGraph::isConnected should be optimized.
+- BipartiteGraph::isConnected should be optimized using boost::graph
 - http://www.boost.org/development/requirements.html#Design_and_Programming
 - Would it be a good idea to cache second-order neighborhoods (delta's) in BipGraph?
 - Would it be a good idea to add the variable label -> index hashmap in FactorGraph,
   to replace the linear searches that are performed every time for findVar()?
   No, a better idea is to avoid calls to findVar() as much as possible.
+- Would it be a good idea to remove the states() caching from VarSet? 
+  Then, we could turn a VarSet into an IndexSet (set of integers).
 - Can the FactorGraph constructors be optimized further?
+- Idea: use a PropertySet as output of a DAIAlg, instead of functions like maxDiff and Iterations().
+- A DAIAlg<T> should not inherit from a FactorGraph/RegionGraph, but should store a reference to it
+- Add max-product to BP
 - Remove x2x?
-- Add iterations (like maxdiff, but counts iterations) to some algorithms
-  Idea: use a PropertySet as output of a DAIAlg
-- A DAIAlg<T> should not inherit from a FactorGraph/RegionGraph, but should store a
-reference to it
-- limit???
+- Remove MR?
 
 TODO FOR RELEASE:
-- Test Visual C++ compatibility
+- Test {Visual C++, cygwin, gcc various version} compatibility; state tested compilers/build environments in README
 - Figure out which libraries are required and document in README
   boost headers, boost::program_options library, boost::graph library,
   boost::math library (under Windows)
@@ -101,21 +87,17 @@ bp.h
 bp.cpp
 jtree.h
 jtree.cpp
-
-FILES IN SVN HEAD THAT ARE STILL RELEVANT:
-ChangeLog
-README
-TODO
-
 hak.h
 hak.cpp
 lc.h
 lc.cpp
-mr.h
-mr.cpp
 treeep.h
 treeep.cpp
-
+mr.h
+mr.cpp
+TODO
+ChangeLog
+README
 
 FILES IN SVN HEAD RELEVANT FOR A LATER RELEASE:
        matlab/dai_potstrength.*
@@ -124,6 +106,9 @@ FILES IN SVN HEAD RELEVANT FOR A LATER RELEASE:
        utils/viewfg
        tests/aliases.conf
        tests/testall
+       tests/test.cpp ("limit")
+       all subdirs not included in git:
+               bugs, experimental, tests/errorbounds, tests/kees, tests/large, tests/ldpc, tests/uai
 
 DOCUMENTATION READY:
 - bipgraph.h, bipgraph.cpp
diff --git a/TODO b/TODO
index a660621..ba34d11 100644 (file)
--- a/TODO
+++ b/TODO
@@ -1,15 +1,90 @@
 IMPORTANT
 
-- Find out which methods need connected factor graphs (at least TreeEP and JunctionTree).
+- Idea: a FactorGraph and a RegionGraph are often equipped with
+extra properties for nodes and edges. The code to initialize those
+is often quite similar; maybe this can be abstracted to classes
+like ExtFactorGraph and ExtRegionGraph (Ext stands for Extended), e.g.
+template <typename NodeProperties, typename EdgeProperties>
+class ExtFactorGraph : public FactorGraph {
+       public:
+               std::vector<NodeProperties>               nodeProps;
+               std::vector<std::vector<EdgeProperties> > edgeProps;
+       // blabla
+}
+A disadvantage of this approach may be worse cachability.
+A problem is if there are nog properties for nodes (type 1), nodes (type 2)
+or edges. Maybe this can be solved using specializations, or using variadac
+template arguments or something similar?
+Idea: you could define a "class Empty {}", and add some code that checks for
+the typeid, comparing it with Empty, and doing something special in that case 
+(e.g., not allocating memory).
+The main disadvantage of this approach seems to be that it leads to even more
+entanglement.
+
+- THIS SMELLS LIKE A GOOD IDEA: instead of polymorphism by inheritance,
+use polymorphism by template parameterization. For example, the real reason
+you introduced the complicated inheritance scheme was for functions like
+    Factor calcMarginal( const InferenceAlgorithm & obj, const VarSet & ns, bool reInit );
+Now instead, you could use a template function:
+template<typename InferenceAlgorithm>
+    Factor calcMarginal( const InferenceAlgorithm &obj, const VarSet &ns, bool reInit );
+This would assume that InferenceAlgorithm supports certain methods, like
+clone(), grm(), ......
+Ideally, you would use concepts to define different classes of inference
+algorithms with different capabilities, for example the ability to calculate logZ,
+the ability to calculate marginals, the ability to calculate bounds, the ability
+to calculate MAP states, etcetera. Then, use traits classes in order to be able to
+query the capabilities of the model. For example, you would be able to query whether
+the inference algorithm supports calculation of logZ.
+
+- If you ever do a rewrite, make sure that the graphical properties are
+not entangled with the probabilistic properties. E.g., a factor graph
+really should be represented as a bipartite graph, with a separate array
+of variable labels and dimensions, and a seperate array of (pointers to)
+factor tables. In this way, each factor could be implemented differently,
+e.g., we could have some sparse factors, some noisy-OR factors, some dense
+factors, some arbitrary precision factors, etc. Or we could make more use
+of templates to have a more generic factor graph. Maybe in the end, 
+HasA relations are better than IsA relations...
+Also, the current setup is stupid: I wrote a new function that works
+on FactorGraphs, and I had to write boiler plate code for it in graphicalmodel.h
+and in regiongraph.h (which is stupid).
+
+- Clean up 
+
+- Use Boost::uBLAS framework to deal with matrices, especially, with
+2D sparse matrices. See http://www.boost.org/libs/numeric/ublas/doc/matrix_sparse.htm
+and tests/errorbounds/errorbounds3
+
+- Introduce naming scheme:
+all Vars and VarSets should be named v_..., e.g. v_i instead of i
+all Factors should be named f_..., e.g. f_I instead of I
+all indices should be named _..., e.g. _k instead of k
+all iterators should be named i_, e.g. i_i is an iterator to i
+all const_iterators should be named ci_, e.g. ci_i is an iterator to i
+
+- Improve weightedgraph or use Boost::Graph
+
+- Iterations and maxDiff are only interesting for iterative inference
+algorithms. Yet, tests/test wants to know these values in a generic
+way. Maybe we have to think of some way (e.g. using a Properties object)
+to return these values from run(). Then we can simply look whether a
+InferenceAlgorithm supports these fields. What other results could
+we return? Time. MaxDiff during initialization for LC methods.
+Or maybe we could use some traits mechanism which we can ask whether the
+object has _iterations and _maxdiff variables.
+
+- Think about whether the cavity initialization belongs to init() or to run().
+
+- Fix LCLin.
+
+- setFactor and setFactors should only change Probs, not Factors.
+
+- Simplify Properties framework: each Property should be a std::string.
+Each DAIAlg should have its own _properties struct and handle conversion.
 
 - Forwardport the Bart patch
 
-- Get rid of multiple inheritance
-
-- Introduce a namespace libDAI
-
-- Investigate the performance penalty from the use of Properties map
-
 - Another design question that needs to be answered: should a BP object own a
   FactorGraph, or only store a reference to a FactorGraph (which can optimize
   memory), or should it be a FactorGraph with additional variables and functions?
@@ -21,6 +96,12 @@ IMPORTANT
   they construct a new object without copying the FactorGraph. Or they can be made
   simply methods of the general InfAlg class.
 
+- Add comments in example.cpp, add documentation.
+
+- Write documentation.
+
+- Improve error handling.
+
 
 DIFFICULT
 
@@ -28,69 +109,49 @@ DIFFICULT
 
 - Bug: TreeEP::logZ() seems to be wrong (?).
 
-- Bug: strange things happen for at least JTree and TreeEP if the factor graph is not connected.
-  Should be fixed by finding out which methods demand connected factor graphs and these
-  methods should check whether the supplied factor graph is connected.
-
 - Kees' trick for preventing NANs in GBP updates:  if quantity smaller than epsilon, replace by epsilon.
 
-- Bug: MF_SEQRANDOM does not work on alarm2.fg (normalization error: Z = 0)
-
 
 OPTIONAL
 
-- Restore brute force EXACT method.
+- Use Expat XML parser for IO and define a XML based fileformat for .fg files?
+
+- Port to GNU autotool chain?
+
+- Another design question that needs to be answered: should a BP object own a
+  FactorGraph, or only store a reference to a FactorGraph (which can optimize
+  memory), or should it be a FactorGraph with additional variables and functions?
+  Probably, the first or second option would be the best. Since FactorGraphs are
+  assumed to be rather static, it would make sense to store *references* to
+  FactorGraphs inside AIAlgs. Or maybe some smart_ptr? It would make sense to 
+  prevent unnecessary copying of FactorGraphs. Also, general marginalization 
+  functions now copy a complete object. Instead, it would make more sense that
+  they construct a new object without copying the FactorGraph. Or they can be made
+  simply methods of the general InfAlg class.
+
+- Forwardport the Bart/Max patch.
 
 - Alternative implementation of undo factor changes: the only things that have to be
 undone are setting a factor to all ones and setting a factor to a Kronecker delta. This
 could also be implemented in the factor itself, which could maintain its state 
 (one/delta/full) and act accordingly.
 
-- Regression test could be enhanced by 
-  * reporting number of iterations
-  * reporting all single node marginal errors
-
 - Think about the _fac2ind variable: a better implementation of this feature is
 probably possible.
 
 - Changed() methods?
 
-- It might be an idea to forbid modification of the structure of a factor graph
-  (except in constructors). Then we could do away with Regenerate and similar methods.
-
-- Think about whether using the L-1 norm as a distance measure between new and
-  old messages/beliefs would be better, since || P(x_A, x_B) - Q(x_A, xB) ||_1 < epsilon
-  implies that || P(x_A) - Q(x_A) ||_1 < epsilon. However, for a product of messages
-  we have the weaker || P1(x_A) P2(x_A) - Q1(x_A) Q2(x_A) ||_1 < 2 epsilon if
-  || P1(x_A) - Q1(x_A) ||_1 < epsilon AND || P2(x_A) - Q1(x_A) ||_1 < epsilon...
-
-- Implement damping.
-
 - Optimize GBP with precalculated indices.
 
+- Optimize all indices as follows: keep a cache of all indices that have been
+computed (use a hash). Then, instead of computing on the fly, use the precomputed
+ones.
+
 - Variable elimination should be implemented generically using a function
   object that tells you which variable to delete.
 
 - Improve general support for graphs and trees.
 
-- What to do in case of divide-by-zero? For HUGIN, the correct way to handle
-  this is to define x/0 = 0. For other methods, this may not be the case...
-  Maybe a good solution would be to make different flavours of Prob which differ
-  in the way that they handle zeroes (two issues here: normalization of an all-zero
-  factor, which appears to be a sign that something is wrong, and division-by-zero
-  in operator/, operator/= and inverse(), which is probably harmless and can be fixed
-  by defining x/0 = 0. I think it is harmless because usually it occurs by calculating
-  some marginal _in the absence of some factor_ as an intermediate step. Afterwards,
-  this factor is always multiplied in again, which sets any value on these positions
-  to zero. Another difference would be how to deal with underflows: GBP often has 
-  numbers that converge to zero, but instead this gives an underflow. A solution is
-  to set everything smaller than some epsilon to zero. Finally, one could decide to
-  add support for sparse potentials. Oh, and whether it can handle negative entries.
+- Add support for sparse potentials.
 
 - Optimize BP_SEQMAX (it should use a better data structure than a vector for the residuals).
-  Also it does strange things in case of negative potentials.
-
-- Investigate the possibility to calculate logZ *during* the inference algorithm
-  by keeping track of all normalizations.
-
-- Write documentation.
index 8853c94..f5b4886 100644 (file)
@@ -1,6 +1,6 @@
 /*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
     Radboud University Nijmegen, The Netherlands
-    
+
     This file is part of libDAI.
 
     libDAI is free software; you can redistribute it and/or modify
@@ -71,16 +71,16 @@ class BP : public DAIAlgFG {
         }
 
         /// Copy constructor
-        BP( const BP & x ) : DAIAlgFG(x), _edges(x._edges), _maxdiff(x._maxdiff), _iters(x._iters), props(x.props) {}
+        BP( const BP &x ) : DAIAlgFG(x), _edges(x._edges), _maxdiff(x._maxdiff), _iters(x._iters), props(x.props) {}
 
         /// Clone *this (virtual copy constructor)
         virtual BP* clone() const { return new BP(*this); }
 
-        /// Create (virtual constructor)
+        /// Create (virtual default constructor)
         virtual BP* create() const { return new BP(); }
 
         /// Assignment operator
-        BP& operator=( const BP & x ) {
+        BP& operator=( const BP &x ) {
             if( this != &x ) {
                 DAIAlgFG::operator=( x );
                 _edges = x._edges;
@@ -91,35 +91,36 @@ class BP : public DAIAlgFG {
             return *this;
         }
 
-        /// Return number of passes over the factorgraph
-        size_t Iterations() const { return _iters; }
-
-        /// Return maximum difference between single node beliefs for two consecutive iterations
-        double maxDiff() const { return _maxdiff; }
-
         /// Identifies itself for logging purposes
-        std::string identify() const;
+        virtual std::string identify() const;
 
         /// Get single node belief
-        Factor belief( const Var &n ) const;
+        virtual Factor belief( const Var &n ) const;
 
         /// Get general belief
-        Factor belief( const VarSet &n ) const;
+        virtual Factor belief( const VarSet &ns ) const;
 
         /// Get all beliefs
-        std::vector<Factor> beliefs() const;
+        virtual std::vector<Factor> beliefs() const;
 
         /// Get log partition sum
-        Real logZ() const;
+        virtual Real logZ() const;
 
         /// Clear messages and beliefs
-        void init();
+        virtual void init();
 
         /// Clear messages and beliefs corresponding to the nodes in ns
         virtual void init( const VarSet &ns );
 
         /// The actual approximate inference algorithm
-        double run();
+        virtual double run();
+
+        /// Return maximum difference between single node beliefs in the last pass
+        virtual double maxDiff() const { return _maxdiff; }
+
+        /// Return number of passes over the factorgraph
+        virtual size_t Iterations() const { return _iters; }
+
 
         Factor beliefV( size_t i ) const;
         Factor beliefF( size_t I ) const;
@@ -136,15 +137,18 @@ class BP : public DAIAlgFG {
 
         void calcNewMessage( size_t i, size_t _I );
         void updateMessage( size_t i, size_t _I ) {
-            if( props.damping == 0.0 )
+            if( props.damping == 0.0 ) {
                 message(i,_I) = newMessage(i,_I);
-            else
+                residual(i,_I) = 0.0;
+            } else {
                 message(i,_I) = (message(i,_I) ^ props.damping) * (newMessage(i,_I) ^ (1.0 - props.damping));
+                residual(i,_I) = dist( newMessage(i,_I), message(i,_I), Prob::DISTLINF );
+            }
         }
         void findMaxResidual( size_t &i, size_t &_I );
 
-        /// Set Props according to the PropertySet opts, where the values can be stored as std::strings or as the type of the corresponding Props member
         void construct();
+        /// Set Props according to the PropertySet opts, where the values can be stored as std::strings or as the type of the corresponding Props member
         void setProperties( const PropertySet &opts );
         PropertySet getProperties() const;
         std::string printProperties() const;
index ba3edee..1147b22 100644 (file)
@@ -38,14 +38,13 @@ namespace dai {
 /// together with an inference algorithm.
 class InfAlg {
     public:
-        /// Clone (virtual copy constructor)
+        /// Clone *this (virtual copy constructor)
         virtual InfAlg* clone() const = 0;
 
-        /// Create (virtual constructor)
+        /// Create (virtual default constructor)
         virtual InfAlg* create() const = 0;
         
-        /// Virtual desctructor
-        // (this is needed because this class contains virtual functions)
+        /// Virtual desctructor (needed because this class contains virtual functions)
         virtual ~InfAlg() {}
         
         /// Identifies itself for logging purposes
@@ -94,8 +93,11 @@ class InfAlg {
         /// Get const reference to underlying FactorGraph
         virtual const FactorGraph &fg() const = 0;
 
-        /// Return maximum difference between beliefs in the last pass
+        /// Return maximum difference between single node beliefs in the last pass
         virtual double maxDiff() const = 0;
+
+        /// Return number of passes over the factorgraph
+        virtual size_t Iterations() const = 0;
 };
 
 
index 77c3a64..246b433 100644 (file)
@@ -1,6 +1,6 @@
 /*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
     Radboud University Nijmegen, The Netherlands
-    
+
     This file is part of libDAI.
 
     libDAI is free software; you can redistribute it and/or modify
@@ -37,6 +37,8 @@ class ExactInf : public DAIAlgFG {
         struct Properties {
             size_t verbose;
         } props;
+        /// Name of this inference method
+        static const char *Name;
 
     private:
         std::vector<Factor> _beliefsV;
@@ -46,14 +48,6 @@ class ExactInf : public DAIAlgFG {
     public:
         /// Default constructor
         ExactInf() : DAIAlgFG(), props(), _beliefsV(), _beliefsF(), _logZ(0) {}
-        
-        /// Copy constructor
-        ExactInf( const ExactInf &x ) : DAIAlgFG(x), props(x.props), _beliefsV(x._beliefsV), _beliefsF(x._beliefsF), _logZ(x._logZ) {}
-
-        /// Clone (virtual copy constructor)
-        virtual ExactInf* clone() const {
-            return new ExactInf(*this);
-        }
 
         /// Construct from FactorGraph fg and PropertySet opts
         ExactInf( const FactorGraph &fg, const PropertySet &opts ) : DAIAlgFG(fg), props(), _beliefsV(), _beliefsF(), _logZ() {
@@ -61,6 +55,15 @@ class ExactInf : public DAIAlgFG {
             construct();
         }
         
+        /// Copy constructor
+        ExactInf( const ExactInf &x ) : DAIAlgFG(x), props(x.props), _beliefsV(x._beliefsV), _beliefsF(x._beliefsF), _logZ(x._logZ) {}
+
+        /// Clone *this (virtual copy constructor)
+        virtual ExactInf* clone() const { return new ExactInf(*this); }
+
+        /// Create (virtual default constructor)
+        virtual ExactInf* create() const { return new ExactInf(); }
+
         /// Assignment operator
         ExactInf& operator=( const ExactInf &x ) {
             if( this != &x ) {
@@ -73,34 +76,20 @@ class ExactInf : public DAIAlgFG {
             return *this;
         }
 
-        /// Create (virtual constructor)
-        virtual ExactInf* create() const { return new ExactInf(); }
-
-        /// Return maximum difference between single node 
-        /// beliefs for two consecutive iterations
-        virtual double maxDiff() const {
-            DAI_THROW(NOT_IMPLEMENTED);
-            return 0.0;
-        }
-        
         /// Identifies itself for logging purposes
         virtual std::string identify() const;
 
         /// Get single node belief
-        virtual Factor belief( const Var &n ) const {
-            return beliefV( findVar( n ) ); 
-        }
+        virtual Factor belief( const Var &n ) const { return beliefV( findVar( n ) ); }
 
         /// Get general belief
-        virtual Factor belief( const VarSet &n ) const;
+        virtual Factor belief( const VarSet &ns ) const;
 
         /// Get all beliefs
         virtual std::vector<Factor> beliefs() const;
 
         /// Get log partition sum
-        virtual Real logZ() const {
-            return _logZ;
-        }
+        virtual Real logZ() const { return _logZ; }
 
         /// Clear messages and beliefs
         virtual void init();
@@ -113,22 +102,25 @@ class ExactInf : public DAIAlgFG {
         /// The actual approximate inference algorithm
         virtual double run();
 
-        /// Name of this inference method
-        static const char *Name;
+        /// Return maximum difference between single node beliefs in the last pass
+        virtual double maxDiff() const {
+            DAI_THROW(NOT_IMPLEMENTED);
+            return 0.0;
+        }
 
+        /// Return number of passes over the factorgraph
+        virtual size_t Iterations() const { 
+            DAI_THROW(NOT_IMPLEMENTED);
+            return 0; 
+        }
+        
         void construct();
-        void restoreFactors( const VarSet &ns ) { FactorGraph::restoreFactors(ns); init(ns); }
         void setProperties( const PropertySet &opts );
         PropertySet getProperties() const;
         std::string printProperties() const;
 
-        Factor beliefV( size_t i ) const { 
-            return _beliefsV[i]; 
-        }
-
-        Factor beliefF( size_t I ) const { 
-            return _beliefsF[I]; 
-        }
+        Factor beliefV( size_t i ) const { return _beliefsV[i]; }
+        Factor beliefF( size_t I ) const { return _beliefsF[I]; }
 };
 
 
index 0260864..7e7f4d9 100644 (file)
@@ -67,12 +67,12 @@ class FactorGraph {
         }
         virtual ~FactorGraph() {}
 
+        /// Clone *this (virtual copy constructor)
+        virtual FactorGraph* clone() const { return new FactorGraph(); }
+
         /// Create (virtual default constructor)
         virtual FactorGraph* create() const { return new FactorGraph(*this); }
 
-        /// Clone (virtual copy constructor)
-        virtual FactorGraph* clone() const { return new FactorGraph(); }
-
         // aliases
         Var & var(size_t i) { return vars[i]; }
         /// Get const reference to i'th variable
index 112c881..159c2f7 100644 (file)
@@ -1,6 +1,6 @@
 /*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
     Radboud University Nijmegen, The Netherlands
-    
+
     This file is part of libDAI.
 
     libDAI is free software; you can redistribute it and/or modify
@@ -40,53 +40,89 @@ class HAK : public DAIAlgRG {
         std::vector<Factor>                _Qb;
         std::vector<std::vector<Factor> >  _muab;
         std::vector<std::vector<Factor> >  _muba;
+        /// Maximum difference encountered so far
+        double _maxdiff;
+        /// Number of iterations needed
+        size_t _iters;
 
     public:
         struct Properties {
             size_t verbose;
             size_t maxiter;
             double tol;
+            double damping;
             DAI_ENUM(ClustersType,MIN,DELTA,LOOP)
             ClustersType clusters;
             bool doubleloop;
             size_t loopdepth;
         } props;
-        double maxdiff;
-        
+        /// Name of this inference method
+        static const char *Name;
+
     public:
         /// Default constructor
-        HAK() : DAIAlgRG(), _Qa(), _Qb(), _muab(), _muba(), props(), maxdiff() {}
+        HAK() : DAIAlgRG(), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {}
 
-        /// Copy constructor
-        HAK(const HAK & x) : DAIAlgRG(x), _Qa(x._Qa), _Qb(x._Qb), _muab(x._muab), _muba(x._muba), props(x.props), maxdiff(x.maxdiff) {}
+        /// Construct from FactorGraph fg and PropertySet opts
+        HAK( const FactorGraph &fg, const PropertySet &opts );
 
-        /// Clone function
-        HAK* clone() const { return new HAK(*this); }
-        
-        /// Create (virtual constructor)
-        virtual HAK* create() const { return new HAK(); }
+        /// Construct from RegionGraph rg and PropertySet opts
+        HAK( const RegionGraph &rg, const PropertySet &opts );
+
+        /// Copy constructor
+        HAK( const HAK &x ) : DAIAlgRG(x), _Qa(x._Qa), _Qb(x._Qb), _muab(x._muab), _muba(x._muba), _maxdiff(x._maxdiff), _iters(x._iters), props(x.props) {}
 
-        /// Construct from RegionGraph
-        HAK(const RegionGraph & rg, const PropertySet &opts);
+        /// Clone *this (virtual copy constructor)
+        virtual HAK* clone() const { return new HAK(*this); }
 
-        /// Construct from RactorGraph using "clusters" option
-        HAK(const FactorGraph & fg, const PropertySet &opts);
+        /// Create (virtual default constructor)
+        virtual HAK* create() const { return new HAK(); }
 
         /// Assignment operator
-        HAK & operator=(const HAK & x) {
+        HAK& operator=( const HAK &x ) {
             if( this != &x ) {
-                DAIAlgRG::operator=(x);
-                _Qa    = x._Qa;
-                _Qb    = x._Qb;
-                _muab  = x._muab;
-                _muba  = x._muba;
-                props  = x.props;
-                maxdiff = x.maxdiff;
+                DAIAlgRG::operator=( x );
+                _Qa      = x._Qa;
+                _Qb      = x._Qb;
+                _muab    = x._muab;
+                _muba    = x._muba;
+                _maxdiff = x._maxdiff;
+                _iters   = x._iters;
+                props    = x.props;
             }
             return *this;
         }
-        
-        static const char *Name;
+
+        /// Identifies itself for logging purposes
+        virtual std::string identify() const;
+
+        /// Get single node belief
+        virtual Factor belief( const Var &n ) const;
+
+        /// Get general belief
+        virtual Factor belief( const VarSet &ns ) const;
+
+        /// Get all beliefs
+        virtual std::vector<Factor> beliefs() const;
+
+        /// Get log partition sum
+        virtual Real logZ() const;
+
+        /// Clear messages and beliefs
+        virtual void init();
+
+        /// Clear messages and beliefs corresponding to the nodes in ns
+        virtual void init( const VarSet &ns );
+
+        /// The actual approximate inference algorithm
+        virtual double run();
+
+        /// Return maximum difference between single node beliefs in the last pass
+        virtual double maxDiff() const { return _maxdiff; }
+
+        /// Return number of passes over the factorgraph
+        virtual size_t Iterations() const { return _iters; }
+
 
         Factor & muab( size_t alpha, size_t _beta ) { return _muab[alpha][_beta]; }
         Factor & muba( size_t alpha, size_t _beta ) { return _muba[alpha][_beta]; }
@@ -95,21 +131,10 @@ class HAK : public DAIAlgRG {
 
         double doGBP();
         double doDoubleLoop();
-        double run();
-        void init();
-        /// Clear messages and beliefs corresponding to the nodes in ns
-        virtual void init( const VarSet &ns );
-        std::string identify() const;
-        Factor belief( const Var &n ) const;
-        Factor belief( const VarSet &ns ) const;
-        std::vector<Factor> beliefs() const;
-        Real logZ () const;
 
-        void restoreFactors( const VarSet &ns ) { RegionGraph::restoreFactors( ns ); init( ns ); }
         void setProperties( const PropertySet &opts );
         PropertySet getProperties() const;
         std::string printProperties() const;
-        double maxDiff() const { return maxdiff; }
 
     private:
         void constructMessages();
index ccb733c..eff98de 100644 (file)
@@ -1,6 +1,6 @@
 /*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
     Radboud University Nijmegen, The Netherlands
-    
+
     This file is part of libDAI.
 
     libDAI is free software; you can redistribute it and/or modify
@@ -59,14 +59,20 @@ class JTree : public DAIAlgRG {
         /// Default constructor
         JTree() : DAIAlgRG(), _RTree(), _Qa(), _Qb(), _mes(), _logZ(), props() {}
 
-        /// Construct JTree object using the specified properties
+        /// Construct from FactorGraph fg and PropertySet opts
         JTree( const FactorGraph &fg, const PropertySet &opts, bool automatic=true );
 
         /// Copy constructor
         JTree( const JTree &x ) : DAIAlgRG(x), _RTree(x._RTree), _Qa(x._Qa), _Qb(x._Qb), _mes(x._mes), _logZ(x._logZ), props(x.props) {}
 
+        /// Clone *this (virtual copy constructor)
+        virtual JTree* clone() const { return new JTree(*this); }
+
+        /// Create (virtual default constructor)
+        virtual JTree* create() const { return new JTree(); }
+
         /// Assignment operator
-        JTree & operator=( const JTree &x ) {
+        JTree& operator=( const JTree &x ) {
             if( this != &x ) {
                 DAIAlgRG::operator=( x );
                 _RTree  = x._RTree;
@@ -79,41 +85,35 @@ class JTree : public DAIAlgRG {
             return *this;
         }
 
-        /// Clone (virtual copy constructor)
-        virtual JTree* clone() const { return new JTree(*this); }
-
-        /// Create (virtual constructor)
-        virtual JTree* create() const { return new JTree(); }
-
-        /// Return number of passes over the factorgraph
-        virtual size_t Iterations() const { return 1UL; }
-
-        /// Return maximum difference between single node beliefs for two consecutive iterations
-        virtual double maxDiff() const { return 0.0; }
-
         /// Identifies itself for logging purposes
-        std::string identify() const;
+        virtual std::string identify() const;
 
         /// Get single node belief
-        Factor belief( const Var &n ) const;
+        virtual Factor belief( const Var &n ) const;
 
         /// Get general belief
-        Factor belief( const VarSet &ns ) const;
+        virtual Factor belief( const VarSet &ns ) const;
 
         /// Get all beliefs
-        std::vector<Factor> beliefs() const;
+        virtual std::vector<Factor> beliefs() const;
 
         /// Get log partition sum
-        Real logZ() const;
+        virtual Real logZ() const;
 
         /// Clear messages and beliefs
-        void init() {}
+        virtual void init() {}
 
         /// Clear messages and beliefs corresponding to the nodes in ns
         virtual void init( const VarSet &/*ns*/ ) {}
 
         /// The actual approximate inference algorithm
-        double run();
+        virtual double run();
+
+        /// Return maximum difference between single node beliefs in the last pass
+        virtual double maxDiff() const { return 0.0; }
+
+        /// Return number of passes over the factorgraph
+        virtual size_t Iterations() const { return 1UL; }
 
 
         void GenerateJT( const std::vector<VarSet> &Cliques );
index 3ec9d30..b49da41 100644 (file)
@@ -44,91 +44,109 @@ class LC : public DAIAlgFG {
         /// Single variable beliefs
         std::vector<Factor>      _beliefs;
 
+        /// Maximum difference encountered so far
+        double                  _maxdiff;
+        /// Number of iterations needed
+        size_t                  _iters;
+
     public:
         struct Properties {
             size_t verbose;
             size_t maxiter;
             double tol;
+            bool reinit;
+            double damping;
             DAI_ENUM(CavityType,FULL,PAIR,PAIR2,UNIFORM)
             CavityType cavity;
             DAI_ENUM(UpdateType,SEQFIX,SEQRND,NONE)
             UpdateType updates;
-            std::string cavainame;
-            PropertySet cavaiopts;
-            bool reinit;
+            std::string cavainame;      // FIXME: needs assignment operator?
+            PropertySet cavaiopts;      // FIXME: needs assignment operator?
         } props;
-        double maxdiff;
+        /// Name of this inference method
+        static const char *Name;
 
     public:
         /// Default constructor
-        LC() : DAIAlgFG(), _pancakes(), _cavitydists(), _phis(), _beliefs(), props(), maxdiff() {}
+        LC() : DAIAlgFG(), _pancakes(), _cavitydists(), _phis(), _beliefs(), _maxdiff(), _iters(), props() {}
+
+        /// Construct from FactorGraph fg and PropertySet opts
+        LC( const FactorGraph &fg, const PropertySet &opts );
+
         /// Copy constructor
-        LC(const LC & x) : DAIAlgFG(x), _pancakes(x._pancakes), _cavitydists(x._cavitydists), _phis(x._phis), _beliefs(x._beliefs), props(x.props), maxdiff(x.maxdiff) {}
-        /// Clone function
-        LC* clone() const { return new LC(*this); }
-        /// Create (virtual constructor)
+        LC( const LC &x ) : DAIAlgFG(x), _pancakes(x._pancakes), _cavitydists(x._cavitydists), _phis(x._phis), _beliefs(x._beliefs), _maxdiff(x._maxdiff), _iters(x._iters), props(x.props) {}
+
+        /// Clone *this (virtual copy constructor)
+        virtual LC* clone() const { return new LC(*this); }
+
+        /// Create (virtual default constructor)
         virtual LC* create() const { return new LC(); }
-        /// Construct LC object from a FactorGraph and parameters
-        LC( const FactorGraph & fg, const PropertySet &opts );
+
         /// Assignment operator
-        LC& operator=(const LC & x) {
+        LC& operator=( const LC &x ) {
             if( this != &x ) {
-                DAIAlgFG::operator=(x);
-                _pancakes       = x._pancakes;
-                _cavitydists    = x._cavitydists;
-                _phis           = x._phis;
-                _beliefs        = x._beliefs;
-                props           = x.props;
-                maxdiff         = x.maxdiff;
+                DAIAlgFG::operator=( x );
+                _pancakes     = x._pancakes;
+                _cavitydists  = x._cavitydists;
+                _phis         = x._phis;
+                _beliefs      = x._beliefs;
+                _maxdiff      = x._maxdiff;
+                _iters        = x._iters;
+                props         = x.props;
             }
             return *this;
         }
 
-        static const char *Name;
-        double CalcCavityDist( size_t i, const std::string &name, const PropertySet &opts );
-        double InitCavityDists( const std::string &name, const PropertySet &opts );
-        long SetCavityDists( std::vector<Factor> &Q );
+        /// Identifies itself for logging purposes
+        virtual std::string identify() const;
 
-        void init();
-        /// Clear messages and beliefs corresponding to the nodes in ns
-        virtual void init( const VarSet &/*ns*/ ) {
-            DAI_THROW(NOT_IMPLEMENTED);
-        }
-        Factor NewPancake (size_t i, size_t _I, bool & hasNaNs);
-        double run();
+        /// Get single node belief
+        virtual Factor belief( const Var &n ) const { return( _beliefs[findVar(n)] ); }
 
-        std::string identify() const;
-        Factor belief (const Var &n) const { return( _beliefs[findVar(n)] ); }
-        Factor belief (const VarSet &/*ns*/) const { 
+        /// Get general belief
+        virtual Factor belief( const VarSet &/*ns*/ ) const {
             DAI_THROW(NOT_IMPLEMENTED);
             return Factor(); 
         }
-        std::vector<Factor> beliefs() const { return _beliefs; }
-        Real logZ() const { 
+
+        /// Get all beliefs
+        virtual std::vector<Factor> beliefs() const { return _beliefs; }
+
+        /// Get log partition sum
+        virtual Real logZ() const { 
             DAI_THROW(NOT_IMPLEMENTED);
             return 0.0; 
         }
+
+        /// Clear messages and beliefs
+        virtual void init();
+
+        /// Clear messages and beliefs corresponding to the nodes in ns
+        virtual void init( const VarSet &/*ns*/ ) { init(); }
+
+        /// The actual approximate inference algorithm
+        virtual double run();
+
+        /// Return maximum difference between single node beliefs in the last pass
+        virtual double maxDiff() const { return _maxdiff; }
+
+        /// Return number of passes over the factorgraph
+        virtual size_t Iterations() const { return _iters; }
+
+        double CalcCavityDist( size_t i, const std::string &name, const PropertySet &opts );
+        double InitCavityDists( const std::string &name, const PropertySet &opts );
+        long SetCavityDists( std::vector<Factor> &Q );
+
+        Factor NewPancake (size_t i, size_t _I, bool & hasNaNs);
+
         void CalcBelief (size_t i);
         const Factor &belief (size_t i) const { return _beliefs[i]; };
         const Factor &pancake (size_t i) const { return _pancakes[i]; };
         const Factor &cavitydist (size_t i) const { return _cavitydists[i]; };
 
-        void clamp( const Var &/*n*/, size_t /*i*/ ) { 
-            DAI_THROW(NOT_IMPLEMENTED);
-        }
-        void restoreFactors( const VarSet &/*ns*/ ) { 
-            DAI_THROW(NOT_IMPLEMENTED);
-        }
-        void backupFactors( const VarSet &/*ns*/ ) { 
-            DAI_THROW(NOT_IMPLEMENTED);
-        }
-        virtual void makeCavity(const Var & /*n*/) { 
-            DAI_THROW(NOT_IMPLEMENTED);
-        }
         void setProperties( const PropertySet &opts );
         PropertySet getProperties() const;
         std::string printProperties() const;
-        double maxDiff() const { return maxdiff; }
 };
 
 
index 60a5b1f..6daef4b 100644 (file)
@@ -1,6 +1,6 @@
 /*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
     Radboud University Nijmegen, The Netherlands
-    
+
     This file is part of libDAI.
 
     libDAI is free software; you can redistribute it and/or modify
@@ -53,7 +53,7 @@ class MF : public DAIAlgFG {
         /// Default constructor
         MF() : DAIAlgFG(), _beliefs(), _maxdiff(0.0), _iters(0U), props() {}
 
-        // construct MF object from FactorGraph
+        /// Construct from FactorGraph fg and PropertySet opts
         MF( const FactorGraph &fg, const PropertySet &opts ) : DAIAlgFG(fg), _beliefs(), _maxdiff(0.0), _iters(0U), props() {
             setProperties( opts );
             construct();
@@ -62,8 +62,14 @@ class MF : public DAIAlgFG {
         /// Copy constructor
         MF( const MF &x ) : DAIAlgFG(x), _beliefs(x._beliefs), _maxdiff(x._maxdiff), _iters(x._iters), props(x.props) {}
 
+        /// Clone *this (virtual copy constructor)
+        virtual MF* clone() const { return new MF(*this); }
+
+        /// Create (virtual default constructor)
+        virtual MF* create() const { return new MF(); }
+
         /// Assignment operator
-        MF & operator=( const MF &x ) {
+        MF& operator=( const MF &x ) {
             if( this != &x ) {
                 DAIAlgFG::operator=( x );
                 _beliefs = x._beliefs;
@@ -74,42 +80,38 @@ class MF : public DAIAlgFG {
             return *this;
         }
 
-        /// Clone *this (virtual copy constructor)
-        virtual MF* clone() const { return new MF(*this); }
-
-        /// Create (virtual constructor)
-        virtual MF* create() const { return new MF(); }
-
-        /// Return number of passes over the factorgraph needed
-        virtual size_t Iterations() const { return _iters; }
-
-        /// Return maximum difference between single node beliefs for two consecutive iterations
-        double maxDiff() const { return _maxdiff; }
-
-        /// Identify *this for logging purposes
-        std::string identify() const;
+        /// Identifies itself for logging purposes
+        virtual std::string identify() const;
 
         /// Get single node belief
-        Factor belief( const Var &n ) const;
+        virtual Factor belief( const Var &n ) const;
 
         /// Get general belief
-        Factor belief( const VarSet &ns ) const;
+        virtual Factor belief( const VarSet &ns ) const;
 
         /// Get all beliefs
-        std::vector<Factor> beliefs() const;
+        virtual std::vector<Factor> beliefs() const;
 
         /// Get log partition sum
-        Real logZ() const;
+        virtual Real logZ() const;
 
-        void construct();
-
-        void init();
+        /// Clear messages and beliefs
+        virtual void init();
 
         /// Clear messages and beliefs corresponding to the nodes in ns
         virtual void init( const VarSet &ns );
 
         /// The actual approximate inference algorithm
-        double run();
+        virtual double run();
+
+        /// Return maximum difference between single node beliefs in the last pass
+        virtual double maxDiff() const { return _maxdiff; }
+
+        /// Return number of passes over the factorgraph
+        virtual size_t Iterations() const { return _iters; }
+
+
+        void construct();
 
         void setProperties( const PropertySet &opts );
         PropertySet getProperties() const;
index 2b55e96..34ec957 100644 (file)
@@ -1,6 +1,6 @@
 /*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
     Radboud University Nijmegen, The Netherlands
-    
+
     This file is part of libDAI.
 
     libDAI is free software; you can redistribute it and/or modify
@@ -56,6 +56,9 @@ class MR : public DAIAlgFG {
 
         std::vector<double> Mag;
 
+        double _maxdiff;
+        size_t _iters;
+
     public:
         struct Properties {
             size_t verbose;
@@ -65,36 +68,92 @@ class MR : public DAIAlgFG {
             UpdateType updates;
             InitType inits;
         } props;
-        double maxdiff;
+        static const char *Name;
 
     public:
-        MR() {}
-        MR( const FactorGraph & fg, const PropertySet &opts );
-        void init(size_t Nin, double *_w, double *_th);
-        void makekindex();
-        void read_files();
-        void init_cor();
-        double init_cor_resp();
-        void solvemcav();
-        void solveM();
-        double run();
-        Factor belief( const Var &n ) const;
-        Factor belief( const VarSet &/*ns*/ ) const { 
+        /// Default constructor
+        MR() : DAIAlgFG(), supported(), con(), nb(), tJ(), theta(), M(), kindex(), cors(), N(), Mag(), _maxdiff(), _iters(), props() {}
+
+        /// Construct from FactorGraph fg and PropertySet opts
+        MR( const FactorGraph &fg, const PropertySet &opts );
+
+        /// Copy constructor
+        MR( const MR &x ) : DAIAlgFG(x), supported(x.supported), con(x.con), nb(x.nb), tJ(x.tJ), theta(x.theta), M(x.M), kindex(x.kindex), cors(x.cors), N(x.N), Mag(x.Mag), _maxdiff(x._maxdiff), _iters(x._iters), props(x.props) {}
+
+        /// Clone *this (virtual copy constructor)
+        virtual MR* clone() const { return new MR(*this); }
+
+        /// Create (virtual default constructor)
+        virtual MR* create() const { return new MR(); }
+
+        /// Assignment operator
+        MR& operator=( const MR &x ) {
+            if( this != &x ) {
+                DAIAlgFG::operator=(x);
+                supported = x.supported;
+                con       = x.con; 
+                nb        = x.nb;
+                tJ        = x.tJ;
+                theta     = x.theta;
+                M         = x.M;
+                kindex    = x.kindex;
+                cors      = x.cors;
+                N         = x.N;
+                Mag       = x.Mag;
+                _maxdiff  = x._maxdiff;
+                _iters    = x._iters;
+                props     = x.props;
+            }
+            return *this;
+        }
+
+        /// Identifies itself for logging purposes
+        virtual std::string identify() const;
+
+        /// Get single node belief
+        virtual Factor belief( const Var &n ) const;
+
+        /// Get general belief
+        virtual Factor belief( const VarSet &/*ns*/ ) const { 
             DAI_THROW(NOT_IMPLEMENTED);
             return Factor(); 
         }
-        std::vector<Factor> beliefs() const;
-        Real logZ() const { 
+
+        /// Get all beliefs
+        virtual std::vector<Factor> beliefs() const;
+
+        /// Get log partition sum
+        virtual Real logZ() const { 
             DAI_THROW(NOT_IMPLEMENTED);
             return 0.0; 
         }
-        void init() {}
+
+        /// Clear messages and beliefs
+        virtual void init() {}
+
         /// Clear messages and beliefs corresponding to the nodes in ns
         virtual void init( const VarSet &/*ns*/ ) {
             DAI_THROW(NOT_IMPLEMENTED);
         }
-        static const char *Name;
-        std::string identify() const;
+
+        /// The actual approximate inference algorithm
+        virtual double run();
+
+        /// Return maximum difference between single node beliefs in the last pass
+        virtual double maxDiff() const { return _maxdiff; }
+
+        /// Return number of passes over the factorgraph
+        virtual size_t Iterations() const { return _iters; }
+
+
+        void init(size_t Nin, double *_w, double *_th);
+        void makekindex();
+        void read_files();
+        void init_cor();
+        double init_cor_resp();
+        void solvemcav();
+        void solveM();
+
         double _tJ(size_t i, sub_nb A);
 
         double Omega(size_t i, size_t _j, size_t _l);
@@ -107,14 +166,10 @@ class MR : public DAIAlgFG {
         void sum_subs(size_t j, sub_nb A, double *sum_even, double *sum_odd);
 
         double sign(double a) { return (a >= 0) ? 1.0 : -1.0; }
-        MR* clone() const { return new MR(*this); }
-        /// Create (virtual constructor)
-        virtual MR* create() const { return new MR(); }
-
+        
         void setProperties( const PropertySet &opts );
         PropertySet getProperties() const;
         std::string printProperties() const;
-        double maxDiff() const { return maxdiff; }
 }; 
 
 
index 051d7b9..9d0d924 100644 (file)
@@ -136,12 +136,12 @@ class RegionGraph : public FactorGraph {
             return *this;
         }
 
+        /// Clone *this (virtual copy constructor)
+        virtual RegionGraph* clone() const { return new RegionGraph(*this); }
+
         /// Create (virtual default constructor)
         virtual RegionGraph* create() const { return new RegionGraph(); }
 
-        /// Clone (virtual copy constructor)
-        virtual RegionGraph* clone() const { return new RegionGraph(*this); }
-
         /// Set the content of the I'th factor and make a backup of its old content if backup == true
         virtual void setFactor( size_t I, const Factor &newFactor, bool backup = false ) {
             FactorGraph::setFactor( I, newFactor, backup ); 
index 9120085..fd039a6 100644 (file)
@@ -1,6 +1,6 @@
 /*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
     Radboud University Nijmegen, The Netherlands
-    
+
     This file is part of libDAI.
 
     libDAI is free software; you can redistribute it and/or modify
 namespace dai {
 
 
-class TreeEPSubTree {
-    protected:
-        std::vector<Factor>  _Qa;
-        std::vector<Factor>  _Qb;
-        DEdgeVec             _RTree;
-        std::vector<size_t>  _a;        // _Qa[alpha]  <->  superTree._Qa[_a[alpha]]
-        std::vector<size_t>  _b;        // _Qb[beta]   <->  superTree._Qb[_b[beta]]
-                                        // _Qb[beta]   <->  _RTree[beta]    
-        const Factor *       _I;
-        VarSet               _ns;
-        VarSet               _nsrem;
-        double               _logZ;
-        
-        
-    public:
-        TreeEPSubTree() : _Qa(), _Qb(), _RTree(), _a(), _b(), _I(NULL), _ns(), _nsrem(), _logZ(0.0) {}
-        TreeEPSubTree( const TreeEPSubTree &x) : _Qa(x._Qa), _Qb(x._Qb), _RTree(x._RTree), _a(x._a), _b(x._b), _I(x._I), _ns(x._ns), _nsrem(x._nsrem), _logZ(x._logZ) {}
-        TreeEPSubTree & operator=( const TreeEPSubTree& x ) {
-            if( this != &x ) {
-                _Qa         = x._Qa;
-                _Qb         = x._Qb;
-                _RTree      = x._RTree;
-                _a          = x._a;
-                _b          = x._b;
-                _I          = x._I;
-                _ns         = x._ns;
-                _nsrem      = x._nsrem;
-                _logZ       = x._logZ;
-            }
-            return *this;
-        }
-
-        TreeEPSubTree( const DEdgeVec &subRTree, const DEdgeVec &jt_RTree, const std::vector<Factor> &jt_Qa, const std::vector<Factor> &jt_Qb, const Factor *I );
-        void init();
-        void InvertAndMultiply( const std::vector<Factor> &Qa, const std::vector<Factor> &Qb );
-        void HUGIN_with_I( std::vector<Factor> &Qa, std::vector<Factor> &Qb );
-        double logZ( const std::vector<Factor> &Qa, const std::vector<Factor> &Qb ) const;
-        const Factor *& I() { return _I; }
-};
-
-
 class TreeEP : public JTree {
     protected:
-        std::map<size_t, TreeEPSubTree>  _Q;
+        /// Maximum difference encountered so far
+        double                  _maxdiff;
+        /// Number of iterations needed
+        size_t                  _iters;
 
     public:
         struct Properties {
@@ -91,52 +53,116 @@ class TreeEP : public JTree {
             double tol;
             DAI_ENUM(TypeType,ORG,ALT)
             TypeType type;
-        } props;
-        double maxdiff;
+        } props; // FIXME: should be props2 because of conflict with JTree::props?
+        /// Name of this inference method
+        static const char *Name;
+
+    protected:
+        class TreeEPSubTree {
+            protected:
+                std::vector<Factor>  _Qa;
+                std::vector<Factor>  _Qb;
+                DEdgeVec             _RTree;
+                std::vector<size_t>  _a;        // _Qa[alpha]  <->  superTree._Qa[_a[alpha]]
+                std::vector<size_t>  _b;        // _Qb[beta]   <->  superTree._Qb[_b[beta]]
+                                                // _Qb[beta]   <->  _RTree[beta]    
+                const Factor *       _I;
+                VarSet               _ns;
+                VarSet               _nsrem;
+                double               _logZ;
+                
+                
+            public:
+                TreeEPSubTree() : _Qa(), _Qb(), _RTree(), _a(), _b(), _I(NULL), _ns(), _nsrem(), _logZ(0.0) {}
+                TreeEPSubTree( const TreeEPSubTree &x) : _Qa(x._Qa), _Qb(x._Qb), _RTree(x._RTree), _a(x._a), _b(x._b), _I(x._I), _ns(x._ns), _nsrem(x._nsrem), _logZ(x._logZ) {}
+                TreeEPSubTree & operator=( const TreeEPSubTree& x ) {
+                    if( this != &x ) {
+                        _Qa         = x._Qa;
+                        _Qb         = x._Qb;
+                        _RTree      = x._RTree;
+                        _a          = x._a;
+                        _b          = x._b;
+                        _I          = x._I;
+                        _ns         = x._ns;
+                        _nsrem      = x._nsrem;
+                        _logZ       = x._logZ;
+                    }
+                    return *this;
+                }
+
+                TreeEPSubTree( const DEdgeVec &subRTree, const DEdgeVec &jt_RTree, const std::vector<Factor> &jt_Qa, const std::vector<Factor> &jt_Qb, const Factor *I );
+                void init();
+                void InvertAndMultiply( const std::vector<Factor> &Qa, const std::vector<Factor> &Qb );
+                void HUGIN_with_I( std::vector<Factor> &Qa, std::vector<Factor> &Qb );
+                double logZ( const std::vector<Factor> &Qa, const std::vector<Factor> &Qb ) const;
+                const Factor *& I() { return _I; }
+        };
+
+        std::map<size_t, TreeEPSubTree>  _Q;
 
     public:
         /// Default constructor
-        TreeEP() : JTree(), _Q(), props(), maxdiff() {};
+        TreeEP() : JTree(), _maxdiff(0.0), _iters(0), props(), _Q() {}
+
+        /// Construct from FactorGraph fg and PropertySet opts
+        TreeEP( const FactorGraph &fg, const PropertySet &opts );
+
         /// Copy constructor
-        TreeEP( const TreeEP& x ) : JTree(x), _Q(x._Q), props(x.props), maxdiff(x.maxdiff) {
+        TreeEP( const TreeEP &x ) : JTree(x), _maxdiff(x._maxdiff), _iters(x._iters), props(x.props), _Q(x._Q) {
             for( size_t I = 0; I < nrFactors(); I++ )
                 if( offtree( I ) )
                     _Q[I].I() = &factor(I);
         }
-        TreeEP* clone() const { return new TreeEP(*this); }
-        /// Create (virtual constructor)
+
+        /// Clone *this (virtual copy constructor)
+        virtual TreeEP* clone() const { return new TreeEP(*this); }
+
+        /// Create (virtual default constructor)
         virtual TreeEP* create() const { return new TreeEP(); }
-        TreeEP & operator=( const TreeEP& x ) {
+
+        /// Assignment operator
+        TreeEP& operator=( const TreeEP &x ) {
             if( this != &x ) {
-                JTree::operator=(x);
-                _Q   = x._Q;
+                JTree::operator=( x );
+                _maxdiff = x._maxdiff;
+                _iters   = x._iters;
+                props    = x.props;
+                _Q       = x._Q;
                 for( size_t I = 0; I < nrFactors(); I++ )
                     if( offtree( I ) )
                         _Q[I].I() = &factor(I);
-                props = x.props;
-                maxdiff = x.maxdiff;
             }
             return *this;
         }
-        TreeEP( const FactorGraph &fg, const PropertySet &opts );
-        void ConstructRG( const DEdgeVec &tree );
 
-        static const char *Name;
-        std::string identify() const;
-        void init();
+        /// Identifies itself for logging purposes
+        virtual std::string identify() const;
+
+        /// Get log partition sum
+        virtual Real logZ() const;
+
+        /// Clear messages and beliefs
+        virtual void init();
+
         /// Clear messages and beliefs corresponding to the nodes in ns
         virtual void init( const VarSet &/*ns*/ ) { init(); }
-        double run();
-        Real logZ() const;
 
-        bool offtree( size_t I ) const { return (fac2OR[I] == -1U); }
+        /// The actual approximate inference algorithm
+        virtual double run();
+
+        /// Return maximum difference between single node beliefs in the last pass
+        virtual double maxDiff() const { return _maxdiff; }
 
-        void restoreFactors( const VarSet &ns ) { RegionGraph::restoreFactors( ns ); init( ns ); }
+        /// Return number of passes over the factorgraph
+        virtual size_t Iterations() const { return _iters; }
+
+
+        void ConstructRG( const DEdgeVec &tree );
+        bool offtree( size_t I ) const { return (fac2OR[I] == -1U); }
 
         void setProperties( const PropertySet &opts );
         PropertySet getProperties() const;
         std::string printProperties() const;
-        double maxDiff() const { return maxdiff; }
 };
 
 
index 8162bf7..fb7220a 100644 (file)
@@ -247,7 +247,6 @@ double BP::run() {
                 size_t i, _I;
                 findMaxResidual( i, _I );
                 updateMessage( i, _I );
-                residual( i, _I ) = 0.0;
 
                 // I->i has been updated, which means that residuals for all
                 // J->j with J in nb[i]\I and j in nb[J]\i have to be updated
@@ -291,7 +290,7 @@ double BP::run() {
         }
 
         if( props.verbose >= 3 )
-            cout << "BP::run:  maxdiff " << diffs.maxDiff() << " after " << _iters+1 << " passes" << endl;
+            cout << Name << "::run:  maxdiff " << diffs.maxDiff() << " after " << _iters+1 << " passes" << endl;
     }
 
     if( diffs.maxDiff() > _maxdiff )
@@ -301,11 +300,11 @@ double BP::run() {
         if( diffs.maxDiff() > props.tol ) {
             if( props.verbose == 1 )
                 cout << endl;
-                cout << "BP::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.maxDiff() << endl;
+                cout << Name << "::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << diffs.maxDiff() << endl;
         } else {
             if( props.verbose >= 3 )
-                cout << "BP::run:  ";
-                cout << "converged in " << _iters << " passes (" << toc() - tic << " clocks)." << endl;
+                cout << Name << "::run:  ";
+                cout << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
         }
     }
 
index 3d2f887..adcc545 100644 (file)
@@ -52,6 +52,10 @@ void HAK::setProperties( const PropertySet &opts ) {
         props.loopdepth = opts.getStringAs<size_t>("loopdepth");
     else
         assert( props.clusters != Properties::ClustersType::LOOP );
+    if( opts.hasKey("damping") )
+        props.damping = opts.getStringAs<double>("damping");
+    else
+        props.damping = 0.0;
 }
 
 
@@ -63,6 +67,7 @@ PropertySet HAK::getProperties() const {
     opts.Set( "doubleloop", props.doubleloop );
     opts.Set( "clusters", props.clusters );
     opts.Set( "loopdepth", props.loopdepth );
+    opts.Set( "damping", props.damping );
     return opts;
 }
 
@@ -75,7 +80,8 @@ string HAK::printProperties() const {
     s << "verbose=" << props.verbose << ",";
     s << "doubleloop=" << props.doubleloop << ",";
     s << "clusters=" << props.clusters << ",";
-    s << "loopdepth=" << props.loopdepth << "]";
+    s << "loopdepth=" << props.loopdepth << ",";
+    s << "damping=" << props.damping << "]";
     return s.str();
 }
 
@@ -111,7 +117,7 @@ void HAK::constructMessages() {
 }
 
 
-HAK::HAK(const RegionGraph & rg, const PropertySet &opts ) : DAIAlgRG(rg) {
+HAK::HAK( const RegionGraph &rg, const PropertySet &opts ) : DAIAlgRG(rg), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
     setProperties( opts );
 
     constructMessages();
@@ -121,7 +127,7 @@ HAK::HAK(const RegionGraph & rg, const PropertySet &opts ) : DAIAlgRG(rg) {
 void HAK::findLoopClusters( const FactorGraph & fg, std::set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars ) {
     for( VarSet::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
         VarSet ind = fg.delta( fg.findVar( *in ) );
-        if( (newcl.size()) >= 2 && (ind >> root) ) {
+        if( (newcl.size()) >= 2 && ind.contains( root ) ) {
             allcl.insert( newcl | *in );
         }
         else if( length > 1 )
@@ -130,7 +136,7 @@ void HAK::findLoopClusters( const FactorGraph & fg, std::set<VarSet> &allcl, Var
 }
 
 
-HAK::HAK(const FactorGraph & fg, const PropertySet &opts) : DAIAlgRG(), props(), maxdiff(0.0) {
+HAK::HAK(const FactorGraph & fg, const PropertySet &opts) : DAIAlgRG(), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
     setProperties( opts );
 
     vector<VarSet> cl;
@@ -150,7 +156,7 @@ HAK::HAK(const FactorGraph & fg, const PropertySet &opts) : DAIAlgRG(), props(),
         for( set<VarSet>::const_iterator c = scl.begin(); c != scl.end(); c++ )
             cl.push_back(*c);
         if( props.verbose >= 3 ) {
-            cout << "HAK uses the following clusters: " << endl;
+            cout << Name << " uses the following clusters: " << endl;
             for( vector<VarSet>::const_iterator cli = cl.begin(); cli != cl.end(); cli++ )
                 cout << *cli << endl;
         }
@@ -162,7 +168,7 @@ HAK::HAK(const FactorGraph & fg, const PropertySet &opts) : DAIAlgRG(), props(),
     constructMessages();
 
     if( props.verbose >= 3 )
-        cout << "HAK regiongraph: " << *this << endl;
+        cout << Name << " regiongraph: " << *this << endl;
 }
 
 
@@ -225,14 +231,24 @@ double HAK::doGBP() {
     // Differences in single node beliefs
     Diffs diffs(nrVars(), 1.0);
 
-    size_t iter = 0;
     // do several passes over the network until maximum number of iterations has
     // been reached or until the maximum belief difference is smaller than tolerance
-    for( iter = 0; iter < props.maxiter && diffs.maxDiff() > props.tol; iter++ ) {
+    for( _iters = 0; _iters < props.maxiter && diffs.maxDiff() > props.tol; _iters++ ) {
         for( size_t beta = 0; beta < nrIRs(); beta++ ) {
             foreach( const Neighbor &alpha, nbIR(beta) ) {
                 size_t _beta = alpha.dual;
                 muab( alpha, _beta ) = _Qa[alpha].marginal(IR(beta)).divided_by( muba(alpha,_beta) );
+                /* TODO: INVESTIGATE THIS PROBLEM
+                 *
+                 * In some cases, the muab's can have very large entries because the muba's have very
+                 * small entries. This may cause NANs later on (e.g., multiplying large quantities may
+                 * result in +inf; normalization then tries to calculate inf / inf which is NAN). 
+                 * A fix of this problem would consist in normalizing the messages muab.
+                 * However, it is not obvious whether this is a real solution, because it has a
+                 * negative performance impact and the NAN's seem to be a symptom of a fundamental
+                 * numerical unstability.
+                 */
+                 muab(alpha,_beta).normalize(); 
             }
 
             Factor Qb_new;
@@ -243,28 +259,49 @@ double HAK::doGBP() {
 
             Qb_new.normalize();
             if( Qb_new.hasNaNs() ) {
-                cout << "HAK::doGBP:  Qb_new has NaNs!" << endl;
+                // TODO: WHAT TO DO IN THIS CASE?
+                cout << Name << "::doGBP:  Qb_new has NaNs!" << endl;
                 return 1.0;
             }
-//          _Qb[beta] = Qb_new.makeZero(1e-100);    // damping?
-            _Qb[beta] = Qb_new;
+            /* TODO: WHAT IS THE PURPOSE OF THE FOLLOWING CODE?
+             *
+             *   _Qb[beta] = Qb_new.makeZero(1e-100);
+             */
+
+            if( props.doubleloop || props.damping == 0.0 )
+                _Qb[beta] = Qb_new; // no damping for double loop
+            else
+                _Qb[beta] = (Qb_new^(1.0 - props.damping)) * (_Qb[beta]^props.damping);
 
             foreach( const Neighbor &alpha, nbIR(beta) ) {
                 size_t _beta = alpha.dual;
-
                 muba(alpha,_beta) = _Qb[beta].divided_by( muab(alpha,_beta) );
 
+                /* TODO: INVESTIGATE WHETHER THIS HACK (INVENTED BY KEES) TO PREVENT NANS MAKES SENSE 
+                 *
+                 *   muba(beta,*alpha).makePositive(1e-100);
+                 *
+                 */
+
                 Factor Qa_new = OR(alpha);
                 foreach( const Neighbor &gamma, nbOR(alpha) )
                     Qa_new *= muba(alpha,gamma.iter);
                 Qa_new ^= (1.0 / OR(alpha).c());
                 Qa_new.normalize();
                 if( Qa_new.hasNaNs() ) {
-                    cout << "HAK::doGBP:  Qa_new has NaNs!" << endl;
+                    cout << Name << "::doGBP:  Qa_new has NaNs!" << endl;
                     return 1.0;
                 }
-//              _Qa[alpha] = Qa_new.makeZero(1e-100); // damping?
-                _Qa[alpha] = Qa_new;
+                /* TODO: WHAT IS THE PURPOSE OF THE FOLLOWING CODE?
+                 *
+                 *   _Qb[beta] = Qb_new.makeZero(1e-100);
+                 */
+
+                if( props.doubleloop || props.damping == 0.0 )
+                    _Qa[alpha] = Qa_new; // no damping for double loop
+                else
+                    // FIXME: GEOMETRIC DAMPING IS SLOW!
+                _Qa[alpha] = (Qa_new^(1.0 - props.damping)) * (_Qa[alpha]^props.damping);
             }
         }
 
@@ -276,21 +313,21 @@ double HAK::doGBP() {
         }
 
         if( props.verbose >= 3 )
-            cout << "HAK::doGBP:  maxdiff " << diffs.maxDiff() << " after " << iter+1 << " passes" << endl;
+            cout << Name << "::doGBP:  maxdiff " << diffs.maxDiff() << " after " << _iters+1 << " passes" << endl;
     }
 
-    if( diffs.maxDiff() > maxdiff )
-        maxdiff = diffs.maxDiff();
+    if( diffs.maxDiff() > _maxdiff )
+        _maxdiff = diffs.maxDiff();
 
     if( props.verbose >= 1 ) {
         if( diffs.maxDiff() > props.tol ) {
             if( props.verbose == 1 )
                 cout << endl;
-            cout << "HAK::doGBP:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.maxDiff() << endl;
+            cout << Name << "::doGBP:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << diffs.maxDiff() << endl;
         } else {
             if( props.verbose >= 2 )
-                cout << "HAK::doGBP:  ";
-            cout << "converged in " << iter << " passes (" << toc() - tic << " clocks)." << endl;
+                cout << Name << "::doGBP:  ";
+            cout << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
         }
     }
 
@@ -326,16 +363,17 @@ double HAK::doDoubleLoop() {
     // Differences in single node beliefs
     Diffs diffs(nrVars(), 1.0);
 
-    size_t  outer_maxiter   = props.maxiter;
-    double  outer_tol       = props.tol;
-    size_t  outer_verbose   = props.verbose;
-    double  org_maxdiff     = maxdiff;
+    size_t outer_maxiter   = props.maxiter;
+    double outer_tol       = props.tol;
+    size_t outer_verbose   = props.verbose;
+    double org_maxdiff     = _maxdiff;
 
     // Set parameters for inner loop
     props.maxiter = 5;
     props.verbose = outer_verbose ? outer_verbose - 1 : 0;
 
     size_t outer_iter = 0;
+    size_t total_iter = 0;
     for( outer_iter = 0; outer_iter < outer_maxiter && diffs.maxDiff() > outer_tol; outer_iter++ ) {
         // Calculate new outer regions
         for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
@@ -355,17 +393,20 @@ double HAK::doDoubleLoop() {
             old_beliefs[i] = new_belief;
         }
 
+        total_iter += Iterations();
+
         if( props.verbose >= 3 )
-            cout << "HAK::doDoubleLoop:  maxdiff " << diffs.maxDiff() << " after " << outer_iter+1 << " passes" << endl;
+            cout << Name << "::doDoubleLoop:  maxdiff " << diffs.maxDiff() << " after " << total_iter << " passes" << endl;
     }
 
     // restore _maxiter, _verbose and _maxdiff
     props.maxiter = outer_maxiter;
     props.verbose = outer_verbose;
-    maxdiff = org_maxdiff;
+    _maxdiff = org_maxdiff;
 
-    if( diffs.maxDiff() > maxdiff )
-        maxdiff = diffs.maxDiff();
+    _iters = total_iter;
+    if( diffs.maxDiff() > _maxdiff )
+        _maxdiff = diffs.maxDiff();
 
     // Restore original outer regions
     ORs = org_ORs;
@@ -378,11 +419,11 @@ double HAK::doDoubleLoop() {
         if( diffs.maxDiff() > props.tol ) {
             if( props.verbose == 1 )
                 cout << endl;
-                cout << "HAK::doDoubleLoop:  WARNING: not converged within " << outer_maxiter << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.maxDiff() << endl;
+                cout << Name << "::doDoubleLoop:  WARNING: not converged within " << outer_maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << diffs.maxDiff() << endl;
             } else {
                 if( props.verbose >= 3 )
-                    cout << "HAK::doDoubleLoop:  ";
-                cout << "converged in " << outer_iter << " passes (" << toc() - tic << " clocks)." << endl;
+                    cout << Name << "::doDoubleLoop:  ";
+                cout << "converged in " << total_iter << " passes (" << toc() - tic << " seconds)." << endl;
             }
         }
 
index 24a5eed..0f28d5b 100644 (file)
@@ -330,22 +330,14 @@ size_t JTree::findEfficientTree( const VarSet& ns, DEdgeVec &Tree, size_t Previo
         }
     }
 
-//  for( size_t e = 0; e < _RTree.size(); e++ )
-//      cout << OR(_RTree[e].n1).vars() << "->" << OR(_RTree[e].n2).vars() << ",  ";
-//  cout << endl;
     // grow new tree
     Graph oldTree;
     for( DEdgeVec::const_iterator e = _RTree.begin(); e != _RTree.end(); e++ )
         oldTree.insert( UEdge(e->n1, e->n2) );
     DEdgeVec newTree = GrowRootedTree( oldTree, maxalpha );
-//  cout << ns << ": ";
-//  for( size_t e = 0; e < newTree.size(); e++ )
-//      cout << OR(newTree[e].n1).vars() << "->" << OR(newTree[e].n2).vars() << ",  ";
-//  cout << endl;
     
     // identify subtree that contains variables of ns which are not in the new root
     VarSet nsrem = ns / OR(maxalpha).vars();
-//  cout << "nsrem:" << nsrem << endl;
     set<DEdge> subTree;
     // for each variable in ns that is not in the root clique
     for( VarSet::const_iterator n = nsrem.begin(); n != nsrem.end(); n++ ) {
@@ -384,10 +376,6 @@ size_t JTree::findEfficientTree( const VarSet& ns, DEdgeVec &Tree, size_t Previo
                 pos = newTree[e-1].n1;
             }
     }
-//  cout << "subTree: " << endl;
-//  for( set<DEdge>::const_iterator sTi = subTree.begin(); sTi != subTree.end(); sTi++ )
-//      cout << OR(sTi->n1).vars() << "->" << OR(sTi->n2).vars() << ",  ";
-//  cout << endl;
 
     // Resulting Tree is a reordered copy of newTree
     // First add edges in subTree to Tree
@@ -395,9 +383,7 @@ size_t JTree::findEfficientTree( const VarSet& ns, DEdgeVec &Tree, size_t Previo
     for( DEdgeVec::const_iterator e = newTree.begin(); e != newTree.end(); e++ )
         if( subTree.count( *e ) ) {
             Tree.push_back( *e );
-//          cout << OR(e->n1).vars() << "->" << OR(e->n2).vars() << ",  ";
         }
-//  cout << endl;
     // Then add edges pointing away from nsrem
     // FIXME
 /*  for( DEdgeVec::const_iterator e = newTree.begin(); e != newTree.end(); e++ )
@@ -406,7 +392,6 @@ size_t JTree::findEfficientTree( const VarSet& ns, DEdgeVec &Tree, size_t Previo
                 if( e->n1 == sTi->n1 || e->n1 == sTi->n2 ||
                     e->n2 == sTi->n1 || e->n2 == sTi->n2 ) {
                     Tree.push_back( *e );
-//                  cout << OR(e->n1).vars() << "->" << OR(e->n2).vars() << ",  ";
                 }
             }*/
     // FIXME
@@ -420,10 +405,8 @@ size_t JTree::findEfficientTree( const VarSet& ns, DEdgeVec &Tree, size_t Previo
                 }
             if( found ) {
                 Tree.push_back( *e );
-                cout << OR(e->n1).vars() << "->" << OR(e->n2).vars() << ",  ";
             }
-        }
-    cout << endl;*/
+        }*/
     size_t subTreeSize = Tree.size();
     // Then add remaining edges
     for( DEdgeVec::const_iterator e = newTree.begin(); e != newTree.end(); e++ )
index 0aec5d9..2c796fe 100644 (file)
@@ -27,7 +27,6 @@
 #include <dai/diffs.h>
 #include <dai/util.h>
 #include <dai/alldai.h>
-#include <dai/x2x.h>
 
 
 namespace dai {
@@ -57,6 +56,10 @@ void LC::setProperties( const PropertySet &opts ) {
         props.cavaiopts = opts.getStringAs<PropertySet>("cavaiopts");
     if( opts.hasKey("reinit") )
         props.reinit = opts.getStringAs<bool>("reinit");
+    if( opts.hasKey("damping") )
+        props.damping = opts.getStringAs<double>("damping");
+    else
+        props.damping = 0.0;
 }
 
 
@@ -70,6 +73,7 @@ PropertySet LC::getProperties() const {
     opts.Set( "cavainame", props.cavainame );
     opts.Set( "cavaiopts", props.cavaiopts );
     opts.Set( "reinit", props.reinit );
+    opts.Set( "damping", props.damping );
     return opts;
 }
 
@@ -84,20 +88,21 @@ string LC::printProperties() const {
     s << "updates=" << props.updates << ",";
     s << "cavainame=" << props.cavainame << ",";
     s << "cavaiopts=" << props.cavaiopts << ",";
-    s << "reinit=" << props.reinit << "]";
+    s << "reinit=" << props.reinit << ",";
+    s << "damping=" << props.damping << "]";
     return s.str();
 }
 
 
-LC::LC( const FactorGraph & fg, const PropertySet &opts ) : DAIAlgFG(fg), _pancakes(), _cavitydists(), _phis(), _beliefs(), props(), maxdiff(0.0) {
+LC::LC( const FactorGraph & fg, const PropertySet &opts ) : DAIAlgFG(fg), _pancakes(), _cavitydists(), _phis(), _beliefs(), _maxdiff(0.0), _iters(0), props() {
     setProperties( opts );
 
     // create pancakes
-    _pancakes.resize(nrVars());
+    _pancakes.resize( nrVars() );
    
     // create cavitydists
     for( size_t i=0; i < nrVars(); i++ )
-        _cavitydists.push_back(Factor(delta(i)));
+        _cavitydists.push_back(Factor( delta(i) ));
 
     // create phis
     _phis.reserve( nrVars() );
@@ -109,6 +114,7 @@ LC::LC( const FactorGraph & fg, const PropertySet &opts ) : DAIAlgFG(fg), _panca
     }
 
     // create beliefs
+    _beliefs.reserve( nrVars() );
     for( size_t i=0; i < nrVars(); i++ )
         _beliefs.push_back(Factor(var(i)));
 }
@@ -160,7 +166,7 @@ double LC::InitCavityDists( const std::string &name, const PropertySet &opts ) {
     double tic = toc();
 
     if( props.verbose >= 1 ) {
-        cout << "LC::InitCavityDists:  ";
+        cout << Name << "::InitCavityDists:  ";
         if( props.cavity == Properties::CavityType::UNIFORM )
             cout << "Using uniform initial cavity distributions" << endl;
         else if( props.cavity == Properties::CavityType::FULL )
@@ -180,7 +186,7 @@ double LC::InitCavityDists( const std::string &name, const PropertySet &opts ) {
     init();
 
     if( props.verbose >= 1 ) {
-        cout << "used " << toc() - tic << " clocks." << endl;
+        cout << Name << "::InitCavityDists used " << toc() - tic << " seconds." << endl;
     }
 
     return maxdiff;
@@ -189,7 +195,7 @@ double LC::InitCavityDists( const std::string &name, const PropertySet &opts ) {
 
 long LC::SetCavityDists( std::vector<Factor> &Q ) {
     if( props.verbose >= 1 ) 
-        cout << "LC::SetCavityDists:  Setting initial cavity distributions" << endl;
+        cout << Name << "::SetCavityDists:  Setting initial cavity distributions" << endl;
     if( Q.size() != nrVars() )
         return -1;
     for( size_t i = 0; i < nrVars(); i++ ) {
@@ -240,6 +246,8 @@ Factor LC::NewPancake (size_t i, size_t _I, bool & hasNaNs) {
         A_I ^= (1.0 / (Ivars.size() - 1));
     Factor A_Ii = (_pancakes[i] * factor(I).inverse() * _phis[i][_I].inverse()).partSum( Ivars / var(i) );
     Factor quot = A_I.divided_by(A_Ii);
+    if( props.damping != 0.0 )
+        quot = (quot^(1.0 - props.damping)) * (_phis[i][_I]^props.damping);
 
     piet *= quot.divided_by( _phis[i][_I] ).normalized();
     _phis[i][_I] = quot.normalized();
@@ -247,7 +255,7 @@ Factor LC::NewPancake (size_t i, size_t _I, bool & hasNaNs) {
     piet.normalize();
 
     if( piet.hasNaNs() ) {
-        cout << "LC::NewPancake(" << i << ", " << _I << "):  has NaNs!" << endl;
+        cout << Name << "::NewPancake(" << i << ", " << _I << "):  has NaNs!" << endl;
         hasNaNs = true;
     }
 
@@ -265,8 +273,8 @@ double LC::run() {
     Diffs diffs(nrVars(), 1.0);
 
     double md = InitCavityDists( props.cavainame, props.cavaiopts );
-    if( md > maxdiff )
-        maxdiff = md;
+    if( md > _maxdiff )
+        _maxdiff = md;
 
     vector<Factor> old_beliefs;
     for(size_t i=0; i < nrVars(); i++ )
@@ -279,8 +287,8 @@ double LC::run() {
             break;
         }
     if( hasNaNs ) {
-        cout << "LC::run:  initial _pancakes has NaNs!" << endl;
-        return -1.0;
+        cout << Name << "::run:  initial _pancakes has NaNs!" << endl;
+        return 1.0;
     }
 
     size_t nredges = nrEdges();
@@ -290,11 +298,9 @@ double LC::run() {
         foreach( const Neighbor &I, nbV(i) )
             update_seq.push_back( Edge( i, I.iter ) );
 
-    size_t iter = 0;
-
     // do several passes over the network until maximum number of iterations has
     // been reached or until the maximum belief difference is smaller than tolerance
-    for( iter=0; iter < props.maxiter && diffs.maxDiff() > props.tol; iter++ ) {
+    for( _iters=0; _iters < props.maxiter && diffs.maxDiff() > props.tol; _iters++ ) {
         // Sequential updates
         if( props.updates == Properties::UpdateType::SEQRND )
             random_shuffle( update_seq.begin(), update_seq.end() );
@@ -304,7 +310,7 @@ double LC::run() {
             size_t _I = update_seq[t].second;
             _pancakes[i] = NewPancake( i, _I, hasNaNs);
             if( hasNaNs )
-                return -1.0;
+                return 1.0;
             CalcBelief( i );
         }
 
@@ -315,21 +321,21 @@ double LC::run() {
         }
 
         if( props.verbose >= 3 )
-            cout << "LC::run:  maxdiff " << diffs.maxDiff() << " after " << iter+1 << " passes" << endl;
+            cout << Name << "::run:  maxdiff " << diffs.maxDiff() << " after " << _iters+1 << " passes" << endl;
     }
 
-    if( diffs.maxDiff() > maxdiff )
-        maxdiff = diffs.maxDiff();
+    if( diffs.maxDiff() > _maxdiff )
+        _maxdiff = diffs.maxDiff();
 
     if( props.verbose >= 1 ) {
         if( diffs.maxDiff() > props.tol ) {
             if( props.verbose == 1 )
                 cout << endl;
-                cout << "LC::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.maxDiff() << endl;
+                cout << Name << "::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << diffs.maxDiff() << endl;
         } else {
             if( props.verbose >= 2 )
-                cout << "LC::run:  ";
-                cout << "converged in " << iter << " passes (" << toc() - tic << " clocks)." << endl;
+                cout << Name << "::run:  ";
+                cout << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
         }
     }
 
index aab55ab..5a45b9f 100644 (file)
@@ -126,7 +126,7 @@ double MF::run() {
         jan.normalize();
 
         if( jan.hasNaNs() ) {
-            cout << "MF::run():  ERROR: jan has NaNs!" << endl;
+            cout << Name << "::run():  ERROR: jan has NaNs!" << endl;
             return 1.0;
         }
 
@@ -145,11 +145,11 @@ double MF::run() {
         if( diffs.maxDiff() > props.tol ) {
             if( props.verbose == 1 )
                 cout << endl;
-            cout << "MF::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.maxDiff() << endl;
+            cout << Name << "::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << diffs.maxDiff() << endl;
         } else {
             if( props.verbose >= 2 )
-                cout << "MF::run:  ";
-            cout << "converged in " << t / pass_size << " passes (" << toc() - tic << " clocks)." << endl;
+                cout << Name << "::run:  ";
+            cout << "converged in " << t / pass_size << " passes (" << toc() - tic << " seconds)." << endl;
         }
     }
 
index 13f3734..76b4ef8 100644 (file)
@@ -466,8 +466,9 @@ void MR::solvemcav() {
         }
     } while((maxdev>props.tol)&&(run<maxruns));
 
-    if( maxdev > maxdiff )
-        maxdiff = maxdev;
+    _iters = run;
+    if( maxdev > _maxdiff )
+        _maxdiff = maxdev;
 
     if(run==maxruns){
         if( props.verbose >= 1 )
@@ -507,7 +508,7 @@ void MR::init_cor() {
     for( size_t i = 0; i < nrVars(); i++ ) {
         vector<Factor> pairq;
         if( props.inits == Properties::InitType::CLAMPING ) {
-            BP bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", string("1e-9"))("maxiter", string("1000UL"))("verbose", string("0UL"))("logdomain", string("0")));
+            BP bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", string("1e-9"))("maxiter", string("10000UL"))("verbose", string("0UL"))("logdomain", string("0")));
             bpcav.makeCavity( i );
             pairq = calcPairBeliefs( bpcav, delta(i), false );
         } else if( props.inits == Properties::InitType::EXACT ) {
@@ -561,8 +562,8 @@ double MR::run() {
 
         if( props.inits == Properties::InitType::RESPPROP ) {
             double md = init_cor_resp();
-            if( md > maxdiff )
-                maxdiff = md;
+            if( md > _maxdiff )
+                _maxdiff = md;
         } else if( props.inits == Properties::InitType::EXACT )
             init_cor(); // FIXME no MaxDiff() calculation
         else if( props.inits == Properties::InitType::CLAMPING )
@@ -574,11 +575,11 @@ double MR::run() {
         solveM();
 
         if( props.verbose >= 1 )
-            cout << "MR needed " << toc() - tic << " clocks." << endl;
+            cout << Name << " needed " << toc() - tic << " seconds." << endl;
 
         return 0.0;
     } else
-        return -1.0;
+        return 1.0;
 }
 
 
@@ -617,7 +618,7 @@ vector<Factor> MR::beliefs() const {
 
 
 
-MR::MR( const FactorGraph &fg, const PropertySet &opts ) : DAIAlgFG(fg), supported(true), maxdiff(0.0) {
+MR::MR( const FactorGraph &fg, const PropertySet &opts ) : DAIAlgFG(fg), supported(true), _maxdiff(0.0), _iters(0) {
     setProperties( opts );
 
     // check whether all vars in fg are binary
index b780d26..ce7d167 100644 (file)
@@ -71,7 +71,7 @@ string TreeEP::printProperties() const {
 }
 
 
-TreeEPSubTree::TreeEPSubTree( const DEdgeVec &subRTree, const DEdgeVec &jt_RTree, const std::vector<Factor> &jt_Qa, const std::vector<Factor> &jt_Qb, const Factor *I ) : _Qa(), _Qb(), _RTree(), _a(), _b(), _I(I), _ns(), _nsrem(), _logZ(0.0) {
+TreeEP::TreeEPSubTree::TreeEPSubTree( const DEdgeVec &subRTree, const DEdgeVec &jt_RTree, const std::vector<Factor> &jt_Qa, const std::vector<Factor> &jt_Qb, const Factor *I ) : _Qa(), _Qb(), _RTree(), _a(), _b(), _I(I), _ns(), _nsrem(), _logZ(0.0) {
     _ns = _I->vars();
 
     // Make _Qa, _Qb, _a and _b corresponding to the subtree
@@ -106,10 +106,10 @@ TreeEPSubTree::TreeEPSubTree( const DEdgeVec &subRTree, const DEdgeVec &jt_RTree
 
     // Find remaining variables (which are not in the new root)
     _nsrem = _ns / _Qa[0].vars();
-};
+}
 
 
-void TreeEPSubTree::init() { 
+void TreeEP::TreeEPSubTree::init() { 
     for( size_t alpha = 0; alpha < _Qa.size(); alpha++ )
         _Qa[alpha].fill( 1.0 );
     for( size_t beta = 0; beta < _Qb.size(); beta++ )
@@ -117,7 +117,7 @@ void TreeEPSubTree::init() {
 }
 
 
-void TreeEPSubTree::InvertAndMultiply( const std::vector<Factor> &Qa, const std::vector<Factor> &Qb ) {
+void TreeEP::TreeEPSubTree::InvertAndMultiply( const std::vector<Factor> &Qa, const std::vector<Factor> &Qb ) {
     for( size_t alpha = 0; alpha < _Qa.size(); alpha++ )
         _Qa[alpha] = Qa[_a[alpha]].divided_by( _Qa[alpha] );
 
@@ -126,7 +126,7 @@ void TreeEPSubTree::InvertAndMultiply( const std::vector<Factor> &Qa, const std:
 }
 
 
-void TreeEPSubTree::HUGIN_with_I( std::vector<Factor> &Qa, std::vector<Factor> &Qb ) {
+void TreeEP::TreeEPSubTree::HUGIN_with_I( std::vector<Factor> &Qa, std::vector<Factor> &Qb ) {
     // Backup _Qa and _Qb
     vector<Factor> _Qa_old(_Qa);
     vector<Factor> _Qb_old(_Qb);
@@ -187,7 +187,7 @@ void TreeEPSubTree::HUGIN_with_I( std::vector<Factor> &Qa, std::vector<Factor> &
 }
 
 
-double TreeEPSubTree::logZ( const std::vector<Factor> &Qa, const std::vector<Factor> &Qb ) const {
+double TreeEP::TreeEPSubTree::logZ( const std::vector<Factor> &Qa, const std::vector<Factor> &Qb ) const {
     double sum = 0.0;
     for( size_t alpha = 0; alpha < _Qa.size(); alpha++ )
         sum += (Qa[_a[alpha]] * _Qa[alpha].log0()).totalSum();
@@ -197,7 +197,7 @@ double TreeEPSubTree::logZ( const std::vector<Factor> &Qa, const std::vector<Fac
 }
 
 
-TreeEP::TreeEP( const FactorGraph &fg, const PropertySet &opts ) : JTree(fg, opts("updates",string("HUGIN")), false), props(), maxdiff(0.0) {
+TreeEP::TreeEP( const FactorGraph &fg, const PropertySet &opts ) : JTree(fg, opts("updates",string("HUGIN")), false), _maxdiff(0.0), _iters(0), props(), _Q() {
     setProperties( opts );
 
     assert( fg.isConnected() );
@@ -205,65 +205,49 @@ TreeEP::TreeEP( const FactorGraph &fg, const PropertySet &opts ) : JTree(fg, opt
     if( opts.hasKey("tree") ) {
         ConstructRG( opts.GetAs<DEdgeVec>("tree") );
     } else {
-        if( props.type == Properties::TypeType::ORG ) {
-            // construct weighted graph with as weights a crude estimate of the
+        if( props.type == Properties::TypeType::ORG || props.type == Properties::TypeType::ALT ) {
+            // ORG: construct weighted graph with as weights a crude estimate of the
             // mutual information between the nodes
-            WeightedGraph<double> wg;
-            for( size_t i = 0; i < nrVars(); ++i ) {
-                Var v_i = var(i);
-                VarSet di = delta(i);
-                for( VarSet::const_iterator j = di.begin(); j != di.end(); j++ )
-                    if( v_i < *j ) {
-                        Factor piet;
-                        for( size_t I = 0; I < nrFactors(); I++ ) {
-                            VarSet Ivars = factor(I).vars();
-                            if( (Ivars == v_i) || (Ivars == *j) )
-                                piet *= factor(I);
-                            else if( Ivars >> (v_i | *j) )
-                                piet *= factor(I).marginal( v_i | *j );
-                        }
-                        if( piet.vars() >> (v_i | *j) ) {
-                            piet = piet.marginal( v_i | *j );
-                            Factor pietf = piet.marginal(v_i) * piet.marginal(*j);
-                            wg[UEdge(i,findVar(*j))] = KL_dist( piet, pietf );
-                        } else
-                            wg[UEdge(i,findVar(*j))] = 0;
-                    }
-            }
-
-            // find maximal spanning tree
-            ConstructRG( MaxSpanningTreePrims( wg ) );
-
-//            cout << "Constructing maximum spanning tree..." << endl;
-//            DEdgeVec MST = MaxSpanningTreePrims( wg );
-//            cout << "Maximum spanning tree:" << endl;
-//            for( DEdgeVec::const_iterator e = MST.begin(); e != MST.end(); e++ )
-//                cout << *e << endl; 
-//            ConstructRG( MST );
-        } else if( props.type == Properties::TypeType::ALT ) {
-            // construct weighted graph with as weights an upper bound on the
+            // ALT: construct weighted graph with as weights an upper bound on the
             // effective interaction strength between pairs of nodes
+            
             WeightedGraph<double> wg;
             for( size_t i = 0; i < nrVars(); ++i ) {
                 Var v_i = var(i);
                 VarSet di = delta(i);
                 for( VarSet::const_iterator j = di.begin(); j != di.end(); j++ )
                     if( v_i < *j ) {
+                        VarSet ij(v_i,*j);
                         Factor piet;
                         for( size_t I = 0; I < nrFactors(); I++ ) {
                             VarSet Ivars = factor(I).vars();
-                            if( Ivars >> (v_i | *j) )
-                                piet *= factor(I);
+                            if( props.type == Properties::TypeType::ORG ) {
+                                if( (Ivars == v_i) || (Ivars == *j) )
+                                    piet *= factor(I);
+                                else if( Ivars >> ij )
+                                    piet *= factor(I).marginal( ij );
+                            } else {
+                                if( Ivars >> ij )
+                                    piet *= factor(I);
+                            }
+                        }
+                        if( props.type == Properties::TypeType::ORG ) {
+                            if( piet.vars() >> ij ) {
+                                piet = piet.marginal( ij );
+                                Factor pietf = piet.marginal(v_i) * piet.marginal(*j);
+                                wg[UEdge(i,findVar(*j))] = KL_dist( piet, pietf );
+                            } else
+                                wg[UEdge(i,findVar(*j))] = 0;
+                        } else {
+                            wg[UEdge(i,findVar(*j))] = piet.strength(v_i, *j);
                         }
-                        wg[UEdge(i,findVar(*j))] = piet.strength(v_i, *j);
                     }
             }
 
             // find maximal spanning tree
             ConstructRG( MaxSpanningTreePrims( wg ) );
-        } else {
+        } else
             DAI_THROW(INTERNAL_ERROR);
-        }
     }
 }
 
@@ -280,7 +264,8 @@ void TreeEP::ConstructRG( const DEdgeVec &tree ) {
     for( size_t i = 0; i < Cliques.size(); i++ )
         for( size_t j = i+1; j < Cliques.size(); j++ ) {
             size_t w = (Cliques[i] & Cliques[j]).size();
-            JuncGraph[UEdge(i,j)] = w;
+            if( w )
+                JuncGraph[UEdge(i,j)] = w;
         }
     
     // Construct maximal spanning tree using Prim's algorithm
@@ -353,26 +338,6 @@ void TreeEP::ConstructRG( const DEdgeVec &tree ) {
             //subTree.resize( subTreeSize );  // FIXME
 //          cout << "subtree " << I << " has size " << subTreeSize << endl;
 
-/*
-            char fn[30];
-            sprintf( fn, "/tmp/subtree_%d.dot", I );
-            std::ofstream dots(fn);
-            dots << "graph G {" << endl;
-            dots << "graph[size=\"9,9\"];" << endl;
-            dots << "node[shape=circle,width=0.4,fixedsize=true];" << endl;
-            for( size_t i = 0; i < nrVars(); i++ )
-                dots << "\tx" << var(i).label() << ((factor(I).vars() >> var(i)) ? "[color=blue];" : ";") << endl;
-            dots << "node[shape=box,style=filled,color=lightgrey,width=0.3,height=0.3,fixedsize=true];" << endl;
-            for( size_t J = 0; J < nrFactors(); J++ )
-                dots << "\tp" << J << ";" << endl;
-            for( size_t iI = 0; iI < FactorGraph::nr_edges(); iI++ )
-                dots << "\tx" << var(FactorGraph::edge(iI).first).label() << " -- p" << FactorGraph::edge(iI).second << ";" << endl;
-            for( size_t a = 0; a < tree.size(); a++ )
-                dots << "\tx" << var(tree[a].n1).label() << " -- x" << var(tree[a].n2).label() << " [color=red];" << endl;
-            dots << "}" << endl;
-            dots.close();
-*/
-            
             TreeEPSubTree QI( subTree, _RTree, _Qa, _Qb, &factor(I) );
             _Q[I] = QI;
         }
@@ -425,11 +390,9 @@ double TreeEP::run() {
     for( size_t i = 0; i < nrVars(); i++ )
         old_beliefs.push_back(belief(var(i)));
 
-    size_t iter = 0;
-    
     // do several passes over the network until maximum number of iterations has
     // been reached or until the maximum belief difference is smaller than tolerance
-    for( iter=0; iter < props.maxiter && diffs.maxDiff() > props.tol; iter++ ) {
+    for( _iters=0; _iters < props.maxiter && diffs.maxDiff() > props.tol; _iters++ ) {
         for( size_t I = 0; I < nrFactors(); I++ )
             if( offtree(I) ) {  
                 _Q[I].InvertAndMultiply( _Qa, _Qb );
@@ -445,21 +408,21 @@ double TreeEP::run() {
         }
 
         if( props.verbose >= 3 )
-            cout << "TreeEP::run:  maxdiff " << diffs.maxDiff() << " after " << iter+1 << " passes" << endl;
+            cout << Name << "::run:  maxdiff " << diffs.maxDiff() << " after " << _iters+1 << " passes" << endl;
     }
 
-    if( diffs.maxDiff() > maxdiff )
-        maxdiff = diffs.maxDiff();
+    if( diffs.maxDiff() > _maxdiff )
+        _maxdiff = diffs.maxDiff();
 
     if( props.verbose >= 1 ) {
         if( diffs.maxDiff() > props.tol ) {
             if( props.verbose == 1 )
                 cout << endl;
-            cout << "TreeEP::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.maxDiff() << endl;
+            cout << Name << "::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << diffs.maxDiff() << endl;
         } else {
             if( props.verbose >= 3 )
-                cout << "TreeEP::run:  ";
-            cout << "converged in " << iter << " passes (" << toc() - tic << " clocks)." << endl;
+                cout << Name << "::run:  ";
+            cout << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
         }
     }
 
index 7d06692..5f24776 100644 (file)
@@ -107,8 +107,12 @@ class TestDAI {
                 } catch( Exception &e ) {
                     has_maxdiff = false;
                 }
-                has_iters = false;
-                iters = 0;
+                try {
+                    iters = obj->Iterations();
+                    has_iters = true;
+                } catch( Exception &e ) {
+                    has_iters = false;
+                }
                 q = allBeliefs();
             };
         }
@@ -340,7 +344,7 @@ int main( int argc, char *argv[] ) {
                 }
                 cout.width( 10 );
                 if( piet.has_iters ) {
-                    cout << piet.iters << "  " << endl;
+                    cout << piet.iters << "  ";
                 } else {
                     cout << "N/A         ";
                 }
index 9587ee1..0f02e9f 100644 (file)
@@ -17,7 +17,7 @@ EXACT[verbose=0]
 # ([13] <0.534838 0.465162 >)
 # ([14] <0.62914 0.37086 >)
 # ([15] <0.135686 0.864314 >)
-JTREE_HUGIN                               1.000e-09   1.000e-09   1.000e-09   1.000e-09   N/A         
+JTREE_HUGIN                               1.000e-09   1.000e-09   1.000e-09   1.000e-09   1           
 # ([0] <3.888e-01 6.112e-01 >)
 # ([1] <5.556e-01 4.444e-01 >)
 # ([2] <4.587e-01 5.413e-01 >)
@@ -34,7 +34,7 @@ JTREE_HUGIN                               1.000e-09   1.000e-09   1.000e-09   1.
 # ([13] <5.348e-01 4.652e-01 >)
 # ([14] <6.291e-01 3.709e-01 >)
 # ([15] <1.357e-01 8.643e-01 >)
-JTREE_SHSH                                1.000e-09   1.000e-09   1.000e-09   1.000e-09   N/A         
+JTREE_SHSH                                1.000e-09   1.000e-09   1.000e-09   1.000e-09   1           
 # ([0] <3.888e-01 6.112e-01 >)
 # ([1] <5.556e-01 4.444e-01 >)
 # ([2] <4.587e-01 5.413e-01 >)
@@ -51,7 +51,7 @@ JTREE_SHSH                                1.000e-09   1.000e-09   1.000e-09   1.
 # ([13] <5.348e-01 4.652e-01 >)
 # ([14] <6.291e-01 3.709e-01 >)
 # ([15] <1.357e-01 8.643e-01 >)
-BP_SEQFIX                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09   N/A         
+BP_SEQFIX                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09   45          
 # ([0] <4.233e-01 5.767e-01 >)
 # ([1] <5.422e-01 4.578e-01 >)
 # ([2] <4.662e-01 5.338e-01 >)
@@ -68,7 +68,7 @@ BP_SEQFIX                                 9.483e-02   3.078e-02   1.737e-02   1.
 # ([13] <5.266e-01 4.734e-01 >)
 # ([14] <6.033e-01 3.967e-01 >)
 # ([15] <1.558e-01 8.442e-01 >)
-BP_SEQRND                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09   N/A         
+BP_SEQRND                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09   43          
 # ([0] <4.233e-01 5.767e-01 >)
 # ([1] <5.422e-01 4.578e-01 >)
 # ([2] <4.662e-01 5.338e-01 >)
@@ -85,7 +85,7 @@ BP_SEQRND                                 9.483e-02   3.078e-02   1.737e-02   1.
 # ([13] <5.266e-01 4.734e-01 >)
 # ([14] <6.033e-01 3.967e-01 >)
 # ([15] <1.558e-01 8.442e-01 >)
-BP_SEQMAX                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09   N/A         
+BP_SEQMAX                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09   13          
 # ([0] <4.233e-01 5.767e-01 >)
 # ([1] <5.422e-01 4.578e-01 >)
 # ([2] <4.662e-01 5.338e-01 >)
@@ -102,7 +102,7 @@ BP_SEQMAX                                 9.483e-02   3.078e-02   1.737e-02   1.
 # ([13] <5.266e-01 4.734e-01 >)
 # ([14] <6.033e-01 3.967e-01 >)
 # ([15] <1.558e-01 8.442e-01 >)
-BP_PARALL                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09   N/A         
+BP_PARALL                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09   71          
 # ([0] <4.233e-01 5.767e-01 >)
 # ([1] <5.422e-01 4.578e-01 >)
 # ([2] <4.662e-01 5.338e-01 >)
@@ -119,7 +119,7 @@ BP_PARALL                                 9.483e-02   3.078e-02   1.737e-02   1.
 # ([13] <5.266e-01 4.734e-01 >)
 # ([14] <6.033e-01 3.967e-01 >)
 # ([15] <1.558e-01 8.442e-01 >)
-BP_SEQFIX_LOG                             9.483e-02   3.078e-02   1.737e-02   1.000e-09   N/A         
+BP_SEQFIX_LOG                             9.483e-02   3.078e-02   1.737e-02   1.000e-09   45          
 # ([0] <4.233e-01 5.767e-01 >)
 # ([1] <5.422e-01 4.578e-01 >)
 # ([2] <4.662e-01 5.338e-01 >)
@@ -136,7 +136,7 @@ BP_SEQFIX_LOG                             9.483e-02   3.078e-02   1.737e-02   1.
 # ([13] <5.266e-01 4.734e-01 >)
 # ([14] <6.033e-01 3.967e-01 >)
 # ([15] <1.558e-01 8.442e-01 >)
-BP_SEQRND_LOG                             9.483e-02   3.078e-02   1.737e-02   1.000e-09   N/A         
+BP_SEQRND_LOG                             9.483e-02   3.078e-02   1.737e-02   1.000e-09   43          
 # ([0] <4.233e-01 5.767e-01 >)
 # ([1] <5.422e-01 4.578e-01 >)
 # ([2] <4.662e-01 5.338e-01 >)
@@ -153,7 +153,7 @@ BP_SEQRND_LOG                             9.483e-02   3.078e-02   1.737e-02   1.
 # ([13] <5.266e-01 4.734e-01 >)
 # ([14] <6.033e-01 3.967e-01 >)
 # ([15] <1.558e-01 8.442e-01 >)
-BP_SEQMAX_LOG                             9.483e-02   3.078e-02   1.737e-02   1.000e-09   N/A         
+BP_SEQMAX_LOG                             9.483e-02   3.078e-02   1.737e-02   1.000e-09   13          
 # ([0] <4.233e-01 5.767e-01 >)
 # ([1] <5.422e-01 4.578e-01 >)
 # ([2] <4.662e-01 5.338e-01 >)
@@ -170,7 +170,7 @@ BP_SEQMAX_LOG                             9.483e-02   3.078e-02   1.737e-02   1.
 # ([13] <5.266e-01 4.734e-01 >)
 # ([14] <6.033e-01 3.967e-01 >)
 # ([15] <1.558e-01 8.442e-01 >)
-BP_PARALL_LOG                             9.483e-02   3.078e-02   1.737e-02   1.000e-09   N/A         
+BP_PARALL_LOG                             9.483e-02   3.078e-02   1.737e-02   1.000e-09   71          
 # ([0] <4.233e-01 5.767e-01 >)
 # ([1] <5.422e-01 4.578e-01 >)
 # ([2] <4.662e-01 5.338e-01 >)
@@ -187,7 +187,7 @@ BP_PARALL_LOG                             9.483e-02   3.078e-02   1.737e-02   1.
 # ([13] <5.266e-01 4.734e-01 >)
 # ([14] <6.033e-01 3.967e-01 >)
 # ([15] <1.558e-01 8.442e-01 >)
-MF_SEQRND                                 3.607e-01   1.904e-01   -9.409e-02  1.000e-09   N/A         
+MF_SEQRND                                 3.607e-01   1.904e-01   -9.409e-02  1.000e-09   58          
 # ([0] <2.053e-01 7.947e-01 >)
 # ([1] <9.163e-01 8.373e-02 >)
 # ([2] <1.579e-01 8.421e-01 >)
@@ -204,7 +204,7 @@ MF_SEQRND                                 3.607e-01   1.904e-01   -9.409e-02  1.
 # ([13] <8.148e-01 1.852e-01 >)
 # ([14] <8.338e-01 1.662e-01 >)
 # ([15] <5.661e-03 9.943e-01 >)
-TREEEP                                    3.268e-02   8.023e-03   6.340e-04   1.000e-09   N/A         
+TREEEP                                    3.268e-02   8.023e-03   6.340e-04   1.000e-09   15          
 # ([0] <3.980e-01 6.020e-01 >)
 # ([1] <5.520e-01 4.480e-01 >)
 # ([2] <4.620e-01 5.380e-01 >)
@@ -221,7 +221,7 @@ TREEEP                                    3.268e-02   8.023e-03   6.340e-04   1.
 # ([13] <5.293e-01 4.707e-01 >)
 # ([14] <6.409e-01 3.591e-01 >)
 # ([15] <1.364e-01 8.636e-01 >)
-TREEEPWC                                  2.356e-02   1.026e-02   4.876e-03   1.000e-09   N/A         
+TREEEPWC                                  2.356e-02   1.026e-02   4.876e-03   1.000e-09   14          
 # ([0] <4.091e-01 5.909e-01 >)
 # ([1] <5.429e-01 4.571e-01 >)
 # ([2] <4.697e-01 5.303e-01 >)
@@ -238,7 +238,7 @@ TREEEPWC                                  2.356e-02   1.026e-02   4.876e-03   1.
 # ([13] <5.298e-01 4.702e-01 >)
 # ([14] <6.299e-01 3.701e-01 >)
 # ([15] <1.384e-01 8.616e-01 >)
-GBP_MIN                                   9.483e-02   3.078e-02   1.737e-02   1.000e-09   N/A         
+GBP_MIN                                   9.483e-02   3.078e-02   1.737e-02   1.000e-09   44          
 # ([0] <4.233e-01 5.767e-01 >)
 # ([1] <5.422e-01 4.578e-01 >)
 # ([2] <4.662e-01 5.338e-01 >)
@@ -255,7 +255,7 @@ GBP_MIN                                   9.483e-02   3.078e-02   1.737e-02   1.
 # ([13] <5.266e-01 4.734e-01 >)
 # ([14] <6.033e-01 3.967e-01 >)
 # ([15] <1.558e-01 8.442e-01 >)
-GBP_DELTA                                 6.291e-01   3.350e-01   -3.303e-01  1.000e-09   N/A         
+GBP_DELTA                                 6.291e-01   3.350e-01   -3.303e-01  1.000e-09   6           
 # ([0] <0.000e+00 1.000e+00 >)
 # ([1] <1.000e+00 0.000e+00 >)
 # ([2] <0.000e+00 1.000e+00 >)
@@ -272,7 +272,7 @@ GBP_DELTA                                 6.291e-01   3.350e-01   -3.303e-01  1.
 # ([13] <0.000e+00 1.000e+00 >)
 # ([14] <0.000e+00 1.000e+00 >)
 # ([15] <0.000e+00 1.000e+00 >)
-GBP_LOOP3                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09   N/A         
+GBP_LOOP3                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09   44          
 # ([0] <4.233e-01 5.767e-01 >)
 # ([1] <5.422e-01 4.578e-01 >)
 # ([2] <4.662e-01 5.338e-01 >)
@@ -289,7 +289,7 @@ GBP_LOOP3                                 9.483e-02   3.078e-02   1.737e-02   1.
 # ([13] <5.266e-01 4.734e-01 >)
 # ([14] <6.033e-01 3.967e-01 >)
 # ([15] <1.558e-01 8.442e-01 >)
-GBP_LOOP4                                 7.893e-01   3.569e-01   -3.781e-01  1.000e-09   N/A         
+GBP_LOOP4                                 7.893e-01   3.408e-01   -2.745e-01  1.000e-09   12          
 # ([0] <0.000e+00 1.000e+00 >)
 # ([1] <1.000e+00 0.000e+00 >)
 # ([2] <0.000e+00 1.000e+00 >)
@@ -303,10 +303,10 @@ GBP_LOOP4                                 7.893e-01   3.569e-01   -3.781e-01  1.
 # ([10] <1.000e+00 0.000e+00 >)
 # ([11] <0.000e+00 1.000e+00 >)
 # ([12] <1.338e-115 1.000e+00 >)
-# ([13] <1.000e+00 3.283e-118 >)
-# ([14] <0.000e+00 1.000e+00 >)
+# ([13] <1.000e+00 3.513e-118 >)
+# ([14] <1.000e+00 0.000e+00 >)
 # ([15] <0.000e+00 1.000e+00 >)
-GBP_LOOP5                                 7.893e-01   3.569e-01   -3.781e-01  1.000e-09   N/A         
+GBP_LOOP5                                 7.893e-01   3.408e-01   -2.745e-01  1.000e-09   12          
 # ([0] <0.000e+00 1.000e+00 >)
 # ([1] <1.000e+00 0.000e+00 >)
 # ([2] <0.000e+00 1.000e+00 >)
@@ -320,10 +320,10 @@ GBP_LOOP5                                 7.893e-01   3.569e-01   -3.781e-01  1.
 # ([10] <1.000e+00 0.000e+00 >)
 # ([11] <0.000e+00 1.000e+00 >)
 # ([12] <1.338e-115 1.000e+00 >)
-# ([13] <1.000e+00 3.283e-118 >)
-# ([14] <0.000e+00 1.000e+00 >)
+# ([13] <1.000e+00 3.513e-118 >)
+# ([14] <1.000e+00 0.000e+00 >)
 # ([15] <0.000e+00 1.000e+00 >)
-GBP_LOOP6                                 7.948e-01   4.458e-01   -4.096e-01  1.000e-09   N/A         
+GBP_LOOP6                                 7.948e-01   4.458e-01   -4.096e-01  1.000e-09   4           
 # ([0] <1.000e+00 0.000e+00 >)
 # ([1] <0.000e+00 1.000e+00 >)
 # ([2] <1.000e+00 0.000e+00 >)
@@ -340,7 +340,7 @@ GBP_LOOP6                                 7.948e-01   4.458e-01   -4.096e-01  1.
 # ([13] <1.000e+00 0.000e+00 >)
 # ([14] <0.000e+00 1.000e+00 >)
 # ([15] <0.000e+00 1.000e+00 >)
-GBP_LOOP7                                 7.948e-01   4.458e-01   -4.096e-01  1.000e-09   N/A         
+GBP_LOOP7                                 7.948e-01   4.458e-01   -4.096e-01  1.000e-09   4           
 # ([0] <1.000e+00 0.000e+00 >)
 # ([1] <0.000e+00 1.000e+00 >)
 # ([2] <1.000e+00 0.000e+00 >)
@@ -357,7 +357,7 @@ GBP_LOOP7                                 7.948e-01   4.458e-01   -4.096e-01  1.
 # ([13] <1.000e+00 0.000e+00 >)
 # ([14] <0.000e+00 1.000e+00 >)
 # ([15] <0.000e+00 1.000e+00 >)
-HAK_MIN                                   9.483e-02   3.078e-02   1.737e-02   1.000e-09   N/A         
+HAK_MIN                                   9.483e-02   3.078e-02   1.737e-02   1.000e-09   618         
 # ([0] <4.233e-01 5.767e-01 >)
 # ([1] <5.422e-01 4.578e-01 >)
 # ([2] <4.662e-01 5.338e-01 >)
@@ -374,7 +374,7 @@ HAK_MIN                                   9.483e-02   3.078e-02   1.737e-02   1.
 # ([13] <5.266e-01 4.734e-01 >)
 # ([14] <6.033e-01 3.967e-01 >)
 # ([15] <1.558e-01 8.442e-01 >)
-HAK_DELTA                                 3.684e-01   1.892e-01   9.675e-01   1.000e-09   N/A         
+HAK_DELTA                                 3.684e-01   1.892e-01   9.675e-01   1.000e-09   319         
 # ([0] <4.902e-01 5.098e-01 >)
 # ([1] <5.098e-01 4.902e-01 >)
 # ([2] <4.902e-01 5.098e-01 >)
@@ -391,7 +391,7 @@ HAK_DELTA                                 3.684e-01   1.892e-01   9.675e-01   1.
 # ([13] <5.098e-01 4.902e-01 >)
 # ([14] <5.098e-01 4.902e-01 >)
 # ([15] <4.902e-01 5.098e-01 >)
-HAK_LOOP3                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09   N/A         
+HAK_LOOP3                                 9.483e-02   3.078e-02   1.737e-02   1.000e-09   618         
 # ([0] <4.233e-01 5.767e-01 >)
 # ([1] <5.422e-01 4.578e-01 >)
 # ([2] <4.662e-01 5.338e-01 >)
@@ -408,7 +408,7 @@ HAK_LOOP3                                 9.483e-02   3.078e-02   1.737e-02   1.
 # ([13] <5.266e-01 4.734e-01 >)
 # ([14] <6.033e-01 3.967e-01 >)
 # ([15] <1.558e-01 8.442e-01 >)
-HAK_LOOP4                                 4.970e-03   1.486e-03   -2.503e-04  1.000e-09   N/A         
+HAK_LOOP4                                 4.970e-03   1.486e-03   -2.503e-04  1.000e-09   1580        
 # ([0] <3.909e-01 6.091e-01 >)
 # ([1] <5.556e-01 4.444e-01 >)
 # ([2] <4.585e-01 5.415e-01 >)
@@ -425,7 +425,7 @@ HAK_LOOP4                                 4.970e-03   1.486e-03   -2.503e-04  1.
 # ([13] <5.351e-01 4.649e-01 >)
 # ([14] <6.319e-01 3.681e-01 >)
 # ([15] <1.307e-01 8.693e-01 >)
-HAK_LOOP5                                 4.970e-03   1.486e-03   -2.503e-04  1.000e-09   N/A         
+HAK_LOOP5                                 4.970e-03   1.486e-03   -2.503e-04  1.000e-09   1580        
 # ([0] <3.909e-01 6.091e-01 >)
 # ([1] <5.556e-01 4.444e-01 >)
 # ([2] <4.585e-01 5.415e-01 >)
@@ -442,7 +442,7 @@ HAK_LOOP5                                 4.970e-03   1.486e-03   -2.503e-04  1.
 # ([13] <5.351e-01 4.649e-01 >)
 # ([14] <6.319e-01 3.681e-01 >)
 # ([15] <1.307e-01 8.693e-01 >)
-MR_RESPPROP_FULL                          1.676e-02   4.933e-03   N/A         1.000e-09   N/A         
+MR_RESPPROP_FULL                          1.676e-02   4.933e-03   N/A         1.000e-09   45          
 # ([0] <3.941e-01 6.059e-01 >)
 # ([1] <5.534e-01 4.466e-01 >)
 # ([2] <4.602e-01 5.398e-01 >)
@@ -459,7 +459,7 @@ MR_RESPPROP_FULL                          1.676e-02   4.933e-03   N/A         1.
 # ([13] <5.407e-01 4.593e-01 >)
 # ([14] <6.252e-01 3.748e-01 >)
 # ([15] <1.346e-01 8.654e-01 >)
-MR_CLAMPING_FULL                          8.359e-03   2.508e-03   N/A         1.000e-09   N/A         
+MR_CLAMPING_FULL                          8.359e-03   2.508e-03   N/A         1.000e-09   45          
 # ([0] <3.900e-01 6.100e-01 >)
 # ([1] <5.560e-01 4.440e-01 >)
 # ([2] <4.589e-01 5.411e-01 >)
@@ -476,7 +476,7 @@ MR_CLAMPING_FULL                          8.359e-03   2.508e-03   N/A         1.
 # ([13] <5.404e-01 4.596e-01 >)
 # ([14] <6.273e-01 3.727e-01 >)
 # ([15] <1.344e-01 8.656e-01 >)
-MR_EXACT_FULL                             3.527e-03   1.038e-03   N/A         1.000e-09   N/A         
+MR_EXACT_FULL                             3.527e-03   1.038e-03   N/A         1.000e-09   44          
 # ([0] <3.867e-01 6.133e-01 >)
 # ([1] <5.572e-01 4.428e-01 >)
 # ([2] <4.585e-01 5.415e-01 >)
@@ -493,7 +493,7 @@ MR_EXACT_FULL                             3.527e-03   1.038e-03   N/A         1.
 # ([13] <5.384e-01 4.616e-01 >)
 # ([14] <6.296e-01 3.704e-01 >)
 # ([15] <1.344e-01 8.656e-01 >)
-MR_RESPPROP_LINEAR                        1.932e-02   5.506e-03   N/A         1.000e-09   N/A         
+MR_RESPPROP_LINEAR                        1.932e-02   5.506e-03   N/A         1.000e-09   45          
 # ([0] <3.923e-01 6.077e-01 >)
 # ([1] <5.529e-01 4.471e-01 >)
 # ([2] <4.608e-01 5.392e-01 >)
@@ -510,7 +510,7 @@ MR_RESPPROP_LINEAR                        1.932e-02   5.506e-03   N/A         1.
 # ([13] <5.388e-01 4.612e-01 >)
 # ([14] <6.239e-01 3.761e-01 >)
 # ([15] <1.369e-01 8.631e-01 >)
-MR_CLAMPING_LINEAR                        1.089e-02   3.076e-03   N/A         1.000e-09   N/A         
+MR_CLAMPING_LINEAR                        1.089e-02   3.076e-03   N/A         1.000e-09   45          
 # ([0] <3.884e-01 6.116e-01 >)
 # ([1] <5.551e-01 4.449e-01 >)
 # ([2] <4.598e-01 5.402e-01 >)
@@ -527,7 +527,7 @@ MR_CLAMPING_LINEAR                        1.089e-02   3.076e-03   N/A         1.
 # ([13] <5.383e-01 4.617e-01 >)
 # ([14] <6.257e-01 3.743e-01 >)
 # ([15] <1.370e-01 8.630e-01 >)
-MR_EXACT_LINEAR                           5.617e-03   1.742e-03   N/A         1.000e-09   N/A         
+MR_EXACT_LINEAR                           5.617e-03   1.742e-03   N/A         1.000e-09   44          
 # ([0] <3.853e-01 6.147e-01 >)
 # ([1] <5.560e-01 4.440e-01 >)
 # ([2] <4.596e-01 5.404e-01 >)
@@ -544,7 +544,7 @@ MR_EXACT_LINEAR                           5.617e-03   1.742e-03   N/A         1.
 # ([13] <5.362e-01 4.638e-01 >)
 # ([14] <6.279e-01 3.721e-01 >)
 # ([15] <1.372e-01 8.628e-01 >)
-LCBP_FULLCAV_SEQFIX                       1.225e-03   5.589e-04   N/A         1.000e-09   N/A         
+LCBP_FULLCAV_SEQFIX                       1.225e-03   5.589e-04   N/A         1.000e-09   39          
 # ([0] <3.888e-01 6.112e-01 >)
 # ([1] <5.559e-01 4.441e-01 >)
 # ([2] <4.583e-01 5.417e-01 >)
@@ -561,7 +561,7 @@ LCBP_FULLCAV_SEQFIX                       1.225e-03   5.589e-04   N/A         1.
 # ([13] <5.358e-01 4.642e-01 >)
 # ([14] <6.285e-01 3.715e-01 >)
 # ([15] <1.359e-01 8.641e-01 >)
-LCBP_FULLCAVin_SEQFIX                     1.225e-03   5.589e-04   N/A         1.000e-09   N/A         
+LCBP_FULLCAVin_SEQFIX                     1.225e-03   5.589e-04   N/A         1.000e-09   39          
 # ([0] <3.888e-01 6.112e-01 >)
 # ([1] <5.559e-01 4.441e-01 >)
 # ([2] <4.583e-01 5.417e-01 >)
@@ -578,7 +578,7 @@ LCBP_FULLCAVin_SEQFIX                     1.225e-03   5.589e-04   N/A         1.
 # ([13] <5.358e-01 4.642e-01 >)
 # ([14] <6.285e-01 3.715e-01 >)
 # ([15] <1.359e-01 8.641e-01 >)
-LCBP_FULLCAV_SEQRND                       1.225e-03   5.589e-04   N/A         1.000e-09   N/A         
+LCBP_FULLCAV_SEQRND                       1.225e-03   5.589e-04   N/A         1.000e-09   36          
 # ([0] <3.888e-01 6.112e-01 >)
 # ([1] <5.559e-01 4.441e-01 >)
 # ([2] <4.583e-01 5.417e-01 >)
@@ -595,7 +595,7 @@ LCBP_FULLCAV_SEQRND                       1.225e-03   5.589e-04   N/A         1.
 # ([13] <5.358e-01 4.642e-01 >)
 # ([14] <6.285e-01 3.715e-01 >)
 # ([15] <1.359e-01 8.641e-01 >)
-LCBP_FULLCAVin_SEQRND                     1.225e-03   5.589e-04   N/A         1.000e-09   N/A         
+LCBP_FULLCAVin_SEQRND                     1.225e-03   5.589e-04   N/A         1.000e-09   35          
 # ([0] <3.888e-01 6.112e-01 >)
 # ([1] <5.559e-01 4.441e-01 >)
 # ([2] <4.583e-01 5.417e-01 >)
@@ -612,7 +612,7 @@ LCBP_FULLCAVin_SEQRND                     1.225e-03   5.589e-04   N/A         1.
 # ([13] <5.358e-01 4.642e-01 >)
 # ([14] <6.285e-01 3.715e-01 >)
 # ([15] <1.359e-01 8.641e-01 >)
-LCBP_FULLCAV_NONE                         1.318e-02   2.644e-03   N/A         1.000e+00   N/A         
+LCBP_FULLCAV_NONE                         1.318e-02   2.644e-03   N/A         1.000e+00   0           
 # ([0] <3.859e-01 6.141e-01 >)
 # ([1] <5.569e-01 4.431e-01 >)
 # ([2] <4.719e-01 5.281e-01 >)
@@ -629,7 +629,7 @@ LCBP_FULLCAV_NONE                         1.318e-02   2.644e-03   N/A         1.
 # ([13] <5.350e-01 4.650e-01 >)
 # ([14] <6.290e-01 3.710e-01 >)
 # ([15] <1.364e-01 8.636e-01 >)
-LCBP_FULLCAVin_NONE                       1.318e-02   2.644e-03   N/A         1.000e+00   N/A         
+LCBP_FULLCAVin_NONE                       1.318e-02   2.644e-03   N/A         1.000e+00   0           
 # ([0] <3.859e-01 6.141e-01 >)
 # ([1] <5.569e-01 4.431e-01 >)
 # ([2] <4.719e-01 5.281e-01 >)
@@ -646,7 +646,7 @@ LCBP_FULLCAVin_NONE                       1.318e-02   2.644e-03   N/A         1.
 # ([13] <5.350e-01 4.650e-01 >)
 # ([14] <6.290e-01 3.710e-01 >)
 # ([15] <1.364e-01 8.636e-01 >)
-LCBP_PAIRCAV_SEQFIX                       1.564e-02   5.284e-03   N/A         1.000e-09   N/A         
+LCBP_PAIRCAV_SEQFIX                       1.564e-02   5.284e-03   N/A         1.000e-09   47          
 # ([0] <3.872e-01 6.128e-01 >)
 # ([1] <5.540e-01 4.460e-01 >)
 # ([2] <4.596e-01 5.404e-01 >)
@@ -663,7 +663,7 @@ LCBP_PAIRCAV_SEQFIX                       1.564e-02   5.284e-03   N/A         1.
 # ([13] <5.282e-01 4.718e-01 >)
 # ([14] <6.240e-01 3.760e-01 >)
 # ([15] <1.409e-01 8.591e-01 >)
-LCBP_PAIRCAVin_SEQFIX                     1.564e-02   5.284e-03   N/A         1.000e-09   N/A         
+LCBP_PAIRCAVin_SEQFIX                     1.564e-02   5.284e-03   N/A         1.000e-09   47          
 # ([0] <3.872e-01 6.128e-01 >)
 # ([1] <5.540e-01 4.460e-01 >)
 # ([2] <4.596e-01 5.404e-01 >)
@@ -680,7 +680,7 @@ LCBP_PAIRCAVin_SEQFIX                     1.564e-02   5.284e-03   N/A         1.
 # ([13] <5.282e-01 4.718e-01 >)
 # ([14] <6.240e-01 3.760e-01 >)
 # ([15] <1.409e-01 8.591e-01 >)
-LCBP_PAIRCAV_SEQRND                       1.564e-02   5.284e-03   N/A         1.000e-09   N/A         
+LCBP_PAIRCAV_SEQRND                       1.564e-02   5.284e-03   N/A         1.000e-09   36          
 # ([0] <3.872e-01 6.128e-01 >)
 # ([1] <5.540e-01 4.460e-01 >)
 # ([2] <4.596e-01 5.404e-01 >)
@@ -697,7 +697,7 @@ LCBP_PAIRCAV_SEQRND                       1.564e-02   5.284e-03   N/A         1.
 # ([13] <5.282e-01 4.718e-01 >)
 # ([14] <6.240e-01 3.760e-01 >)
 # ([15] <1.409e-01 8.591e-01 >)
-LCBP_PAIRCAVin_SEQRND                     1.564e-02   5.284e-03   N/A         1.000e-09   N/A         
+LCBP_PAIRCAVin_SEQRND                     1.564e-02   5.284e-03   N/A         1.000e-09   38          
 # ([0] <3.872e-01 6.128e-01 >)
 # ([1] <5.540e-01 4.460e-01 >)
 # ([2] <4.596e-01 5.404e-01 >)
@@ -714,7 +714,7 @@ LCBP_PAIRCAVin_SEQRND                     1.564e-02   5.284e-03   N/A         1.
 # ([13] <5.282e-01 4.718e-01 >)
 # ([14] <6.240e-01 3.760e-01 >)
 # ([15] <1.409e-01 8.591e-01 >)
-LCBP_PAIRCAV_NONE                         1.869e-01   6.816e-02   N/A         1.000e+00   N/A         
+LCBP_PAIRCAV_NONE                         1.869e-01   6.816e-02   N/A         1.000e+00   0           
 # ([0] <3.558e-01 6.442e-01 >)
 # ([1] <5.320e-01 4.680e-01 >)
 # ([2] <4.581e-01 5.419e-01 >)
@@ -731,7 +731,7 @@ LCBP_PAIRCAV_NONE                         1.869e-01   6.816e-02   N/A         1.
 # ([13] <6.323e-01 3.677e-01 >)
 # ([14] <8.161e-01 1.839e-01 >)
 # ([15] <1.676e-01 8.324e-01 >)
-LCBP_PAIRCAVin_NONE                       1.869e-01   6.816e-02   N/A         1.000e+00   N/A         
+LCBP_PAIRCAVin_NONE                       1.869e-01   6.816e-02   N/A         1.000e+00   0           
 # ([0] <3.558e-01 6.442e-01 >)
 # ([1] <5.320e-01 4.680e-01 >)
 # ([2] <4.581e-01 5.419e-01 >)
@@ -748,7 +748,7 @@ LCBP_PAIRCAVin_NONE                       1.869e-01   6.816e-02   N/A         1.
 # ([13] <6.323e-01 3.677e-01 >)
 # ([14] <8.161e-01 1.839e-01 >)
 # ([15] <1.676e-01 8.324e-01 >)
-LCBP_PAIR2CAV_SEQFIX                      1.535e-02   4.445e-03   N/A         1.000e-09   N/A         
+LCBP_PAIR2CAV_SEQFIX                      1.535e-02   4.445e-03   N/A         1.000e-09   47          
 # ([0] <3.844e-01 6.156e-01 >)
 # ([1] <5.557e-01 4.443e-01 >)
 # ([2] <4.588e-01 5.412e-01 >)
@@ -765,7 +765,7 @@ LCBP_PAIR2CAV_SEQFIX                      1.535e-02   4.445e-03   N/A         1.
 # ([13] <5.280e-01 4.720e-01 >)
 # ([14] <6.250e-01 3.750e-01 >)
 # ([15] <1.411e-01 8.589e-01 >)
-LCBP_PAIR2CAVin_SEQFIX                    1.535e-02   4.445e-03   N/A         1.000e-09   N/A         
+LCBP_PAIR2CAVin_SEQFIX                    1.535e-02   4.445e-03   N/A         1.000e-09   47          
 # ([0] <3.844e-01 6.156e-01 >)
 # ([1] <5.557e-01 4.443e-01 >)
 # ([2] <4.588e-01 5.412e-01 >)
@@ -782,7 +782,7 @@ LCBP_PAIR2CAVin_SEQFIX                    1.535e-02   4.445e-03   N/A         1.
 # ([13] <5.280e-01 4.720e-01 >)
 # ([14] <6.250e-01 3.750e-01 >)
 # ([15] <1.411e-01 8.589e-01 >)
-LCBP_PAIR2CAV_SEQRND                      1.535e-02   4.445e-03   N/A         1.000e-09   N/A         
+LCBP_PAIR2CAV_SEQRND                      1.535e-02   4.445e-03   N/A         1.000e-09   39          
 # ([0] <3.844e-01 6.156e-01 >)
 # ([1] <5.557e-01 4.443e-01 >)
 # ([2] <4.588e-01 5.412e-01 >)
@@ -799,7 +799,7 @@ LCBP_PAIR2CAV_SEQRND                      1.535e-02   4.445e-03   N/A         1.
 # ([13] <5.280e-01 4.720e-01 >)
 # ([14] <6.250e-01 3.750e-01 >)
 # ([15] <1.411e-01 8.589e-01 >)
-LCBP_PAIR2CAVin_SEQRND                    1.535e-02   4.445e-03   N/A         1.000e-09   N/A         
+LCBP_PAIR2CAVin_SEQRND                    1.535e-02   4.445e-03   N/A         1.000e-09   35          
 # ([0] <3.844e-01 6.156e-01 >)
 # ([1] <5.557e-01 4.443e-01 >)
 # ([2] <4.588e-01 5.412e-01 >)
@@ -816,7 +816,7 @@ LCBP_PAIR2CAVin_SEQRND                    1.535e-02   4.445e-03   N/A         1.
 # ([13] <5.280e-01 4.720e-01 >)
 # ([14] <6.250e-01 3.750e-01 >)
 # ([15] <1.411e-01 8.589e-01 >)
-LCBP_PAIR2CAV_NONE                        1.894e-01   7.196e-02   N/A         1.000e+00   N/A         
+LCBP_PAIR2CAV_NONE                        1.894e-01   7.196e-02   N/A         1.000e+00   0           
 # ([0] <3.525e-01 6.475e-01 >)
 # ([1] <5.395e-01 4.605e-01 >)
 # ([2] <4.567e-01 5.433e-01 >)
@@ -833,7 +833,7 @@ LCBP_PAIR2CAV_NONE                        1.894e-01   7.196e-02   N/A         1.
 # ([13] <6.444e-01 3.556e-01 >)
 # ([14] <8.185e-01 1.815e-01 >)
 # ([15] <1.841e-01 8.159e-01 >)
-LCBP_PAIR2CAVin_NONE                      1.894e-01   7.196e-02   N/A         1.000e+00   N/A         
+LCBP_PAIR2CAVin_NONE                      1.894e-01   7.196e-02   N/A         1.000e+00   0           
 # ([0] <3.525e-01 6.475e-01 >)
 # ([1] <5.395e-01 4.605e-01 >)
 # ([2] <4.567e-01 5.433e-01 >)
@@ -850,7 +850,7 @@ LCBP_PAIR2CAVin_NONE                      1.894e-01   7.196e-02   N/A         1.
 # ([13] <6.444e-01 3.556e-01 >)
 # ([14] <8.185e-01 1.815e-01 >)
 # ([15] <1.841e-01 8.159e-01 >)
-LCBP_UNICAV_SEQFIX                        9.483e-02   3.078e-02   N/A         1.000e-09   N/A         
+LCBP_UNICAV_SEQFIX                        9.483e-02   3.078e-02   N/A         1.000e-09   43          
 # ([0] <4.233e-01 5.767e-01 >)
 # ([1] <5.422e-01 4.578e-01 >)
 # ([2] <4.662e-01 5.338e-01 >)
@@ -867,7 +867,7 @@ LCBP_UNICAV_SEQFIX                        9.483e-02   3.078e-02   N/A         1.
 # ([13] <5.266e-01 4.734e-01 >)
 # ([14] <6.033e-01 3.967e-01 >)
 # ([15] <1.558e-01 8.442e-01 >)
-LCBP_UNICAV_SEQRND                        9.483e-02   3.078e-02   N/A         1.000e-09   N/A         
+LCBP_UNICAV_SEQRND                        9.483e-02   3.078e-02   N/A         1.000e-09   43          
 # ([0] <4.233e-01 5.767e-01 >)
 # ([1] <5.422e-01 4.578e-01 >)
 # ([2] <4.662e-01 5.338e-01 >)