1
2 #include <R.h>
3 #include <Rinternals.h>
4
5
6 // // //
7 // // Compute M3 %*% (U %x% U) or M4 %*% (U %x% U %x% U)
8 // // //
9
M3HOOIiterator(SEXP XX,SEXP UU,SEXP PP,SEXP KK)10 SEXP M3HOOIiterator(SEXP XX, SEXP UU, SEXP PP, SEXP KK){
11 /*
12 arguments
13 XX : numeric vector with unique elements of a coskewness matrix
14 UU : numeric vector (vectorized matrix) with the loadings matrix
15 PP : integer, dimension of XX
16 KK : integer, number of columns of UU
17
18 computes Phi %*% (U %x% U)
19
20 Written by Dries Cornilly
21 */
22
23 // declare pointers for the vectors
24 double *X, *U;
25
26 // use REAL() to access the C array inside the numeric vector passed in from R
27 X = REAL(XX);
28 U = REAL(UU);
29 int P = asInteger(PP);
30 int K = asInteger(KK);
31
32 // allocate and compute the result matrix
33 SEXP Z3 = PROTECT(allocMatrix(REALSXP, P, K * K));
34 double *rZ3 = REAL(Z3);
35 for (int ii = 0; ii < P * K * K; ii++) {
36 rZ3[ii] = 0.0;
37 }
38
39 // loop over the coskenwess elements
40 int iter = 0;
41 for (int ii = 0; ii < P; ii++) {
42 for (int jj = ii; jj < P; jj++) {
43 for (int kk = jj; kk < P; kk++) {
44 if (ii == jj) {
45 if (jj == kk) {
46 // phi_iii
47 for (int uu = 0; uu < K; uu++) {
48 for (int vv = 0; vv < K; vv++) {
49 rZ3[(vv * K + uu) * P + ii] += X[iter] * U[uu * P + ii] * U[vv * P + ii];
50 }
51 }
52 } else {
53 // phi_iik
54 for (int uu = 0; uu < K; uu++) {
55 int uuP = uu * P;
56 for (int vv = 0; vv < K; vv++) {
57 int vvP = vv * P;
58 rZ3[(vv * K + uu) * P + ii] += X[iter] * (U[uuP + ii] * U[vvP + kk] + U[uuP + kk] * U[vvP + ii]);
59 rZ3[(vv * K + uu) * P + kk] += X[iter] * U[uuP + ii] * U[vvP + ii];
60 }
61 }
62 }
63 } else {
64 if (jj == kk) {
65 // phi_ijj
66 for (int uu = 0; uu < K; uu++) {
67 int uuP = uu * P;
68 for (int vv = 0; vv < K; vv++) {
69 int vvP = vv * P;
70 rZ3[(vv * K + uu) * P + ii] += X[iter] * U[uuP + jj] * U[vvP + jj];
71 rZ3[(vv * K + uu) * P + jj] += X[iter] * (U[uuP + ii] * U[vvP + jj] + U[uuP + jj] * U[vvP + ii]);
72 }
73 }
74 } else {
75 //phi_ijk
76 for (int uu = 0; uu < K; uu++) {
77 int uuP = uu * P;
78 for (int vv = 0; vv < K; vv++) {
79 int vvP = vv * P;
80 int vvKuuP = (vv * K + uu) * P;
81 rZ3[vvKuuP + ii] += X[iter] * (U[uuP + jj] * U[vvP + kk] + U[uuP + kk] * U[vvP + jj]);
82 rZ3[vvKuuP + jj] += X[iter] * (U[uuP + ii] * U[vvP + kk] + U[uuP + kk] * U[vvP + ii]);
83 rZ3[vvKuuP + kk] += X[iter] * (U[uuP + ii] * U[vvP + jj] + U[uuP + jj] * U[vvP + ii]);
84 }
85 }
86 }
87 }
88 iter++;
89 } // loop kk
90 } // loop jj
91 } // loop ii
92
93 UNPROTECT(1);
94 return Z3;
95 }
96
97
M4HOOIiterator(SEXP XX,SEXP UU,SEXP PP,SEXP KK)98 SEXP M4HOOIiterator(SEXP XX, SEXP UU, SEXP PP, SEXP KK){
99 /*
100 arguments
101 XX : numeric vector with unique elements of a cokurtosis matrix
102 UU : numeric vector (vectorized matrix) with the loadings matrix
103 PP : integer, dimension of XX
104 KK : integer, number of columns of UU
105
106 computes Psi %*% (U %x% U %x% U)
107
108 Written by Dries Cornilly
109 */
110
111 // declare pointers for the vectors
112 double *X, *U;
113
114 // use REAL() to access the C array inside the numeric vector passed in from R
115 X = REAL(XX);
116 U = REAL(UU);
117 int P = asInteger(PP);
118 int K = asInteger(KK);
119 int K2 = K * K;
120
121 // allocate and compute the result matrix
122 SEXP Z4 = PROTECT(allocMatrix(REALSXP, P, K * K * K));
123 double *rZ4 = REAL(Z4);
124 for (int ii = 0; ii < P * K * K * K; ii++) {
125 rZ4[ii] = 0.0;
126 }
127
128 // loop over the cokurtosis elements
129 int iter = 0;
130 for (int ii = 0; ii < P; ii++) {
131 for (int jj = ii; jj < P; jj++) {
132 for (int kk = jj; kk < P; kk++) {
133 for (int ll = kk; ll < P; ll++) {
134 if (ii == jj) {
135 if (jj == kk) {
136 if (kk == ll) {
137 // psi_iiii
138 for (int uu = 0; uu < K; uu++) {
139 for (int vv = 0; vv < K; vv++) {
140 for (int ww = 0; ww < K; ww++) {
141 rZ4[(ww * K2 + vv * K + uu) * P + ii] += X[iter] * U[uu * P + ii] * U[vv * P + ii] * U[ww * P + ii];
142 }
143 }
144 }
145 } else {
146 // psi_iiil
147 for (int uu = 0; uu < K; uu++) {
148 int uuP = uu * P;
149 for (int vv = 0; vv < K; vv++) {
150 int vvP = vv * P;
151 for (int ww = 0; ww < K; ww++) {
152 int wwP = ww * P;
153 rZ4[(ww * K2 + vv * K + uu) * P + ii] += X[iter] * (U[uuP + ii] * U[vvP + ii] * U[wwP + ll] +
154 U[uuP + ii] * U[vvP + ll] * U[wwP + ii] + U[uuP + ll] * U[vvP + ii] * U[wwP + ii]);
155 rZ4[(ww * K2 + vv * K + uu) * P + ll] += X[iter] * U[uuP + ii] * U[vvP + ii] * U[wwP + ii];
156 }
157 }
158 }
159 }
160 } else {
161 if (kk == ll) {
162 // psi_iikk
163 for (int uu = 0; uu < K; uu++) {
164 int uuP = uu * P;
165 for (int vv = 0; vv < K; vv++) {
166 int vvP = vv * P;
167 for (int ww = 0; ww < K; ww++) {
168 int wwP = ww * P;
169 rZ4[(ww * K2 + vv * K + uu) * P + ii] += X[iter] * (U[uuP + ii] * U[vvP + kk] * U[wwP + kk] +
170 U[uuP + kk] * U[vvP + ii] * U[wwP + kk] + U[uuP + kk] * U[vvP + kk] * U[wwP + ii]);
171 rZ4[(ww * K2 + vv * K + uu) * P + kk] += X[iter] * (U[uuP + ii] * U[vvP + ii] * U[wwP + kk] +
172 U[uuP + ii] * U[vvP + kk] * U[wwP + ii] + U[uuP + kk] * U[vvP + ii] * U[wwP + ii]);
173 }
174 }
175 }
176 } else {
177 // psi_iikl
178 for (int uu = 0; uu < K; uu++) {
179 int uuP = uu * P;
180 for (int vv = 0; vv < K; vv++) {
181 int vvP = vv * P;
182 for (int ww = 0; ww < K; ww++) {
183 int wwP = ww * P;
184 rZ4[(ww * K2 + vv * K + uu) * P + ii] += X[iter] * (U[uuP + ii] * U[vvP + kk] * U[wwP + ll] +
185 U[uuP + ii] * U[vvP + ll] * U[wwP + kk] + U[uuP + kk] * U[vvP + ii] * U[wwP + ll] +
186 U[uuP + kk] * U[vvP + ll] * U[wwP + ii] + U[uuP + ll] * U[vvP + ii] * U[wwP + kk] +
187 U[uuP + ll] * U[vvP + kk] * U[wwP + ii]);
188 rZ4[(ww * K2 + vv * K + uu) * P + kk] += X[iter] * (U[uuP + ii] * U[vvP + ii] * U[wwP + ll] +
189 U[uuP + ii] * U[vvP + ll] * U[wwP + ii] + U[uuP + ll] * U[vvP + ii] * U[wwP + ii]);
190 rZ4[(ww * K2 + vv * K + uu) * P + ll] += X[iter] * (U[uuP + ii] * U[vvP + ii] * U[wwP + kk] +
191 U[uuP + ii] * U[vvP + kk] * U[wwP + ii] + U[uuP + kk] * U[vvP + ii] * U[wwP + ii]);
192 }
193 }
194 }
195 }
196 }
197 } else {
198 if (jj == kk) {
199 if (kk == ll) {
200 // psi_ijjj
201 for (int uu = 0; uu < K; uu++) {
202 int uuP = uu * P;
203 for (int vv = 0; vv < K; vv++) {
204 int vvP = vv * P;
205 for (int ww = 0; ww < K; ww++) {
206 int wwP = ww * P;
207 rZ4[(ww * K2 + vv * K + uu) * P + ii] += X[iter] * U[uuP + jj] * U[vvP + jj] * U[wwP + jj];
208 rZ4[(ww * K2 + vv * K + uu) * P + jj] += X[iter] * (U[uuP + ii] * U[vvP + jj] * U[wwP + jj] +
209 U[uuP + jj] * U[vvP + ii] * U[wwP + jj] + U[uuP + jj] * U[vvP + jj] * U[wwP + ii]);
210 }
211 }
212 }
213 } else {
214 // psi_ijjl
215 for (int uu = 0; uu < K; uu++) {
216 int uuP = uu * P;
217 for (int vv = 0; vv < K; vv++) {
218 int vvP = vv * P;
219 for (int ww = 0; ww < K; ww++) {
220 int wwP = ww * P;
221 rZ4[(ww * K2 + vv * K + uu) * P + ii] += X[iter] * (U[uuP + jj] * U[vvP + jj] * U[wwP + ll] +
222 U[uuP + jj] * U[vvP + ll] * U[wwP + jj] + U[uuP + ll] * U[vvP + jj] * U[wwP + jj]);
223 rZ4[(ww * K2 + vv * K + uu) * P + jj] += X[iter] * (U[uuP + ii] * U[vvP + jj] * U[wwP + ll] +
224 U[uuP + ii] * U[vvP + ll] * U[wwP + jj] + U[uuP + jj] * U[vvP + ii] * U[wwP + ll] +
225 U[uuP + jj] * U[vvP + ll] * U[wwP + ii] + U[uuP + ll] * U[vvP + ii] * U[wwP + jj] +
226 U[uuP + ll] * U[vvP + jj] * U[wwP + ii]);
227 rZ4[(ww * K2 + vv * K + uu) * P + ll] += X[iter] * (U[uuP + ii] * U[vvP + jj] * U[wwP + jj] +
228 U[uuP + jj] * U[vvP + ii] * U[wwP + jj] + U[uuP + jj] * U[vvP + jj] * U[wwP + ii]);
229 }
230 }
231 }
232 }
233 } else {
234 if (kk == ll) {
235 // psi_ijkk
236 for (int uu = 0; uu < K; uu++) {
237 int uuP = uu * P;
238 for (int vv = 0; vv < K; vv++) {
239 int vvP = vv * P;
240 for (int ww = 0; ww < K; ww++) {
241 int wwP = ww * P;
242 rZ4[(ww * K2 + vv * K + uu) * P + ii] += X[iter] * (U[uuP + jj] * U[vvP + kk] * U[wwP + kk] +
243 U[uuP + kk] * U[vvP + jj] * U[wwP + kk] + U[uuP + kk] * U[vvP + kk] * U[wwP + jj]);
244 rZ4[(ww * K2 + vv * K + uu) * P + jj] += X[iter] * (U[uuP + ii] * U[vvP + kk] * U[wwP + kk] +
245 U[uuP + kk] * U[vvP + ii] * U[wwP + kk] + U[uuP + kk] * U[vvP + kk] * U[wwP + ii]);
246 rZ4[(ww * K2 + vv * K + uu) * P + kk] += X[iter] * (U[uuP + ii] * U[vvP + jj] * U[wwP + kk] +
247 U[uuP + ii] * U[vvP + kk] * U[wwP + jj] + U[uuP + jj] * U[vvP + ii] * U[wwP + kk] +
248 U[uuP + jj] * U[vvP + kk] * U[wwP + ii] + U[uuP + kk] * U[vvP + ii] * U[wwP + jj] +
249 U[uuP + kk] * U[vvP + jj] * U[wwP + ii]);
250 }
251 }
252 }
253 } else {
254 // psi_ijkl
255 for (int uu = 0; uu < K; uu++) {
256 int uuP = uu * P;
257 for (int vv = 0; vv < K; vv++) {
258 int vvP = vv * P;
259 for (int ww = 0; ww < K; ww++) {
260 int wwP = ww * P;
261 rZ4[(ww * K2 + vv * K + uu) * P + ii] += X[iter] * (U[uuP + jj] * U[vvP + kk] * U[wwP + ll] +
262 U[uuP + jj] * U[vvP + ll] * U[wwP + kk] + U[uuP + kk] * U[vvP + jj] * U[wwP + ll] +
263 U[uuP + kk] * U[vvP + ll] * U[wwP + jj] + U[uuP + ll] * U[vvP + jj] * U[wwP + kk] +
264 U[uuP + ll] * U[vvP + kk] * U[wwP + jj]);
265 rZ4[(ww * K2 + vv * K + uu) * P + jj] += X[iter] * (U[uuP + ii] * U[vvP + kk] * U[wwP + ll] +
266 U[uuP + ii] * U[vvP + ll] * U[wwP + kk] + U[uuP + kk] * U[vvP + ii] * U[wwP + ll] +
267 U[uuP + kk] * U[vvP + ll] * U[wwP + ii] + U[uuP + ll] * U[vvP + ii] * U[wwP + kk] +
268 U[uuP + ll] * U[vvP + kk] * U[wwP + ii]);
269 rZ4[(ww * K2 + vv * K + uu) * P + kk] += X[iter] * (U[uuP + ii] * U[vvP + jj] * U[wwP + ll] +
270 U[uuP + ii] * U[vvP + ll] * U[wwP + jj] + U[uuP + jj] * U[vvP + ii] * U[wwP + ll] +
271 U[uuP + jj] * U[vvP + ll] * U[wwP + ii] + U[uuP + ll] * U[vvP + ii] * U[wwP + jj] +
272 U[uuP + ll] * U[vvP + jj] * U[wwP + ii]);
273 rZ4[(ww * K2 + vv * K + uu) * P + ll] += X[iter] * (U[uuP + ii] * U[vvP + jj] * U[wwP + kk] +
274 U[uuP + ii] * U[vvP + kk] * U[wwP + jj] + U[uuP + jj] * U[vvP + ii] * U[wwP + kk] +
275 U[uuP + jj] * U[vvP + kk] * U[wwP + ii] + U[uuP + kk] * U[vvP + ii] * U[wwP + jj] +
276 U[uuP + kk] * U[vvP + jj] * U[wwP + ii]);
277 }
278 }
279 }
280 }
281 }
282 }
283 iter++;
284 } // loop ll
285 } // loop kk
286 } // loop jj
287 } // loop ii
288
289 UNPROTECT(1);
290 return Z4;
291 }
292