Mercurial > hg > camir-aes2014
diff toolboxes/FullBNT-1.0.7/bnt/inference/static/@quickscore_inf_engine/private/C_quickscore.c @ 0:e9a9cd732c1e tip
first hg version after svn
author | wolffd |
---|---|
date | Tue, 10 Feb 2015 15:05:51 +0000 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/FullBNT-1.0.7/bnt/inference/static/@quickscore_inf_engine/private/C_quickscore.c Tue Feb 10 15:05:51 2015 +0000 @@ -0,0 +1,164 @@ +/* To compile, type "mex C_quickscore.c" */ + +#include <stdio.h> +#include "nrutil.h" +#include "nrutil.c" +#include <math.h> +#include "mex.h" + +#define MAX(X,Y) (X)>(Y)?(X):(Y) + +int two_to_the(int n) +{ + return 1 << n; +} + +void int2bin(int num, int nbits, int bits[]) +{ + int i, mask; + mask = 1 << (nbits-1); /* mask = 0010...0 , where the 1 is in col nbits (rightmost = col 1) */ + for (i = 0; i < nbits; i++) { + bits[i] = ((num & mask) == 0) ? 0 : 1; + num <<= 1; + } +} + + +void quickscore(int ndiseases, int nfindings, const double *fpos, int npos, const double *fneg, int nneg, + const double *inhibit, const double *prior, const double *leak, double *prob) +{ + double *Pon, *Poff, **Uon, **Uoff, **post, *pterm, *ptermOff, *ptermOn, temp, p, myp; + int *bits, nsubsets, *fmask; + int f, d, i, j, si, size_subset, sign; + + Pon = dvector(0, ndiseases); + Poff = dvector(0, ndiseases); + Pon[0] = 1; + Poff[0] = 0; + for (i=1; i <= ndiseases; i++) { + Pon[i] = prior[i-1]; + Poff[i] = 1-Pon[i]; + } + + Uon = dmatrix(0, nfindings-1, 0, ndiseases); + Uoff = dmatrix(0, nfindings-1, 0, ndiseases); + d = 0; + for (f=0; f < nfindings; f++) { + Uon[f][d] = leak[f]; + Uoff[f][d] = leak[f]; + } + for (f=0; f < nfindings; f++) { + for (d=1; d <= ndiseases; d++) { + Uon[f][d] = inhibit[f + nfindings*(d-1)]; + Uoff[f][d] = 1; + } + } + + post = dmatrix(0, ndiseases, 0, 1); + for (d = 0; d <= ndiseases; d++) { + post[d][0] = 0; + post[d][1] = 0; + } + + bits = ivector(0, npos-1); + fmask = ivector(0, nfindings-1); + pterm = dvector(0, ndiseases); + ptermOff = dvector(0, ndiseases); + ptermOn = dvector(0, ndiseases); + + nsubsets = two_to_the(npos); + + for (si = 0; si < nsubsets; si++) { + int2bin(si, npos, bits); + for (i=0; i < nfindings; i++) fmask[i] = 0; + for (i=0; i < nneg; i++) fmask[(int)fneg[i]-1] = 1; + size_subset = 0; + for (i=0; i < npos; i++) { + if (bits[i]) { + size_subset++; + fmask[(int)fpos[i]-1] = 1; + } + } + p = 1; + for (d=0; d <= ndiseases; d++) { + temp = 1; + for (j = 0; j < nfindings; j++) { + if (fmask[j]) temp *= Uoff[j][d]; + } + ptermOff[d] = temp; + + temp = 1; + for (j = 0; j < nfindings; j++) { + if (fmask[j]) temp *= Uon[j][d]; + } + ptermOn[d] = temp; + + pterm[d] = Poff[d]*ptermOff[d] + Pon[d]*ptermOn[d]; + p *= pterm[d]; + } + sign = (int) pow(-1, size_subset); + for (d=0; d <= ndiseases; d++) { + myp = p / pterm[d]; + post[d][0] += sign*(myp * ptermOff[d]); + post[d][1] += sign*(myp * ptermOn[d]); + } + } /* next si */ + + + for (d=0; d <= ndiseases; d++) { + post[d][0] *= Poff[d]; + post[d][1] *= Pon[d]; + } + for (d=0; d <= ndiseases; d++) { + temp = post[d][0] + post[d][1]; + post[d][0] /= temp; + post[d][1] /= temp; + if (d>0) { prob[d-1] = post[d][1]; } + } + + + free_dvector(Pon, 0, ndiseases); + free_dvector(Poff, 0, ndiseases); + free_dmatrix(Uon, 0, nfindings-1, 0, ndiseases); + free_dmatrix(Uoff, 0, nfindings-1, 0, ndiseases); + free_dmatrix(post, 0, ndiseases, 0, 1); + free_ivector(bits, 0, npos-1); + free_ivector(fmask, 0, nfindings-1); + free_dvector(pterm, 0, ndiseases); + free_dvector(ptermOff, 0, ndiseases); + free_dvector(ptermOn, 0, ndiseases); +} + + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *fpos, *fneg, *inhibit, *prior, *leak, *prob; + int npos, nneg, ndiseases, nfindings; + double *p; + + /* read the input args */ + fpos = mxGetPr(prhs[0]); + npos = MAX(mxGetM(prhs[0]), mxGetN(prhs[0])); + + fneg = mxGetPr(prhs[1]); + nneg = MAX(mxGetM(prhs[1]), mxGetN(prhs[1])); + + inhibit = mxGetPr(prhs[2]); /* inhibit(finding, disease) */ + nfindings = mxGetM(prhs[2]); + ndiseases = mxGetN(prhs[2]); + + prior = mxGetPr(prhs[3]); + + leak = mxGetPr(prhs[4]); + + + /* set the output pointers */ + plhs[0] = mxCreateDoubleMatrix(1, ndiseases, mxREAL); + prob = mxGetPr(plhs[0]); + + quickscore(ndiseases, nfindings, fpos, npos, fneg, nneg, inhibit, prior, leak, prob); +} +