1 #define USE_FC_LEN_T 1
2 #include <stddef.h>
3 #include <R.h>
4 #include <math.h>
5 #include <stdio.h>
6 #include <stdlib.h>
7 #include <limits.h>
8
9 static int min_JM (int *, int *);
10 static int max_JM (int *, int *);
11
12 static void rowcentre_JM (double *, int, int);
13 static void colstandard_JM (double *, int, int);
14 static void rowstd_JM (double *, int, int, int);
15 static void transpose_mat_JM (double *, int *, int *, double *);
16 static void mmult_JM (double *, int, int, double *, int, int, double *);
17 static void orthog_mat_JM (double *, int, double *);
18 static void gramsch_JM (double *, int, int, int);
19 static void svd_JM (double *, int *, int *, double *, double *, double *);
20
21 static void Symm_logcosh_JM (double *, int, double *, int, int, double, double *, double *);
22 static void Symm_exp_JM (double *, int, double *, int, int, double, double *, double *);
23 static void Def_logcosh_JM (double *, int, double *, int, int, double, double *);
24 static void Def_exp_JM (double *, int, double *, int, int, double, double *);
25 static void calc_A_JM(double*, double*, double*, int*, int*, int*, double*, double*);
26 static void calc_K_JM(double*, int*, int*, double*);
27
28 #include <R_ext/Lapack.h>
29
30 static void
rowcentre_JM(double * ans,int n,int p)31 rowcentre_JM (double *ans, int n, int p)
32 {
33 /* mean centres nxp matrix ans */
34 double tmp;
35 int i, j;
36 for (i = 0; i < n; i++) {
37 tmp = 0;
38 for (j = 0; j < p; j++) {
39 tmp = tmp + ((double) ans[p * i + j]) / p;
40 }
41 for (j = 0; j < p; j++) {
42 ans[p * i + j] -= (double) tmp;
43 }
44 }
45 }
46
47 static void
colstandard_JM(double * ans,int n,int p)48 colstandard_JM (double *ans, int n, int p)
49 {
50 /* transform columns of nxp matrix ans to have zero mean and unit variance */
51 double tmp[2];
52 double tmp1;
53 int i, j;
54 for (i = 0; i < p; i++) {
55 tmp[0] = 0;
56 tmp[1] = 0;
57
58 for (j = 0; j < n; j++) {
59 tmp[0] += (double) ans[p * j + i];
60 tmp[1] += ((double) ans[p * j + i]) * ((double) ans[p * j + i]);
61 }
62
63 tmp[0] = tmp[0] / n;
64 tmp1 = (tmp[1] - n * (tmp[0]) * (tmp[0])) / (n - 1);
65
66 tmp[1] = sqrt (tmp1);
67 for (j = 0; j < n; j++) {
68 ans[p * j + i] =
69 (double) ((((double) ans[p * j + i]) - tmp[0]) / tmp[1]);
70 }
71 }
72 }
73
74
75 static void
svd_JM(double * mat,int * n,int * p,double * u,double * d,double * v)76 svd_JM (double *mat, int *n, int *p, double *u, double *d, double *v)
77 {
78
79 /* calculates svd decomposition of nxp matrix mat */
80 /* mat is a pointer to an nxp array of doubles */
81 /* n is a pointer to an integer specifying the no. of rows of mat */
82 /* p is a pointer to an integer specifying the no. of cols of mat */
83 /* u is a pointer to a double array of dimension (n,n) */
84 /* d is a pointer to a double array of dimension min(n,p) */
85 /* v is a pointer to a double array of dimension (p,p) */
86
87
88 int info, *iwork, lwork, a, b;
89 size_t iwork_size, ilwork, nn = *n, pp = *p, mm;
90 double *work, *mat1, *u1, *v1;
91 char jobz = 'A';
92
93 mm = min_JM(n,p);
94 iwork_size = 8 * (size_t) mm;
95
96 a = max_JM(n,p);
97 b = 4 * mm * mm + 4 * mm;
98 ilwork = 3 * mm * mm + max_JM(&a, &b);
99 if (ilwork > INT_MAX)
100 error("svd on %d x %d exceeds Fortran indexing limits");
101
102 work = Calloc (ilwork, double);
103 iwork = Calloc (iwork_size, int);
104 mat1 = Calloc (nn * pp, double);
105 u1 = Calloc (nn * nn, double);
106 v1 = Calloc (pp * pp, double);
107
108 transpose_mat_JM (mat, n, p, mat1);
109
110 lwork = ilwork;
111 F77_CALL (dgesdd) (&jobz, n, p, mat1, n, d, u1, n, v1, p, work,
112 &lwork, iwork, &info FCONE);
113
114 transpose_mat_JM (u1, n, n, u);
115 transpose_mat_JM (v1, p, p, v);
116
117 Free (mat1);
118 Free (u1);
119 Free (v1);
120 Free (work);
121 Free (iwork);
122 }
123
124
125 static void
transpose_mat_JM(double * mat,int * n,int * p,double * ans)126 transpose_mat_JM (double *mat, int *n, int *p, double *ans)
127 {
128 /* transpose nxp matrix mat */
129 int i, j;
130
131 for (i = 0; i < *n; i++) {
132 for (j = 0; j < *p; j++) {
133 *(ans + j * (*n) + i) = *(mat + i * (*p) + j);
134 }
135 }
136 }
137
138
min_JM(int * a,int * b)139 static int min_JM (int *a, int *b)
140 {
141 /* find minimum of a and b */
142 int ans;
143
144 ans = *b;
145 if (*a < *b) ans = *a;
146
147 return ans;
148 }
149
max_JM(int * a,int * b)150 static int max_JM (int *a, int *b)
151 {
152 /* find maximum of a and b */
153
154 int ans;
155
156 ans = *b;
157 if (*a > *b) ans = *a;
158
159 return ans;
160 }
161
162
163 static void
mmult_JM(double * A,int n,int p,double * B,int q,int r,double * C)164 mmult_JM (double *A, int n, int p, double *B, int q, int r, double *C)
165 {
166 /* matrix multiplication using FORTRAN BLAS routine SGEMM */
167 /* A is (n*p) and B is (q*r), A*B returned to C */
168
169 double alpha = 1.0, beta = 0.0;
170 int M, K, N;
171 char transA = 'N', transB = 'N';
172
173 if (p != q) {
174 error ("Error, matrices not suitable\nfor multiplication");
175 } else {
176 M = n;
177 K = p;
178 N = r;
179 F77_CALL (dgemm) (&transA, &transB, &N, &M, &K, &alpha, B, &N,
180 A, &K, &beta, C, &N FCONE FCONE);
181 }
182 }
183
184 static void
orthog_mat_JM(double * mat,int e,double * orthog)185 orthog_mat_JM (double *mat, int e, double *orthog)
186 {
187 /* take Wmat, (e*e), and return orthogonalized version to orthog_W */
188 double *u, *v, *d, *temp;
189 int i;
190 size_t ee = e;
191
192
193 u = Calloc (ee * ee, double);
194 d = Calloc (ee, double);
195 v = Calloc (ee * ee, double);
196 temp = Calloc (ee * ee, double);
197
198 svd_JM (mat, &e, &e, u, d, v);
199 for (i = 0; i < e; i++) {
200 temp[i * e + i] = 1 / (d[i]);
201 }
202
203 mmult_JM (u, e, e, temp, e, e, v);
204 transpose_mat_JM (u, &e, &e, temp);
205 mmult_JM (v, e, e, temp, e, e, u);
206 mmult_JM (u, e, e, mat, e, e, orthog);
207
208
209 Free (u);
210 Free (v);
211 Free (d);
212 Free (temp);
213
214 }
215
216 static void
Symm_logcosh_JM(double * w_init,int e,double * data,int f,int p,double alpha,double * w_final,double * Tol)217 Symm_logcosh_JM (double *w_init, int e, double *data, int f, int p, double alpha, double *w_final, double *Tol)
218 {
219
220 /* Function that carries out Symmetric ICA using a logcosh approximation to the neg. entropy function */
221
222 double *mat1, *mat2, *mat3, *mat4, *mat5, *mat6;
223 int i, j;
224 double mean;
225
226 if (e != f) {
227 error ("error in Symm_logcosh_JM, dims dont match");
228 }
229 else {
230 size_t es = (size_t)e * (size_t)e;
231 size_t ep = (size_t)e * (size_t)p;
232 mat1 = Calloc (ep, double);
233 mat2 = Calloc (ep, double);
234 mat3 = Calloc (es, double);
235 mat4 = Calloc (es, double);
236 mat5 = Calloc (es, double);
237 mat6 = Calloc (es, double);
238
239 mmult_JM (w_init, e, e, data, e, p, mat1);
240
241
242 for (i = 0; i < e; i++) {
243 for (j = 0; j < p; j++) {
244 mat1[i * p + j] = tanh (alpha * mat1[i * p + j]);
245 }
246 }
247 transpose_mat_JM (data, &e, &p, mat2);
248 for (i = 0; i < e; i++) {
249 for (j = 0; j < p; j++) {
250 mat2[i * p + j] = (mat2[i * p + j]) / p;
251 }
252 }
253 mmult_JM (mat1, e, p, mat2, p, e, mat3);
254 for (i = 0; i < e; i++) {
255 for (j = 0; j < p; j++) {
256 mat1[i * p + j] =
257 (alpha * (1 - (mat1[i * p + j]) * (mat1[i * p + j])));
258 }
259 }
260
261 for (i = 0; i < e; i++) {
262 mean = 0;
263 for (j = 0; j < p; j++) {
264 mean += ((mat1[i * p + j]) / p);
265 }
266 mat4[i * e + i] = mean;
267 }
268 mmult_JM (mat4, e, e, w_init, e, e, mat5);
269 for (i = 0; i < e; i++) {
270 for (j = 0; j < e; j++) {
271 mat4[i * e + j] = (mat3[i * e + j] - mat5[i * e + j]);
272 }
273 }
274
275 transpose_mat_JM (w_init, &e, &e, mat6);
276 orthog_mat_JM (mat4, e, w_final);
277
278
279 mmult_JM (w_final, e, e, mat6, e, e, mat5);
280 mean = 0;
281 for (i = 0; i < e; i++) {
282 if (fabs (1 - fabs (mat5[i * e + i])) > mean) {
283 mean = (fabs (1 - fabs (mat5[i * e + i])));
284 }
285 }
286 *Tol = mean;
287 Free (mat1);
288 Free (mat2);
289 Free (mat3);
290 Free (mat4);
291 Free (mat5);
292 Free (mat6);
293 }
294 }
295
296 static void
Def_logcosh_JM(double * w_init,int e,double * data,int f,int p,double alpha,double * w_final)297 Def_logcosh_JM (double *w_init, int e, double *data, int f, int p, double alpha, double *w_final)
298 {
299 /* Function that carries out Deflation ICA using an logcosh approximation to the neg. entropy function */
300
301 double *mat1, *mat2, *mat3, *mat4;
302 int i, j;
303 double mean;
304
305 if (e != f) {
306 error ("error in Def_logcosh_JM, dims dont match");
307 }
308 else {
309 mat1 = Calloc (p, double);
310 mat2 = Calloc ((size_t)e * (size_t)p, double);
311 mat3 = Calloc (e, double);
312 mat4 = Calloc (e, double);
313
314 mmult_JM (w_init, 1, e, data, e, p, mat1);
315
316
317 for (i = 0; i < p; i++) {
318 mat1[i] = tanh (alpha * mat1[i]);
319 }
320 transpose_mat_JM (data, &e, &p, mat2);
321 for (i = 0; i < e; i++) {
322 for (j = 0; j < p; j++) {
323 mat2[i * p + j] = (mat2[i * p + j]) / p;
324 }
325 }
326
327 mmult_JM (mat1, 1, p, mat2, p, e, mat3);
328 for (i = 0; i < p; i++) {
329 mat1[i] = (alpha * (1 - (mat1[i]) * (mat1[i])));
330 }
331
332 mean = 0;
333 for (j = 0; j < p; j++) {
334 mean += ((mat1[j]) / p);
335 }
336 for (i = 0; i < e; i++) {
337 mat4[i] = (w_init[i]) * mean;
338 }
339 for (i = 0; i < e; i++) {
340 w_final[i] = (mat3[i] - mat4[i]);
341 }
342
343 Free (mat1);
344 Free (mat2);
345 Free (mat3);
346 Free (mat4);
347
348 }
349 }
350
351 static void
Symm_exp_JM(double * w_init,int e,double * data,int f,int p,double alpha,double * w_final,double * Tol)352 Symm_exp_JM (double *w_init, int e, double *data, int f, int p, double alpha, double *w_final, double *Tol)
353 {
354 /* Function that carries out Symmetric ICA using a exponential approximation to the neg. entropy function */
355
356 double *mat1, *mat2, *mat3, *mat4, *mat5, *mat0, *mat6;
357 int i, j;
358 double mean;
359
360 if (e != f) {
361 error ("error in Symm_exp_JM, dims dont match");
362 }
363 else {
364 size_t ep = (size_t)e * (size_t)p;
365 size_t ee = (size_t)e * (size_t)e;
366 mat0 = Calloc (ep, double);
367 mat1 = Calloc (ep, double);
368 mat2 = Calloc (ep, double);
369 mat3 = Calloc (ee, double);
370 mat4 = Calloc (ee, double);
371 mat5 = Calloc (ee, double);
372 mat6 = Calloc (ee, double);
373 mmult_JM (w_init, e, e, data, e, p, mat1);
374 for (i = 0; i < e; i++) {
375 for (j = 0; j < p; j++) {
376 mat0[i * p + j] =
377 (mat1[i * p + j]) * exp (-0.5 * (mat1[i * p + j]) *
378 (mat1[i * p + j]));
379 }
380 }
381 transpose_mat_JM (data, &e, &p, mat2);
382 for (i = 0; i < e; i++) {
383 for (j = 0; j < p; j++) {
384 mat2[i * p + j] = (mat2[i * p + j]) / p;
385 }
386 }
387 mmult_JM (mat0, e, p, mat2, p, e, mat3);
388 for (i = 0; i < e; i++) {
389 for (j = 0; j < p; j++) {
390 mat1[i * p + j] =
391 ((1 - (mat1[i * p + j]) * (mat1[i * p + j])) *
392 exp (-0.5 * (mat1 [i * p + j]) * (mat1 [i * p + j])));
393 }
394 }
395
396 for (i = 0; i < e; i++) {
397 mean = 0;
398 for (j = 0; j < p; j++) {
399 mean += ((mat1[i * p + j]) / p);
400 }
401 mat4[i * e + i] = mean;
402 }
403 mmult_JM (mat4, e, e, w_init, e, e, mat5);
404 for (i = 0; i < e; i++) {
405 for (j = 0; j < e; j++) {
406 mat4[i * e + j] = (mat3[i * e + j] - mat5[i * e + j]);
407 }
408 }
409
410 transpose_mat_JM (w_init, &e, &e, mat6);
411 orthog_mat_JM (mat4, e, w_final);
412
413 mmult_JM (w_final, e, e, mat6, e, e, mat5);
414 mean = 0;
415 for (i = 0; i < e; i++) {
416 if (fabs (1 - fabs (mat5[i * e + i])) > mean) {
417 mean = (fabs (1 - fabs (mat5[i * e + i])));
418 }
419 }
420 *Tol = mean;
421 Free (mat1);
422 Free (mat2);
423 Free (mat3);
424 Free (mat4);
425 Free (mat5);
426 Free (mat0);
427 Free (mat6);
428 }
429 }
430
431 static void
Def_exp_JM(double * w_init,int e,double * data,int f,int p,double alpha,double * w_final)432 Def_exp_JM (double *w_init, int e, double *data, int f, int p, double alpha, double *w_final)
433 {
434 /* Function that carries out Deflation ICA using an exponential approximation to the neg. entropy function */
435
436 double *mat1, *mat2, *mat3, *mat4;
437 int i, j;
438 double mean;
439
440 if (e != f) {
441 error ("error in Def_exp_JM, dims dont match");
442 }
443 else {
444 mat1 = Calloc (p, double);
445 mat2 = Calloc ((size_t)e * (size_t)p, double);
446 mat3 = Calloc (e, double);
447 mat4 = Calloc (e, double);
448
449 mmult_JM (w_init, 1, e, data, e, p, mat1);
450
451 for (i = 0; i < p; i++) {
452 mat1[i] = ((mat1[i]) * exp (-0.5 * (mat1[i]) * (mat1[i])));
453 }
454
455 transpose_mat_JM (data, &e, &p, mat2);
456 for (i = 0; i < e; i++) {
457 for (j = 0; j < p; j++) {
458 mat2[i * p + j] = (mat2[i * p + j]) / p;
459 }
460 }
461
462 mmult_JM (mat1, 1, p, mat2, p, e, mat3);
463
464 mmult_JM (w_init, 1, e, data, e, p, mat1);
465 for (i = 0; i < p; i++) {
466 mat1[i] =
467 ((1 -
468 (mat1[i]) * (mat1[i])) * exp (-.5 * (mat1[i]) * (mat1[i])));
469 }
470 mean = 0;
471 for (j = 0; j < p; j++) {
472 mean += ((mat1[j]) / p);
473 }
474 for (i = 0; i < e; i++) {
475 mat4[i] = (w_init[i]) * mean;
476 }
477 for (i = 0; i < e; i++) {
478 w_final[i] = (mat3[i] - mat4[i]);
479 }
480
481
482 Free (mat1);
483 Free (mat2);
484 Free (mat3);
485 Free (mat4);
486
487 }
488 }
489
490 static void
gramsch_JM(double * ww,int n,int m,int k)491 gramsch_JM (double *ww, int n, int m, int k)
492 {
493 int ip, jp;
494 double tmp;
495 /* do Gram-Schmidt on row k of (n*m) matrix ww */
496 k -= 1;
497 if (k > n) {
498 error ("Error in gramsch");
499 }
500 else {
501 for (ip = 0; ip < k; ip++) {
502 tmp = 0;
503 for (jp = 0; jp < m; jp++) {
504 tmp += ((ww[m * ip + jp]) * (ww[m * k + jp]));
505 }
506 for (jp = 0; jp < m; jp++) {
507 ww[m * k + jp] = (ww[m * k + jp] - ((ww[m * ip + jp]) * tmp));
508 }
509 }
510 }
511 }
512
513 static void
rowstd_JM(double * ww,int n,int m,int k)514 rowstd_JM (double *ww, int n, int m, int k)
515 {
516 /* for ww (n*m), make ||ww[k, ]|| equal 1 */
517 double tmp = 0;
518 int i;
519 k -= 1;
520 if (k > n) {
521 error ("Error in rowstd");
522 }
523 else {
524 for (i = 0; i < m; i++) {
525 tmp += ((ww[k * m + i]) * (ww[k * m + i]));
526 }
527 tmp = sqrt (tmp);
528 for (i = 0; i < m; i++) {
529 ww[k * m + i] = ((ww[k * m + i]) / tmp);
530 }
531 }
532 }
533
534
535 static void
calc_K_JM(double * x,int * n,int * p,double * K)536 calc_K_JM(double *x, int *n, int *p, double *K)
537 {
538 int i, j;
539 double *xxt, *xt, *u, *d, *v, *temp1, *temp2;
540 size_t nn = *n, pp = *p;
541
542 xxt = Calloc (nn * nn, double);
543 xt = Calloc (nn * pp, double);
544
545 /* transpose x matrix */
546 transpose_mat_JM (x, n, p, xt);
547
548 /* calculate sample covariance matrix xxt */
549 mmult_JM (x, *n, *p, xt, *p, *n, xxt);
550 for (i = 0; i < *n; i++) {
551 for (j = 0; j < *n; j++) {
552 xxt[*n * i + j] = xxt[*n * i + j] / *p;
553 }
554 }
555 Free (xt);
556
557 /* calculate svd decomposition of xxt */
558 u = Calloc (nn * nn, double);
559 d = Calloc (nn, double);
560 v = Calloc (nn * nn, double);
561
562 svd_JM (xxt, n, n, u, d, v);
563
564 /* calculate K matrix*/
565 temp1 = Calloc (nn * nn, double);
566 temp2 = Calloc (nn * nn, double);
567
568 for (i = 0; i < *n; i++) {
569 temp1[*n * i + i] = 1 / sqrt (d[i]);
570 }
571
572 transpose_mat_JM (u, n, n, temp2);
573 mmult_JM (temp1, *n, *n, temp2, *n, *n, K);
574
575 Free (temp1);
576 Free (temp2);
577 Free(xxt);
578 Free(u);
579 Free(d);
580 Free(v);
581 }
582
583 static void
calc_A_JM(double * w,double * k,double * data,int * e,int * n,int * p,double * A,double * unmixed_data)584 calc_A_JM(double *w, double *k, double *data,
585 int *e, int *n, int *p, double *A, double *unmixed_data)
586 {
587 /* calculate un-mixing matrix A */
588 int i;
589 double *um, *umt, *umumt, *uu, *dd, *vv, *temp1, *temp2, *temp3;
590 size_t nn = *n, ee = *e;
591
592 um = Calloc (ee * nn, double);
593 umt = Calloc (nn * ee, double);
594
595 mmult_JM (w, *e, *e, k, *e, *n, um);
596 mmult_JM (um, *e, *n, data, *n, *p, unmixed_data);
597 transpose_mat_JM (um, e, n, umt);
598
599 umumt = Calloc (ee * ee, double);
600 mmult_JM (um, *e, *n, umt, *n, *e, umumt);
601
602 uu = Calloc (ee * ee, double);
603 dd = Calloc (ee, double);
604 vv = Calloc (ee * ee, double);
605 svd_JM (umumt, e, e, uu, dd, vv);
606
607 temp1 = Calloc (ee * ee, double);
608 for (i = 0; i < *e; i++) {
609 temp1[*e * i + i] = 1 / (dd[i]);
610 }
611
612 temp2 = Calloc (ee * ee, double);
613 temp3 = Calloc (ee * ee, double);
614 transpose_mat_JM (vv, e, e, temp3);
615 mmult_JM (temp3, *e, *e, temp1, *e, *e, temp2);
616 transpose_mat_JM (uu, e, e, vv);
617 mmult_JM (temp2, *e, *e, vv, *e, *e, uu);
618
619 mmult_JM (umt, *n, *e, uu, *e, *e, A);
620
621 Free(um);
622 Free(umt);
623 Free(umumt);
624 Free(uu);
625 Free(dd);
626 Free(vv);
627 Free(temp1);
628 Free(temp2);
629 Free(temp3);
630
631 }
632
633 static void
icainc_JM(double * data_matrix,double * w_matrix,int * nn,int * pp,int * ee,double * alpha,int * rowflag,int * colflag,int * funflag,int * maxit,double * lim,int * defflag,int * verbose,double * data_pre,double * Kmat1,double * w_final,double * ansa,double * ansx2)634 icainc_JM (double *data_matrix, double *w_matrix, int *nn, int *pp, int *ee,
635 double *alpha, int *rowflag, int *colflag, int *funflag, int *maxit,
636 double *lim, int *defflag, int *verbose, double *data_pre, double *Kmat1,
637 double *w_final, double *ansa, double *ansx2)
638 {
639
640 /* main ICA function */
641
642 int i, j, k;
643 size_t n = *nn, p = *pp, e = *ee;
644 double tol;
645 double *temp_w1, *temp_w2;
646 double *data1, *Kmat, *temp1, *w_init;
647
648 /* make a copy of the data matrix*/
649 data1 = Calloc (n * p, double);
650 for (i = 0; i < n; i++) {
651 for (j = 0; j < p; j++) {
652 data_pre[i * p + j] = data_matrix[i * p + j];
653 }
654 }
655
656 /* row center data matrix if required*/
657 if (*rowflag == 1) {
658 rowcentre_JM (data_pre, n, p);
659 if (*verbose == 1)
660 Rprintf ("Centering\n");
661 }
662
663 /* standardize columns of data matrix if required*/
664 if (*colflag == 1) {
665 colstandard_JM (data_pre, n, p);
666 Rprintf("colstandard\n");
667 }
668
669 /* calculate pre-whitening matrix Kmat */
670 if (*verbose == 1) Rprintf ("Whitening\n");
671 Kmat = Calloc (n * n, double);
672 calc_K_JM(data_pre, nn, pp, Kmat);
673
674 /* pre-whiten data and reduce dimension from size n to size e */
675
676 for (i = 0; i < e; i++) {
677 for (j = 0; j < n; j++) {
678 Kmat1[i * n + j] = Kmat[i * n + j];
679 }
680 }
681 mmult_JM (Kmat1, e, n, data_pre, n, p, data1);
682
683 /* calculate initial (orthogonal) unmixing matrix w */
684 temp1 = Calloc (e * e, double);
685 w_init = Calloc (e * e, double);
686 for (i = 0; i < e; i++) {
687 for (j = 0; j < e; j++) {
688 temp1[i * e + j] = w_matrix[i * e + j];
689 }
690 }
691 orthog_mat_JM (temp1, e, w_init);
692
693 /* Main ICA code */
694 if (*defflag == 0) {
695 if (*funflag == 1) {
696 if (*verbose == 1)
697 Rprintf("Symmetric FastICA using logcosh approx. to neg-entropy function\n");
698
699 i = 1;
700 Symm_logcosh_JM (w_init, e, data1, e, p, *alpha, w_final, &tol);
701 if (*verbose == 1)
702 Rprintf ("Iteration %d tol=%f\n", i, tol);
703 i = 2;
704
705 while ((tol > (*lim)) && (i < (*maxit))) {
706 Symm_logcosh_JM (w_final, e, data1, e, p, *alpha, w_final, &tol);
707 if (*verbose == 1)
708 Rprintf ("Iteration %d tol=%f\n", i, tol);
709 i += 1;
710 }
711 }
712
713 if (*funflag == 2) {
714 if (*verbose == 1)
715 Rprintf("Symmetric FastICA using exponential approx. to neg-entropy function\n");
716
717 i = 1;
718 Symm_exp_JM (w_init, e, data1, e, p, *alpha, w_final, &tol);
719 if (*verbose == 1) Rprintf ("Iteration %d tol=%f\n", i, tol);
720
721 i = 2;
722 while ((tol > (*lim)) && (i < (*maxit))) {
723 Symm_exp_JM (w_final, e, data1, e, p, *alpha, w_final, &tol);
724 if (*verbose == 1) Rprintf ("Iteration %d tol=%f\n", i, tol);
725 i += 1;
726 }
727 }
728 }
729
730 if (*defflag == 1) {
731 temp_w1 = Calloc (e, double);
732 temp_w2 = Calloc (e, double);
733
734 if (*funflag == 1) {
735 if (*verbose == 1)
736 Rprintf ("Deflation FastICA using logcosh approx. to neg-entropy function\n");
737
738 for (i = 0; i < e; i++) {
739 k = 0;
740 gramsch_JM (w_init, e, e, i + 1);
741 rowstd_JM (w_init, e, e, i + 1);
742 tol = 1;
743
744 while ((tol > (*lim)) && (k < (*maxit))) {
745 for (j = 0; j < e; j++) {
746 temp_w1[j] = w_init[i * e + j];
747 }
748 Def_logcosh_JM (temp_w1, e, data1, e, p, *alpha, temp_w2);
749 for (j = 0; j < e; j++) {
750 w_init[i * e + j] = temp_w2[j];
751 }
752 gramsch_JM (w_init, e, e, i + 1);
753 rowstd_JM (w_init, e, e, i + 1);
754 tol = 0;
755 for (j = 0; j < e; j++) {
756 tol += ((temp_w1[j]) * (w_init[i * e + j]));
757 }
758 tol = (fabs (fabs (tol) - 1));
759 k += 1;
760 }
761
762 if (*verbose == 1)
763 Rprintf ("Component %d needed %d iterations tol=%f\n",
764 i + 1, k, tol);
765
766 }
767 }
768 if (*funflag == 2) {
769
770 if (*verbose == 1)
771 Rprintf ("Deflation FastICA using exponential approx. to neg-entropy function\n");
772
773 for (i = 0; i < e; i++) {
774 k = 0;
775 gramsch_JM (w_init, e, e, i + 1);
776 rowstd_JM (w_init, e, e, i + 1);
777 tol = 1;
778
779 while ((tol > (*lim)) && (k < (*maxit))) {
780 for (j = 0; j < e; j++) {
781 temp_w1[j] = w_init[i * e + j];
782 }
783 Def_exp_JM (temp_w1, e, data1, e, p, *alpha, temp_w2);
784 for (j = 0; j < e; j++) {
785 w_init[i * e + j] = temp_w2[j];
786 }
787 gramsch_JM (w_init, e, e, i + 1);
788 rowstd_JM (w_init, e, e, i + 1);
789 tol = 0;
790 for (j = 0; j < e; j++) {
791 tol += ((temp_w1[j]) * (w_init[i * e + j]));
792 }
793 tol = (fabs (fabs (tol) - 1));
794 k += 1;
795 }
796
797 if (*verbose == 1)
798 Rprintf ("Component %d needed %d iterations tol=%f\n",
799 i + 1, k, tol);
800
801 }
802 }
803 for (i = 0; i < e; i++) {
804 for (j = 0; j < e; j++) {
805 w_final[i * e + j] = w_init[i * e + j];
806 }
807 }
808 Free (temp_w1);
809 Free (temp_w2);
810 }
811
812 /* calculate mixing matrix ansa */
813 calc_A_JM(w_final, Kmat1, data_pre, ee, nn, pp, ansa, ansx2);
814
815 Free (data1);
816 Free (Kmat);
817 Free (temp1);
818 Free (w_init);
819 }
820
821 #include <R_ext/Rdynload.h>
822
823 static const R_CMethodDef CEntries[] = {
824 {"icainc_JM", (DL_FUNC) &icainc_JM, 18},
825 {NULL, NULL, 0}
826 };
827
828
829 void
R_init_fastICA(DllInfo * dll)830 R_init_fastICA(DllInfo *dll)
831 {
832 R_registerRoutines(dll, CEntries, NULL, NULL, NULL);
833 R_useDynamicSymbols(dll, FALSE);
834 }
835