84 static void selcol(
const mat oldMatrix,
const vec maskVector, mat & newMatrix);
85 static void pcamat(
const mat vectors,
const int numOfIC,
int firstEig,
int lastEig, mat & Es, vec & Ds);
86 static void remmean(mat inVectors, mat & outVectors, vec & meanValue);
87 static void whitenv(
const mat vectors,
const mat E,
const mat D, mat & newVectors, mat & whiteningMatrix, mat & dewhiteningMatrix);
88 static mat
orth(
const mat A);
89 static mat
mpower(
const mat A,
const double y);
90 static ivec
getSamples(
const int max,
const double percentage);
91 static vec
sumcol(
const mat A);
92 static void 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);
111 stabilization =
false;
112 maxNumIterations = 100000;
116 mixedSig = ma_mixedSig;
118 lastEig = mixedSig.rows();
119 numOfIC = mixedSig.rows();
136 guess =
zeros(Dim, Dim);
138 guess = mat(initGuess);
140 VecPr =
zeros(mixedSig.rows(), numOfIC);
142 icasig =
zeros(numOfIC, mixedSig.cols());
144 remmean(mixedSig, mixedSigC, mixedMean);
146 pcamat(mixedSigC, numOfIC, firstEig, lastEig, E, D);
148 whitenv(mixedSigC, E,
diag(D), whitesig, whiteningMatrix, dewhiteningMatrix);
153 for (
int i = 0; i < NcFirst.size(); i++) {
156 NcVp(NcFirst(i)) = 0.0;
157 VecPr.set_col(i, dewhiteningMatrix.get_col(i));
161 if (PCAonly ==
false) {
163 Dim = whitesig.rows();
165 if (numOfIC > Dim) numOfIC = Dim;
167 fpica(whitesig, whiteningMatrix, dewhiteningMatrix, approach, numOfIC, g, finetune, a1, a2, mu, stabilization, epsilon, maxNumIterations, maxFineTune, initState, guess, sampleSize, A, W);
169 icasig = W * mixedSig;
210 initGuess = ma_initGuess;
233 static void selcol(
const mat oldMatrix,
const vec maskVector, mat & newMatrix)
238 for (
int i = 0; i <
size(maskVector); i++)
if (maskVector(i) == 1) numTaken++;
240 newMatrix =
zeros(oldMatrix.rows(), numTaken);
244 for (
int i = 0; i <
size(maskVector); i++) {
246 if (maskVector(i) == 1) {
248 newMatrix.set_col(numTaken, oldMatrix.get_col(i));
258 static void pcamat(
const mat vectors,
const int numOfIC,
int firstEig,
int lastEig, mat & Es, vec & Ds)
265 double lowerLimitValue = 0.0,
266 higherLimitValue = 0.0;
268 int oldDimension = vectors.rows();
272 eig_sym(covarianceMatrix, Dt, Et);
277 for (
int i = 0; i < Dt.length(); i++)
if (Dt(i) >
FICA_TOL) maxLastEig++;
280 if (maxLastEig > numOfIC) maxLastEig = numOfIC;
291 for (
int i = 0; i <
size(Dt); i++) eigenvalues(i) = eigenvalues2(
size(Dt) - i - 1);
293 if (lastEig > maxLastEig) lastEig = maxLastEig;
295 if (lastEig < oldDimension) lowerLimitValue = (eigenvalues(lastEig - 1) + eigenvalues(lastEig)) / 2;
296 else lowerLimitValue = eigenvalues(oldDimension - 1) - 1;
298 for (
int i = 0; i < size(Dt); i++) if (Dt(i) > lowerLimitValue) lowerColumns(i) = 1;
300 if (firstEig > 1) higherLimitValue = (eigenvalues(firstEig - 2) + eigenvalues(firstEig - 1)) / 2;
301 else higherLimitValue = eigenvalues(0) + 1;
304 for (
int i = 0; i <
size(Dt); i++)
if (Dt(i) < higherLimitValue) higherColumns(i) = 1;
307 for (
int i = 0; i <
size(Dt); i++) selectedColumns(i) = (lowerColumns(i) == 1 && higherColumns(i) == 1) ? 1 : 0;
309 selcol(Et, selectedColumns, Es);
313 for (
int i = 0; i <
size(selectedColumns); i++)
if (selectedColumns(i) == 1) numTaken++;
315 Ds =
zeros(numTaken);
319 for (
int i = 0; i <
size(Dt); i++)
320 if (selectedColumns(i) == 1) {
321 Ds(numTaken) = Dt(i);
330 static void remmean(mat inVectors, mat & outVectors, vec & meanValue)
333 outVectors =
zeros(inVectors.rows(), inVectors.cols());
334 meanValue =
zeros(inVectors.rows());
336 for (
int i = 0; i < inVectors.rows(); i++) {
338 meanValue(i) =
mean(inVectors.get_row(i));
340 for (
int j = 0; j < inVectors.cols(); j++) outVectors(i, j) = inVectors(i, j) - meanValue(i);
346 static void whitenv(
const mat vectors,
const mat E,
const mat D, mat & newVectors, mat & whiteningMatrix, mat & dewhiteningMatrix)
349 whiteningMatrix =
zeros(E.cols(), E.rows());
350 dewhiteningMatrix =
zeros(E.rows(), E.cols());
352 for (
int i = 0; i < D.cols(); i++) {
354 dewhiteningMatrix.set_col(i,
std::sqrt(D(i, i))*E.get_col(i));
357 newVectors = whiteningMatrix * vectors;
369 double eps = 2.2e-16;
375 if (A.rows() > A.cols()) {
377 U = U(0, U.rows() - 1, 0, A.cols() - 1);
378 S = S(0, A.cols() - 1);
381 mmax = (A.rows() > A.cols()) ? A.rows() : A.cols();
383 tol = mmax * eps *
max(S);
385 for (
int i = 0; i < size(S); i++) if (S(i) > tol) r++;
387 Q = U(0, U.rows() - 1, 0, r - 1);
392 static mat
mpower(
const mat A,
const double y)
395 mat T =
zeros(A.rows(), A.cols());
396 mat dd =
zeros(A.rows(), A.cols());
397 vec d =
zeros(A.rows());
398 vec dOut =
zeros(A.rows());
406 for (
int i = 0; i < T.cols(); i++) T.set_col(i, T.get_col(i) /
norm(T.get_col(i)));
420 for (
int i = 0; i <
max; i++)
if (rd(i) < percentage) { sV.add_elem(sZ, i); sZ++; }
431 vec out =
zeros(A.cols());
433 for (
int i = 0; i < A.cols(); i++) { out(i) =
sum(A.get_col(i)); }
439 static void 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)
442 int vectorSize = X.rows();
443 int numSamples = X.cols();
445 int gFine = finetune + 1;
446 double myyOrig = myy;
448 int failureLimit = 5;
449 int usedNlinearity = 0;
453 int initialStateMode = initState;
454 double minAbsCos = 0.0, minAbsCos2 = 0.0;
456 if (sampleSize * numSamples < 1000) sampleSize = (1000 / (double)numSamples < 1.0) ? 1000 / (double)numSamples : 1.0;
458 if (sampleSize != 1.0) gOrig += 2;
459 if (myy != 1.0) gOrig += 1;
461 int fineTuningEnabled = 1;
464 if (myy != 1.0) gFine = gOrig;
465 else gFine = gOrig + 1;
466 fineTuningEnabled = 0;
469 int stabilizationEnabled = stabilization;
471 if (!stabilization && myy != 1.0) stabilizationEnabled =
true;
473 usedNlinearity = gOrig;
475 if (initState ==
FICA_INIT_GUESS && guess.rows() != whiteningMatrix.cols()) {
476 initialStateMode = 0;
479 else if (guess.cols() < numOfIC) {
481 mat guess2 =
randu(guess.rows(), numOfIC - guess.cols()) - 0.5;
484 else if (guess.cols() > numOfIC) guess = guess(0, guess.rows() - 1, 0, numOfIC - 1);
488 usedNlinearity = gOrig;
493 A =
zeros(vectorSize, numOfIC);
494 mat B =
zeros(vectorSize, numOfIC);
496 if (initialStateMode == 0) B =
orth(
randu(vectorSize, numOfIC) - 0.5);
497 else B = whiteningMatrix * guess;
499 mat BOld =
zeros(B.rows(), B.cols());
500 mat BOld2 =
zeros(B.rows(), B.cols());
504 if (
round == maxNumIterations - 1) {
510 A = dewhiteningMatrix * B;
521 if (1 - minAbsCos < epsilon) {
523 if (fineTuningEnabled && notFine) {
526 usedNlinearity = gFine;
527 myy = myyK * myyOrig;
528 BOld =
zeros(B.rows(), B.cols());
529 BOld2 =
zeros(B.rows(), B.cols());
535 A = dewhiteningMatrix * B;
542 else if (stabilizationEnabled) {
543 if (!stroke && (1 - minAbsCos2 < epsilon)) {
547 if (
mod(usedNlinearity, 2) == 0) usedNlinearity += 1 ;
554 if (myy == 1 &&
mod(usedNlinearity, 2) != 0) usedNlinearity -= 1;
557 else if (!loong && (
round > maxNumIterations / 2)) {
561 if (
mod(usedNlinearity, 2) == 0) usedNlinearity += 1;
570 switch (usedNlinearity) {
574 B = (X *
pow(
transpose(X) * B, 3)) / numSamples - 3 * B;
579 mat Gpow3 =
pow(Y, 3);
581 mat D =
diag(
pow(Beta - 3 * numSamples , -1));
586 mat Xsub = X.get_cols(
getSamples(numSamples, sampleSize));
587 B = (Xsub *
pow(
transpose(Xsub) * B, 3)) / Xsub.cols() - 3 * B;
592 mat Gpow3 =
pow(Ysub, 3);
594 mat D =
diag(
pow(Beta - 3 * Ysub.rows() , -1));
595 B = B + myy * B * (
transpose(Ysub) * Gpow3 -
diag(Beta)) * D;
607 mat hypTan =
tanh(a1 * Y);
610 mat D =
diag(
pow(Beta - a1 * Beta2 , -1));
615 mat Xsub = X.get_cols(
getSamples(numSamples, sampleSize));
622 mat hypTan =
tanh(a1 * Ysub);
625 mat D =
diag(
pow(Beta - a1 * Beta2 , -1));
626 B = B + myy * B * (
transpose(Ysub) * hypTan -
diag(Beta)) * D;
633 mat Usquared =
pow(U, 2);
634 mat ex =
exp(-a2 * Usquared / 2);
636 mat dGauss =
elem_mult(1 - a2 * Usquared, ex);
642 mat ex =
exp(-a2 *
pow(Y, 2) / 2);
646 mat D =
diag(
pow(Beta - Beta2 , -1));
651 mat Xsub = X.get_cols(
getSamples(numSamples, sampleSize));
653 mat Usquared =
pow(U, 2);
654 mat ex =
exp(-a2 * Usquared / 2);
656 mat dGauss =
elem_mult(1 - a2 * Usquared, ex);
662 mat ex =
exp(-a2 *
pow(Ysub, 2) / 2);
666 mat D =
diag(
pow(Beta - Beta2 , -1));
667 B = B + myy * B * (
transpose(Ysub) * gauss -
diag(Beta)) * D;
678 mat Gskew =
pow(Y, 2);
685 mat Xsub = X.get_cols(
getSamples(numSamples, sampleSize));
686 B = (Xsub * (
pow(
transpose(Xsub) * B, 2))) / Xsub.cols();
691 mat Gskew =
pow(Ysub, 2);
694 B = B + myy * B * (
transpose(Ysub) * Gskew -
diag(Beta)) * D;
712 A =
zeros(whiteningMatrix.cols(), numOfIC);
714 mat B =
zeros(vectorSize, numOfIC);
719 while (round <= numOfIC) {
723 usedNlinearity = gOrig;
728 int endFinetuning = 0;
730 vec w =
zeros(vectorSize);
732 if (initialStateMode == 0)
734 w =
randu(vectorSize) - 0.5;
736 else w = whiteningMatrix * guess.get_col(round);
742 vec wOld =
zeros(vectorSize);
743 vec wOld2 =
zeros(vectorSize);
748 while (i <= maxNumIterations + gabba) {
756 if (i == maxNumIterations + 1) {
762 if (numFailures > failureLimit) {
766 A = dewhiteningMatrix * B;
781 else if (i >= endFinetuning) wOld = w;
783 if (
norm(w - wOld) < epsilon ||
norm(w + wOld) < epsilon) {
785 if (fineTuningEnabled && notFine) {
789 wOld =
zeros(vectorSize);
790 wOld2 =
zeros(vectorSize);
791 usedNlinearity = gFine;
792 myy = myyK * myyOrig;
793 endFinetuning = maxFinetune + i;
801 B.set_col(round - 1, w);
803 A.set_col(round - 1, dewhiteningMatrix*w);
805 W.set_row(round - 1,
transpose(whiteningMatrix)*w);
813 else if (stabilizationEnabled) {
815 if (stroke == 0.0 && (
norm(w - wOld2) < epsilon ||
norm(w + wOld2) < epsilon)) {
820 if (
mod(usedNlinearity, 2) == 0) {
828 else if (stroke != 0.0) {
833 if (myy == 1 &&
mod(usedNlinearity, 2) != 0) {
839 else if (notFine && !loong && i > maxNumIterations / 2) {
844 if (
mod(usedNlinearity, 2) == 0) {
858 switch (usedNlinearity) {
862 w = (X *
pow(
transpose(X) * w, 3)) / numSamples - 3 * w;
867 vec Gpow3 = X *
pow(Y, 3) / numSamples;
868 double Beta =
dot(w, Gpow3);
869 w = w - myy * (Gpow3 - Beta * w) / (3 - Beta);
873 mat Xsub = X.get_cols(
getSamples(numSamples, sampleSize));
874 w = (Xsub *
pow(
transpose(Xsub) * w, 3)) / Xsub.cols() - 3 * w;
878 mat Xsub = X.get_cols(
getSamples(numSamples, sampleSize));
879 vec Gpow3 = Xsub *
pow(
transpose(Xsub) * w, 3) / (Xsub.cols());
880 double Beta =
dot(w, Gpow3);
881 w = w - myy * (Gpow3 - Beta * w) / (3 - Beta);
888 w = (X * hypTan - a1 *
sum(1 -
pow(hypTan, 2)) * w) / numSamples;
893 vec hypTan =
tanh(a1 * Y);
894 double Beta =
dot(w, X * hypTan);
895 w = w - myy * ((X * hypTan - Beta * w) / (a1 *
sum(1 -
pow(hypTan, 2)) - Beta));
899 mat Xsub = X.get_cols(
getSamples(numSamples, sampleSize));
901 w = (Xsub * hypTan - a1 *
sum(1 -
pow(hypTan, 2)) * w) / Xsub.cols();
905 mat Xsub = X.get_cols(
getSamples(numSamples, sampleSize));
907 double Beta =
dot(w, Xsub * hypTan);
908 w = w - myy * ((Xsub * hypTan - Beta * w) / (a1 *
sum(1 -
pow(hypTan, 2)) - Beta));
915 vec Usquared =
pow(u, 2);
916 vec ex =
exp(-a2 * Usquared / 2);
918 vec dGauss =
elem_mult(1 - a2 * Usquared, ex);
919 w = (X * gauss -
sum(dGauss) * w) / numSamples;
924 vec Usquared =
pow(u, 2);
926 vec ex =
exp(-a2 * Usquared / 2);
928 vec dGauss =
elem_mult(1 - a2 * Usquared, ex);
929 double Beta =
dot(w, X * gauss);
930 w = w - myy * ((X * gauss - Beta * w) / (
sum(dGauss) - Beta));
934 mat Xsub = X.get_cols(
getSamples(numSamples, sampleSize));
936 vec Usquared =
pow(u, 2);
937 vec ex =
exp(-a2 * Usquared / 2);
939 vec dGauss =
elem_mult(1 - a2 * Usquared, ex);
940 w = (Xsub * gauss -
sum(dGauss) * w) / Xsub.cols();
944 mat Xsub = X.get_cols(
getSamples(numSamples, sampleSize));
946 vec Usquared =
pow(u, 2);
947 vec ex =
exp(-a2 * Usquared / 2);
949 vec dGauss =
elem_mult(1 - a2 * Usquared, ex);
950 double Beta =
dot(w, Xsub * gauss);
951 w = w - myy * ((Xsub * gauss - Beta * w) / (
sum(dGauss) - Beta));
962 vec Gskew = X *
pow(Y, 2) / numSamples;
963 double Beta =
dot(w, Gskew);
964 w = w - myy * (Gskew - Beta * w / (-Beta));
968 mat Xsub = X.get_cols(
getSamples(numSamples, sampleSize));
969 w = (Xsub * (
pow(
transpose(Xsub) * w, 2))) / Xsub.cols();
973 mat Xsub = X.get_cols(
getSamples(numSamples, sampleSize));
974 vec Gskew = Xsub *
pow(
transpose(Xsub) * w, 2) / Xsub.cols();
975 double Beta =
dot(w, Gskew);
976 w = w - myy * (Gskew - Beta * w) / (-Beta);