AM.xs

#define PERL_NO_GET_CONTEXT
#define NO_XSLOCKS
#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"

#include "ppport.h"

#include <assert.h>

#define NUM_LATTICES 4

/*
 * This program must deal with integers that are too big to be
 * represented by 32 bits.
 *
 * They are represented by AM_BIG_INT, which is typedef'd to
 *
 * unsigned long a[8]
 *
 * where each a[i] < 2*16.  Such an array represents the integer
 *
 * a[0] + a[1] * 2^16 + ... + a[7] * 2^(7*16).
 *
 * We only use 16 bits of the unsigned long instead of 32, so that
 * when we add or multiply two large integers, we have room for overflow.
 * After any addition or multiplication, the result is carried so that
 * each element of the array is again < 2*16.
 *
 * Someday I may rewrite this in assembler.
 *
 */
typedef unsigned short AM_SHORT;
typedef unsigned long AM_LONG;
typedef AM_LONG AM_BIG_INT[8];

#define high_bits(x) x >> 16
#define low_bits(x) x & 0xffff

#define carry(var, ind) \
  var[ind + 1] += high_bits(var[ind]); \
  var[ind] = low_bits(var[ind])

/* carry macros for math using AM_BIG_INT */
#define carry_pointer(p) \
  *(p + 1) += high_bits(*(p)); \
  *(p) = low_bits(*(p))

#define carry_replace(var, ind) \
  var[ind + 1] = high_bits(var[ind]); \
  var[ind] = low_bits(var[ind])

#define hash_pointer_from_stack(ind) \
  (HV *) SvRV(ST(ind))

#define array_pointer_from_stack(ind) \
  AvARRAY((AV *)SvRV(ST(ind)))

#define unsigned_int_from_stack(ind) \
  SvUVX(ST(ind))

/* AM_SUPRAs form a linked list; using for(iter_supra(x, supra)) loops over the list members using the temp variable x */
#define iter_supras(loop_var, supra_ptr) \
  loop_var = supra_ptr + supra_ptr->next; loop_var != supra_ptr; loop_var = supra_ptr + loop_var->next

#define sublist_top(supra) \
  supra->data + supra->data[0] + 1

/*
 * structure for the supracontexts
 */

typedef struct AM_supra {
  /* list of subcontexts in this supracontext
   *
   * data[0] is the number of subcontexts in
   * the array;
   *
   * data[1] is always 0 (useful for finding
   * intersections; see below)
   *
   * data[i] is not an actually subcontext
   * label; instead, all the subcontext labels
   * are kept in an array called subcontext
   * (bad choice of name?)  created in
   * function _fillandcount().  Thus, the
   * actual subcontexts in the supracontext
   * are subcontext[data[2]], ...
   *
   * data[i] < data[i+1] if i > 1 and
   * i < data[0].
   *
   * Using an array of increasing positive
   * integers makes it easy to take
   * intersections (see lattice.pod).
   */
  AM_SHORT *data;

  /* number of supracontexts that contain
   * precisely these subcontexts;
   *
   * According to the AM algorithm, we're
   * supposed to look at all the homogeneous
   * supracontexts to compute the analogical
   * set.  Instead of traversing the
   * supracontextual lattice to find them, we
   * can instead traverse the list of AM_SUPRA
   * with count > 0 and use the value of count
   * to do our computing.
   *
   * Since we're actually traversing four
   * small lattices and taking intersections,
   * we'll be multiplying the four values of
   * count to get what we want.
   *
   */
  AM_SHORT count;

  /*
   * used to implement two linked lists
   *
   * One linked list contains all the nonempty
   * supracontexts (i.e., data[0] is not 0).
   * This linked list is in fact circular.
   *
   * One linked list contains all the unused
   * memory that can be used for new
   * supracontexts.
   */
  AM_SHORT next;

  /*
   * used during the filling of the
   * supracontextual lattice (see below)
   */
  unsigned char touched;
} AM_SUPRA;

/*
 * There is quite a bit of data that must pass between AM.pm and
 * AM.xs.  Instead of repeatedly passing it back and forth on
 * the argument stack, AM.pm sends references to the variables
 * holding this shared data, by calling _xs_initialize() (defined later
 * on).  These pointers are then stored in the following structure,
 * which is put into the magic part of $self (since $self is an HV,
 * it is perforce an SvPVMG as well).
 *
 * Note that for arrays, we store a pointer to the array data itself,
 * not the AV*.  That means that in AM.pm, we have to be careful
 * how we make assignments to array variables; a reassignment such as
 *
 * @sum = (pack "L!8", 0, 0, 0, 0, 0, 0, 0, 0) x @sum;
 *
 * breaks everything because the pointer stored here then won't point
 * to the actual data anymore.  That's why the appropriate line in
 * AM.pm is
 *
 * foreach (@sum) {
 *   $_ = pack "L!8", 0, 0, 0, 0, 0, 0, 0, 0;
 * }
 *
 * Most of the identifiers in the struct have the same names as the
 * variables created in AM.pm and are documented there.  Those
 * that don't are documented below.
 *
 * This trick of storing pointers like this is borrowed from the
 * source code of Perl/Tk.  Thanks, Nick!
 *
 */

typedef struct AM_guts {

  /*
   * Let i be an integer from 0 to 3; this represents which of the
   * four sublattices we are considering.
   *
   * Let lattice = lattice_list[i] and supralist = supra_list[i]; then lattice and
   * supralist taken together tell us which subcontexts are in a
   * particular supracontext.  If s is the label of a supracontext,
   * then it contains the subcontexts listed in
   * supralist[lattice[s]].data[].
   *
   */

  AM_SHORT *lattice_list[NUM_LATTICES];
  AM_SUPRA *supra_list[NUM_LATTICES];

  /* array ref containing number of active features in
   * each lattice (currently we us four lattices)
   */
  SV **lattice_sizes;
  /* array ref containing class labels for whole data set;
   * array index is data item index in data set.
   */
  SV **classes;
  /* TODO: ??? */
  SV **itemcontextchain;
  /* TODO: ??? */
  HV *itemcontextchainhead;
  /* Maps subcontext binary labels to class indices (or 0 if subcontext is heterogeneous) */
  HV *context_to_class;
  /* Maps subcontext binary labels to the number of training items
   * contained in that subcontext
   */
  HV *context_size;
  /* Maps binary context labels to the number of pointers to each,
   * or to the number of pointers to each class label if heterogenous.
   * The key "grand_total" maps to the total number of pointers.
   */
  HV *pointers;
  /* Maps binary context labels to the size of the gang effect of
   * that subcontext. A gang effect is the number of pointers in
   * the given subcontext multiplied by the number of training items
   * contained in the context.
   */
  HV *raw_gang;
  /* number of pointers to each class label;
   * keys are class indices and values are numbers
   * of pointers (AM_BIG_INT).
   */
  SV **sum;
  /*
   * contains the total number of possible class labels;
   * used for computing gang effects.
   */
  IV num_classes;
} AM_GUTS;

/*
 * A function and a vtable necessary for the use of Perl magic
 * TODO: explain the necessity
 */

static int AMguts_mgFree(pTHX_ SV *sv, MAGIC *magic) {
  int i;
  AM_GUTS *guts = (AM_GUTS *) SvPVX(magic->mg_obj);
  for (i = 0; i < NUM_LATTICES; ++i) {
    Safefree(guts->lattice_list[i]);
    Safefree(guts->supra_list[i][0].data);
    Safefree(guts->supra_list[i]);
  }
  return 0;
}

MGVTBL AMguts_vtab = {
  NULL,
  NULL,
  NULL,
  NULL,
  AMguts_mgFree
};

/*
 * arrays used in the change-of-base portion of normalize(SV *s)
 * they are initialized in BOOT
 */

AM_LONG tens[16]; /* 10, 10*2, 10*4, ... */
AM_LONG ones[16]; /*  1,  1*2,  1*4, ... */

/*
 * function: normalize(SV *s)
 *
 * s is an SvPV whose PV* is an unsigned long array representing a very
 * large integer
 *
 * this function modifies s so that its NV is the floating point
 * representation of the very large integer value, while its PV* is
 * the decimal representation of the very large integer value in ASCII
 * (cool, a double-valued scalar)
 *
 * computing the NV is straightforward
 *
 * computing the PV is done using the old change-of-base algorithm:
 * repeatedly divide by 10, and use the remainders to construct the
 * ASCII digits from least to most significant
 *
 */
const unsigned int ASCII_0 = 0x30;
const unsigned int DIVIDE_SPACE = 10;
const int OUTSPACE_SIZE = 55;
void normalize(pTHX_ SV *s) {
  AM_LONG *p = (AM_LONG *)SvPVX(s);

  AM_LONG dspace[DIVIDE_SPACE];
  AM_LONG qspace[DIVIDE_SPACE];
  AM_LONG *dividend, *quotient, *dptr, *qptr;

  STRLEN length = SvCUR(s) / sizeof(AM_LONG);
  /* length indexes into dspace and qspace */
  assert(length <= DIVIDE_SPACE);

  /*
   * outptr iterates outspace from end to beginning, and an ASCII digit is inserted at each location.
   * No need to 0-terminate, since we track the final string length in outlength and pass it to sv_setpvn.
   */
  char outspace[OUTSPACE_SIZE];
  char *outptr;
  outptr = outspace + (OUTSPACE_SIZE - 1);
  unsigned int outlength = 0;

  /* TODO: is this required to be a certain number of bits? */
  long double nn = 0;

  /* nn will be assigned to the NV */
  for (int j = 8; j; --j) {
    /*   2^16    * nn +           p[j-1] */
    nn = 65536.0 * nn + (double) *(p + j - 1);
  }

  dividend = &dspace[0];
  quotient = &qspace[0];
  Copy(p, dividend, length, AM_LONG);

  while (1) {
    while (length && (*(dividend + length - 1) == 0)) {
      --length;
    }
    if (length == 0) {
      sv_setpvn(s, outptr, outlength);
      break;
    }
    dptr = dividend + length - 1;
    qptr = quotient + length - 1;
    AM_LONG carry = 0;
    while (dptr >= dividend) {
      unsigned int i;
      *dptr += carry << 16;
      *qptr = 0;
      for (i = 16; i; ) {
        --i;
        if (tens[i] <= *dptr) {
          *dptr -= tens[i];
          *qptr += ones[i];
        }
      }
      carry = *dptr;
      --dptr;
      --qptr;
    }
    --outptr;
    *outptr = (char)(ASCII_0 + *dividend) & 0x00ff;
    ++outlength;
    AM_LONG *temp = dividend;
    dividend = quotient;
    quotient = temp;
  }

  SvNVX(s) = nn;
  SvNOK_on(s);
}

 /* Given 2 lists of training item indices sorted in descending order,
  * fill a third list with the intersection of items in these lists.
  * This is a simple intersection, and no check for heterogeneity is
  * performed.
  * Return the next empty (available) index address in the third list.
  * If the two lists have no intersection, then the return value is
  * just the same as the third input.
  */
unsigned short *intersect_supras(
    AM_SHORT *intersection_list_top, AM_SHORT *subcontext_list_top, AM_SHORT *k){
  while (1) {
    while (*intersection_list_top > *subcontext_list_top) {
      --intersection_list_top;
    }
    if (*intersection_list_top == 0) {
      break;
    }
    if (*intersection_list_top < *subcontext_list_top) {
      AM_SHORT *temp = intersection_list_top;
      intersection_list_top = subcontext_list_top;
      subcontext_list_top = temp;
      continue;
    }
    *k = *intersection_list_top;
    --intersection_list_top;
    --subcontext_list_top;
    --k;
  }
  return k;
}
 /* The first three inputs are the same as for intersect_supra above,
  * and the fourth paramater should be a list containing the class
  * index for all of the training items. In addition to combining
  * the first two lists into the third via intersection, the final
  * list is checked for heterogeneity and the non-deterministic
  * heterogeneous supracontexts are removed.
  * The return value is the number of items contained in the resulting
  * list.
  */
AM_SHORT intersect_supras_final(
    AM_SHORT *intersection_list_top, AM_SHORT *subcontext_list_top,
    AM_SHORT *intersect, AM_SHORT *subcontext_class){
  AM_SHORT class = 0;
  AM_SHORT length = 0;
  while (1) {
    while (*intersection_list_top > *subcontext_list_top) {
      --intersection_list_top;
    }
    if (*intersection_list_top == 0) {
      break;
    }
    if (*intersection_list_top < *subcontext_list_top) {
      AM_SHORT *temp = intersection_list_top;
      intersection_list_top = subcontext_list_top;
      subcontext_list_top = temp;
      continue;
    }
    *intersect = *intersection_list_top;
    ++intersect;
    ++length;

     /* is it heterogeneous? */
    if (class == 0) {
      /* is it not deterministic? */
      if (length > 1) {
        length = 0;
        break;
      } else {
        class = subcontext_class[*intersection_list_top];
      }
    } else {
      /* Do the classes not match? */
      if (class != subcontext_class[*intersection_list_top]) {
        length = 0;
        break;
      }
    }
    --intersection_list_top;
    --subcontext_list_top;
  }
  return length;
}

/* clear out the supracontexts */
void clear_supras(AM_SUPRA **supra_list, int supras_length)
{
  AM_SUPRA *p;
  for (int i = 0; i < supras_length; i++)
  {
    for (iter_supras(p, supra_list[i]))
    {
      Safefree(p->data);
    }
  }
}

MODULE = Algorithm::AM PACKAGE = Algorithm::AM

PROTOTYPES: DISABLE

BOOT:
  {
    AM_LONG ten = 10;
    AM_LONG one = 1;
    AM_LONG *tensptr = &tens[0];
    AM_LONG *onesptr = &ones[0];
    unsigned int i;
    for (i = 16; i; i--) {
      *tensptr = ten;
      *onesptr = one;
      ++tensptr;
      ++onesptr;
      ten <<= 1;
      one <<= 1;
    }
  }

 /*
  * This function is called by from AM.pm right after creating
  * a blessed reference to Algorithm::AM. It stores the necessary
  * pointers in the AM_GUTS structure and attaches it to the magic
  * part of the reference.
  *
  */

void
_xs_initialize(...)
 PPCODE:
  /* NOT A POINTER THIS TIME! (let memory allocate automatically) */
  AM_GUTS guts;
  /* 9 arguments are passed to the _xs_initialize method: */
  /* $self, the AM object */
  HV *self = hash_pointer_from_stack(0);
  /* For explanations on these, see the comments on AM_guts */
  SV **lattice_sizes = array_pointer_from_stack(1);
  guts.classes = array_pointer_from_stack(2);
  guts.itemcontextchain = array_pointer_from_stack(3);
  guts.itemcontextchainhead = hash_pointer_from_stack(4);
  guts.context_to_class = hash_pointer_from_stack(5);
  guts.context_size = hash_pointer_from_stack(6);
  guts.pointers = hash_pointer_from_stack(7);
  guts.raw_gang = hash_pointer_from_stack(8);
  guts.sum = array_pointer_from_stack(9);
  /* Length of guts.sum */
  guts.num_classes = av_len((AV *) SvRV(ST(9)));

  /*
   * Since the sublattices are small, we just take a chunk of memory
   * here that will be large enough for our purposes and do the actual
   * memory allocation within the code; this reduces the overhead of
   * repeated system calls.
   *
   */

  for (int i = 0; i < NUM_LATTICES; ++i) {
    UV v = SvUVX(lattice_sizes[i]);
    Newxz(guts.lattice_list[i], 1 << v, AM_SHORT);
    Newxz(guts.supra_list[i], 1 << (v + 1), AM_SUPRA); /* CHANGED */ /* TODO: what changed? */
    Newxz(guts.supra_list[i][0].data, 2, AM_SHORT);
  }

  /* Perl magic invoked here */

  SV *svguts = newSVpv((char *)&guts, sizeof(AM_GUTS));
  sv_magic((SV *) self, svguts, PERL_MAGIC_ext, NULL, 0);
  SvRMAGICAL_off((SV *) self);
  MAGIC *magic = mg_find((SV *)self, PERL_MAGIC_ext);
  magic->mg_virtual = &AMguts_vtab;
  mg_magical((SV *) self);

void
_fillandcount(...)
 PPCODE:
  /* Input args are the AM object ($self), number of features in each
   * lattice, and a flag to indicate whether to count occurrences
   * (true) or pointers (false), also known as linear/quadratic.
   */
  HV *self = hash_pointer_from_stack(0);
  SV **lattice_sizes_input = array_pointer_from_stack(1);
  UV linear_flag = unsigned_int_from_stack(2);
  MAGIC *magic = mg_find((SV *)self, PERL_MAGIC_ext);
  AM_GUTS *guts = (AM_GUTS *)SvPVX(magic->mg_obj);

  /*
   * We initialize the memory for the sublattices, including setting up the
   * linked lists.
   */

  AM_SHORT **lattice_list = guts->lattice_list;
  AM_SUPRA **supra_list = guts->supra_list;
  /* this helps us manage the free list in supra_list[i] */
  AM_SHORT nptr[NUM_LATTICES];
  AM_SHORT lattice_sizes[NUM_LATTICES];
  for (int sublattice_index = 0; sublattice_index < NUM_LATTICES; ++sublattice_index) {
    /* Extract numeric values for the specified lattice_sizes */
    lattice_sizes[sublattice_index] = (AM_SHORT) SvUVX(lattice_sizes_input[sublattice_index]);
    /* TODO: explain the lines below */
    Zero(lattice_list[sublattice_index], 1 << lattice_sizes[sublattice_index], AM_SHORT);
    supra_list[sublattice_index][0].next = 0;
    nptr[sublattice_index] = 1;
    for (int i = 1; i < 1 << (lattice_sizes[sublattice_index] + 1); ++i) {/* CHANGED (TODO: changed what?) */
      supra_list[sublattice_index][i].next = (AM_SHORT) i + 1;
    }
  }

  /*
   * Instead of adding subcontext labels directly to the supracontexts,
   * we store all of these labels in an array called subcontext.  We
   * then store the array indices of the subcontext labels in the
   * supracontexts.  That means the list of subcontexts in the
   * supracontexts is an increasing sequence of positive integers, handy
   * for taking intersections (see lattice.pod).
   *
   * The index into the array is called subcontextnumber.
   *
   * The array of matching classes is called subcontext_class.
   *
   */

  HV *context_to_class = guts->context_to_class;
  AM_SHORT subcontextnumber = (AM_SHORT)HvUSEDKEYS(context_to_class);
  AM_SHORT *subcontext;
  Newxz(subcontext, NUM_LATTICES *(subcontextnumber + 1), AM_SHORT);
  subcontext += NUM_LATTICES * subcontextnumber;
  AM_SHORT *subcontext_class;
  Newxz(subcontext_class, subcontextnumber + 1, AM_SHORT);
  subcontext_class += subcontextnumber;

  AM_SHORT *intersectlist, *intersectlist2, *intersectlist3;
  AM_SHORT *ilist2top, *ilist3top;
  Newxz(intersectlist, subcontextnumber + 1, AM_SHORT);
  Newxz(intersectlist2, subcontextnumber + 1, AM_SHORT);
  ilist2top = intersectlist2 + subcontextnumber;
  Newxz(intersectlist3, subcontextnumber + 1, AM_SHORT);
  ilist3top = intersectlist3 + subcontextnumber;

  hv_iterinit(context_to_class);
  HE *context_to_class_entry;
  while ((context_to_class_entry = hv_iternext(context_to_class))) {
    AM_SHORT *contextptr = (AM_SHORT *) HeKEY(context_to_class_entry);
    AM_SHORT class = (AM_SHORT) SvUVX(HeVAL(context_to_class_entry));
    for (int sublattice_index = 0; sublattice_index < NUM_LATTICES; ++sublattice_index, ++contextptr) {
      AM_SHORT active = lattice_sizes[sublattice_index];
      AM_SHORT *lattice = lattice_list[sublattice_index];
      AM_SUPRA *supralist = supra_list[sublattice_index];
      AM_SHORT nextsupra = nptr[sublattice_index];
      AM_SHORT context = *contextptr;

      /* We want to add subcontextnumber to the appropriate
       * supracontexts in the four smaller lattices.
       *
       * Suppose we want to add subcontextnumber to the supracontext
       * labeled by d.  supralist[lattice[d]] is an AM_SUPRA which
       * reflects the current state of the supracontext.  Suppose this
       * state is
       *
       * data:    2 0 x y (i.e., currently contains two subcontexts)
       * count:   5
       * next:    7
       * touched: 0
       *
       * Then we pluck an unused AM_SUPRA off of the free list;
       * suppose that it's located at supralist[9] (the variable
       * nextsupra tells us where).  Then supralist[lattice[d]] will
       * change to
       *
       * data:    2 0 x y
       * count:   4 (decrease by 1)
       * next:    9
       * touched: 1
       *
       * and supralist[9] will become
       *
       * data:    3 0 subcontextnumber x y (now has three subcontexts)
       * count:   1
       * next:    7
       * touched: 0
       *
       * (note: the entries in data[] are added in decreasing order)
       *
       *
       * If, on the other hand, if supralist[lattice[d]] looks like
       *
       * data:    2 0 x y
       * count:   8
       * next:    11
       * touched: 1
       *
       * that means that supralist[11] must look something like
       *
       * data:    3 0 subcontextnumber x y
       * count:   4
       * next:    2
       * touched: 0
       *
       * There already exists a supracontext with subcontextnumber
       * added in!  So we change supralist[lattice[d]] to
       *
       * data:    2 0 x y
       * count:   7 (decrease by 1)
       * next:    11
       * touched: 1
       *
       * change supralist[11] to
       *
       * data:    3 0 subcontextnumber x y
       * count:   5 (increase by 1)
       * next:    2
       * touched: 0
       *
       * and set lattice[d] = 11.
       */

      subcontext[sublattice_index] = context;
      AM_SHORT gaps[16];
      if (context == 0) {
        AM_SUPRA *p;
        for (iter_supras(p, supralist)) {
          AM_SHORT *data;
          Newxz(data, p->data[0] + 3, AM_SHORT);
          Copy(p->data + 2, data + 3, p->data[0], AM_SHORT);
          data[2] = subcontextnumber;
          data[0] = p->data[0] + 1;
          Safefree(p->data);
          p->data = data;
        }
        if (lattice[context] == 0) {

          /* in this case, the subcontext will be
           * added to all supracontexts, so there's
           * no need to hassle with a Gray code and
           * move pointers
           */

          AM_SHORT count = 0;
          AM_SHORT ci = nptr[sublattice_index];
          nptr[sublattice_index] = supralist[ci].next;
          AM_SUPRA *c = supralist + ci;
          c->next = supralist->next;
          supralist->next = ci;
          Newxz(c->data, 3, AM_SHORT);
          c->data[2] = subcontextnumber;
          c->data[0] = 1;
          for (int i = 0; i < (1 << active); ++i) {
            if (lattice[i] == 0) {
              lattice[i] = ci;
              ++count;
            }
          }
          c->count = count;
        }
        continue;
      }

      /* set up traversal using Gray code */
      AM_SHORT d = context;
      AM_SHORT numgaps = 0;
      for (int i = 1 << (active - 1); i; i >>= 1) {
        if (!(i & context)) {
          gaps[numgaps++] = i;
        }
      }
      AM_SHORT t = 1 << numgaps;

      AM_SHORT pi = lattice[context];
      AM_SUPRA *p = supralist + pi;
      if (pi) {
        --(p->count);
      }
      AM_SHORT ci = nextsupra;
      nextsupra = supralist[ci].next;
      p->touched = 1;
      AM_SUPRA *c = supralist + ci;
      c->touched = 0;
      c->next = p->next;
      p->next = ci;
      c->count = 1;
      Newxz(c->data, p->data[0] + 3, AM_SHORT);
      Copy(p->data + 2, c->data + 3, p->data[0], AM_SHORT);
      c->data[2] = subcontextnumber;
      c->data[0] = p->data[0] + 1;
      lattice[context] = ci;

      /* traverse */
      while (--t) {
        AM_SHORT tt;
        int i;
        /* find the rightmost 1 in t; from HAKMEM, I believe */
        for (i = 0, tt = ~t & (t - 1); tt; tt >>= 1, ++i) {
          ;
        }
        d ^= gaps[i];

        p = supralist + (pi = lattice[d]);
        if (pi) {
          --(p->count);
        }
        switch (p->touched) {
          case 1:
            ++supralist[lattice[d] = p->next].count;
            break;
          case 0:
            ci = nextsupra;
            nextsupra = supralist[ci].next;
            p->touched = 1;
            c = supralist + ci;
            c->touched = 0;
            c->next = p->next;
            p->next = ci;
            c->count = 1;
            Newxz(c->data, p->data[0] + 3, AM_SHORT);
            Copy(p->data + 2, c->data + 3, p->data[0], AM_SHORT);
            c->data[2] = subcontextnumber;
            c->data[0] = p->data[0] + 1;
            lattice[d] = ci;
        }
      }

      /* Here we return all AM_SUPRA with count 0 back to the free
       * list and set touched = 0 for all remaining.
       */

      p = supralist;
      p->touched = 0;
      int i;
      do {
        if (supralist[i = p->next].count == 0) {
          Safefree(supralist[i].data);
          p->next = supralist[i].next;
          supralist[i].next = nextsupra;
          nextsupra = (AM_SHORT) i;
        } else {
          p = supralist + p->next;
          p->touched = 0;
        }
      } while (p->next);
      nptr[sublattice_index] = nextsupra;
    } /*end for(sublattice_index = 0...*/
    subcontext -= NUM_LATTICES;
    *subcontext_class = class;
    --subcontext_class;
    --subcontextnumber;
  } /*end while (context_to_class_entry = hv_iternext(...*/

  HV *context_size = guts->context_size;
  HV *pointers = guts->pointers;

  /*
   * The code is in three parts:
   *
   * 1. We successively take one nonempty supracontext from each of the
   *    four small lattices and take their intersection to find a
   *    supracontext of the big lattice.  If at any point we get the
   *    empty set, we move on.
   *
   * 2. We determine if the supracontext so found is heterogeneous; if
   *    so, we skip it.
   *
   * 3. Otherwise, we count the pointers or occurrences.
   *
   */
  {
    /* find intersections */
    AM_SUPRA * p0;
    for (iter_supras(p0, supra_list[0])) {
      AM_SUPRA *p1;
      for (iter_supras(p1, supra_list[1])) {
        /* Find intersection between p0 and p1 */
        AM_SHORT *k = intersect_supras(
          sublist_top(p0),
          sublist_top(p1),
          ilist2top
        );
        /* If k has not been increased then intersection was empty */
        if (k == ilist2top) {
          continue;
        }
        *k = 0;

        AM_SUPRA *p2;
        for (iter_supras(p2, supra_list[2])) {

          /*Find intersection between previous intersection and p2*/
          k = intersect_supras(
            ilist2top,
            sublist_top(p2),
            ilist3top
          );
          /* If k has not been increased then intersection was empty */
          if (k == ilist3top) {
            continue;
          }
          *k = 0;

          AM_SUPRA *p3;
          for (iter_supras(p3, supra_list[3])) {

            /* Find intersection between previous intersection and p3;
             * check for disqualified supras this time.
             */
            AM_SHORT length = intersect_supras_final(
              ilist3top,
              sublist_top(p3),
              intersectlist,
              subcontext_class
            );

            /* count occurrences */
            if (length) {
              AM_BIG_INT count = {0, 0, 0, 0, 0, 0, 0, 0};

              count[0]  = p0->count;

              count[0] *= p1->count;
              carry(count, 0);

              count[0] *= p2->count;
              count[1] *= p2->count;
              carry(count, 0);
              carry(count, 1);

              count[0] *= p3->count;
              count[1] *= p3->count;
              count[2] *= p3->count;
              carry(count, 0);
              carry(count, 1);
              carry(count, 2);
              if(!linear_flag){
                /* If scoring is pointers (quadratic) instead of linear*/
                AM_LONG pointercount = 0;
                for (int i = 0; i < length; ++i) {
                  pointercount += (AM_LONG) SvUV(*hv_fetch(context_size,
                      (char *) (subcontext + (NUM_LATTICES * intersectlist[i])), 8, 0));
                }
                if (pointercount & 0xffff0000) {
                  AM_SHORT pchi = (AM_SHORT) (high_bits(pointercount));
                  AM_SHORT pclo = (AM_SHORT) (low_bits(pointercount));
                  AM_LONG hiprod[6];
                  hiprod[1] = pchi * count[0];
                  hiprod[2] = pchi * count[1];
                  hiprod[3] = pchi * count[2];
                  hiprod[4] = pchi * count[3];
                  count[0] *= pclo;
                  count[1] *= pclo;
                  count[2] *= pclo;
                  count[3] *= pclo;
                  carry(count, 0);
                  carry(count, 1);
                  carry(count, 2);
                  carry(count, 3);

                  count[1] += hiprod[1];
                  count[2] += hiprod[2];
                  count[3] += hiprod[3];
                  count[4] += hiprod[4];
                  carry(count, 1);
                  carry(count, 2);
                  carry(count, 3);
                  carry(count, 4);
                  } else {
                    count[0] *= pointercount;
                    count[1] *= pointercount;
                    count[2] *= pointercount;
                    count[3] *= pointercount;
                    carry(count, 0);
                    carry(count, 1);
                    carry(count, 2);
                    carry(count, 3);
                }
              }
              for (int i = 0; i < length; ++i) {
                SV *final_pointers_sv = *hv_fetch(pointers,
                    (char *) (subcontext + (NUM_LATTICES * intersectlist[i])), 8, 1);
                if (!SvPOK(final_pointers_sv)) {
                  SvUPGRADE(final_pointers_sv, SVt_PVNV);
                  SvGROW(final_pointers_sv, 8 * sizeof(AM_LONG) + 1);
                  Zero(SvPVX(final_pointers_sv), 8, AM_LONG);
                  SvCUR_set(final_pointers_sv, 8 * sizeof(AM_LONG));
                  SvPOK_on(final_pointers_sv);
                }
                AM_LONG *final_pointers = (AM_LONG *) SvPVX(final_pointers_sv);
                for (int j = 0; j < 7; ++j) {
                  *(final_pointers + j) += count[j];
                  carry_pointer(final_pointers + j);
                }
              } /* end for (i = 0;... */
            } /* end if (length) */
          } /* end for (iter_supras(p3... */
        } /* end  for (iter_supras(p2... */
      } /* end  for (iter_supras(p1... */
    } /* end  for (iter_supras(p0... */

    clear_supras(supra_list, 4);

    /*
     * compute analogical set and raw gang effects
     *
     * Technically, we don't compute the analogical set; instead, we
     * compute how many pointers/occurrences there are for each of the
     * data items in a particular subcontext, and associate that number
     * with the subcontext label, not directly with the data item.  We can
     * do this because if two data items are in the same subcontext, they
     * will have the same number of pointers/occurrences.
     *
     * If the user wants the detailed analogical set, it will be created
     * in Result.pm.
     *
     */

    HV *raw_gang = guts->raw_gang;
    SV **classes = guts->classes;
    SV **itemcontextchain = guts->itemcontextchain;
    HV *itemcontextchainhead = guts->itemcontextchainhead;
    SV **sum = guts->sum;
    IV num_classes = guts->num_classes;
    AM_BIG_INT grand_total = {0, 0, 0, 0, 0, 0, 0, 0};
    hv_iterinit(pointers);
    HE * pointers_entry;
    while ((pointers_entry = hv_iternext(pointers))) {
      AM_BIG_INT p;
      Copy(SvPVX(HeVAL(pointers_entry)), p, 8, AM_LONG);

      SV *num_examplars = *hv_fetch(context_size, HeKEY(pointers_entry), NUM_LATTICES * sizeof(AM_SHORT), 0);
      AM_LONG count = (AM_LONG)SvUVX(num_examplars);
      AM_SHORT counthi = (AM_SHORT)(high_bits(count));
      AM_SHORT countlo = (AM_SHORT)(low_bits(count));

      /* initialize 0 because it won't be overwritten */
      /*
       * TODO: multiply through p[7] into gangcount[7]
       * and warn if there's potential overflow
       */
      AM_BIG_INT gangcount;
      gangcount[0] = 0;
      for (int i = 0; i < 7; ++i) {
        gangcount[i] += countlo * p[i];
        carry_replace(gangcount, i);
      }
      gangcount[7] += countlo * p[7];

      /* TODO: why is element 0 not considered here? */
      if (counthi) {
        for (int i = 0; i < 6; ++i) {
          gangcount[i + 1] += counthi * p[i];
          carry(gangcount, i + 1);
        }
      }
      for (int i = 0; i < 7; ++i) {
        grand_total[i] += gangcount[i];
        carry(grand_total, i);
      }
      grand_total[7] += gangcount[7];

      normalize(aTHX_ HeVAL(pointers_entry));

      SV* gang_pointers = *hv_fetch(raw_gang, HeKEY(pointers_entry), NUM_LATTICES * sizeof(AM_SHORT), 1);
      SvUPGRADE(gang_pointers, SVt_PVNV);
      sv_setpvn(gang_pointers, (char *) gangcount, 8 * sizeof(AM_LONG));
      normalize(aTHX_ gang_pointers);

      SV* this_class_sv = *hv_fetch(context_to_class, HeKEY(pointers_entry), NUM_LATTICES * sizeof(AM_SHORT), 0);
      AM_SHORT this_class = (AM_SHORT) SvUVX(this_class_sv);
      if (this_class) {
        AM_LONG *s = (AM_LONG *) SvPVX(sum[this_class]);
        for (int i = 0; i < 7; ++i) {
          *(s + i) += gangcount[i];
          carry_pointer(s + i);
        }
      } else {
      SV *exemplar = *hv_fetch(itemcontextchainhead, HeKEY(pointers_entry), NUM_LATTICES * sizeof(AM_SHORT), 0);
        while (SvIOK(exemplar)) {
          IV datanum = SvIVX(exemplar);
          IV ocnum = SvIVX(classes[datanum]);
          AM_LONG *s = (AM_LONG *) SvPVX(sum[ocnum]);
          for (int i = 0; i < 7; ++i) {
            *(s + i) += p[i];
            carry_pointer(s + i);
            exemplar = itemcontextchain[datanum];
          }
        }
      }
    }
    for (int i = 1; i <= num_classes; ++i) {
      normalize(aTHX_ sum[i]);
    }

    SV *grand_total_entry = *hv_fetch(pointers, "grand_total", 11, 1);
    SvUPGRADE(grand_total_entry, SVt_PVNV);
    sv_setpvn(grand_total_entry, (char *) grand_total, 8 * sizeof(AM_LONG));
    normalize(aTHX_ grand_total_entry);

    Safefree(subcontext);
    Safefree(subcontext_class);
    Safefree(intersectlist);
    Safefree(intersectlist2);
    Safefree(intersectlist3);
  }