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