fisher.c

Go to the documentation of this file.

/*
Copyright (C) 2013 David Novotny and Andrea Vedaldi.
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 "fisher.h"
#include "gmm.h"
#include "mathop.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#ifdef VL_FISHER_INSTANTIATING

static vl_size
VL_XCAT(_vl_fisher_encode_, SFX)
(TYPE * enc,
 TYPE const * means, vl_size dimension, vl_size numClusters,
 TYPE const * covariances,
 TYPE const * priors,
 TYPE const * data, vl_size numData,
 int flags)
{
  vl_size dim;
  vl_index i_cl, i_d;
  vl_size numTerms = 0 ;
  TYPE * posteriors ;
  TYPE * sqrtInvSigma;

  assert(numClusters >= 1) ;
  assert(dimension >= 1) ;

  posteriors = vl_malloc(sizeof(TYPE) * numClusters * numData);
  sqrtInvSigma = vl_malloc(sizeof(TYPE) * dimension * numClusters);

  memset(enc, 0, sizeof(TYPE) * 2 * dimension * numClusters) ;

  for (i_cl = 0 ; i_cl < (signed)numClusters ; ++i_cl) {
    for(dim = 0; dim < dimension; dim++) {
      sqrtInvSigma[i_cl*dimension + dim] = sqrt(1.0 / covariances[i_cl*dimension + dim]);
    }
  }

  VL_XCAT(vl_get_gmm_data_posteriors_, SFX)(posteriors, numClusters, numData,
                                            priors,
                                            means, dimension,
                                            covariances,
                                            data) ;

  /* sparsify posterior assignments with the FAST option */
  if (flags & VL_FISHER_FLAG_FAST) {
    for(i_d = 0; i_d < (signed)numData; i_d++) {
      /* find largest posterior assignment for datum i_d */
      vl_index best = 0 ;
      TYPE bestValue = posteriors[i_d * numClusters] ;
      for (i_cl = 1 ; i_cl < (signed)numClusters; ++ i_cl) {
        TYPE p = posteriors[i_cl + i_d * numClusters] ;
        if (p > bestValue) {
          bestValue = p ;
          best = i_cl ;
        }
      }
      /* make all posterior assignments zero but the best one */
      for (i_cl = 0 ; i_cl < (signed)numClusters; ++ i_cl) {
        posteriors[i_cl + i_d * numClusters] =
        (TYPE)(i_cl == best) ;
      }
    }
  }

#if defined(_OPENMP)
#pragma omp parallel for default(shared) private(i_cl, i_d, dim) num_threads(vl_get_max_threads()) reduction(+:numTerms)
#endif
  for(i_cl = 0; i_cl < (signed)numClusters; ++ i_cl) {
    TYPE uprefix;
    TYPE vprefix;

    TYPE * uk = enc + i_cl*dimension ;
    TYPE * vk = enc + i_cl*dimension + numClusters * dimension ;

    /*
    If the GMM component is degenerate and has a null prior, then it
    must have null posterior as well. Hence it is safe to skip it.  In
    practice, we skip over it even if the prior is very small; if by
    any chance a feature is assigned to such a mode, then its weight
    would be very high due to the division by priors[i_cl] below.
    */
    if (priors[i_cl] < 1e-6) { continue ; }

    for(i_d = 0; i_d < (signed)numData; i_d++) {
      TYPE p = posteriors[i_cl + i_d * numClusters] ;
      if (p < 1e-6) continue ;
      numTerms += 1;
      for(dim = 0; dim < dimension; dim++) {
        TYPE diff = data[i_d*dimension + dim] - means[i_cl*dimension + dim] ;
        diff *= sqrtInvSigma[i_cl*dimension + dim] ;
        *(uk + dim) += p * diff ;
        *(vk + dim) += p * (diff * diff - 1);
      }
    }

    if (numData > 0) {
      uprefix = 1/(numData*sqrt(priors[i_cl]));
      vprefix = 1/(numData*sqrt(2*priors[i_cl]));
      for(dim = 0; dim < dimension; dim++) {
        *(uk + dim) = *(uk + dim) * uprefix;
        *(vk + dim) = *(vk + dim) * vprefix;
      }
    }
  }

  vl_free(posteriors);
  vl_free(sqrtInvSigma) ;

  if (flags & VL_FISHER_FLAG_SQUARE_ROOT) {
    for(dim = 0; dim < 2 * dimension * numClusters ; dim++) {
      TYPE z = enc [dim] ;
      if (z >= 0) {
        enc[dim] = VL_XCAT(vl_sqrt_, SFX)(z) ;
      } else {
        enc[dim] = - VL_XCAT(vl_sqrt_, SFX)(- z) ;
      }
    }
  }

  if (flags & VL_FISHER_FLAG_NORMALIZED) {
    TYPE n = 0 ;
    for(dim = 0 ; dim < 2 * dimension * numClusters ; dim++) {
      TYPE z = enc [dim] ;
      n += z * z ;
    }
    n = VL_XCAT(vl_sqrt_, SFX)(n) ;
    n = VL_MAX(n, 1e-12) ;
    for(dim = 0 ; dim < 2 * dimension * numClusters ; dim++) {
      enc[dim] /= n ;
    }
  }

  return numTerms ;
}

#else
/* not VL_FISHER_INSTANTIATING */

#ifndef __DOXYGEN__
#define FLT VL_TYPE_FLOAT
#define TYPE float
#define SFX f
#define VL_FISHER_INSTANTIATING
#include "fisher.c"

#define FLT VL_TYPE_DOUBLE
#define TYPE double
#define SFX d
#define VL_FISHER_INSTANTIATING
#include "fisher.c"
#endif

/* not VL_FISHER_INSTANTIATING */
#endif

/* ================================================================ */
#ifndef VL_FISHER_INSTANTIATING

VL_EXPORT vl_size
vl_fisher_encode
(void * enc, vl_type dataType,
 void const * means, vl_size dimension, vl_size numClusters,
 void const * covariances,
 void const * priors,
 void const * data,  vl_size numData,
 int flags
)
{
  switch(dataType) {
    case VL_TYPE_FLOAT:
      return _vl_fisher_encode_f
      ((float *) enc,
       (float const *) means, dimension, numClusters,
       (float const *) covariances,
       (float const *) priors,
       (float const *) data, numData,
       flags);
    case VL_TYPE_DOUBLE:
      return _vl_fisher_encode_d
      ((double *) enc,
       (double const *) means, dimension, numClusters,
       (double const *) covariances,
       (double const *) priors,
       (double const *) data, numData,
       flags);
      break;
    default:
      abort();
  }
}
/* not VL_FISHER_INSTANTIATING */
#endif

#undef SFX
#undef TYPE
#undef FLT
#undef VL_FISHER_INSTANTIATING