Algorithm-AM
view release on metacpan or search on metacpan
/*
* 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
--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);
}
}
( run in 0.837 second using v1.01-cache-2.11-cpan-140bd7fdf52 )