gammaFunc.c

Go to the documentation of this file.

/* Copyright 2013 Jorge Merlino

   This file is part of Context.

   Context is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   Context is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with Context.  If not, see <http://www.gnu.org/licenses/>.
*/
 
#include "gammaTable.h"
#include "gammaFunc.h"
#include "debug.h"
#include "reset.h"
#include <math.h> /* for log */

#ifdef DEBUG

static int gammaHits=0;

static int gammaMisses=0;

#endif

#define LNPI_2 M_LNPI / 2

static double alphasizeLog = 0;

static double alphaEntropy = 0;

static double evalLogGamma(int n)
{
  int index = n-2;

  if (index >= 0 && index < TBL_SIZE) {
    DEBUGCODE(gammaHits++);
    return logGammaTbl[index];
  }
  else if (n == 1) {
    DEBUGCODE(gammaHits++);
    return 0;
  }
  else {
    DEBUGCODE(gammaMisses++);
    return gsl_sf_lngamma(n);
  }
}


static double evalLogGamma2(int n)
{
  int index = n-1;

  if (index >= 0 && index < TBL_SIZE) {
    DEBUGCODE(gammaHits++);
    return logGammaTbl2[index];
  }
  else if (n == 0) {
    DEBUGCODE(gammaHits++);
    return LNPI_2; /* Log (gamma (1/2)) */
  }
  else {
    DEBUGCODE(gammaMisses++);
    return gsl_sf_lngamma(n + 0.5);
  }
}


double kt (statistics_t stats) {
  Uint i, ns = 0, alpha_2 = alphasize / 2; /* integer division */
  double sum = 0, gamma1, gamma2, gamma3;
  
  for(i=0; i<alphasize; i++) {
    sum += evalLogGamma2(stats->count[i]); 
    ns += stats->count[i];
  }

  if (alphasize % 2 == 0) {
    gamma1 = evalLogGamma(ns + alpha_2);
    gamma2 = alpha_2 * M_LNPI;
    gamma3 = evalLogGamma(alpha_2);
  }
  else {
    gamma1 = evalLogGamma2(ns + alpha_2);
    gamma2 = (alpha_2 + 0.5) * M_LNPI;
    gamma3 = evalLogGamma2(alpha_2);
  }

  return (gamma1 + gamma2 - gamma3 - sum) / M_LN2;
}

double nodeCost (statistics_t stats) {
  Uint ns = 0, i;
  double sum = 0;
  double loghalf = -0.6931471805599452862;
  double x;

  for (i=0; i<stats->symbolCount; i++) {
    ns += stats->count[stats->symbols[i]];
    if (stats->count[stats->symbols[i]] > 0) {
      sum += evalLogGamma2(stats->count[stats->symbols[i]]-1);
    }
  }
  x = (evalLogGamma(ns) - ((stats->symbolCount - 1) * loghalf) - evalLogGamma(stats->symbolCount) - sum + (stats->symbolCount * 0.5 * M_LNPI)) / M_LN2;
  return x;
}

double escapeCost (statistics_t stats, Uint * distinct) {
  Uint ns = 0, i;
  double sum = 0;
  double loghalf = -0.6931471805599452862;

  for (i=0; i<stats->symbolCount; i++) {
    ns += distinct[stats->symbols[i]];
    if (distinct[stats->symbols[i]] > 0) {
      sum += evalLogGamma2(distinct[stats->symbols[i]]-1);
    }
  }
  return (evalLogGamma(ns) - ((stats->symbolCount - 1) * loghalf) - evalLogGamma(stats->symbolCount) - sum + (stats->symbolCount * 0.5 * M_LNPI)) / M_LN2;
}

double log2Alpha () {
  if (alphasizeLog == 0) {
    alphasizeLog = gsl_sf_log(alphasize) / M_LN2;
  }
  return alphasizeLog;
}


double hAlpha () {
  if (alphaEntropy == 0) {
    double invAlpha = (double)1 / alphasize;
    alphaEntropy = (invAlpha * log2Alpha()) - ((1 - invAlpha) * gsl_sf_log(1 - invAlpha) / M_LN2);
    DEBUGCODE(printf(">>>> %e\n", alphaEntropy));
  }
  return alphaEntropy;
}

#ifdef DEBUG

int getHits() {
  return gammaHits;
}


int getMisses() {
  return gammaMisses;
}

#endif