Browse code

add packages to the repository

edge/ pwOmics/ EMDomics/



git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/edge@102444 bc3139a8-67e5-0310-9ffc-ced21a209358

Sonali Arora authored on 14/04/2015 17:06:48
Showing1 changed files
1 1
new file mode 100644
... ...
@@ -0,0 +1,195 @@
1
+#include "edgeKLODP.h"
2
+
3
+/********************************************************************************************
4
+functions for KLODP:
5
+  odpScoreCluster: compute sum of normal densities to be used as numerator or denominator in score
6
+  with c.member.
7
+********************************************************************************************/
8
+
9
+void odpScoreCluster(double *sumDat, double *mu, double *sigma, int *m, int *n, int *p, int *null, int *cluster, double *scr) {
10
+  int i, j, g;
11
+  double *first, *middle;
12
+
13
+  /* if alternative component, set up a couple of vectors */
14
+  /* allocate memory */
15
+   first = vector(0, *m - 1);
16
+
17
+  /* initialize to zero */
18
+	for(i = 0; i < *m; i++)
19
+      first[i] = 0.0;
20
+
21
+  if(*null == 0) {
22
+    /* allocate memory */
23
+    middle = vector(0, *p - 1);
24
+
25
+    /* initialize to zero */
26
+	for(i = 0; i < *p; i++) {
27
+      middle[i] = 0.0;
28
+    }
29
+  }
30
+
31
+  for(i = 0; i < *m; i++) {
32
+	for(j=0; j< *n ; j++){
33
+		  first[i] += sumDat[j + i * *n]*sumDat[j + i * *n];
34
+	}
35
+  }
36
+
37
+  for(i = 0; i < *m; i++) {
38
+    scr[i] = 0.0;
39
+
40
+    for(g = 0; g < *p; g++) {			/* g scans genes */
41
+      /* alternative component */
42
+      if(*null == 0) {
43
+         /* middle[j] += 2 * sumDat[i + (l + 1) * *m] * mu[g + l * *m];*/
44
+		    for(j=0; j< *n ; j++){
45
+			    middle[g] += 2 * sumDat[j + i * *n]*sumDat[j + g * *n + *n * *m];
46
+		    }
47
+  		  /*last[g] += nGrp[l] * mu[g + l * *m] * mu[g + l * *m];*/
48
+		    scr[i] += pow(1 / sigma[g], *n) * exp(-0.5 / sigma[g] / sigma[g] * (first[i] - middle[g] + mu[g])) * cluster[g];
49
+      } else /* null component */
50
+        scr[i] += pow(1 / sigma[g], *n) * exp(-0.5 / sigma[g] / sigma[g] * first[i]) * cluster[g];
51
+    }
52
+	    /* reset vectors to zero, if necessary */
53
+    if(*null == 0) {
54
+      for(g = 0; g < *p; g++) {
55
+        middle[g] = 0.0;
56
+      }
57
+    }
58
+
59
+  }
60
+
61
+  /* free memory, if necessary */
62
+    free_vector(first, 0, *m - 1);
63
+
64
+  if(*null == 0) {
65
+    free_vector(middle, 0, *p - 1);
66
+  }
67
+}
68
+
69
+void kldistance(double *centerFit, double *centerVar, double *fit, double *var, int *m, int *nc, int *n, double *kldd) {
70
+  int i, j, l;
71
+  double sum;
72
+
73
+  for(i = 0; i < *m; i++) {
74
+    for(j = 0; j < *nc; j++) {			/* l scans clusters */
75
+	kldd[j + i* *nc] = 0.0;
76
+	sum = 0.0;
77
+	for(l=0; l< *n ; l++){
78
+		sum += pow((centerFit[l + j* *n]-fit[l + i* *n]),2);
79
+	}
80
+		kldd[j + i* *nc] =  (sum * (1 / centerVar[j] + 1 / var[i]))/2 + *n * (centerVar[j] / var[i] + var[i] / centerVar[j])/2 - *n;
81
+    }
82
+  }
83
+}
84
+
85
+
86
+
87
+/* quicksort routine */
88
+void sortQK(int low, int high, int n, double *w) {
89
+  if(low < high) {
90
+    int lo = low, hi = high + 1;
91
+    double elem = w[low];
92
+    for (;;) {
93
+      while ((lo < n) && (w[++lo] < elem));
94
+      while ((hi >= 0) && (w[--hi] > elem));
95
+      if (lo < hi) swapQK(lo, hi, w);
96
+      else break;
97
+    }
98
+
99
+    swapQK(low, hi, w);
100
+    sortQK(low, hi - 1, n, w);
101
+    sortQK(hi + 1, high, n, w);
102
+  }
103
+}
104
+
105
+/* swap function for use with sortQK() */
106
+void swapQK(int i, int j, double *w) {
107
+  double tmp = w[i];
108
+
109
+  w[i] = w[j];
110
+  w[j] = tmp;
111
+}
112
+
113
+/* allocate a int vector with subscript range v[nl...nh] */
114
+int *ivector(int nl, int nh) {
115
+  int *v;
116
+
117
+  v = (int *) malloc((size_t)((nh - nl + 1 + NR_END) * sizeof(int)));
118
+  if(!v) Rprintf("\n allocation failure in ivector()\n");
119
+  return v - nl + NR_END;
120
+}
121
+
122
+/* free a int vector allocated with ivector() */
123
+void free_ivector(int *v, int nl, int nh) {
124
+  free((FREE_ARG) (v + nl - NR_END));
125
+}
126
+
127
+/* allocate a int matrix with subscript ranges m[nrl...nrh][ncl...nch] */
128
+int **imatrix(int nrl, int nrh, int ncl, int nch) {
129
+  int i, nrow = nrh - nrl + 1, ncol = nch - ncl + 1;
130
+  int **m;
131
+
132
+  /* allocate pointers to rows */
133
+  m = (int **) malloc((size_t)((nrow + NR_END) * sizeof(int*)));
134
+  if(!m) Rprintf("%s", "allocation fialure\n");
135
+
136
+  m += NR_END;
137
+  m -= nrl;
138
+
139
+  /* set pointer to rows */
140
+  m[nrl] = (int *) malloc((size_t)((nrow * ncol + NR_END) * sizeof(int)));
141
+  if(!m[nrl]) Rprintf("%s", "allocation fialure\n");
142
+  m[nrl] += NR_END;
143
+  m[nrl] -= ncl;
144
+
145
+  for(i = nrl + 1; i <= nrh; i++) m[i] = m[i - 1] + ncol;
146
+  return m;
147
+}
148
+
149
+/* free int matrix allocated with imatrix() */
150
+void free_imatrix(int **m, int nrl, int nrh, int ncl, int nch) {
151
+  free((FREE_ARG) (m[nrl] + ncl - NR_END));
152
+  free((FREE_ARG) (m + nrl - NR_END));
153
+}
154
+
155
+/* allocate a double matrix with subscript ranges m[nrl...nrh][ncl...nch] */
156
+double **matrix(int nrl, int nrh, int ncl, int nch) {
157
+  int i, nrow = nrh - nrl + 1, ncol = nch - ncl + 1;
158
+  double **m;
159
+
160
+  /* allocate pointers to rows */
161
+  m = (double **) malloc((size_t)((nrow + NR_END) * sizeof(double*)));
162
+  if(!m) Rprintf("%s", "allocation fialure\n");
163
+
164
+  m += NR_END;
165
+  m -= nrl;
166
+
167
+  /* set pointer to rows */
168
+  m[nrl] = (double *) malloc((size_t)((nrow * ncol + NR_END) * sizeof(double)));
169
+  if(!m[nrl]) Rprintf("%s", "allocation fialure\n");
170
+  m[nrl] += NR_END;
171
+  m[nrl] -= ncl;
172
+
173
+  for(i = nrl + 1; i <= nrh; i++) m[i] = m[i - 1] + ncol;
174
+  return m;
175
+}
176
+
177
+/* free double matrix allocated with matrix() */
178
+void free_matrix(double **m, int nrl, int nrh, int ncl, int nch) {
179
+  free((FREE_ARG) (m[nrl] + ncl - NR_END));
180
+  free((FREE_ARG) (m + nrl - NR_END));
181
+}
182
+
183
+/* allocate a double vector with subscript range v[nl...nh] */
184
+double *vector(int nl, int nh) {
185
+  double *v;
186
+
187
+  v = (double *) malloc((size_t) ((nh - nl + 1 + NR_END) * sizeof(double)));
188
+  if(!v) Rprintf("\n allocation failure in vector()\n");
189
+  return v - nl + NR_END;
190
+}
191
+
192
+/* free double vector allocated with vector() */
193
+void free_vector(double *v, int nl, int nh) {
194
+  free((FREE_ARG) (v + nl - NR_END));
195
+}