IMPLEMENTATION DETAILS


From tests on toy models we have drawn the following conclusions regarding
optimization of BP:

- The log-representation of probabilities turned out to be a factor 2 or 3
  slower than the direct representation.

- Out of the four BP update schemes that we tried (sequential-edge,
  random-edge, sequential-factor, random-factor) the -edge schemes were
  significantly faster than the -factor schemes, and the sequential- schemes were
  faster than the random- schemes. The fastest combination was sequential-edge,
  which usually offered a reduction of about one order of magnitude of the number
  of iterations needed for convergence. 

- Precalculating the indices yields an additional improvement of about a factor 2.

Hence, the fastest BP implementation uses direct representation, precalculated
indices and a sequential-edge update scheme.  I have added a parallel update
scheme for completeness. 

I have also implemented RBP by Koller et al. This seems to yield speed
improvements particularly in difficult cases.

Although it is not the simplest way, the best way to calculate a marginal is
not to add one large delta factor to the network, but to pre-multiply all the
factors with the appropriate delta functions. Even better would be to perform
surgery on the factor graph to simplify its structure.


Precalculating the indices actually yields more improvement than the factor 2
if you also use precalculated indices for calculating the beliefs (and not only
for the messages). Furthermore, when using precalculated indices it is not
necessary to keep tracks of the underlying set<Var> of a probability
distribution; replacing Factors by vector<double> yields an additional speed
improvement. It turns out that the construction, copying and iteration over
set<Var> takes relatively much time. Finally, the decision to use set<size_t>
for the neighbour sets of nodes in the factor graph (FactorGraph::_nb1 and
_nb2) turns out to be not so good: it is significantly faster to iterate over a
vector<size_t> than a set<size_t>. Combining these optimizations yields
approximately a factor 5 improvement in performance for a 20x20 binary Ising
grid. (Thanks go to Vicenc for these lessons).


The main conclusions: don't use set<...> in inner loops when you can also use
vector<...>; and (slightly less general) use a vector<double> instead of a
Factor when you can.


As an alternative, one might consider to use the BOOST Graph library instead of
my own BipartiteGraph (and descendants FactorGraph, RegionGraph)
implementation.

- Copying an array into a vector can be done either using the reserve-push_back approach:
    vector<T> a;
    a.reserve( N );
    for( size_t i = 0; i < N; i++ )
        a.push_back( b[i] );
  or using the resize-[] approach:
    vector<T> a;
    a.resize( N );
    for( size_t i = 0; i < N; i++ )
        a[i] = b[i];
  The first approach is almost twice as fast.


A Factor HAS a Prob and HAS a VarSet. One might think that it would be a good idea to
use multiple inheritance instead, so that a Factor inherits from Prob and VarSet. Then
we would not need to redefine things like Prob::entropy and Prob::operator^ because it
would be automatic. It turns out that this is the case for things like Prob::entropy
(which return a double), but not for operators and other methods that return an object
of the same class, because they then return an object of the ancestor class instead of
of the child class (e.g. operator^ would return a Prob instead of a Factor).


Often, one wants to multiply together some factors, apart from one (where the excluded
factor varies, the totality keeps fixed). In this case, it is generally faster to
cache the product of all factors, and to divide by the excluded factor. This can lead
to problems if entries of some factors contain zeroes. But, since usually (in BP and
LCBP, probably also HAK) one later multiplies with the missing factor anyway (to 
calculate the beliefs e.g.), the places with zeroes do not really matter. A solution that
seems to work is to simply define x/0 = 0 for factor entries. In this way, one can use
the faster divide-out approach instead of the slower multiply-except-one approach.


# Using NOINIT for LCBP_FULLCAV_SEQFIX multiplies total running time by a factor  0.285416
# Using NOINIT for LCBPI_PAIRCAV_SEQFIX multiplies total running time by a factor  0.972613
# Using NOINIT for LCBP_PAIRCAV_SEQFIX multiplies total running time by a factor  0.814879
(cumulative timings on the large set of PROMEDAS potentials); the bigger the network,
the larger the improvement.
