1 /* Copyright 2010 by Roger S. Bivand. */
2
3 #include "spdep.h"
4
5 /** */
6 static int c__1 = 1;
7
8 typedef struct opt_error_sse {
9 double *y;
10 double *x;
11 double *yl;
12 double *wy1;
13 double *xlq;
14 double *wx1;
15 double *qy;
16 double *xlqyl;
17 double *work;
18 double *qraux;
19 int *jpvt;
20 int set;
21 } OPT_ERROR_SSE;
22
23 typedef struct hess_error_sse {
24 double *y;
25 double *x;
26 double *yl;
27 double *wy1;
28 double *xl;
29 double *wx1;
30 double *beta1;
31 double *xlb;
32 int set;
33 } HESS_ERROR_SSE;
34
35 typedef struct hess_lag_sse {
36 double *y;
37 double *x;
38 double *yl;
39 double *wy1;
40 double *beta1;
41 double *xb;
42 int set;
43 } HESS_LAG_SSE;
44
opt_error_init()45 SEXP opt_error_init() {
46
47 OPT_ERROR_SSE *pt;
48 SEXP ptr;
49
50 pt = Calloc(1, OPT_ERROR_SSE);
51 pt->set = FALSE;
52
53 PROTECT(ptr = R_MakeExternalPtr(pt, R_NilValue, R_NilValue));
54
55 UNPROTECT(1);
56 return(ptr);
57
58 }
59
60
opt_error_set(SEXP env)61 void opt_error_set(SEXP env) {
62
63 OPT_ERROR_SSE *pt;
64 SEXP y, x, wy, WX;
65 int i, n, p, np, pc=0;
66
67 n = INTEGER_POINTER(findVarInFrame(env, install("n")))[0];
68 p = INTEGER_POINTER(findVarInFrame(env, install("p")))[0];
69 np = n*p;
70
71 pt = (OPT_ERROR_SSE *) R_ExternalPtrAddr(findVarInFrame(env,
72 install("ptr")));
73 if (pt->set) error("opt_error_set: function called out of order");
74
75 PROTECT(y = findVarInFrame(env, install("y"))); pc++;
76 PROTECT(x = findVarInFrame(env, install("x"))); pc++;
77 PROTECT(wy = findVarInFrame(env, install("wy"))); pc++;
78 PROTECT(WX = findVarInFrame(env, install("WX"))); pc++;
79
80 pt->y = Calloc(n, double);
81 pt->x = Calloc(np, double);
82 pt->yl = Calloc(n, double);
83 pt->wy1 = Calloc(n, double);
84 pt->xlq = Calloc(np, double);
85 pt->wx1 = Calloc(np, double);
86 pt->qy = Calloc(np, double);
87 pt->xlqyl = Calloc(p, double);
88 pt->jpvt = Calloc(p, int);
89 pt->work = Calloc(p*2, double);
90 /* pt->work = Calloc(p, double); */
91 pt->qraux = Calloc(p, double);
92
93 for (i=0; i<n; i++) {
94 pt->y[i] = NUMERIC_POINTER(y)[i];
95 pt->wy1[i] = NUMERIC_POINTER(wy)[i];
96 }
97 for (i=0; i<np; i++) {
98 pt->x[i] = NUMERIC_POINTER(x)[i];
99 pt->wx1[i] = NUMERIC_POINTER(WX)[i];
100 }
101 pt->set = TRUE;
102 UNPROTECT(pc);
103
104 return;
105 }
106
107
opt_error_free(SEXP ptr)108 SEXP opt_error_free(SEXP ptr) {
109
110 OPT_ERROR_SSE *pt;
111
112 pt = (OPT_ERROR_SSE *) R_ExternalPtrAddr(ptr);
113
114 Free(pt->qraux);
115 Free(pt->work);
116 Free(pt->jpvt);
117 Free(pt->xlqyl);
118 Free(pt->qy);
119 Free(pt->wx1);
120 Free(pt->xlq);
121 Free(pt->wy1);
122 Free(pt->yl);
123 Free(pt->x);
124 Free(pt->y);
125
126 Free(pt);
127 R_ClearExternalPtr(ptr);
128 return(R_NilValue);
129 }
130
hess_error_init()131 SEXP hess_error_init() {
132
133 HESS_ERROR_SSE *pt;
134 SEXP ptr;
135
136 pt = Calloc(1, HESS_ERROR_SSE);
137 pt->set = FALSE;
138
139 PROTECT(ptr = R_MakeExternalPtr(pt, R_NilValue, R_NilValue));
140
141 UNPROTECT(1);
142 return(ptr);
143
144 }
145
146
hess_error_set(SEXP env)147 void hess_error_set(SEXP env) {
148
149 HESS_ERROR_SSE *pt;
150 SEXP y, x, wy, WX;
151 int i, n, p, np, pc=0;
152
153 n = INTEGER_POINTER(findVarInFrame(env, install("n")))[0];
154 p = INTEGER_POINTER(findVarInFrame(env, install("p")))[0];
155 np = n*p;
156
157 pt = (HESS_ERROR_SSE *) R_ExternalPtrAddr(findVarInFrame(env,
158 install("ptr")));
159 if (pt->set) error("hess_error_set: function called out of order");
160
161 PROTECT(y = findVarInFrame(env, install("y"))); pc++;
162 PROTECT(x = findVarInFrame(env, install("x"))); pc++;
163 PROTECT(wy = findVarInFrame(env, install("wy"))); pc++;
164 PROTECT(WX = findVarInFrame(env, install("WX"))); pc++;
165
166 pt->y = Calloc(n, double);
167 pt->x = Calloc(np, double);
168 pt->yl = Calloc(n, double);
169 pt->wy1 = Calloc(n, double);
170 pt->xl = Calloc(np, double);
171 pt->wx1 = Calloc(np, double);
172 pt->beta1 = Calloc(p, double);
173 pt->xlb = Calloc(n, double);
174
175 for (i=0; i<n; i++) {
176 pt->y[i] = NUMERIC_POINTER(y)[i];
177 pt->wy1[i] = NUMERIC_POINTER(wy)[i];
178 }
179 for (i=0; i<np; i++) {
180 pt->x[i] = NUMERIC_POINTER(x)[i];
181 pt->wx1[i] = NUMERIC_POINTER(WX)[i];
182 }
183 pt->set = TRUE;
184 UNPROTECT(pc);
185
186 return;
187 }
188
189
hess_error_free(SEXP ptr)190 SEXP hess_error_free(SEXP ptr) {
191
192 HESS_ERROR_SSE *pt;
193
194 pt = (HESS_ERROR_SSE *) R_ExternalPtrAddr(ptr);
195
196 Free(pt->xlb);
197 Free(pt->beta1);
198 Free(pt->wx1);
199 Free(pt->xl);
200 Free(pt->wy1);
201 Free(pt->yl);
202 Free(pt->x);
203 Free(pt->y);
204
205 Free(pt);
206 R_ClearExternalPtr(ptr);
207 return(R_NilValue);
208 }
209
hess_lag_init()210 SEXP hess_lag_init() {
211
212 HESS_LAG_SSE *pt;
213 SEXP ptr;
214
215 pt = Calloc(1, HESS_LAG_SSE);
216 pt->set = FALSE;
217
218 PROTECT(ptr = R_MakeExternalPtr(pt, R_NilValue, R_NilValue));
219
220 UNPROTECT(1);
221 return(ptr);
222
223 }
224
225
hess_lag_set(SEXP env)226 void hess_lag_set(SEXP env) {
227
228 HESS_LAG_SSE *pt;
229 SEXP y, x, wy;
230 int i, n, p, np, pc=0;
231
232 n = INTEGER_POINTER(findVarInFrame(env, install("n")))[0];
233 p = INTEGER_POINTER(findVarInFrame(env, install("m")))[0];
234 np = n*p;
235
236 pt = (HESS_LAG_SSE *) R_ExternalPtrAddr(findVarInFrame(env,
237 install("ptr")));
238 if (pt->set) error("hess_lag_set: function called out of order");
239
240 PROTECT(y = findVarInFrame(env, install("y"))); pc++;
241 PROTECT(x = findVarInFrame(env, install("x"))); pc++;
242 PROTECT(wy = findVarInFrame(env, install("wy"))); pc++;
243
244 pt->y = Calloc(n, double);
245 pt->x = Calloc(np, double);
246 pt->yl = Calloc(n, double);
247 pt->wy1 = Calloc(n, double);
248 pt->beta1 = Calloc(p, double);
249 pt->xb = Calloc(n, double);
250
251 for (i=0; i<n; i++) {
252 pt->y[i] = NUMERIC_POINTER(y)[i];
253 pt->wy1[i] = NUMERIC_POINTER(wy)[i];
254 }
255 for (i=0; i<np; i++) pt->x[i] = NUMERIC_POINTER(x)[i];
256 pt->set = TRUE;
257 UNPROTECT(pc);
258
259 return;
260 }
261
262
hess_lag_free(SEXP ptr)263 SEXP hess_lag_free(SEXP ptr) {
264
265 HESS_LAG_SSE *pt;
266
267 pt = (HESS_LAG_SSE *) R_ExternalPtrAddr(ptr);
268
269 Free(pt->xb);
270 Free(pt->beta1);
271 Free(pt->wy1);
272 Free(pt->yl);
273 Free(pt->x);
274 Free(pt->y);
275
276 Free(pt);
277 R_ClearExternalPtr(ptr);
278 return(R_NilValue);
279 }
280
281 /**
282 * Calculate the sum of squared errors term for spatial regression
283 * using an environment to hold data
284 *
285 * @param env pointer to an SEXP environment
286 * @param coef current value of coefficient being optimzed
287 *
288 * @return double, value of SSE for current coef
289 *
290 */
R_ml_sse_env(SEXP env,SEXP coef)291 SEXP R_ml_sse_env(SEXP env, SEXP coef) {
292
293 SEXP res;
294 // SEXP y, x, wy, WX;
295 int i, k, n, p, np;
296 double tol=1e-7, cyl, cxlqyl, sse;
297 char *trans = "T";
298 double one = 1.0, zero = 0.0;
299 double m_lambda = - NUMERIC_POINTER(coef)[0];
300 int pc=0, first_time;
301 OPT_ERROR_SSE *pt;
302
303 first_time = LOGICAL_POINTER(findVarInFrame(env, install("first_time")))[0];
304 if (first_time) {
305 opt_error_set(env);
306 }
307
308 n = INTEGER_POINTER(findVarInFrame(env, install("n")))[0];
309 p = INTEGER_POINTER(findVarInFrame(env, install("p")))[0];
310 np = n*p;
311 pt = (OPT_ERROR_SSE *) R_ExternalPtrAddr(findVarInFrame(env,
312 install("ptr")));
313
314 for (i=0; i<n; i++) pt->yl[i] = pt->y[i];
315 for (i=0; i<np; i++) pt->xlq[i] = pt->x[i];
316
317 F77_CALL(daxpy)(&n, &m_lambda, pt->wy1, &c__1, pt->yl, &c__1);
318
319 F77_CALL(daxpy)(&np, &m_lambda, pt->wx1, &c__1, pt->xlq, &c__1);
320
321 F77_CALL(dqrdc2)(pt->xlq, &n, &n, &p, &tol, &k, pt->qraux, pt->jpvt,
322 pt->work);
323 if (p != k) warning("Q looses full rank");
324 /* k = 0;
325 F77_CALL(dqrdc)(pt->xlq, &n, &n, &p, pt->qraux, pt->jpvt, pt->work, &k);*/
326
327 for (i=0; i<n*k; i++) pt->qy[i] = 0.0;
328 for (i=0; i<k; i++) pt->qy[(i +(n*i))] = 1.0;
329
330 F77_CALL(dqrqy)(pt->xlq, &n, &k, pt->qraux, pt->qy, &k, pt->qy);
331
332 F77_CALL(dgemv)(trans, &n, &k, &one, pt->qy, &n, pt->yl, &c__1, &zero,
333 pt->xlqyl, &c__1 FCONE);
334
335 cyl = F77_CALL(ddot)(&n, pt->yl, &c__1, pt->yl, &c__1);
336
337 cxlqyl = F77_CALL(ddot)(&k, pt->xlqyl, &c__1, pt->xlqyl, &c__1);
338
339 sse = cyl - cxlqyl;
340
341 PROTECT(res=NEW_NUMERIC(1)); pc++;
342 NUMERIC_POINTER(res)[0] = sse;
343 UNPROTECT(pc);
344
345 return(res);
346
347 }
348
R_ml1_sse_env(SEXP env,SEXP lambda,SEXP beta)349 SEXP R_ml1_sse_env(SEXP env, SEXP lambda, SEXP beta) {
350
351 SEXP res;
352 int i, n, p, np;
353 double sse;
354 char *trans = "N";
355 double one = 1.0, zero = 0.0, m_one = -1.0;
356 double m_lambda = - NUMERIC_POINTER(lambda)[0];
357 int pc=0, first_time;
358 HESS_ERROR_SSE *pt;
359
360 first_time = LOGICAL_POINTER(findVarInFrame(env, install("first_time")))[0];
361 if (first_time) {
362 hess_error_set(env);
363 }
364
365 n = INTEGER_POINTER(findVarInFrame(env, install("n")))[0];
366 p = INTEGER_POINTER(findVarInFrame(env, install("p")))[0];
367 np = n*p;
368 pt = (HESS_ERROR_SSE *) R_ExternalPtrAddr(findVarInFrame(env,
369 install("ptr")));
370
371 for (i=0; i<n; i++) pt->yl[i] = pt->y[i];
372 for (i=0; i<np; i++) pt->xl[i] = pt->x[i];
373
374 for (i=0; i<p; i++) pt->beta1[i] = NUMERIC_POINTER(beta)[i];
375
376 F77_CALL(daxpy)(&n, &m_lambda, pt->wy1, &c__1, pt->yl, &c__1);
377
378 F77_CALL(daxpy)(&np, &m_lambda, pt->wx1, &c__1, pt->xl, &c__1);
379
380 F77_CALL(dgemv)(trans, &n, &p, &one, pt->xl, &n, pt->beta1, &c__1, &zero,
381 pt->xlb, &c__1 FCONE);
382
383 F77_CALL(daxpy)(&n, &m_one, pt->xlb, &c__1, pt->yl, &c__1);
384
385 sse = F77_CALL(ddot)(&n, pt->yl, &c__1, pt->yl, &c__1);
386
387 PROTECT(res=NEW_NUMERIC(1)); pc++;
388 NUMERIC_POINTER(res)[0] = sse;
389 UNPROTECT(pc);
390
391 return(res);
392
393 }
394
R_ml2_sse_env(SEXP env,SEXP rho,SEXP beta)395 SEXP R_ml2_sse_env(SEXP env, SEXP rho, SEXP beta) {
396
397 SEXP res;
398 int i, n, p;
399 double sse;
400 char *trans = "N";
401 double one = 1.0, zero = 0.0, m_one = -1.0;
402 double m_rho1 = - NUMERIC_POINTER(rho)[0];
403 int pc=0, first_time;
404 HESS_LAG_SSE *pt;
405
406 first_time = LOGICAL_POINTER(findVarInFrame(env, install("first_time")))[0];
407 if (first_time) {
408 hess_lag_set(env);
409 }
410
411 n = INTEGER_POINTER(findVarInFrame(env, install("n")))[0];
412 p = INTEGER_POINTER(findVarInFrame(env, install("m")))[0];
413 pt = (HESS_LAG_SSE *) R_ExternalPtrAddr(findVarInFrame(env,
414 install("ptr")));
415
416 for (i=0; i<n; i++) pt->yl[i] = pt->y[i];
417 for (i=0; i<p; i++) pt->beta1[i] = NUMERIC_POINTER(beta)[i];
418
419 F77_CALL(daxpy)(&n, &m_rho1, pt->wy1, &c__1, pt->yl, &c__1);
420
421 F77_CALL(dgemv)(trans, &n, &p, &one, pt->x, &n, pt->beta1, &c__1, &zero,
422 pt->xb, &c__1 FCONE);
423
424 F77_CALL(daxpy)(&n, &m_one, pt->xb, &c__1, pt->yl, &c__1);
425
426 sse = F77_CALL(ddot)(&n, pt->yl, &c__1, pt->yl, &c__1);
427
428 PROTECT(res=NEW_NUMERIC(1)); pc++;
429 NUMERIC_POINTER(res)[0] = sse;
430 UNPROTECT(pc);
431
432 return(res);
433
434 }
435
436