Browse code

more "R_rlm_interfaces.c" from affyPLM to preprocessCore

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

Ben Bolstad authored on 12/08/2007 18:55:18
Showing1 changed files

1 1
new file mode 100644
... ...
@@ -0,0 +1,221 @@
1
+/*********************************************************************
2
+ **
3
+ ** file: R_rlm_interfaces.c
4
+ **
5
+ ** Aim: Code which provides .Call() interfaces to the rlm code.
6
+ **
7
+ ** Copyright (C) 2006 Ben Bolstad
8
+ **
9
+ ** created by: B. M. Bolstad <bmb@bmbolstad.com>
10
+ ** 
11
+ ** created on: Aug 16, 2006
12
+ **
13
+ ** History
14
+ ** Aug 16, 2006 Initial version
15
+ ** Nov 1, 2006 - add SE to output of function
16
+ **
17
+ **
18
+ **
19
+ *********************************************************************/
20
+
21
+#include <R.h>
22
+#include <Rdefines.h>
23
+#include <Rmath.h>
24
+#include <Rinternals.h>
25
+#include "rlm.h"
26
+#include "rlm_se.h"
27
+
28
+#include "psi_fns.h"
29
+
30
+
31
+
32
+
33
+/**********************************************************************************
34
+ **
35
+ ** SEXP R_rlm_rma_default_model(SEXP Y, SEXP PsiCode, SEXP transform)
36
+ **
37
+ ** 
38
+ ** SEXP Y - A matrix with probes in rows and arrays in columns
39
+ ** SEXP PsiCode - An integer code corresponding to the function that should be used to determine
40
+ **                how outliers are down weighted.
41
+ ** SEXP PsiK - a parameter for weighting algorithm.
42
+ **
43
+ ** Returns 
44
+ ** parameter estimates. weights, residuals, Standard error estimates
45
+ **
46
+ *********************************************************************/
47
+
48
+
49
+
50
+
51
+
52
+SEXP R_rlm_rma_default_model(SEXP Y, SEXP PsiCode, SEXP PsiK){
53
+
54
+
55
+  SEXP R_return_value;
56
+  SEXP R_weights;
57
+  SEXP R_residuals;
58
+  SEXP R_beta;
59
+  SEXP R_SE;
60
+  
61
+  SEXP R_return_value_names;
62
+
63
+  SEXP dim1;
64
+
65
+  double *beta;
66
+  double *residuals;
67
+  double *weights;
68
+  double *se;
69
+
70
+  double residSE;
71
+
72
+  double *Ymat;
73
+
74
+  int rows;
75
+  int cols;
76
+  
77
+  PROTECT(dim1 = getAttrib(Y,R_DimSymbol));
78
+  rows = INTEGER(dim1)[0];
79
+  cols = INTEGER(dim1)[1];
80
+  UNPROTECT(1);
81
+
82
+  PROTECT(R_return_value = allocVector(VECSXP,4));
83
+  PROTECT(R_beta = allocVector(REALSXP, rows + cols));
84
+  PROTECT(R_weights = allocMatrix(REALSXP,rows,cols));
85
+  PROTECT(R_residuals = allocMatrix(REALSXP,rows,cols));
86
+  PROTECT(R_SE = allocVector(REALSXP,rows+cols));
87
+
88
+  SET_VECTOR_ELT(R_return_value,0,R_beta);
89
+  SET_VECTOR_ELT(R_return_value,1,R_weights);
90
+  SET_VECTOR_ELT(R_return_value,2,R_residuals);
91
+  SET_VECTOR_ELT(R_return_value,3,R_SE);
92
+
93
+  UNPROTECT(4);
94
+
95
+  beta = NUMERIC_POINTER(R_beta);
96
+  residuals = NUMERIC_POINTER(R_residuals);
97
+  weights = NUMERIC_POINTER(R_weights);
98
+  se = NUMERIC_POINTER(R_SE);
99
+
100
+  Ymat = NUMERIC_POINTER(Y);
101
+  
102
+  
103
+  
104
+
105
+  
106
+  rlm_fit_anova(Ymat, rows, cols, beta, residuals, weights, PsiFunc(asInteger(PsiCode)),asReal(PsiK), 20, 0);
107
+  
108
+  rlm_compute_se_anova(Ymat, rows, cols, beta, residuals, weights,se, (double *)NULL, &residSE, 4, PsiFunc(asInteger(PsiCode)),asReal(PsiK));
109
+
110
+
111
+  
112
+
113
+
114
+
115
+
116
+
117
+
118
+
119
+
120
+
121
+  PROTECT(R_return_value_names= allocVector(STRSXP,4));
122
+  SET_VECTOR_ELT(R_return_value_names,0,mkChar("Estimates"));
123
+  SET_VECTOR_ELT(R_return_value_names,1,mkChar("Weights"));
124
+  SET_VECTOR_ELT(R_return_value_names,2,mkChar("Residuals"));
125
+  SET_VECTOR_ELT(R_return_value_names,3,mkChar("StdErrors"));
126
+  setAttrib(R_return_value, R_NamesSymbol,R_return_value_names);
127
+  UNPROTECT(2);
128
+  return R_return_value;
129
+
130
+}
131
+
132
+
133
+
134
+
135
+
136
+
137
+
138
+
139
+
140
+
141
+
142
+SEXP R_wrlm_rma_default_model(SEXP Y, SEXP PsiCode, SEXP PsiK, SEXP Weights){
143
+
144
+
145
+  SEXP R_return_value;
146
+  SEXP R_weights;
147
+  SEXP R_residuals;
148
+  SEXP R_beta;
149
+  SEXP R_SE;
150
+  
151
+  SEXP R_return_value_names;
152
+
153
+  SEXP dim1;
154
+
155
+  double *beta;
156
+  double *residuals;
157
+  double *weights;
158
+  double *se;
159
+
160
+  double residSE;
161
+
162
+  double *Ymat;
163
+  double *w;
164
+
165
+  int rows;
166
+  int cols;
167
+  
168
+  PROTECT(dim1 = getAttrib(Y,R_DimSymbol));
169
+  rows = INTEGER(dim1)[0];
170
+  cols = INTEGER(dim1)[1];
171
+  UNPROTECT(1);
172
+
173
+  PROTECT(R_return_value = allocVector(VECSXP,4));
174
+  PROTECT(R_beta = allocVector(REALSXP, rows + cols));
175
+  PROTECT(R_weights = allocMatrix(REALSXP,rows,cols));
176
+  PROTECT(R_residuals = allocMatrix(REALSXP,rows,cols));
177
+  PROTECT(R_SE = allocVector(REALSXP,rows+cols));
178
+
179
+  SET_VECTOR_ELT(R_return_value,0,R_beta);
180
+  SET_VECTOR_ELT(R_return_value,1,R_weights);
181
+  SET_VECTOR_ELT(R_return_value,2,R_residuals);
182
+  SET_VECTOR_ELT(R_return_value,3,R_SE);
183
+
184
+  UNPROTECT(4);
185
+
186
+  beta = NUMERIC_POINTER(R_beta);
187
+  residuals = NUMERIC_POINTER(R_residuals);
188
+  weights = NUMERIC_POINTER(R_weights);
189
+  se = NUMERIC_POINTER(R_SE);
190
+
191
+  Ymat = NUMERIC_POINTER(Y);
192
+  
193
+  w = NUMERIC_POINTER(Weights);
194
+  
195
+
196
+  
197
+  rlm_wfit_anova(Ymat, rows, cols, w, beta, residuals, weights, PsiFunc(asInteger(PsiCode)),asReal(PsiK), 20, 0);
198
+  
199
+  rlm_compute_se_anova(Ymat, rows, cols, beta, residuals, weights,se, (double *)NULL, &residSE, 4, PsiFunc(asInteger(PsiCode)),asReal(PsiK));
200
+
201
+
202
+  
203
+
204
+
205
+
206
+
207
+
208
+
209
+
210
+
211
+
212
+  PROTECT(R_return_value_names= allocVector(STRSXP,4));
213
+  SET_VECTOR_ELT(R_return_value_names,0,mkChar("Estimates"));
214
+  SET_VECTOR_ELT(R_return_value_names,1,mkChar("Weights"));
215
+  SET_VECTOR_ELT(R_return_value_names,2,mkChar("Residuals"));
216
+  SET_VECTOR_ELT(R_return_value_names,3,mkChar("StdErrors"));
217
+  setAttrib(R_return_value, R_NamesSymbol,R_return_value_names);
218
+  UNPROTECT(2);
219
+  return R_return_value;
220
+
221
+}