wolffd@0: /* To compile, type "mex C_quickscore.c" */
wolffd@0: 
wolffd@0: #include <stdio.h>
wolffd@0: #include "nrutil.h"
wolffd@0: #include "nrutil.c"
wolffd@0: #include <math.h>
wolffd@0: #include "mex.h"
wolffd@0: 
wolffd@0: #define MAX(X,Y) (X)>(Y)?(X):(Y)
wolffd@0: 
wolffd@0: int two_to_the(int n)
wolffd@0: {
wolffd@0:   return 1 << n;
wolffd@0: }
wolffd@0: 
wolffd@0: void int2bin(int num, int nbits, int bits[])
wolffd@0: {
wolffd@0:   int i, mask;
wolffd@0:   mask = 1 << (nbits-1); /* mask = 0010...0 , where the 1 is in col nbits (rightmost = col 1) */
wolffd@0:   for (i = 0; i < nbits; i++) {
wolffd@0:     bits[i] = ((num & mask) == 0) ? 0 : 1;
wolffd@0:     num <<= 1;
wolffd@0:   }
wolffd@0: }
wolffd@0: 
wolffd@0: 
wolffd@0: void quickscore(int ndiseases, int nfindings, const double *fpos, int npos, const double *fneg, int nneg,
wolffd@0:                   const double *inhibit, const double *prior, const double *leak, double *prob)
wolffd@0: {
wolffd@0:   double *Pon, *Poff, **Uon, **Uoff, **post, *pterm, *ptermOff, *ptermOn, temp, p, myp;
wolffd@0:   int *bits, nsubsets, *fmask;
wolffd@0:   int f, d, i, j, si, size_subset, sign;
wolffd@0: 
wolffd@0:   Pon = dvector(0, ndiseases);
wolffd@0:   Poff = dvector(0, ndiseases);
wolffd@0:   Pon[0] = 1;
wolffd@0:   Poff[0] = 0;
wolffd@0:   for (i=1; i <= ndiseases; i++) {
wolffd@0:     Pon[i] = prior[i-1];
wolffd@0:     Poff[i] = 1-Pon[i];
wolffd@0:   }
wolffd@0: 
wolffd@0:   Uon = dmatrix(0, nfindings-1, 0, ndiseases);
wolffd@0:   Uoff = dmatrix(0, nfindings-1, 0, ndiseases);
wolffd@0:   d = 0;
wolffd@0:   for (f=0; f < nfindings; f++) {
wolffd@0:     Uon[f][d] = leak[f];
wolffd@0:     Uoff[f][d] = leak[f];
wolffd@0:   }
wolffd@0:   for (f=0; f < nfindings; f++) {
wolffd@0:     for (d=1; d <= ndiseases; d++) {
wolffd@0:       Uon[f][d] = inhibit[f + nfindings*(d-1)];
wolffd@0:       Uoff[f][d] = 1;
wolffd@0:     }
wolffd@0:   }
wolffd@0:   
wolffd@0:   post = dmatrix(0, ndiseases, 0, 1);
wolffd@0:   for (d = 0; d <= ndiseases; d++) {
wolffd@0:     post[d][0] = 0;
wolffd@0:     post[d][1] = 0;
wolffd@0:   }
wolffd@0:   
wolffd@0:   bits = ivector(0, npos-1);
wolffd@0:   fmask = ivector(0, nfindings-1);
wolffd@0:   pterm = dvector(0, ndiseases);
wolffd@0:   ptermOff = dvector(0, ndiseases);
wolffd@0:   ptermOn = dvector(0, ndiseases);
wolffd@0: 
wolffd@0:   nsubsets = two_to_the(npos);
wolffd@0: 
wolffd@0:   for (si = 0; si < nsubsets; si++) {
wolffd@0:     int2bin(si, npos, bits);
wolffd@0:     for (i=0; i < nfindings; i++) fmask[i] = 0;
wolffd@0:     for (i=0; i < nneg; i++) fmask[(int)fneg[i]-1] = 1;
wolffd@0:     size_subset = 0;
wolffd@0:     for (i=0; i < npos; i++) {
wolffd@0:       if (bits[i]) {
wolffd@0: 	size_subset++;
wolffd@0: 	fmask[(int)fpos[i]-1] = 1;
wolffd@0:       }
wolffd@0:     }
wolffd@0:     p = 1;
wolffd@0:     for (d=0; d <= ndiseases; d++) {
wolffd@0:       temp = 1;
wolffd@0:       for (j = 0; j < nfindings; j++) {
wolffd@0: 	if (fmask[j]) temp *= Uoff[j][d];
wolffd@0:       }
wolffd@0:       ptermOff[d] = temp;
wolffd@0: 
wolffd@0:       temp = 1;
wolffd@0:       for (j = 0; j < nfindings; j++) {
wolffd@0: 	if (fmask[j]) temp *= Uon[j][d];
wolffd@0:       }
wolffd@0:       ptermOn[d] = temp;
wolffd@0: 
wolffd@0:       pterm[d] = Poff[d]*ptermOff[d] + Pon[d]*ptermOn[d];
wolffd@0:       p *= pterm[d];
wolffd@0:     }
wolffd@0:     sign = (int) pow(-1, size_subset);
wolffd@0:     for (d=0; d <= ndiseases; d++) {
wolffd@0:       myp = p / pterm[d];
wolffd@0:       post[d][0] += sign*(myp * ptermOff[d]);
wolffd@0:       post[d][1] += sign*(myp * ptermOn[d]);
wolffd@0:     }
wolffd@0:   } /* next si */
wolffd@0: 
wolffd@0:   
wolffd@0:   for (d=0; d <= ndiseases; d++) {
wolffd@0:     post[d][0] *= Poff[d];
wolffd@0:     post[d][1] *= Pon[d];
wolffd@0:   }
wolffd@0:   for (d=0; d <= ndiseases; d++) {
wolffd@0:     temp = post[d][0] + post[d][1];
wolffd@0:     post[d][0] /= temp;
wolffd@0:     post[d][1] /= temp;
wolffd@0:     if (d>0) { prob[d-1] = post[d][1]; }
wolffd@0:   }
wolffd@0: 
wolffd@0:   
wolffd@0:   free_dvector(Pon, 0, ndiseases);
wolffd@0:   free_dvector(Poff, 0, ndiseases);
wolffd@0:   free_dmatrix(Uon, 0, nfindings-1, 0, ndiseases);
wolffd@0:   free_dmatrix(Uoff, 0, nfindings-1, 0, ndiseases);
wolffd@0:   free_dmatrix(post, 0, ndiseases, 0, 1);
wolffd@0:   free_ivector(bits, 0, npos-1);
wolffd@0:   free_ivector(fmask, 0, nfindings-1);
wolffd@0:   free_dvector(pterm, 0, ndiseases);
wolffd@0:   free_dvector(ptermOff, 0, ndiseases);
wolffd@0:   free_dvector(ptermOn, 0, ndiseases);
wolffd@0: }
wolffd@0: 
wolffd@0: 
wolffd@0: void mexFunction(
wolffd@0:                  int nlhs,       mxArray *plhs[],
wolffd@0:                  int nrhs, const mxArray *prhs[]
wolffd@0:                  )
wolffd@0: {
wolffd@0:   double *fpos, *fneg, *inhibit, *prior, *leak, *prob;
wolffd@0:   int npos, nneg, ndiseases, nfindings;
wolffd@0:   double *p;
wolffd@0: 
wolffd@0:   /* read the input args */
wolffd@0:   fpos = mxGetPr(prhs[0]);
wolffd@0:   npos = MAX(mxGetM(prhs[0]), mxGetN(prhs[0]));
wolffd@0: 
wolffd@0:   fneg = mxGetPr(prhs[1]);
wolffd@0:   nneg = MAX(mxGetM(prhs[1]), mxGetN(prhs[1]));
wolffd@0: 
wolffd@0:   inhibit = mxGetPr(prhs[2]); /* inhibit(finding, disease) */
wolffd@0:   nfindings = mxGetM(prhs[2]);
wolffd@0:   ndiseases = mxGetN(prhs[2]);
wolffd@0: 
wolffd@0:   prior = mxGetPr(prhs[3]);
wolffd@0: 
wolffd@0:   leak = mxGetPr(prhs[4]);
wolffd@0: 
wolffd@0: 
wolffd@0:  /* set the output pointers */
wolffd@0:   plhs[0] = mxCreateDoubleMatrix(1, ndiseases, mxREAL);
wolffd@0:   prob = mxGetPr(plhs[0]);
wolffd@0: 
wolffd@0:   quickscore(ndiseases, nfindings, fpos, npos, fneg, nneg, inhibit, prior, leak, prob);
wolffd@0: }
wolffd@0: