wotd.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 <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_math.h>
#include <sys/types.h>
#include "types.h"
#include "debug.h"
#include "spacedef.h"
#include "fsmTree.h"
#include "stack.h"
#include "alpha.h"
#include "gammaFunc.h"
#include "statistics.h"
#include "text.h"

#define INTWORDSIZE (UintConst(1) << LOGWORDSIZE)    

#define FIRSTBIT (UintConst(1) << (INTWORDSIZE-1)) 

#define SECONDBIT (FIRSTBIT >> 1)                  

#define BRANCHWIDTH UintConst(3) 

#define MAX_HEIGHT 30

#define NODEINDEX(N) ((Uint) ((N) - streetab))

#define LEAFBIT FIRSTBIT 

#define RIGHTMOSTCHILDBIT SECONDBIT 

#define UNEVALUATEDBIT FIRSTBIT 

#define ISLEAF(P) ((*(P)) & LEAFBIT)

#define ISRIGHTMOSTCHILD(P) ((*(P)) & RIGHTMOSTCHILDBIT)

#define ISUNEVALUATED(P) ((*((P)+1)) & UNEVALUATEDBIT)

#define GETLP(P) ((*(P)) & ~(LEAFBIT | RIGHTMOSTCHILDBIT))

#define GETFIRSTCHILD(P) (*((P)+1))

#define SETLP(P,LP) *(P) = (*(P) & RIGHTMOSTCHILDBIT) | (LP)

#define SETFIRSTCHILD(P,C) *((P)+1) = C

#define SETLEAF(P,L)        *(P) = (L) | LEAFBIT

#define GETSTATS(P)         (*((P)+2))

#define SETSTATS(P,C)       *((P)+2) = C

#define UNDEFREFERENCE      (UINT_MAX)   

#define SUFFIXNUMBER(L)     ((Uint) (*(L) - text))  

#define STOREBOUNDARIES(P,L,R) *(P) = (Uint) ((L) - suffixbase);\
                               *((P)+1) = ((R) - suffixbase) | UNEVALUATEDBIT

#define GETLEFTBOUNDARY(P)  (suffixbase + *(P))

#define GETRIGHTBOUNDARY(P) (suffixbase + ((*((P)+1)) & ~UNEVALUATEDBIT))

#define UNDEFINEDSUCC  (UINT_MAX)    

Uchar *sentinel,               
      **suffixes,              
      **suffixbase,            
      **sbuffer,               
      **sbufferspace = NULL,   
      **bound[UCHAR_MAX+1];    
Uint  occurrence[UCHAR_MAX+1], 
      *streetab = NULL,        
      streetabsize,            
      *nextfreeentry,          
      sbufferwidth,            
      maxsbufferwidth,         
      suffixessize,            
      maxunusedsuffixes,       
      rootchildtab[UCHAR_MAX+1]; 
BOOL  rootevaluated;   


static Uchar **getsbufferspaceeager(Uchar **left,Uchar **right)
{
  Uint width = (Uint) (right-left+1);

  if(width > (Uint) (left-suffixes))
  {
    if(width > sbufferwidth)
    {
      sbufferwidth = width;
      REALLOC(sbufferspace,sbufferspace,Uchar *,sbufferwidth);
    }
    return sbufferspace;
  }
  return left - width;
}


#define MAXSUCCSPACE (BRANCHWIDTH * alphasize + 1)


static void allocstreetab()
{
  Uint tmpindex = NODEINDEX(nextfreeentry); 
  if (tmpindex + MAXSUCCSPACE >= streetabsize) {
    REALLOC(streetab,streetab,Uint,streetabsize + MAXSUCCSPACE);
    streetabsize += MAXSUCCSPACE;
    /* update necessary, since streetab may have been moved. */
    nextfreeentry = streetab + tmpindex;
  }
}


static void sortByChar(Uchar **left, Uchar **right, Uint prefixlen)
{
  Uchar **i, **j, **nextFree = sbuffer;
  Uint a;

  if(*right + prefixlen == sentinel)  /* shortest suffix is sentinel: skip */
  {
    *right += prefixlen;
    right--;
  }
  for(i=left; i<=right; i++) /* determine size for each group */
  {
    *i += prefixlen;         /* drop the common prefix */
    occurrence[(Uint) **i]++;
  }
  for(i=left; i<=right; i++) /* determine right bound for each group */
  {
    a = (Uint) **i;
    if(occurrence[a] > 0)
    {
      bound[a] = nextFree+occurrence[a]-1;
      nextFree = bound[a]+1;
      occurrence[a] = 0;
    }
  }
  for(i=right; i>=left; i--) /* insert suffixes into buffer */
  {
    *(bound[(Uint) **i]--) = *i;
  }
  for(i=left,j=sbuffer; i<=right; i++,j++) /* copy grouped suffixes back */
  {
    *i = *j;
  }
}


static void sortByChar0()
{
  Uchar *cptr, **nextFree = suffixes;
  Uint a;

  for(cptr=text; cptr < text+textlen; cptr++) /* determine size for each group */
  {
    occurrence[(Uint) *cptr]++;
  }
  for(cptr=characters; cptr < characters+alphasize; cptr++)
  {
    a = (Uint) *cptr;
    bound[a] = nextFree+occurrence[a]-1;
    nextFree = bound[a]+1;
    occurrence[a] = 0;
  }
  for(cptr=text+textlen-1; cptr>=text; cptr--) /* insert suffixes into array */
  {
   *(bound[(Uint) *cptr]--) = cptr;
  }
  suffixes[textlen] = sentinel;  /* suffix $ is the largest suffix */
}


static Uint grouplcp(Uchar **left,Uchar **right)
{
  Uchar cmpchar, **i;
  Uint j;

  for(j=UintConst(1); /* nothing */; j++)
  {
    if(*right+j == sentinel)
    {
      return j;
    }
    cmpchar = *(*left+j);
    for(i=left+1; i<=right; i++)
    {
      if(*(*i+j) != cmpchar)
      {
        return j;
      }
    }
  }
}


static Uint evalsuccedges(Uchar **left,Uchar **right)
{
  Uchar firstchar, **r, **l;
  Uint leafnum, firstbranch = UNDEFREFERENCE, *previousnode = NULL;
  BOOL sentineledge = False;

  allocstreetab();
  if(*right == sentinel)
  {
    right--;  /* skip the smallest suffix */
    sentineledge = True;
  }
  for(l=left; l<=right; l=r+1)
  {
    for(firstchar=**l,r=l; r<right && **(r+1)==firstchar; r++)
    {
      /* nothing */ ;
    }
    previousnode = nextfreeentry;
    if(r > l) /* create branching node */
    {
      if(firstbranch == UNDEFREFERENCE)
      {
        firstbranch = NODEINDEX(nextfreeentry);
      }
      STOREBOUNDARIES(nextfreeentry,l,r);
      /* store l and r. resume later with this unevaluated node */
      nextfreeentry += BRANCHWIDTH;
    } else /* create leaf */
    {
      leafnum = SUFFIXNUMBER(l);
      SETLEAF(nextfreeentry,leafnum);
      nextfreeentry++;
    }
  }
  if(sentineledge)
  {
    leafnum = SUFFIXNUMBER(right+1);
    SETLEAF(nextfreeentry,leafnum);
    previousnode = nextfreeentry++;
  }
  assert(previousnode != NULL);
  *previousnode |= RIGHTMOSTCHILDBIT;
  return firstbranch;
}

static Uint evalrootsuccedges(Uchar **left,Uchar **right)
{
  Uchar firstchar, **r, **l;
  Uint *rptr, leafnum, firstbranch = UNDEFREFERENCE;

  for(rptr = rootchildtab; rptr <= rootchildtab + UCHAR_MAX; rptr++)
  {
    *rptr = UNDEFINEDSUCC;
  }
  for(l=left; l<=right; l=r+1) /* first phase */
  {
    for(firstchar=**l,r=l; r<right && **(r+1)==firstchar; r++)
    {
      /* nothing */ ;
    }
    if(r > l) /* create branching node */
    {
      if(firstbranch == UNDEFREFERENCE)
      {
        firstbranch = NODEINDEX(nextfreeentry);
      }
      STOREBOUNDARIES(nextfreeentry,l,r);
      /* store l and r. resume later with this unevaluated branch node */
      rootchildtab[firstchar] = NODEINDEX(nextfreeentry);
      nextfreeentry += BRANCHWIDTH;
    } else /* create leaf */
    {
      leafnum = SUFFIXNUMBER(l);
      SETLEAF(nextfreeentry,leafnum);
      rootchildtab[firstchar] = leafnum | LEAFBIT;
      nextfreeentry++;
    }
  }
  SETLEAF(nextfreeentry,textlen | RIGHTMOSTCHILDBIT);
  nextfreeentry++;
  return firstbranch;
}


static Uint evaluatenodeeager(Uint node, Uint *length)
{
  Uint prefixlen, *nodeptr, unusedsuffixes;
  Uchar **left, **right;

  nodeptr = streetab + node;
  left = GETLEFTBOUNDARY(nodeptr);
  right = GETRIGHTBOUNDARY(nodeptr);
  SETLP(nodeptr,SUFFIXNUMBER(left));
  SETFIRSTCHILD(nodeptr,NODEINDEX(nextfreeentry));

  unusedsuffixes = (Uint) (left - suffixes);
  if(suffixessize > UintConst(10000) && unusedsuffixes > maxunusedsuffixes)
  {
    Uint tmpdiff, width = (Uint) (right - left + 1);
    Uchar **i, **j;
    for(i=left, j=suffixes; i<suffixes+suffixessize; i++, j++)
    {
      *j = *i;  /* move remaining suffixes to the left */
    }
    suffixessize -= unusedsuffixes;
    maxunusedsuffixes = suffixessize >> 1;
    tmpdiff = (Uint) (suffixes - suffixbase);
    REALLOC(suffixes,suffixes,Uchar *,suffixessize);
    suffixbase = suffixes - (tmpdiff + unusedsuffixes);
    left = suffixes;
    right = suffixes + width - 1;
  }
  sbuffer = getsbufferspaceeager(left,right);
  prefixlen = grouplcp(left,right);
  *length = prefixlen;
  sortByChar(left,right,prefixlen);
  return evalsuccedges(left,right);
}


static Uint getnextbranch(Uint previousbranch)
{
  Uint *nodeptr = streetab + previousbranch;

  if(ISRIGHTMOSTCHILD(nodeptr))
  {
    return UNDEFREFERENCE;
  }
  nodeptr += BRANCHWIDTH;
  while(True)
  {
    if(ISLEAF(nodeptr))
    {
      if(ISRIGHTMOSTCHILD(nodeptr))
      {
        return UNDEFREFERENCE;
      }
      nodeptr++;
    } else
    {
      return NODEINDEX(nodeptr);
    }
  }
}

static void pruneRoot() {
  Uint *nodeptr;
  Uint pos, i, idx;
  Uint end;
  statistics_t stats, childStats; 
  double est, auxx;

  Uint * distinct;
  CALLOC(distinct, Uint, alphasize);

  nodeptr = streetab;

  /* counters */
  stats = getStatistics();

  do {
    if (ISLEAF(nodeptr)) {
      if (GETLP(nodeptr) > 0) { /* if false the previous character is outside of string */
    pos = GETLP(nodeptr) - 1; /* position where leaf starts minus previous length */
    idx = alphaindex[*(text+pos)];
    if (stats->count[idx] == 0) {
      stats->symbols[stats->symbolCount++] = idx;
    }
    stats->count[idx]++;
    distinct[idx]++;
      }
      /*stats->cost += log2Alpha();*/
    }
    else {
      childStats = (statistics_t)GETSTATS(nodeptr);
      for (i=0; i<childStats->symbolCount; i++) {
    if (stats->count[childStats->symbols[i]] == 0) {
      stats->symbols[stats->symbolCount++] = childStats->symbols[i];
    } 
    stats->count[childStats->symbols[i]] += childStats->count[childStats->symbols[i]];
    distinct[childStats->symbols[i]]++;
      }
      stats->cost += childStats->cost;
      returnStatistics(childStats);
    }

    end = ISRIGHTMOSTCHILD(nodeptr);
    if (!end) {
      if (ISLEAF(nodeptr)) {
    nodeptr++;
      }
      else {
    nodeptr += BRANCHWIDTH;
      }
    }
  } while (!end);

  stats->cost += hAlpha() * alphasize;

  auxx = escapeCost(stats, distinct);
  stats->cost += auxx;
  assert(auxx >= 0);

  free(distinct);

  est = nodeCost(stats);

  if (est <= stats->cost) { /* pruning needed */
    nextfreeentry = 0; /* indicates an empty tree */
  }
  freeStatistics(stats);
}


static void prune(Uint node, Uint length, Uint branchLength) {
  Uint *nodeptr, *nodeptrPar;
  Uint pos, i, idx;
  statistics_t stats=NULL, childStats; 
  Uint end;
  double est, auxx;

  Uint * distinct;
  CALLOC(distinct, Uint, alphasize);

  nodeptrPar = streetab + node;
  nodeptr = streetab + GETFIRSTCHILD(nodeptrPar);

  do {
    if (ISLEAF(nodeptr)) {
      if (stats == NULL) {
    stats = getStatistics();
      }
      if (GETLP(nodeptr) > length) { /* if false the previous character is outside the string */
    pos = GETLP(nodeptr) - length - 1; /* position where the leaf begins minus previous length */
    idx = alphaindex[*(text+pos)];
    if (stats->count[idx] == 0) {
      stats->symbols[stats->symbolCount++] = idx;
    }
    stats->count[idx]++;
    distinct[idx]++;
      }
    }
    else {
      childStats = (statistics_t)GETSTATS(nodeptr);

      if (stats == NULL) {
    stats = childStats;
    for (i=0; i<stats->symbolCount; i++) {
      distinct[stats->symbols[i]]++;
    } 
      }
      else {
    for (i=0; i<childStats->symbolCount; i++) {
      if (stats->count[childStats->symbols[i]] == 0) {
        stats->symbols[stats->symbolCount++] = childStats->symbols[i];
      } 
      stats->count[childStats->symbols[i]] += childStats->count[childStats->symbols[i]];
      distinct[childStats->symbols[i]]++;
    }
    stats->cost += childStats->cost;
    returnStatistics(childStats);
      }
    }

    end = ISRIGHTMOSTCHILD(nodeptr);
    if (!end) {
      if (ISLEAF(nodeptr)) {
    nodeptr++;
      }
      else {
    nodeptr += BRANCHWIDTH;
      }
    }
  } while (!end);

  stats->cost += hAlpha() * alphasize * branchLength;
  
  auxx = escapeCost(stats, distinct);
  stats->cost += auxx;
  assert(auxx >= 0);

  free(distinct);
  if ((length-branchLength) < MAX_HEIGHT) {
    est = nodeCost(stats);
    assert(est >= 0);

    if (est <= stats->cost) { /* pruning needed */
      nextfreeentry = streetab + GETFIRSTCHILD(nodeptrPar);
      SETFIRSTCHILD(nodeptrPar, UNDEFREFERENCE);
      stats->cost = est;
    }
  }
  else {   
    stats->cost = 0xFFFFFFFF;
    nextfreeentry = streetab + GETFIRSTCHILD(nodeptrPar);
    SETFIRSTCHILD(nodeptrPar, UNDEFREFERENCE);
  }

  SETSTATS(nodeptrPar, (Uint)stats);
}

static void prune2 (Uint node, Uint length) {
  Uint *nodeptr, idx;
  Uchar **left, **right, **i;
  statistics_t stats;

  stats = getStatistics();

  nodeptr = streetab + node;
  left = GETLEFTBOUNDARY(nodeptr);
  right = GETRIGHTBOUNDARY(nodeptr);

  for (i = left; i <= right; i++) {
    if (SUFFIXNUMBER(i) > length) {
      idx = alphaindex[text[SUFFIXNUMBER(i) - length - 1]];
      if (stats->count[idx] == 0) {
    stats->symbols[stats->symbolCount++] = idx;
      }
      stats->count[idx]++;
    }
  }
  stats->cost = 0xFFFFFFFF;
  SETSTATS(nodeptr, (Uint)stats);
  SETFIRSTCHILD(nodeptr, UNDEFREFERENCE);
}

static void evaluateeager() {
  Uint firstbranch, nextbranch, node, newNode, length, branchLength, stacktop=0, stackalloc=0, *stack = NULL;

  sortByChar0();
  firstbranch = evalrootsuccedges(suffixes,suffixes+textlen-1);

  if(firstbranch != UNDEFREFERENCE)
  {
    PUSHNODE(firstbranch);
    PUSHNODE(0);
    PUSHNODE(0);
    while(NOTSTACKEMPTY)
    {
      POPNODE(branchLength);
      POPNODE(length);
      POPNODE(node);
      while(node != UNDEFREFERENCE)
      {
    if (!ISUNEVALUATED(streetab+node)) { /* the node has already been evaluated */
      prune(node, length, branchLength);
      node = UNDEFREFERENCE; /* to continue the while */
    }
    else {
      if((nextbranch = getnextbranch(node)) != UNDEFREFERENCE) {
        PUSHNODE(nextbranch);
        PUSHNODE(length); /* my parent's length */
        PUSHNODE(0); /* does not matter, will be branchLenth when evaluated */
      }

      if (length < MAX_HEIGHT) {
        newNode = evaluatenodeeager(node, &branchLength);
        if (newNode == UNDEFREFERENCE) { /* all children are leaves */
          prune(node, length + branchLength, branchLength);
        }
        else {
          PUSHNODE(node);
          PUSHNODE(branchLength + length);
          PUSHNODE(branchLength);
        }
        node = newNode;
        length += branchLength; /* is the length of the father */
      }
      else {
        prune2(node, length);  
        node = UNDEFREFERENCE; 
      }
    }
      }
    }
    FREE(stack);
  }
}


static void inittree()
{
  Uint i;

  sentinel = text+textlen;
  REALLOC(streetab,streetab,Uint,BRANCHWIDTH + MAXSUCCSPACE);
  streetabsize = BRANCHWIDTH + MAXSUCCSPACE;
  nextfreeentry = streetab;

  /* no partition */
  suffixessize = textlen+1;
  maxunusedsuffixes = suffixessize >> 1;
  CALLOC(suffixes, Uchar *, suffixessize);
  suffixbase = suffixes;

  sbufferwidth = 0;
  maxsbufferwidth = textlen >> 8;
  rootevaluated = False;
  for(i=0; i<=UCHAR_MAX; i++)
  {
    occurrence[i] = 0;
  }
}


static void wotd()
{
  inittree();
  evaluateeager();
  FREE(suffixes);
  FREE(sbufferspace);  
}


Uint getnextsibling(Uint node) {
  Uint *nodeptr = streetab + node;

  if(ISRIGHTMOSTCHILD(nodeptr))
  {
    return UNDEFREFERENCE;
  }
  else {
    if(ISLEAF(nodeptr)) {
      return NODEINDEX(nodeptr + 1);
    }
    else {
      return NODEINDEX(nodeptr + BRANCHWIDTH);
    }
  }  
}


static fsmTree_t buildTree () {
  Uint stacktop=0, stackalloc=0, *stack = NULL, sibling, child, pos, node, *nodeptr, parentptr;
  fsmTree_t ret = initFsmTree(), parent;

  if (nextfreeentry == 0) { /* only root */
    return ret;
  }

  PUSHNODE(0);
  PUSHNODE((Uint)ret);

  while(NOTSTACKEMPTY) {
    POPNODE(parentptr);
    POPNODE(node);
    parent = (fsmTree_t)parentptr;
    nodeptr = streetab + node;

    if (GETLP(nodeptr) != textlen) { /* ignore $ leaves */
      pos = GETINDEX(GETLP(nodeptr));
      parent->children[pos] = initFsmTree();
      parent->children[pos]->parent = parent;
      parent->children[pos]->left = GETLP(nodeptr); 

      /* push the next sibling */
      sibling = getnextsibling(node);
      if (sibling != UNDEFREFERENCE) {
    PUSHNODE(sibling);
    PUSHNODE(parentptr); 
      }

      if (parent->left != -1) { /* != ROOT */
    parent->children[pos]->length = parent->length + parent->right - parent->left + 1;
      }
      if (ISLEAF(nodeptr) || GETFIRSTCHILD(nodeptr) == UNDEFREFERENCE) { /* is a leaf or has been pruned */
    parent->children[pos]->right = GETLP(nodeptr); /* shorten leaf */
      }
      else {
    child = GETFIRSTCHILD(nodeptr);
    parent->children[pos]->right = GETLP(streetab + child) - 1;

    /* push this child */
    PUSHNODE(child);
    PUSHNODE((Uint)parent->children[pos]);
      }
    }
  }

  FREE(stack);
  FREE(streetab);
  return ret;
}


fsmTree_t buildSTree ()
{
  wotd();
  /* as the representation has no root it has to be done specifically */
  pruneRoot();
  freeBuffer();

  return buildTree();
}