IT++ Logo
fastica.cpp
Go to the documentation of this file.
1
62#include <itpp/signal/fastica.h>
63#include <itpp/signal/sigfun.h>
68#include <itpp/base/matfunc.h>
69#include <itpp/base/random.h>
70#include <itpp/base/sort.h>
71#include <itpp/base/specmat.h>
72#include <itpp/base/svec.h>
74#include <itpp/stat/misc_stat.h>
75
76
77using namespace itpp;
78
79
84static void selcol(const mat oldMatrix, const vec maskVector, mat & newMatrix);
85static int pcamat(const mat vectors, const int numOfIC, int firstEig, int lastEig, mat & Es, vec & Ds);
86static void remmean(mat inVectors, mat & outVectors, vec & meanValue);
87static void whitenv(const mat vectors, const mat E, const mat D, mat & newVectors, mat & whiteningMatrix, mat & dewhiteningMatrix);
88static mat orth(const mat A);
89static mat mpower(const mat A, const double y);
90static ivec getSamples(const int max, const double percentage);
91static vec sumcol(const mat A);
92static bool fpica(const mat X, const mat whiteningMatrix, const mat dewhiteningMatrix, const int approach, const int numOfIC, const int g, const int finetune, const double a1, const double a2, double myy, const int stabilization, const double epsilon, const int maxNumIterations, const int maxFinetune, const int initState, mat guess, double sampleSize, mat & A, mat & W);
95namespace itpp
96{
97
98// Constructor, init default values
100{
101
102 // Init default values
103 approach = FICA_APPROACH_SYMM;
105 finetune = true;
106 a1 = 1.0;
107 a2 = 1.0;
108 mu = 1.0;
109 epsilon = 0.0001;
110 sampleSize = 1.0;
111 stabilization = false;
112 maxNumIterations = 100000;
113 maxFineTune = 100;
114 firstEig = 1;
115
116 mixedSig = ma_mixedSig;
117
118 lastEig = mixedSig.rows();
119 numOfIC = mixedSig.rows();
120 PCAonly = false;
121 initState = FICA_INIT_RAND;
122
123}
124
125// Call main function
127{
128
129 int Dim = numOfIC;
130
131 mat mixedSigC;
132 vec mixedMean;
133
134 mat guess;
135 if (initState == FICA_INIT_RAND)
136 guess = zeros(Dim, Dim);
137 else
138 guess = mat(initGuess);
139
140 VecPr = zeros(mixedSig.rows(), numOfIC);
141
142 icasig = zeros(numOfIC, mixedSig.cols());
143
144 remmean(mixedSig, mixedSigC, mixedMean);
145
146 if (pcamat(mixedSigC, numOfIC, firstEig, lastEig, E, D) < 1) {
147 // no principal components could be found (e.g. all-zero data): return the unchanged input
148 icasig = mixedSig;
149 return false;
150 }
151
152 whitenv(mixedSigC, E, diag(D), whitesig, whiteningMatrix, dewhiteningMatrix);
153
154 Dim = whitesig.rows();
155 if (numOfIC > Dim) numOfIC = Dim;
156
157 ivec NcFirst = to_ivec(zeros(numOfIC));
158 vec NcVp = D;
159 for (int i = 0; i < NcFirst.size(); i++) {
160
162 NcVp(NcFirst(i)) = 0.0;
163 VecPr.set_col(i, dewhiteningMatrix.get_col(i));
164
165 }
166
167 bool result = true;
168 if (PCAonly == false) {
169
170 result = fpica(whitesig, whiteningMatrix, dewhiteningMatrix, approach, numOfIC, g, finetune, a1, a2, mu, stabilization, epsilon, maxNumIterations, maxFineTune, initState, guess, sampleSize, A, W);
171
172 icasig = W * mixedSig;
173
174 }
175
176 else { // PCA only : returns E as IcaSig
177 icasig = VecPr;
178 }
179 return result;
180}
181
182void Fast_ICA::set_approach(int in_approach) { approach = in_approach; if (approach == FICA_APPROACH_DEFL) finetune = true; }
183
185
187
189
190void Fast_ICA::set_a1(double fl_a1) { a1 = fl_a1; }
191
192void Fast_ICA::set_a2(double fl_a2) { a2 = fl_a2; }
193
194void Fast_ICA::set_mu(double fl_mu) { mu = fl_mu; }
195
196void Fast_ICA::set_epsilon(double fl_epsilon) { epsilon = fl_epsilon; }
197
199
201
203
205
207
209
211
213{
214 initGuess = ma_initGuess;
215 initState = FICA_INIT_GUESS;
216}
217
218mat Fast_ICA::get_mixing_matrix() { if (PCAonly) { it_warning("No ICA performed."); return (zeros(1, 1));} else return A; }
219
220mat Fast_ICA::get_separating_matrix() { if (PCAonly) { it_warning("No ICA performed."); return(zeros(1, 1)); } else return W; }
221
222mat Fast_ICA::get_independent_components() { if (PCAonly) { it_warning("No ICA performed."); return(zeros(1, 1)); } else return icasig; }
223
225
227
228mat Fast_ICA::get_whitening_matrix() { return whiteningMatrix; }
229
230mat Fast_ICA::get_dewhitening_matrix() { return dewhiteningMatrix; }
231
232mat Fast_ICA::get_white_sig() { return whitesig; }
233
234} // namespace itpp
235
236
237static void selcol(const mat oldMatrix, const vec maskVector, mat & newMatrix)
238{
239
240 int numTaken = 0;
241
242 for (int i = 0; i < size(maskVector); i++) if (maskVector(i) == 1) numTaken++;
243
245
246 numTaken = 0;
247
248 for (int i = 0; i < size(maskVector); i++) {
249
250 if (maskVector(i) == 1) {
251
252 newMatrix.set_col(numTaken, oldMatrix.get_col(i));
253 numTaken++;
254
255 }
256 }
257
258 return;
259
260}
261
262static int pcamat(const mat vectors, const int numOfIC, int firstEig, int lastEig, mat & Es, vec & Ds)
263{
264
265 mat Et;
266 vec Dt;
267 cmat Ec;
268 cvec Dc;
269 double lowerLimitValue = 0.0,
270 higherLimitValue = 0.0;
271
272 int oldDimension = vectors.rows();
273
275
277
278 int maxLastEig = 0;
279
280 // Compute rank
281 for (int i = 0; i < Dt.length(); i++) if (Dt(i) > FICA_TOL) maxLastEig++;
282 if (maxLastEig < 1) return 0;
283
284 // Force numOfIC components
285 if (maxLastEig > numOfIC) maxLastEig = numOfIC;
286
287 vec eigenvalues = zeros(size(Dt));
288 vec eigenvalues2 = zeros(size(Dt));
289
291
292 sort(eigenvalues2);
293
294 vec lowerColumns = zeros(size(Dt));
295
296 for (int i = 0; i < size(Dt); i++) eigenvalues(i) = eigenvalues2(size(Dt) - i - 1);
297
298 if (lastEig > maxLastEig) lastEig = maxLastEig;
299
300 if (lastEig < oldDimension) lowerLimitValue = (eigenvalues(lastEig - 1) + eigenvalues(lastEig)) / 2;
302
303 for (int i = 0; i < size(Dt); i++) if (Dt(i) > lowerLimitValue) lowerColumns(i) = 1;
304
305 if (firstEig > 1) higherLimitValue = (eigenvalues(firstEig - 2) + eigenvalues(firstEig - 1)) / 2;
306 else higherLimitValue = eigenvalues(0) + 1;
307
308 vec higherColumns = zeros(size(Dt));
309 for (int i = 0; i < size(Dt); i++) if (Dt(i) < higherLimitValue) higherColumns(i) = 1;
310
312 for (int i = 0; i < size(Dt); i++) selectedColumns(i) = (lowerColumns(i) == 1 && higherColumns(i) == 1) ? 1 : 0;
313
315
316 int numTaken = 0;
317
318 for (int i = 0; i < size(selectedColumns); i++) if (selectedColumns(i) == 1) numTaken++;
319
320 Ds = zeros(numTaken);
321
322 numTaken = 0;
323
324 for (int i = 0; i < size(Dt); i++)
325 if (selectedColumns(i) == 1) {
326 Ds(numTaken) = Dt(i);
327 numTaken++;
328 }
329
330 return lastEig;
331
332}
333
334
335static void remmean(mat inVectors, mat & outVectors, vec & meanValue)
336{
337
338 outVectors = zeros(inVectors.rows(), inVectors.cols());
339 meanValue = zeros(inVectors.rows());
340
341 for (int i = 0; i < inVectors.rows(); i++) {
342
343 meanValue(i) = mean(inVectors.get_row(i));
344
345 for (int j = 0; j < inVectors.cols(); j++) outVectors(i, j) = inVectors(i, j) - meanValue(i);
346
347 }
348
349}
350
351static void whitenv(const mat vectors, const mat E, const mat D, mat & newVectors, mat & whiteningMatrix, mat & dewhiteningMatrix)
352{
353
354 whiteningMatrix = zeros(E.cols(), E.rows());
355 dewhiteningMatrix = zeros(E.rows(), E.cols());
356
357 for (int i = 0; i < D.cols(); i++) {
358 whiteningMatrix.set_row(i, std::pow(std::sqrt(D(i, i)), -1)*E.get_col(i));
359 dewhiteningMatrix.set_col(i, std::sqrt(D(i, i))*E.get_col(i));
360 }
361
362 newVectors = whiteningMatrix * vectors;
363
364 return;
365
366}
367
368static mat orth(const mat A)
369{
370
371 mat Q;
372 mat U, V;
373 vec S;
374 double eps = 2.2e-16;
375 double tol = 0.0;
376 int mmax = 0;
377 int r = 0;
378
379 svd(A, U, S, V);
380 if (A.rows() > A.cols()) {
381
382 U = U(0, U.rows() - 1, 0, A.cols() - 1);
383 S = S(0, A.cols() - 1);
384 }
385
386 mmax = (A.rows() > A.cols()) ? A.rows() : A.cols();
387
388 tol = mmax * eps * max(S);
389
390 for (int i = 0; i < size(S); i++) if (S(i) > tol) r++;
391
392 Q = U(0, U.rows() - 1, 0, r - 1);
393
394 return (Q);
395}
396
397static mat mpower(const mat A, const double y)
398{
399
400 mat T = zeros(A.rows(), A.cols());
401 mat dd = zeros(A.rows(), A.cols());
402 vec d = zeros(A.rows());
403 vec dOut = zeros(A.rows());
404
405 eig_sym(A, d, T);
406
407 dOut = pow(d, y);
408
409 diag(dOut, dd);
410
411 for (int i = 0; i < T.cols(); i++) T.set_col(i, T.get_col(i) / norm(T.get_col(i)));
412
413 return (T*dd*transpose(T));
414
415}
416
417static ivec getSamples(const int max, const double percentage)
418{
419
420 vec rd = randu(max);
421 sparse_vec sV;
422 ivec out;
423 int sZ = 0;
424
425 for (int i = 0; i < max; i++) if (rd(i) < percentage) { sV.add_elem(sZ, i); sZ++; }
426
427 out = to_ivec(full(sV));
428
429 return (out);
430
431}
432
433static vec sumcol(const mat A)
434{
435
436 vec out = zeros(A.cols());
437
438 for (int i = 0; i < A.cols(); i++) { out(i) = sum(A.get_col(i)); }
439
440 return (out);
441
442}
443
444static bool fpica(const mat X, const mat whiteningMatrix, const mat dewhiteningMatrix, const int approach, const int numOfIC, const int g, const int finetune, const double a1, const double a2, double myy, const int stabilization, const double epsilon, const int maxNumIterations, const int maxFinetune, const int initState, mat guess, double sampleSize, mat & A, mat & W)
445{
446
447 int vectorSize = X.rows();
448 int numSamples = X.cols();
449 int gOrig = g;
450 int gFine = finetune + 1;
451 double myyOrig = myy;
452 double myyK = 0.01;
453 int failureLimit = 5;
454 int usedNlinearity = 0;
455 double stroke = 0.0;
456 int notFine = 1;
457 int loong = 0;
458 int initialStateMode = initState;
459 double minAbsCos = 0.0, minAbsCos2 = 0.0;
460
461 if (sampleSize * numSamples < 1000) sampleSize = (1000 / (double)numSamples < 1.0) ? 1000 / (double)numSamples : 1.0;
462
463 if (sampleSize != 1.0) gOrig += 2;
464 if (myy != 1.0) gOrig += 1;
465
466 int fineTuningEnabled = 1;
467
468 if (!finetune) {
469 if (myy != 1.0) gFine = gOrig;
470 else gFine = gOrig + 1;
472 }
473
474 int stabilizationEnabled = stabilization;
475
476 if (!stabilization && myy != 1.0) stabilizationEnabled = true;
477
479
480 if (initState == FICA_INIT_GUESS && guess.rows() != whiteningMatrix.cols()) {
482
483 }
484 else if (guess.cols() < numOfIC) {
485
486 mat guess2 = randu(guess.rows(), numOfIC - guess.cols()) - 0.5;
488 }
489 else if (guess.cols() > numOfIC) guess = guess(0, guess.rows() - 1, 0, numOfIC - 1);
490
491 if (approach == FICA_APPROACH_SYMM) {
492
494 stroke = 0;
495 notFine = 1;
496 loong = 0;
497
498 A = zeros(vectorSize, numOfIC);
499 mat B = zeros(vectorSize, numOfIC);
500
501 if (initialStateMode == 0) B = orth(randu(vectorSize, numOfIC) - 0.5);
502 else B = whiteningMatrix * guess;
503
504 mat BOld = zeros(B.rows(), B.cols());
505 mat BOld2 = zeros(B.rows(), B.cols());
506
507 for (int round = 0; round < maxNumIterations; round++) {
508
509 if (round == maxNumIterations - 1) {
510
511 // If there is a convergence problem,
512 // we still want ot return something.
513 // This is difference with original
514 // Matlab implementation.
515 A = dewhiteningMatrix * B;
516 W = transpose(B) * whiteningMatrix;
517
518 return false;
519 }
520
521 B = B * mpower(transpose(B) * B , -0.5);
522
525
526 if (1 - minAbsCos < epsilon) {
527
528 if (fineTuningEnabled && notFine) {
529
530 notFine = 0;
532 myy = myyK * myyOrig;
533 BOld = zeros(B.rows(), B.cols());
534 BOld2 = zeros(B.rows(), B.cols());
535
536 }
537
538 else {
539
540 A = dewhiteningMatrix * B;
541 break;
542
543 }
544
545 } // IF epsilon
546
547 else if (stabilizationEnabled) {
548 if (!stroke && (1 - minAbsCos2 < epsilon)) {
549
550 stroke = myy;
551 myy /= 2;
552 if (mod(usedNlinearity, 2) == 0) usedNlinearity += 1 ;
553
554 }
555 else if (stroke) {
556
557 myy = stroke;
558 stroke = 0;
559 if (myy == 1 && mod(usedNlinearity, 2) != 0) usedNlinearity -= 1;
560
561 }
562 else if (!loong && (round > maxNumIterations / 2)) {
563
564 loong = 1;
565 myy /= 2;
566 if (mod(usedNlinearity, 2) == 0) usedNlinearity += 1;
567
568 }
569
570 } // stabilizationEnabled
571
572 BOld2 = BOld;
573 BOld = B;
574
575 switch (usedNlinearity) {
576
577 // pow3
578 case FICA_NONLIN_POW3 : {
579 B = (X * pow(transpose(X) * B, 3)) / numSamples - 3 * B;
580 break;
581 }
582 case(FICA_NONLIN_POW3+1) : {
583 mat Y = transpose(X) * B;
584 mat Gpow3 = pow(Y, 3);
585 vec Beta = sumcol(pow(Y, 4));
586 mat D = diag(pow(Beta - 3 * numSamples , -1));
587 B = B + myy * B * (transpose(Y) * Gpow3 - diag(Beta)) * D;
588 break;
589 }
590 case(FICA_NONLIN_POW3+2) : {
591 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
592 B = (Xsub * pow(transpose(Xsub) * B, 3)) / Xsub.cols() - 3 * B;
593 break;
594 }
595 case(FICA_NONLIN_POW3+3) : {
596 mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B;
597 mat Gpow3 = pow(Ysub, 3);
598 vec Beta = sumcol(pow(Ysub, 4));
599 mat D = diag(pow(Beta - 3 * Ysub.rows() , -1));
600 B = B + myy * B * (transpose(Ysub) * Gpow3 - diag(Beta)) * D;
601 break;
602 }
603
604 // TANH
605 case FICA_NONLIN_TANH : {
606 mat hypTan = tanh(a1 * transpose(X) * B);
607 B = (X * hypTan) / numSamples - elem_mult(reshape(repeat(sumcol(1 - pow(hypTan, 2)), B.rows()), B.rows(), B.cols()), B) / numSamples * a1;
608 break;
609 }
610 case(FICA_NONLIN_TANH+1) : {
611 mat Y = transpose(X) * B;
612 mat hypTan = tanh(a1 * Y);
613 vec Beta = sumcol(elem_mult(Y, hypTan));
614 vec Beta2 = sumcol(1 - pow(hypTan, 2));
615 mat D = diag(pow(Beta - a1 * Beta2 , -1));
616 B = B + myy * B * (transpose(Y) * hypTan - diag(Beta)) * D;
617 break;
618 }
619 case(FICA_NONLIN_TANH+2) : {
620 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
621 mat hypTan = tanh(a1 * transpose(Xsub) * B);
622 B = (Xsub * hypTan) / Xsub.cols() - elem_mult(reshape(repeat(sumcol(1 - pow(hypTan, 2)), B.rows()), B.rows(), B.cols()), B) / Xsub.cols() * a1;
623 break;
624 }
625 case(FICA_NONLIN_TANH+3) : {
626 mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B;
627 mat hypTan = tanh(a1 * Ysub);
629 vec Beta2 = sumcol(1 - pow(hypTan, 2));
630 mat D = diag(pow(Beta - a1 * Beta2 , -1));
631 B = B + myy * B * (transpose(Ysub) * hypTan - diag(Beta)) * D;
632 break;
633 }
634
635 // GAUSS
636 case FICA_NONLIN_GAUSS : {
637 mat U = transpose(X) * B;
638 mat Usquared = pow(U, 2);
639 mat ex = exp(-a2 * Usquared / 2);
640 mat gauss = elem_mult(U, ex);
641 mat dGauss = elem_mult(1 - a2 * Usquared, ex);
642 B = (X * gauss) / numSamples - elem_mult(reshape(repeat(sumcol(dGauss), B.rows()), B.rows(), B.cols()), B) / numSamples;
643 break;
644 }
645 case(FICA_NONLIN_GAUSS+1) : {
646 mat Y = transpose(X) * B;
647 mat ex = exp(-a2 * pow(Y, 2) / 2);
648 mat gauss = elem_mult(Y, ex);
649 vec Beta = sumcol(elem_mult(Y, gauss));
650 vec Beta2 = sumcol(elem_mult(1 - a2 * pow(Y, 2), ex));
651 mat D = diag(pow(Beta - Beta2 , -1));
652 B = B + myy * B * (transpose(Y) * gauss - diag(Beta)) * D;
653 break;
654 }
655 case(FICA_NONLIN_GAUSS+2) : {
656 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
657 mat U = transpose(Xsub) * B;
658 mat Usquared = pow(U, 2);
659 mat ex = exp(-a2 * Usquared / 2);
660 mat gauss = elem_mult(U, ex);
661 mat dGauss = elem_mult(1 - a2 * Usquared, ex);
662 B = (Xsub * gauss) / Xsub.cols() - elem_mult(reshape(repeat(sumcol(dGauss), B.rows()), B.rows(), B.cols()), B) / Xsub.cols();
663 break;
664 }
665 case(FICA_NONLIN_GAUSS+3) : {
666 mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B;
667 mat ex = exp(-a2 * pow(Ysub, 2) / 2);
668 mat gauss = elem_mult(Ysub, ex);
669 vec Beta = sumcol(elem_mult(Ysub, gauss));
670 vec Beta2 = sumcol(elem_mult(1 - a2 * pow(Ysub, 2), ex));
671 mat D = diag(pow(Beta - Beta2 , -1));
672 B = B + myy * B * (transpose(Ysub) * gauss - diag(Beta)) * D;
673 break;
674 }
675
676 // SKEW
677 case FICA_NONLIN_SKEW : {
678 B = (X * (pow(transpose(X) * B, 2))) / numSamples;
679 break;
680 }
681 case(FICA_NONLIN_SKEW+1) : {
682 mat Y = transpose(X) * B;
683 mat Gskew = pow(Y, 2);
684 vec Beta = sumcol(elem_mult(Y, Gskew));
685 mat D = diag(pow(Beta , -1));
686 B = B + myy * B * (transpose(Y) * Gskew - diag(Beta)) * D;
687 break;
688 }
689 case(FICA_NONLIN_SKEW+2) : {
690 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
691 B = (Xsub * (pow(transpose(Xsub) * B, 2))) / Xsub.cols();
692 break;
693 }
694 case(FICA_NONLIN_SKEW+3) : {
695 mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B;
696 mat Gskew = pow(Ysub, 2);
697 vec Beta = sumcol(elem_mult(Ysub, Gskew));
698 mat D = diag(pow(Beta , -1));
699 B = B + myy * B * (transpose(Ysub) * Gskew - diag(Beta)) * D;
700 break;
701 }
702
703
704 } // SWITCH usedNlinearity
705
706 } // FOR maxIterations
707
708 W = transpose(B) * whiteningMatrix;
709
710
711 } // IF FICA_APPROACH_SYMM APPROACH
712
713 // DEFLATION
714 else {
715
716 // FC 01/12/05
717 A = zeros(whiteningMatrix.cols(), numOfIC);
718 // A = zeros( vectorSize, numOfIC );
719 mat B = zeros(vectorSize, numOfIC);
720 W = transpose(B) * whiteningMatrix;
721 int round = 1;
722 int numFailures = 0;
723
724 while (round <= numOfIC) {
725
726 myy = myyOrig;
727
729 stroke = 0;
730
731 notFine = 1;
732 loong = 0;
733 int endFinetuning = 0;
734
735 vec w = zeros(vectorSize);
736
737 if (initialStateMode == 0)
738
739 w = randu(vectorSize) - 0.5;
740
741 else w = whiteningMatrix * guess.get_col(round);
742
743 w = w - B * transpose(B) * w;
744
745 w /= norm(w);
746
747 vec wOld = zeros(vectorSize);
748 vec wOld2 = zeros(vectorSize);
749
750 int i = 1;
751 int gabba = 1;
752
753 while (i <= maxNumIterations + gabba) {
754
755 w = w - B * transpose(B) * w;
756
757 w /= norm(w);
758
759 if (notFine) {
760
761 if (i == maxNumIterations + 1) {
762
763 round--;
764
765 numFailures++;
766
768
769 if (round == 0) {
770
771 A = dewhiteningMatrix * B;
772 W = transpose(B) * whiteningMatrix;
773
774 } // IF round
775
776 return false;
777
778 } // IF numFailures > failureLimit
779
780 break;
781
782 } // IF i == maxNumIterations+1
783
784 } // IF NOTFINE
785
786 else if (i >= endFinetuning) wOld = w;
787
788 if (norm(w - wOld) < epsilon || norm(w + wOld) < epsilon) {
789
790 if (fineTuningEnabled && notFine) {
791
792 notFine = 0;
797 myy = myyK * myyOrig;
799
800 } // IF finetuning
801
802 else {
803
804 numFailures = 0;
805
806 B.set_col(round - 1, w);
807
808 A.set_col(round - 1, dewhiteningMatrix*w);
809
810 W.set_row(round - 1, transpose(whiteningMatrix)*w);
811
812 break;
813
814 } // ELSE finetuning
815
816 } // IF epsilon
817
818 else if (stabilizationEnabled) {
819
820 if (stroke == 0.0 && (norm(w - wOld2) < epsilon || norm(w + wOld2) < epsilon)) {
821
822 stroke = myy;
823 myy /= 2.0 ;
824
825 if (mod(usedNlinearity, 2) == 0) {
826
828
829 } // IF MOD
830
831 }// IF !stroke
832
833 else if (stroke != 0.0) {
834
835 myy = stroke;
836 stroke = 0.0;
837
838 if (myy == 1 && mod(usedNlinearity, 2) != 0) {
840 }
841
842 } // IF Stroke
843
844 else if (notFine && !loong && i > maxNumIterations / 2) {
845
846 loong = 1;
847 myy /= 2.0;
848
849 if (mod(usedNlinearity, 2) == 0) {
850
852
853 } // IF MOD
854
855 } // IF notFine
856
857 } // IF stabilization
858
859
860 wOld2 = wOld;
861 wOld = w;
862
863 switch (usedNlinearity) {
864
865 // pow3
866 case FICA_NONLIN_POW3 : {
867 w = (X * pow(transpose(X) * w, 3)) / numSamples - 3 * w;
868 break;
869 }
870 case(FICA_NONLIN_POW3+1) : {
871 vec Y = transpose(X) * w;
872 vec Gpow3 = X * pow(Y, 3) / numSamples;
873 double Beta = dot(w, Gpow3);
874 w = w - myy * (Gpow3 - Beta * w) / (3 - Beta);
875 break;
876 }
877 case(FICA_NONLIN_POW3+2) : {
878 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
879 w = (Xsub * pow(transpose(Xsub) * w, 3)) / Xsub.cols() - 3 * w;
880 break;
881 }
882 case(FICA_NONLIN_POW3+3) : {
883 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
884 vec Gpow3 = Xsub * pow(transpose(Xsub) * w, 3) / (Xsub.cols());
885 double Beta = dot(w, Gpow3);
886 w = w - myy * (Gpow3 - Beta * w) / (3 - Beta);
887 break;
888 }
889
890 // TANH
891 case FICA_NONLIN_TANH : {
892 vec hypTan = tanh(a1 * transpose(X) * w);
893 w = (X * hypTan - a1 * sum(1 - pow(hypTan, 2)) * w) / numSamples;
894 break;
895 }
896 case(FICA_NONLIN_TANH+1) : {
897 vec Y = transpose(X) * w;
898 vec hypTan = tanh(a1 * Y);
899 double Beta = dot(w, X * hypTan);
900 w = w - myy * ((X * hypTan - Beta * w) / (a1 * sum(1 - pow(hypTan, 2)) - Beta));
901 break;
902 }
903 case(FICA_NONLIN_TANH+2) : {
904 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
905 vec hypTan = tanh(a1 * transpose(Xsub) * w);
906 w = (Xsub * hypTan - a1 * sum(1 - pow(hypTan, 2)) * w) / Xsub.cols();
907 break;
908 }
909 case(FICA_NONLIN_TANH+3) : {
910 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
911 vec hypTan = tanh(a1 * transpose(Xsub) * w);
912 double Beta = dot(w, Xsub * hypTan);
913 w = w - myy * ((Xsub * hypTan - Beta * w) / (a1 * sum(1 - pow(hypTan, 2)) - Beta));
914 break;
915 }
916
917 // GAUSS
918 case FICA_NONLIN_GAUSS : {
919 vec u = transpose(X) * w;
920 vec Usquared = pow(u, 2);
921 vec ex = exp(-a2 * Usquared / 2);
922 vec gauss = elem_mult(u, ex);
923 vec dGauss = elem_mult(1 - a2 * Usquared, ex);
924 w = (X * gauss - sum(dGauss) * w) / numSamples;
925 break;
926 }
927 case(FICA_NONLIN_GAUSS+1) : {
928 vec u = transpose(X) * w;
929 vec Usquared = pow(u, 2);
930
931 vec ex = exp(-a2 * Usquared / 2);
932 vec gauss = elem_mult(u, ex);
933 vec dGauss = elem_mult(1 - a2 * Usquared, ex);
934 double Beta = dot(w, X * gauss);
935 w = w - myy * ((X * gauss - Beta * w) / (sum(dGauss) - Beta));
936 break;
937 }
938 case(FICA_NONLIN_GAUSS+2) : {
939 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
940 vec u = transpose(Xsub) * w;
941 vec Usquared = pow(u, 2);
942 vec ex = exp(-a2 * Usquared / 2);
943 vec gauss = elem_mult(u, ex);
944 vec dGauss = elem_mult(1 - a2 * Usquared, ex);
945 w = (Xsub * gauss - sum(dGauss) * w) / Xsub.cols();
946 break;
947 }
948 case(FICA_NONLIN_GAUSS+3) : {
949 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
950 vec u = transpose(Xsub) * w;
951 vec Usquared = pow(u, 2);
952 vec ex = exp(-a2 * Usquared / 2);
953 vec gauss = elem_mult(u, ex);
954 vec dGauss = elem_mult(1 - a2 * Usquared, ex);
955 double Beta = dot(w, Xsub * gauss);
956 w = w - myy * ((Xsub * gauss - Beta * w) / (sum(dGauss) - Beta));
957 break;
958 }
959
960 // SKEW
961 case FICA_NONLIN_SKEW : {
962 w = (X * (pow(transpose(X) * w, 2))) / numSamples;
963 break;
964 }
965 case(FICA_NONLIN_SKEW+1) : {
966 vec Y = transpose(X) * w;
967 vec Gskew = X * pow(Y, 2) / numSamples;
968 double Beta = dot(w, Gskew);
969 w = w - myy * (Gskew - Beta * w / (-Beta));
970 break;
971 }
972 case(FICA_NONLIN_SKEW+2) : {
973 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
974 w = (Xsub * (pow(transpose(Xsub) * w, 2))) / Xsub.cols();
975 break;
976 }
977 case(FICA_NONLIN_SKEW+3) : {
978 mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
979 vec Gskew = Xsub * pow(transpose(Xsub) * w, 2) / Xsub.cols();
980 double Beta = dot(w, Gskew);
981 w = w - myy * (Gskew - Beta * w) / (-Beta);
982 break;
983 }
984
985 } // SWITCH nonLinearity
986
987 w /= norm(w);
988 i++;
989
990 } // WHILE i<= maxNumIterations+gabba
991
992 round++;
993
994 } // While round <= numOfIC
995
996 } // ELSE Deflation
997 return true;
998} // FPICA
General array class.
Definition array.h:105
int size() const
Returns the number of data elements in the array object.
Definition array.h:155
int length() const
Returns the number of data elements in the array object.
Definition array.h:157
mat get_white_sig()
Get whitened signals.
Definition fastica.cpp:232
Fast_ICA(mat ma_mixed_sig)
Constructor.
Definition fastica.cpp:99
void set_init_guess(mat ma_initGuess)
Set initial guess matrix instead of random (default)
Definition fastica.cpp:212
mat get_dewhitening_matrix()
Get the de-whitening matrix.
Definition fastica.cpp:230
void set_mu(double fl_mu)
Set parameter.
Definition fastica.cpp:194
void set_a1(double fl_a1)
Set parameter.
Definition fastica.cpp:190
void set_approach(int in_approach)
Set approach : FICA_APPROACH_DEFL or FICA_APPROACH_SYMM (default)
Definition fastica.cpp:182
void set_max_fine_tune(int in_maxFineTune)
Set maximum number of iterations for fine tuning.
Definition fastica.cpp:204
void set_pca_only(bool in_PCAonly)
If true, only perform Principal Component Analysis (default = false)
Definition fastica.cpp:210
void set_a2(double fl_a2)
Set parameter.
Definition fastica.cpp:192
bool separate(void)
Explicit launch of main FastICA function.
Definition fastica.cpp:126
mat get_independent_components()
Get separated signals.
Definition fastica.cpp:222
void set_first_eig(int in_firstEig)
Set first eigenvalue index to take into account.
Definition fastica.cpp:206
void set_sample_size(double fl_sampleSize)
Set sample size.
Definition fastica.cpp:198
void set_non_linearity(int in_g)
Set non-linearity.
Definition fastica.cpp:186
mat get_mixing_matrix()
Get mixing matrix.
Definition fastica.cpp:218
void set_epsilon(double fl_epsilon)
Set convergence parameter .
Definition fastica.cpp:196
mat get_whitening_matrix()
Get the whitening matrix.
Definition fastica.cpp:228
void set_max_num_iterations(int in_maxNumIterations)
Set maximum number of iterations.
Definition fastica.cpp:202
void set_last_eig(int in_lastEig)
Set last eigenvalue index to take into account.
Definition fastica.cpp:208
void set_stabilization(bool in_stabilization)
Set stabilization mode true or off.
Definition fastica.cpp:200
mat get_separating_matrix()
Get separating matrix.
Definition fastica.cpp:220
mat get_principal_eigenvectors()
Get nrIC first columns of the de-whitening matrix.
Definition fastica.cpp:226
void set_fine_tune(bool in_finetune)
Set fine tuning.
Definition fastica.cpp:188
int get_nrof_independent_components()
Get number of independent components.
Definition fastica.cpp:224
void set_nrof_independent_components(int in_nrIC)
Set number of independent components to separate.
Definition fastica.cpp:184
Definitions of eigenvalue decomposition functions.
static bool fpica(const mat X, const mat whiteningMatrix, const mat dewhiteningMatrix, const int approach, const int numOfIC, const int g, const int finetune, const double a1, const double a2, double myy, const int stabilization, const double epsilon, const int maxNumIterations, const int maxFinetune, const int initState, mat guess, double sampleSize, mat &A, mat &W)
Local functions for FastICA.
Definition fastica.cpp:444
static mat mpower(const mat A, const double y)
Local functions for FastICA.
Definition fastica.cpp:397
static vec sumcol(const mat A)
Local functions for FastICA.
Definition fastica.cpp:433
static int pcamat(const mat vectors, const int numOfIC, int firstEig, int lastEig, mat &Es, vec &Ds)
Local functions for FastICA.
Definition fastica.cpp:262
static void whitenv(const mat vectors, const mat E, const mat D, mat &newVectors, mat &whiteningMatrix, mat &dewhiteningMatrix)
Local functions for FastICA.
Definition fastica.cpp:351
static void selcol(const mat oldMatrix, const vec maskVector, mat &newMatrix)
Local functions for FastICA.
Definition fastica.cpp:237
static void remmean(mat inVectors, mat &outVectors, vec &meanValue)
Local functions for FastICA.
Definition fastica.cpp:335
static mat orth(const mat A)
Local functions for FastICA.
Definition fastica.cpp:368
static ivec getSamples(const int max, const double percentage)
Local functions for FastICA.
Definition fastica.cpp:417
Definition of FastICA (Independent Component Analysis) for IT++.
#define FICA_INIT_GUESS
Set predefined start for Fast_ICA.
Definition fastica.h:85
#define FICA_NONLIN_POW3
Use x^3 non-linearity.
Definition fastica.h:74
#define FICA_INIT_RAND
Set random start for Fast_ICA.
Definition fastica.h:83
#define FICA_NONLIN_GAUSS
Use Gaussian non-linearity.
Definition fastica.h:78
#define FICA_APPROACH_SYMM
Use symmetric approach : compute all ICs at a time.
Definition fastica.h:71
#define FICA_APPROACH_DEFL
Use deflation approach : compute IC one-by-one in a Gram-Schmidt-like fashion.
Definition fastica.h:69
#define FICA_TOL
Eigenvalues of the covariance matrix lower than FICA_TOL are discarded for analysis.
Definition fastica.h:88
#define FICA_NONLIN_SKEW
Use skew non-linearity.
Definition fastica.h:80
#define FICA_NONLIN_TANH
Use tanh(x) non-linearity.
Definition fastica.h:76
Mat< T > diag(const Vec< T > &v, const int K=0)
Create a diagonal matrix using vector v as its diagonal.
Definition matfunc.h:557
#define it_warning(s)
Display a warning message.
Definition itassert.h:173
vec tanh(const vec &x)
Tan hyperbolic function.
Definition trig_hyp.h:97
vec pow(const double x, const vec &y)
Calculates x to the power of y (x^y)
Definition log_exp.h:176
vec exp(const vec &x)
Exp of the elements of a vector x.
Definition log_exp.h:155
int size(const Vec< T > &v)
Length of vector.
Definition matfunc.h:55
T sum(const Vec< T > &v)
Sum of all elements in the vector.
Definition matfunc.h:59
void transpose(const Mat< T > &m, Mat< T > &out)
Transposition of the matrix m returning the transposed matrix in out.
Definition matfunc.h:308
bool svd(const mat &A, vec &S)
Get singular values s of a real matrix A using SVD.
Definition svd.cpp:185
bool eig_sym(const mat &A, vec &d, mat &V)
Calculates the eigenvalues and eigenvectors of a symmetric real matrix.
Definition eigen.cpp:252
int mod(int k, int n)
Calculates the modulus, i.e. the signed reminder after division.
Definition elem_math.h:166
T min(const Vec< T > &in)
Minimum value of vector.
Definition min_max.h:125
T max(const Vec< T > &v)
Maximum value of vector.
Definition min_max.h:45
int max_index(const Vec< T > &in)
Return the postion of the maximum element in the vector.
Definition min_max.h:208
double randu(void)
Generates a random uniform (0,1) number.
Definition random.h:804
Mat< T > reshape(const Mat< T > &m, int rows, int cols)
Reshape the matrix into an rows*cols matrix.
Definition matfunc.h:822
mat cov(const mat &X, bool is_zero_mean)
Covariance matrix calculation.
Definition sigfun.cpp:226
ITPP_EXPORT vec zeros(int size)
A Double vector of zeros.
double norm(const cvec &v)
Calculate the 2-norm: norm(v)=sqrt(sum(abs(v).^2))
Definition misc_stat.cpp:77
double mean(const vec &v)
The mean value.
Definition misc_stat.cpp:36
Various functions on vectors and matrices - header file.
Minimum and maximum functions on vectors and matrices.
Miscellaneous statistics functions and classes - header file.
itpp namespace
Definition itmex.h:37
Mat< T > full(const Sparse_Mat< T > &s)
Convert a sparse matrix s into its dense representation.
Definition smat.h:998
Mat< Num_T > concat_horizontal(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Horizontal concatenation of two matrices.
Definition mat.h:1194
ITPP_EXPORT double round(double x)
Round to nearest integer, return result in double.
Num_T dot(const Vec< Num_T > &v1, const Vec< Num_T > &v2)
Inner (dot) product of two vectors v1 and v2.
Definition vec.h:1005
const double eps
Constant eps.
Definition misc.h:109
ivec to_ivec(const Vec< T > &v)
Converts a Vec<T> to ivec.
Definition converters.h:79
Mat< Num_T > elem_mult(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Element wise multiplication of two matrices.
Definition mat.h:1582
bin abs(const bin &inbin)
absolute value of bin
Definition binary.h:174
Definition of classes for random number generators.
Resampling functions - header file.
Definitions of signal processing functions.
Sorting functions.
Definitions of special vectors and matrices.
Definitions of Singular Value Decompositions.
Sparse Vector Class definitions.
Trigonometric and hyperbolic functions - header file.

Generated on Tue Aug 17 2021 10:59:15 for IT++ by Doxygen 1.9.8