GRASS GIS 8 Programmer's Manual 8.2.1RC1(2022)-exported
c_reg.c
Go to the documentation of this file.
1#include <math.h>
2
3#include <grass/gis.h>
4#include <grass/raster.h>
5
6enum {
11 REGRESSION_P = 4
12};
13
14static 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) {
66 *result = numer / denom;
67 break;
69 *result = ybar - xbar * numer / denom;
70 break;
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
87void c_reg_m(DCELL * result, DCELL * values, int n, const void *closure)
88{
89 regression(result, values, n, REGRESSION_SLOPE);
90}
91
92void c_reg_c(DCELL * result, DCELL * values, int n, const void *closure)
93{
94 regression(result, values, n, REGRESSION_OFFSET);
95}
96
97void c_reg_r2(DCELL * result, DCELL * values, int n, const void *closure)
98{
99 regression(result, values, n, REGRESSION_COEFF_DET);
100}
101
102void c_reg_t(DCELL * result, DCELL * values, int n, const void *closure)
103{
104 regression(result, values, n, REGRESSION_T);
105}
106
107static 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;
163 *result = ybar - xbar * numer / denom;
164 break;
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
181void w_reg_m(DCELL * result, DCELL(*values)[2], int n, const void *closure)
182{
183 regression_w(result, values, n, REGRESSION_SLOPE);
184}
185
186void w_reg_c(DCELL * result, DCELL(*values)[2], int n, const void *closure)
187{
188 regression_w(result, values, n, REGRESSION_OFFSET);
189}
190
191void 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
196void w_reg_t(DCELL * result, DCELL(*values)[2], int n, const void *closure)
197{
198 regression_w(result, values, n, REGRESSION_T);
199}
void w_reg_m(DCELL *result, DCELL(*values)[2], int n, const void *closure)
Definition: c_reg.c:181
void w_reg_c(DCELL *result, DCELL(*values)[2], int n, const void *closure)
Definition: c_reg.c:186
void c_reg_m(DCELL *result, DCELL *values, int n, const void *closure)
Definition: c_reg.c:87
void w_reg_r2(DCELL *result, DCELL(*values)[2], int n, const void *closure)
Definition: c_reg.c:191
void c_reg_t(DCELL *result, DCELL *values, int n, const void *closure)
Definition: c_reg.c:102
@ REGRESSION_P
Definition: c_reg.c:11
@ REGRESSION_COEFF_DET
Definition: c_reg.c:9
@ REGRESSION_OFFSET
Definition: c_reg.c:8
@ REGRESSION_T
Definition: c_reg.c:10
@ REGRESSION_SLOPE
Definition: c_reg.c:7
void c_reg_c(DCELL *result, DCELL *values, int n, const void *closure)
Definition: c_reg.c:92
void w_reg_t(DCELL *result, DCELL(*values)[2], int n, const void *closure)
Definition: c_reg.c:196
void c_reg_r2(DCELL *result, DCELL *values, int n, const void *closure)
Definition: c_reg.c:97
int count