1 /* Copyright 2015 by Roger S. Bivand. */
2
3 #include "spdep.h"
4
5
lmin21(SEXP nb,SEXP y,SEXP cy,SEXP card)6 SEXP lmin21(SEXP nb, SEXP y, SEXP cy, SEXP card) {
7 int i, j, k, nswitch=0, n=length(card), pc=0;
8 SEXP ans;
9 double t1, t2, ytemp;
10 double *Y, *CY;
11
12 Y = (double *) R_alloc((size_t) n, sizeof(double));
13 CY = (double *) R_alloc((size_t) n, sizeof(double));
14
15 for (i=0; i<n; i++) {
16 Y[i] = NUMERIC_POINTER(y)[i];
17 CY[i] = NUMERIC_POINTER(cy)[i];
18 }
19
20 PROTECT(ans = NEW_LIST(2)); pc++;
21 SET_VECTOR_ELT(ans, 0, NEW_NUMERIC(n));
22 SET_VECTOR_ELT(ans, 1, NEW_INTEGER(1));
23
24 for (i=0; i<n; i++) {
25 if (INTEGER_POINTER(card)[i] > 0) {
26 t1 = fabs(Y[i] - CY[i]);
27 t2 = fabs(-2*CY[i]);
28 for (j=0; j<INTEGER_POINTER(card)[i]; j++) {
29 k = INTEGER_POINTER(VECTOR_ELT(nb, i))[j]-ROFFSET;
30 t1 = t1 + fabs(Y[k] - CY[k]);
31 t2 = t2 + fabs(Y[k] - (CY[k] - Y[i] - CY[i]));
32 }
33 if (t1 <= t2) {
34 nswitch++;
35 ytemp = Y[i];
36 Y[i] = -CY[i];
37 for (j=0; j<INTEGER_POINTER(card)[i]; j++) {
38 k = INTEGER_POINTER(VECTOR_ELT(nb, i))[j]-ROFFSET;
39 CY[k] = CY[k] - ytemp + Y[i];
40 }
41 }
42 }
43 }
44
45 for (i=0; i<n; i++) {
46 NUMERIC_POINTER(VECTOR_ELT(ans, 0))[i] = Y[i];
47 }
48
49 INTEGER_POINTER(VECTOR_ELT(ans, 1))[0] = nswitch;
50 UNPROTECT(pc); /* ans */
51 return(ans);
52 }
53
lmin22(SEXP nb,SEXP y,SEXP cy,SEXP card,SEXP beta)54 SEXP lmin22(SEXP nb, SEXP y, SEXP cy, SEXP card, SEXP beta) {
55 int i, j, k, nswitch=0, n=length(card), pc=0;
56 SEXP ans;
57 double t1, t2, ytemp, yhat;
58 double *Y, *CY, *B;
59
60 Y = (double *) R_alloc((size_t) n, sizeof(double));
61 CY = (double *) R_alloc((size_t) n, sizeof(double));
62 B = (double *) R_alloc((size_t) length(beta), sizeof(double));
63
64 for (i=0; i<n; i++) {
65 Y[i] = NUMERIC_POINTER(y)[i];
66 CY[i] = NUMERIC_POINTER(cy)[i];
67 }
68 for (i=0; i<length(beta); i++) {
69 B[i] = NUMERIC_POINTER(beta)[i];
70 }
71
72 PROTECT(ans = NEW_LIST(2)); pc++;
73 SET_VECTOR_ELT(ans, 0, NEW_NUMERIC(n));
74 SET_VECTOR_ELT(ans, 1, NEW_INTEGER(1));
75
76 for (i=0; i<n; i++) {
77 if (INTEGER_POINTER(card)[i] > 0) {
78 t1 = fabs(Y[i] - CY[i]);
79 yhat = B[0] + B[1]*CY[i];
80 t2 = fabs(yhat - CY[i]);
81 for (j=0; j<INTEGER_POINTER(card)[i]; j++) {
82 k = INTEGER_POINTER(VECTOR_ELT(nb, i))[j]-ROFFSET;
83 t1 = t1 + fabs(Y[k] - CY[k]);
84 t2 = t2 + fabs(Y[k] - (CY[k] - Y[i] + B[0] + B[1]*CY[i]));
85 }
86 if (t1 <= t2) {
87 nswitch++;
88 ytemp = Y[i];
89 Y[i] = yhat;
90 for (j=0; j<INTEGER_POINTER(card)[i]; j++) {
91 k = INTEGER_POINTER(VECTOR_ELT(nb, i))[j]-ROFFSET;
92 CY[k] = CY[k] - ytemp + Y[i];
93 }
94 }
95 }
96 }
97
98 for (i=0; i<n; i++) {
99 NUMERIC_POINTER(VECTOR_ELT(ans, 0))[i] = Y[i];
100 }
101
102 INTEGER_POINTER(VECTOR_ELT(ans, 1))[0] = nswitch;
103 UNPROTECT(pc); /* ans */
104 return(ans);
105 }
106
lmin23(SEXP nb,SEXP y,SEXP cy,SEXP card,SEXP beta,SEXP tol)107 SEXP lmin23(SEXP nb, SEXP y, SEXP cy, SEXP card, SEXP beta, SEXP tol) {
108 int i, j, k, nswitch=0, n=length(card), pc=0;
109 SEXP ans;
110 double tmp, var, yhat;
111 double *Y, *CY, *B;
112
113 Y = (double *) R_alloc((size_t) n, sizeof(double));
114 CY = (double *) R_alloc((size_t) n, sizeof(double));
115 B = (double *) R_alloc((size_t) length(beta), sizeof(double));
116
117 for (i=0; i<n; i++) {
118 Y[i] = NUMERIC_POINTER(y)[i];
119 CY[i] = NUMERIC_POINTER(cy)[i];
120 }
121 for (i=0; i<length(beta); i++) {
122 B[i] = NUMERIC_POINTER(beta)[i];
123 }
124 PROTECT(ans = NEW_LIST(2)); pc++;
125 SET_VECTOR_ELT(ans, 0, NEW_NUMERIC(n));
126 SET_VECTOR_ELT(ans, 1, NEW_INTEGER(1));
127
128 for (i=0; i<n; i++) {
129 if (INTEGER_POINTER(card)[i] > 0) {
130 yhat = B[0] + B[1]*CY[i];
131 var = fabs(Y[i] - yhat);
132 if (var > NUMERIC_POINTER(tol)[0]) {
133 nswitch++;
134 tmp = Y[i];
135 Y[i] = yhat;
136 for (j=0; j<INTEGER_POINTER(card)[i]; j++) {
137 k = INTEGER_POINTER(VECTOR_ELT(nb, i))[j]-ROFFSET;
138 CY[k] = CY[k] - tmp + Y[i];
139 }
140 }
141 }
142 }
143
144 for (i=0; i<n; i++) {
145 NUMERIC_POINTER(VECTOR_ELT(ans, 0))[i] = Y[i];
146 }
147
148 INTEGER_POINTER(VECTOR_ELT(ans, 1))[0] = nswitch;
149 UNPROTECT(pc); /* ans */
150 return(ans);
151 }
152
lmin3(SEXP nb,SEXP ev1,SEXP ev1_lag,SEXP n_nei,SEXP beta,SEXP tol)153 SEXP lmin3(SEXP nb, SEXP ev1, SEXP ev1_lag, SEXP n_nei, SEXP beta, SEXP tol) {
154 int i, j, k, nswitch=0, n=length(n_nei), pc=0;
155 SEXP ans;
156 double tmp, var, yhat, ntmp;
157 double *Y, *CY, *B;
158
159 Y = (double *) R_alloc((size_t) n, sizeof(double));
160 CY = (double *) R_alloc((size_t) n, sizeof(double));
161 B = (double *) R_alloc((size_t) length(beta), sizeof(double));
162
163 for (i=0; i<n; i++) {
164 Y[i] = NUMERIC_POINTER(ev1)[i];
165 CY[i] = NUMERIC_POINTER(ev1_lag)[i];
166 }
167 for (i=0; i<length(beta); i++) {
168 B[i] = NUMERIC_POINTER(beta)[i];
169 }
170 PROTECT(ans = NEW_LIST(2)); pc++;
171 SET_VECTOR_ELT(ans, 0, NEW_NUMERIC(n));
172 SET_VECTOR_ELT(ans, 1, NEW_INTEGER(1));
173
174 for (i=0; i<n; i++) {
175 if (INTEGER_POINTER(n_nei)[i] > 0) {
176 yhat = B[0] + B[1]*CY[i];
177 var = fabs(Y[i] - yhat);
178 if (var > NUMERIC_POINTER(tol)[0]) {
179 nswitch++;
180 tmp = Y[i];
181 Y[i] = yhat;
182 for (j=0; j<INTEGER_POINTER(n_nei)[i]; j++) {
183 k = INTEGER_POINTER(VECTOR_ELT(nb, i))[j]-ROFFSET;
184 ntmp = sqrt(INTEGER_POINTER(n_nei)[i] *
185 INTEGER_POINTER(n_nei)[k]);
186 CY[k] = CY[k] - (tmp/ntmp) + (Y[i]/ntmp);
187 }
188 }
189 }
190 }
191
192 for (i=0; i<n; i++) {
193 NUMERIC_POINTER(VECTOR_ELT(ans, 0))[i] = Y[i];
194 }
195
196 INTEGER_POINTER(VECTOR_ELT(ans, 1))[0] = nswitch;
197 UNPROTECT(pc); /* ans */
198 return(ans);
199 }
200
201
lmin3S(SEXP nb,SEXP ev1,SEXP ev1_lag,SEXP n_nei,SEXP card,SEXP beta,SEXP tol)202 SEXP lmin3S(SEXP nb, SEXP ev1, SEXP ev1_lag, SEXP n_nei, SEXP card, SEXP beta, SEXP tol) {
203 int i, j, k, nswitch=0, n=length(card), pc=0;
204 SEXP ans;
205 double tmp, var, yhat, ntmp;
206 double *Y, *CY, *B;
207
208 Y = (double *) R_alloc((size_t) n, sizeof(double));
209 CY = (double *) R_alloc((size_t) n, sizeof(double));
210 B = (double *) R_alloc((size_t) length(beta), sizeof(double));
211
212 for (i=0; i<n; i++) {
213 Y[i] = NUMERIC_POINTER(ev1)[i];
214 CY[i] = NUMERIC_POINTER(ev1_lag)[i];
215 }
216 for (i=0; i<length(beta); i++) {
217 B[i] = NUMERIC_POINTER(beta)[i];
218 }
219 PROTECT(ans = NEW_LIST(2)); pc++;
220 SET_VECTOR_ELT(ans, 0, NEW_NUMERIC(n));
221 SET_VECTOR_ELT(ans, 1, NEW_INTEGER(1));
222
223 for (i=0; i<n; i++) {
224 if (INTEGER_POINTER(card)[i] > 0) {
225 yhat = B[0] + B[1]*CY[i];
226 var = fabs(Y[i] - yhat);
227 if (var > NUMERIC_POINTER(tol)[0]) {
228 nswitch++;
229 tmp = Y[i];
230 Y[i] = yhat;
231 for (j=0; j<INTEGER_POINTER(card)[i]; j++) {
232 k = INTEGER_POINTER(VECTOR_ELT(nb, i))[j]-ROFFSET;
233 ntmp = sqrt(NUMERIC_POINTER(n_nei)[i] *
234 NUMERIC_POINTER(n_nei)[k]);
235 CY[k] = CY[k] - (tmp/ntmp) + (Y[i]/ntmp);
236 }
237 }
238 }
239 }
240
241 for (i=0; i<n; i++) {
242 NUMERIC_POINTER(VECTOR_ELT(ans, 0))[i] = Y[i];
243 }
244
245 INTEGER_POINTER(VECTOR_ELT(ans, 1))[0] = nswitch;
246 UNPROTECT(pc); /* ans */
247 return(ans);
248 }
249