Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/inference/static/@quickscore_inf_engine/private/C_quickscore.c @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 /* To compile, type "mex C_quickscore.c" */
wolffd@0 2
wolffd@0 3 #include <stdio.h>
wolffd@0 4 #include "nrutil.h"
wolffd@0 5 #include "nrutil.c"
wolffd@0 6 #include <math.h>
wolffd@0 7 #include "mex.h"
wolffd@0 8
wolffd@0 9 #define MAX(X,Y) (X)>(Y)?(X):(Y)
wolffd@0 10
wolffd@0 11 int two_to_the(int n)
wolffd@0 12 {
wolffd@0 13 return 1 << n;
wolffd@0 14 }
wolffd@0 15
wolffd@0 16 void int2bin(int num, int nbits, int bits[])
wolffd@0 17 {
wolffd@0 18 int i, mask;
wolffd@0 19 mask = 1 << (nbits-1); /* mask = 0010...0 , where the 1 is in col nbits (rightmost = col 1) */
wolffd@0 20 for (i = 0; i < nbits; i++) {
wolffd@0 21 bits[i] = ((num & mask) == 0) ? 0 : 1;
wolffd@0 22 num <<= 1;
wolffd@0 23 }
wolffd@0 24 }
wolffd@0 25
wolffd@0 26
wolffd@0 27 void quickscore(int ndiseases, int nfindings, const double *fpos, int npos, const double *fneg, int nneg,
wolffd@0 28 const double *inhibit, const double *prior, const double *leak, double *prob)
wolffd@0 29 {
wolffd@0 30 double *Pon, *Poff, **Uon, **Uoff, **post, *pterm, *ptermOff, *ptermOn, temp, p, myp;
wolffd@0 31 int *bits, nsubsets, *fmask;
wolffd@0 32 int f, d, i, j, si, size_subset, sign;
wolffd@0 33
wolffd@0 34 Pon = dvector(0, ndiseases);
wolffd@0 35 Poff = dvector(0, ndiseases);
wolffd@0 36 Pon[0] = 1;
wolffd@0 37 Poff[0] = 0;
wolffd@0 38 for (i=1; i <= ndiseases; i++) {
wolffd@0 39 Pon[i] = prior[i-1];
wolffd@0 40 Poff[i] = 1-Pon[i];
wolffd@0 41 }
wolffd@0 42
wolffd@0 43 Uon = dmatrix(0, nfindings-1, 0, ndiseases);
wolffd@0 44 Uoff = dmatrix(0, nfindings-1, 0, ndiseases);
wolffd@0 45 d = 0;
wolffd@0 46 for (f=0; f < nfindings; f++) {
wolffd@0 47 Uon[f][d] = leak[f];
wolffd@0 48 Uoff[f][d] = leak[f];
wolffd@0 49 }
wolffd@0 50 for (f=0; f < nfindings; f++) {
wolffd@0 51 for (d=1; d <= ndiseases; d++) {
wolffd@0 52 Uon[f][d] = inhibit[f + nfindings*(d-1)];
wolffd@0 53 Uoff[f][d] = 1;
wolffd@0 54 }
wolffd@0 55 }
wolffd@0 56
wolffd@0 57 post = dmatrix(0, ndiseases, 0, 1);
wolffd@0 58 for (d = 0; d <= ndiseases; d++) {
wolffd@0 59 post[d][0] = 0;
wolffd@0 60 post[d][1] = 0;
wolffd@0 61 }
wolffd@0 62
wolffd@0 63 bits = ivector(0, npos-1);
wolffd@0 64 fmask = ivector(0, nfindings-1);
wolffd@0 65 pterm = dvector(0, ndiseases);
wolffd@0 66 ptermOff = dvector(0, ndiseases);
wolffd@0 67 ptermOn = dvector(0, ndiseases);
wolffd@0 68
wolffd@0 69 nsubsets = two_to_the(npos);
wolffd@0 70
wolffd@0 71 for (si = 0; si < nsubsets; si++) {
wolffd@0 72 int2bin(si, npos, bits);
wolffd@0 73 for (i=0; i < nfindings; i++) fmask[i] = 0;
wolffd@0 74 for (i=0; i < nneg; i++) fmask[(int)fneg[i]-1] = 1;
wolffd@0 75 size_subset = 0;
wolffd@0 76 for (i=0; i < npos; i++) {
wolffd@0 77 if (bits[i]) {
wolffd@0 78 size_subset++;
wolffd@0 79 fmask[(int)fpos[i]-1] = 1;
wolffd@0 80 }
wolffd@0 81 }
wolffd@0 82 p = 1;
wolffd@0 83 for (d=0; d <= ndiseases; d++) {
wolffd@0 84 temp = 1;
wolffd@0 85 for (j = 0; j < nfindings; j++) {
wolffd@0 86 if (fmask[j]) temp *= Uoff[j][d];
wolffd@0 87 }
wolffd@0 88 ptermOff[d] = temp;
wolffd@0 89
wolffd@0 90 temp = 1;
wolffd@0 91 for (j = 0; j < nfindings; j++) {
wolffd@0 92 if (fmask[j]) temp *= Uon[j][d];
wolffd@0 93 }
wolffd@0 94 ptermOn[d] = temp;
wolffd@0 95
wolffd@0 96 pterm[d] = Poff[d]*ptermOff[d] + Pon[d]*ptermOn[d];
wolffd@0 97 p *= pterm[d];
wolffd@0 98 }
wolffd@0 99 sign = (int) pow(-1, size_subset);
wolffd@0 100 for (d=0; d <= ndiseases; d++) {
wolffd@0 101 myp = p / pterm[d];
wolffd@0 102 post[d][0] += sign*(myp * ptermOff[d]);
wolffd@0 103 post[d][1] += sign*(myp * ptermOn[d]);
wolffd@0 104 }
wolffd@0 105 } /* next si */
wolffd@0 106
wolffd@0 107
wolffd@0 108 for (d=0; d <= ndiseases; d++) {
wolffd@0 109 post[d][0] *= Poff[d];
wolffd@0 110 post[d][1] *= Pon[d];
wolffd@0 111 }
wolffd@0 112 for (d=0; d <= ndiseases; d++) {
wolffd@0 113 temp = post[d][0] + post[d][1];
wolffd@0 114 post[d][0] /= temp;
wolffd@0 115 post[d][1] /= temp;
wolffd@0 116 if (d>0) { prob[d-1] = post[d][1]; }
wolffd@0 117 }
wolffd@0 118
wolffd@0 119
wolffd@0 120 free_dvector(Pon, 0, ndiseases);
wolffd@0 121 free_dvector(Poff, 0, ndiseases);
wolffd@0 122 free_dmatrix(Uon, 0, nfindings-1, 0, ndiseases);
wolffd@0 123 free_dmatrix(Uoff, 0, nfindings-1, 0, ndiseases);
wolffd@0 124 free_dmatrix(post, 0, ndiseases, 0, 1);
wolffd@0 125 free_ivector(bits, 0, npos-1);
wolffd@0 126 free_ivector(fmask, 0, nfindings-1);
wolffd@0 127 free_dvector(pterm, 0, ndiseases);
wolffd@0 128 free_dvector(ptermOff, 0, ndiseases);
wolffd@0 129 free_dvector(ptermOn, 0, ndiseases);
wolffd@0 130 }
wolffd@0 131
wolffd@0 132
wolffd@0 133 void mexFunction(
wolffd@0 134 int nlhs, mxArray *plhs[],
wolffd@0 135 int nrhs, const mxArray *prhs[]
wolffd@0 136 )
wolffd@0 137 {
wolffd@0 138 double *fpos, *fneg, *inhibit, *prior, *leak, *prob;
wolffd@0 139 int npos, nneg, ndiseases, nfindings;
wolffd@0 140 double *p;
wolffd@0 141
wolffd@0 142 /* read the input args */
wolffd@0 143 fpos = mxGetPr(prhs[0]);
wolffd@0 144 npos = MAX(mxGetM(prhs[0]), mxGetN(prhs[0]));
wolffd@0 145
wolffd@0 146 fneg = mxGetPr(prhs[1]);
wolffd@0 147 nneg = MAX(mxGetM(prhs[1]), mxGetN(prhs[1]));
wolffd@0 148
wolffd@0 149 inhibit = mxGetPr(prhs[2]); /* inhibit(finding, disease) */
wolffd@0 150 nfindings = mxGetM(prhs[2]);
wolffd@0 151 ndiseases = mxGetN(prhs[2]);
wolffd@0 152
wolffd@0 153 prior = mxGetPr(prhs[3]);
wolffd@0 154
wolffd@0 155 leak = mxGetPr(prhs[4]);
wolffd@0 156
wolffd@0 157
wolffd@0 158 /* set the output pointers */
wolffd@0 159 plhs[0] = mxCreateDoubleMatrix(1, ndiseases, mxREAL);
wolffd@0 160 prob = mxGetPr(plhs[0]);
wolffd@0 161
wolffd@0 162 quickscore(ndiseases, nfindings, fpos, npos, fneg, nneg, inhibit, prior, leak, prob);
wolffd@0 163 }
wolffd@0 164