1 #include <math.h>
2
3 #include <grass/gis.h>
4 #include <grass/raster.h>
5
6 enum {
7 REGRESSION_SLOPE = 0,
8 REGRESSION_OFFSET = 1,
9 REGRESSION_COEFF_DET = 2,
10 REGRESSION_T = 3,
11 REGRESSION_P = 4
12 };
13
regression(DCELL * result,DCELL * values,int n,int which)14 static void regression(DCELL * result, DCELL * values, int n, int which)
15 {
16 DCELL xsum, ysum;
17 DCELL xbar, ybar;
18 DCELL numer, denom, denom2;
19 DCELL Rsq;
20 int count;
21 int i;
22
23 xsum = ysum = 0.0;
24 count = 0;
25
26 for (i = 0; i < n; i++) {
27 if (Rast_is_d_null_value(&values[i]))
28 continue;
29
30 xsum += i;
31 ysum += values[i];
32 count++;
33 }
34
35 if (count < 2) {
36 Rast_set_d_null_value(result, 1);
37 return;
38 }
39
40 xbar = xsum / count;
41 ybar = ysum / count;
42
43 numer = 0.0;
44 for (i = 0; i < n; i++)
45 if (!Rast_is_d_null_value(&values[i]))
46 numer += i * values[i];
47 numer -= count * xbar * ybar;
48
49 denom = 0.0;
50 for (i = 0; i < n; i++)
51 if (!Rast_is_d_null_value(&values[i]))
52 denom += (DCELL) i * i;
53 denom -= count * xbar * xbar;
54
55 if (which >= REGRESSION_COEFF_DET || which == REGRESSION_T) {
56 denom2 = 0.0;
57 for (i = 0; i < n; i++)
58 if (!Rast_is_d_null_value(&values[i]))
59 denom2 += values[i] * values[i];
60 denom2 -= count * ybar * ybar;
61 Rsq = (numer * numer) / (denom * denom2);
62 }
63
64 switch (which) {
65 case REGRESSION_SLOPE:
66 *result = numer / denom;
67 break;
68 case REGRESSION_OFFSET:
69 *result = ybar - xbar * numer / denom;
70 break;
71 case REGRESSION_COEFF_DET:
72 *result = Rsq;
73 break;
74 case REGRESSION_T:
75 *result = sqrt(Rsq * (count - 2) / (1 - Rsq));
76 break;
77 default:
78 Rast_set_d_null_value(result, 1);
79 break;
80 }
81
82 /* Check for NaN */
83 if (*result != *result)
84 Rast_set_d_null_value(result, 1);
85 }
86
c_reg_m(DCELL * result,DCELL * values,int n,const void * closure)87 void c_reg_m(DCELL * result, DCELL * values, int n, const void *closure)
88 {
89 regression(result, values, n, REGRESSION_SLOPE);
90 }
91
c_reg_c(DCELL * result,DCELL * values,int n,const void * closure)92 void c_reg_c(DCELL * result, DCELL * values, int n, const void *closure)
93 {
94 regression(result, values, n, REGRESSION_OFFSET);
95 }
96
c_reg_r2(DCELL * result,DCELL * values,int n,const void * closure)97 void c_reg_r2(DCELL * result, DCELL * values, int n, const void *closure)
98 {
99 regression(result, values, n, REGRESSION_COEFF_DET);
100 }
101
c_reg_t(DCELL * result,DCELL * values,int n,const void * closure)102 void c_reg_t(DCELL * result, DCELL * values, int n, const void *closure)
103 {
104 regression(result, values, n, REGRESSION_T);
105 }
106
regression_w(DCELL * result,DCELL (* values)[2],int n,int which)107 static void regression_w(DCELL * result, DCELL(*values)[2], int n, int which)
108 {
109 DCELL xsum, ysum;
110 DCELL xbar, ybar;
111 DCELL numer, denom, denom2;
112 DCELL Rsq;
113 int count;
114 int i;
115
116 xsum = ysum = 0.0;
117 count = 0;
118
119 for (i = 0; i < n; i++) {
120 if (Rast_is_d_null_value(&values[i][0]))
121 continue;
122
123 xsum += i * values[i][1];
124 ysum += values[i][0] * values[i][1];
125 count += values[i][1];
126 }
127
128 if (count < 2) {
129 Rast_set_d_null_value(result, 1);
130 return;
131 }
132
133 xbar = xsum / count;
134 ybar = ysum / count;
135
136 numer = 0.0;
137 for (i = 0; i < n; i++)
138 if (!Rast_is_d_null_value(&values[i][0]))
139 numer += i * values[i][0] * values[i][1];
140 numer -= count * xbar * ybar;
141
142 denom = 0.0;
143 for (i = 0; i < n; i++)
144 if (!Rast_is_d_null_value(&values[i][0]))
145 denom += (DCELL) i * i * values[i][1];
146
147 denom -= count * xbar * xbar;
148
149 if (which == REGRESSION_COEFF_DET || which == REGRESSION_T) {
150 denom2 = 0.0;
151 for (i = 0; i < n; i++)
152 if (!Rast_is_d_null_value(&values[i][0]))
153 denom2 += values[i][0] * values[i][0] * values[i][1];
154 denom2 -= count * ybar * ybar;
155 Rsq = (numer * numer) / (denom * denom2);
156 }
157
158 switch (which) {
159 case REGRESSION_SLOPE:
160 *result = numer / denom;
161 break;
162 case REGRESSION_OFFSET:
163 *result = ybar - xbar * numer / denom;
164 break;
165 case REGRESSION_COEFF_DET:
166 *result = Rsq;
167 break;
168 case REGRESSION_T:
169 *result = sqrt(Rsq * (count - 2) / (1 - Rsq));
170 break;
171 default:
172 Rast_set_d_null_value(result, 1);
173 break;
174 }
175
176 /* Check for NaN */
177 if (*result != *result)
178 Rast_set_d_null_value(result, 1);
179 }
180
w_reg_m(DCELL * result,DCELL (* values)[2],int n,const void * closure)181 void w_reg_m(DCELL * result, DCELL(*values)[2], int n, const void *closure)
182 {
183 regression_w(result, values, n, REGRESSION_SLOPE);
184 }
185
w_reg_c(DCELL * result,DCELL (* values)[2],int n,const void * closure)186 void w_reg_c(DCELL * result, DCELL(*values)[2], int n, const void *closure)
187 {
188 regression_w(result, values, n, REGRESSION_OFFSET);
189 }
190
w_reg_r2(DCELL * result,DCELL (* values)[2],int n,const void * closure)191 void w_reg_r2(DCELL * result, DCELL(*values)[2], int n, const void *closure)
192 {
193 regression_w(result, values, n, REGRESSION_COEFF_DET);
194 }
195
w_reg_t(DCELL * result,DCELL (* values)[2],int n,const void * closure)196 void w_reg_t(DCELL * result, DCELL(*values)[2], int n, const void *closure)
197 {
198 regression_w(result, values, n, REGRESSION_T);
199 }
200