aib.c

Go to the documentation of this file.

/*
Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
All rights reserved.

This file is part of the VLFeat library and is made available under
the terms of the BSD license (see the COPYING file).
*/

#include "aib.h"
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>

/* The maximum value which beta may take */
#define BETA_MAX DBL_MAX

void vl_aib_normalize_P (double * P, vl_uint nelem)
{
    vl_uint i;
    double sum = 0;
    for(i=0; i<nelem; i++)
        sum += P[i];
    for(i=0; i<nelem; i++)
        P[i] /= sum;
}

vl_uint *vl_aib_new_nodelist (vl_uint nentries)
{
    vl_uint * nodelist = vl_malloc(sizeof(vl_uint)*nentries);
    vl_uint n;
    for(n=0; n<nentries; n++)
        nodelist[n] = n;

    return nodelist;
}

double * vl_aib_new_Px(double * Pcx, vl_uint nvalues, vl_uint nlabels)
{
    double * Px = vl_malloc(sizeof(double)*nvalues);
    vl_uint r,c;
    for(r=0; r<nvalues; r++)
    {
        double sum = 0;
        for(c=0; c<nlabels; c++)
            sum += Pcx[r*nlabels+c];
        Px[r] = sum;
    }
    return Px;
}

double * vl_aib_new_Pc(double * Pcx, vl_uint nvalues, vl_uint nlabels)
{
    double * Pc = vl_malloc(sizeof(double)*nlabels);
    vl_uint r, c;
    for(c=0; c<nlabels; c++)
    {
        double sum = 0;
        for(r=0; r<nvalues; r++)
            sum += Pcx[r*nlabels+c];
        Pc[c] = sum;
    }
    return Pc;
}

void vl_aib_min_beta
(VlAIB * aib, vl_uint * besti, vl_uint * bestj, double * minbeta)
{
    vl_uint i;
    *minbeta = aib->beta[0];
    *besti   = 0;
    *bestj   = aib->bidx[0];

    for(i=0; i<aib->nentries; i++)
    {
        if(aib->beta[i] < *minbeta)
        {
            *minbeta = aib->beta[i];
            *besti = i;
            *bestj = aib->bidx[i];
        }
    }
}

void
vl_aib_merge_nodes (VlAIB * aib, vl_uint i, vl_uint j, vl_uint new)
{
  vl_uint last_entry = aib->nentries - 1 ;
  vl_uint c, n ;

  /* clear the list of nodes to update */
  aib->nwhich = 0;

  /* make sure that i is smaller than j */
  if(i > j) { vl_uint tmp = j; j = i; i = tmp; }

  /* -----------------------------------------------------------------
   *                    Merge entries i and j, storing the result in i
   * -------------------------------------------------------------- */

  aib-> Px   [i] += aib->Px[j] ;
  aib-> beta [i]  = BETA_MAX ;
  aib-> nodes[i]  = new ;

  for (c = 0; c < aib->nlabels; c++)
    aib-> Pcx [i*aib->nlabels + c] += aib-> Pcx [j*aib->nlabels + c] ;

  /* -----------------------------------------------------------------
   *                                              Move last entry to j
   * -------------------------------------------------------------- */

  aib-> Px    [j]  = aib-> Px    [last_entry];
  aib-> beta  [j]  = aib-> beta  [last_entry];
  aib-> bidx  [j]  = aib-> bidx  [last_entry];
  aib-> nodes [j]  = aib-> nodes [last_entry];

  for (c = 0 ;  c < aib->nlabels ; c++)
    aib-> Pcx[j*aib->nlabels + c] = aib-> Pcx [last_entry*aib->nlabels + c] ;

  /* delete last entry */
  aib-> nentries -- ;

  /* -----------------------------------------------------------------
   *                                        Scan for entries to update
   * -------------------------------------------------------------- */

  /*
   * After mergin entries i and j, we need to update all other entries
   * that had one of these two as closest match. We also need to
   * update the renewend entry i. This is added by the loop below
   * since bidx [i] = j exactly because i was merged.
   *
   * Additionaly, since we moved the last entry back to the entry j,
   * we need to adjust the valeus of bidx to reflect this.
   */

  for (n = 0 ; n < aib->nentries; n++) {
    if(aib->bidx[n] == i || aib->bidx[n] == j) {
        aib->bidx  [n] = 0;
        aib->beta  [n] = BETA_MAX;
        aib->which [aib->nwhich++] = n ;
      }
    else if(aib->bidx[n] == last_entry) {
      aib->bidx[n] = j ;
    }
  }
}

void
vl_aib_update_beta (VlAIB * aib)
{

#define PLOGP(x) ((x)*log((x)))

  vl_uint i;
  double * Px  = aib->Px;
  double * Pcx = aib->Pcx;
  double * tmp = vl_malloc(sizeof(double)*aib->nentries);
  vl_uint a, b, c ;

  /*
   * T1 = I(x,c) - I([x]_ij) = A + B - C
   *
   * A  = \sum_c p(xa,c)           \log ( p(xa,c)          /  p(xa)       )
   * B  = \sum_c p(xb,c)           \log ( p(xb,c)          /  p(xb)       )
   * C  = \sum_c (p(xa,c)+p(xb,c)) \log ((p(xa,c)+p(xb,c)) / (p(xa)+p(xb)))
   *
   * C  = C1 + C2
   * C1 = \sum_c (p(xa,c)+p(xb,c)) \log (p(xa,c)+p(xb,c))
   * C2 = - (p(xa)+p(xb) \log (p(xa)+p(xb))
   */

  /* precalculate A and B */
  for (a = 0; a < aib->nentries; a++) {
    tmp[a] = 0;
    for (c = 0; c < aib->nlabels; c++) {
        double Pac = Pcx [a*aib->nlabels + c] ;
        if(Pac != 0) tmp[a] += Pac * log (Pac / Px[a]) ;
    }
  }

  /* for each entry listed in which */
  for (i = 0 ; i < aib->nwhich; i++) {
    a = aib->which[i];

    /* for each other entry */
    for(b = 0 ; b < aib->nentries ; b++) {
      double T1 = 0 ;

      if (a == b || Px [a] == 0 || Px [b] == 0) continue ;


      T1 = PLOGP ((Px[a] + Px[b])) ;                  /* - C2 */
      T1 += tmp[a] + tmp[b] ;                         /* + A + B */

      for (c = 0 ; c < aib->nlabels; ++ c) {
        double Pac = Pcx [a*aib->nlabels + c] ;
        double Pbc = Pcx [b*aib->nlabels + c] ;
        if (Pac == 0 && Pbc == 0) continue;
        T1 += - PLOGP ((Pac + Pbc)) ;                 /* - C1 */
      }

      /*
       * Now we have beta(a,b). We check wether this is the best beta
       * for entries a and b.
       */
      {
        double beta = T1 ;

        if (beta < aib->beta[a])
          {
            aib->beta[a] = beta;
            aib->bidx[a] = b;
          }
        if (beta < aib->beta[b])
          {
            aib->beta[b] = beta;
            aib->bidx[b] = a;
          }
      }
    }
  }
  vl_free(tmp);
}

void vl_aib_calculate_information(VlAIB * aib, double * I, double * H)
{
    vl_uint r, c;
    *H = 0;
    *I = 0;

    /*
     * H(x)   = - sum_x p(x)    \ log p(x)
     * I(x,c) =   sum_xc p(x,c) \ log (p(x,c) / p(x)p(c))
     */

    /* for each entry */
    for(r = 0 ; r< aib->nentries ; r++) {

      if (aib->Px[r] == 0) continue ;
      *H += -log(aib->Px[r]) * aib->Px[r] ;

      for(c=0; c<aib->nlabels; c++) {
        if (aib->Pcx[r*aib->nlabels+c] == 0) continue;
        *I += aib->Pcx[r*aib->nlabels+c] *
          log (aib->Pcx[r*aib->nlabels+c] / (aib->Px[r]*aib->Pc[c])) ;
      }
    }
}

VlAIB * vl_aib_new(double * Pcx, vl_uint nvalues, vl_uint nlabels)
{
    VlAIB * aib = vl_malloc(sizeof(VlAIB));
    vl_uint i ;

    aib->verbosity = 0 ;
    aib->Pcx   = Pcx ;
    aib->nvalues = nvalues ;
    aib->nlabels = nlabels ;

    vl_aib_normalize_P (aib->Pcx, aib->nvalues * aib->nlabels) ;

    aib->Px = vl_aib_new_Px (aib->Pcx, aib->nvalues, aib->nlabels) ;
    aib->Pc = vl_aib_new_Pc (aib->Pcx, aib->nvalues, aib->nlabels) ;

    aib->nentries = aib->nvalues ;
    aib->nodes    = vl_aib_new_nodelist(aib->nentries) ;
    aib->beta     = vl_malloc(sizeof(double) * aib->nentries) ;
    aib->bidx     = vl_malloc(sizeof(vl_uint)   * aib->nentries) ;

    for(i = 0 ; i < aib->nentries ; i++)
      aib->beta [i] = BETA_MAX ;

    /* Initially we must consider all nodes */
    aib->nwhich = aib->nvalues;
    aib->which  = vl_aib_new_nodelist (aib->nwhich) ;

    aib->parents = vl_malloc(sizeof(vl_uint)*(aib->nvalues*2-1));
    /* Initially, all parents point to a nonexistent node */
    for (i = 0 ; i < 2 * aib->nvalues - 1 ; i++)
      aib->parents [i] = 2 * aib->nvalues ;

    /* Allocate cost output vector */
    aib->costs = vl_malloc (sizeof(double) * (aib->nvalues - 1 + 1)) ;


    return aib ;
}

void
vl_aib_delete (VlAIB * aib)
{
  if (aib) {
    if (aib-> nodes)   vl_free (aib-> nodes);
    if (aib-> beta)    vl_free (aib-> beta);
    if (aib-> bidx)    vl_free (aib-> bidx);
    if (aib-> which)   vl_free (aib-> which);
    if (aib-> Px)      vl_free (aib-> Px);
    if (aib-> Pc)      vl_free (aib-> Pc);
    if (aib-> parents) vl_free (aib-> parents);
    if (aib-> costs)   vl_free (aib-> costs);

    vl_free (aib) ;
  }
}

VL_EXPORT
void vl_aib_process(VlAIB *aib)
{
    vl_uint i, besti, bestj, newnode, nodei, nodej;
    double I, H;
    double minbeta;

    /* Calculate initial value of cost function */
    vl_aib_calculate_information (aib, &I, &H) ;
    aib->costs[0] = I;

    /* Initially which = all */

    /* For each merge */
    for(i = 0 ; i < aib->nvalues - 1 ; i++) {

      /* update entries in aib-> which */
      vl_aib_update_beta(aib);

      /* find best pair of nodes to merge */
      vl_aib_min_beta (aib, &besti, &bestj, &minbeta);

      if(minbeta == BETA_MAX)
        /* only null-probability entries remain */
        break;

      /* Add the parent pointers for the new node */
      newnode = aib->nvalues + i ;
      nodei = aib->nodes[besti];
      nodej = aib->nodes[bestj];

      aib->parents [nodei] = newnode ;
      aib->parents [nodej] = newnode ;
      aib->parents [newnode] = 0 ;

      /* Merge the nodes which produced the minimum beta */
      vl_aib_merge_nodes (aib, besti, bestj, newnode) ;
      vl_aib_calculate_information (aib, &I, &H) ;

      aib->costs[i+1] = I;

      if (aib->verbosity > 0) {
        VL_PRINTF ("aib: (%5d,%5d)=%5d dE: %10.3g I: %6.4g H: %6.4g updt: %5d\n",
                   nodei,
                   nodej,
                   newnode,
                   minbeta,
                   I,
                   H,
                   aib->nwhich) ;
      }
    }

    /* fill ignored entries with NaNs */
    for(; i < aib->nvalues - 1 ; i++)
        aib->costs[i+1] = VL_NAN_D ;
}