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
|