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