84 static void selcol(
const mat oldMatrix,
const vec maskVector, mat & newMatrix);
85 static int 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 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);
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 if (
pcamat(mixedSigC, numOfIC, firstEig, lastEig, E, D) < 1) {
152 whitenv(mixedSigC, E,
diag(D), whitesig, whiteningMatrix, dewhiteningMatrix);
154 Dim = whitesig.rows();
155 if (numOfIC > Dim) numOfIC = Dim;
159 for (
int i = 0; i < NcFirst.size(); i++) {
162 NcVp(NcFirst(i)) = 0.0;
163 VecPr.set_col(i, dewhiteningMatrix.get_col(i));
168 if (PCAonly ==
false) {
170 result =
fpica(whitesig, whiteningMatrix, dewhiteningMatrix, approach, numOfIC, g, finetune, a1, a2, mu, stabilization, epsilon, maxNumIterations, maxFineTune, initState, guess, sampleSize, A, W);
172 icasig = W * mixedSig;
214 initGuess = ma_initGuess;
242 for (
int i = 0; i <
size(maskVector); i++)
if (maskVector(i) == 1) numTaken++;
244 newMatrix =
zeros(oldMatrix.rows(), numTaken);
248 for (
int i = 0; i <
size(maskVector); i++) {
250 if (maskVector(i) == 1) {
252 newMatrix.set_col(numTaken, oldMatrix.get_col(i));
262 static int pcamat(
const mat vectors,
const int numOfIC,
int firstEig,
int lastEig,
mat & Es,
vec & Ds)
269 double lowerLimitValue = 0.0,
270 higherLimitValue = 0.0;
272 int oldDimension = vectors.rows();
276 eig_sym(covarianceMatrix, Dt, Et);
281 for (
int i = 0; i < Dt.length(); i++)
if (Dt(i) >
FICA_TOL) maxLastEig++;
282 if (maxLastEig < 1)
return 0;
285 if (maxLastEig > numOfIC) maxLastEig = numOfIC;
296 for (
int i = 0; i <
size(Dt); i++) eigenvalues(i) = eigenvalues2(
size(Dt) - i - 1);
298 if (lastEig > maxLastEig) lastEig = maxLastEig;
300 if (lastEig < oldDimension) lowerLimitValue = (eigenvalues(lastEig - 1) + eigenvalues(lastEig)) / 2;
301 else lowerLimitValue = eigenvalues(oldDimension - 1) - 1;
303 for (
int i = 0; i <
size(Dt); i++)
if (Dt(i) > lowerLimitValue) lowerColumns(i) = 1;
305 if (firstEig > 1) higherLimitValue = (eigenvalues(firstEig - 2) + eigenvalues(firstEig - 1)) / 2;
306 else higherLimitValue = eigenvalues(0) + 1;
309 for (
int i = 0; i <
size(Dt); i++)
if (Dt(i) < higherLimitValue) higherColumns(i) = 1;
312 for (
int i = 0; i <
size(Dt); i++) selectedColumns(i) = (lowerColumns(i) == 1 && higherColumns(i) == 1) ? 1 : 0;
314 selcol(Et, selectedColumns, Es);
318 for (
int i = 0; i <
size(selectedColumns); i++)
if (selectedColumns(i) == 1) numTaken++;
320 Ds =
zeros(numTaken);
324 for (
int i = 0; i <
size(Dt); i++)
325 if (selectedColumns(i) == 1) {
326 Ds(numTaken) = Dt(i);
338 outVectors =
zeros(inVectors.rows(), inVectors.cols());
339 meanValue =
zeros(inVectors.rows());
341 for (
int i = 0; i < inVectors.rows(); i++) {
343 meanValue(i) =
mean(inVectors.get_row(i));
345 for (
int j = 0; j < inVectors.cols(); j++) outVectors(i, j) = inVectors(i, j) - meanValue(i);
354 whiteningMatrix =
zeros(E.cols(), E.rows());
355 dewhiteningMatrix =
zeros(E.rows(), E.cols());
357 for (
int i = 0; i < D.cols(); i++) {
359 dewhiteningMatrix.set_col(i,
std::sqrt(D(i, i))*E.get_col(i));
362 newVectors = whiteningMatrix * vectors;
374 double eps = 2.2e-16;
380 if (A.rows() > A.cols()) {
382 U = U(0, U.rows() - 1, 0, A.cols() - 1);
383 S = S(0, A.cols() - 1);
386 mmax = (A.rows() > A.cols()) ? A.rows() : A.cols();
388 tol = mmax *
eps *
max(S);
390 for (
int i = 0; i <
size(S); i++)
if (S(i) > tol) r++;
392 Q = U(0, U.rows() - 1, 0, r - 1);
411 for (
int i = 0; i < T.cols(); i++) T.set_col(i, T.get_col(i) /
norm(T.get_col(i)));
425 for (
int i = 0; i <
max; i++)
if (rd(i) < percentage) { sV.add_elem(sZ, i); sZ++; }
438 for (
int i = 0; i < A.cols(); i++) { out(i) =
sum(A.get_col(i)); }
444 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)
447 int vectorSize = X.rows();
448 int numSamples = X.cols();
450 int gFine = finetune + 1;
451 double myyOrig = myy;
453 int failureLimit = 5;
454 int usedNlinearity = 0;
458 int initialStateMode = initState;
459 double minAbsCos = 0.0, minAbsCos2 = 0.0;
461 if (sampleSize * numSamples < 1000) sampleSize = (1000 / (double)numSamples < 1.0) ? 1000 / (double)numSamples : 1.0;
463 if (sampleSize != 1.0) gOrig += 2;
464 if (myy != 1.0) gOrig += 1;
466 int fineTuningEnabled = 1;
469 if (myy != 1.0) gFine = gOrig;
470 else gFine = gOrig + 1;
471 fineTuningEnabled = 0;
474 int stabilizationEnabled = stabilization;
476 if (!stabilization && myy != 1.0) stabilizationEnabled =
true;
478 usedNlinearity = gOrig;
480 if (initState ==
FICA_INIT_GUESS && guess.rows() != whiteningMatrix.cols()) {
481 initialStateMode = 0;
484 else if (guess.cols() < numOfIC) {
486 mat guess2 =
randu(guess.rows(), numOfIC - guess.cols()) - 0.5;
489 else if (guess.cols() > numOfIC) guess = guess(0, guess.rows() - 1, 0, numOfIC - 1);
493 usedNlinearity = gOrig;
498 A =
zeros(vectorSize, numOfIC);
501 if (initialStateMode == 0) B =
orth(
randu(vectorSize, numOfIC) - 0.5);
502 else B = whiteningMatrix * guess;
504 mat BOld =
zeros(B.rows(), B.cols());
505 mat BOld2 =
zeros(B.rows(), B.cols());
509 if (
round == maxNumIterations - 1) {
515 A = dewhiteningMatrix * B;
526 if (1 - minAbsCos < epsilon) {
528 if (fineTuningEnabled && notFine) {
531 usedNlinearity = gFine;
532 myy = myyK * myyOrig;
533 BOld =
zeros(B.rows(), B.cols());
534 BOld2 =
zeros(B.rows(), B.cols());
540 A = dewhiteningMatrix * B;
547 else if (stabilizationEnabled) {
548 if (!stroke && (1 - minAbsCos2 < epsilon)) {
552 if (
mod(usedNlinearity, 2) == 0) usedNlinearity += 1 ;
559 if (myy == 1 &&
mod(usedNlinearity, 2) != 0) usedNlinearity -= 1;
562 else if (!loong && (
round > maxNumIterations / 2)) {
566 if (
mod(usedNlinearity, 2) == 0) usedNlinearity += 1;
575 switch (usedNlinearity) {
579 B = (X *
pow(
transpose(X) * B, 3)) / numSamples - 3 * B;
586 mat D =
diag(
pow(Beta - 3 * numSamples , -1));
592 B = (Xsub *
pow(
transpose(Xsub) * B, 3)) / Xsub.cols() - 3 * B;
599 mat D =
diag(
pow(Beta - 3 * Ysub.rows() , -1));
600 B = B + myy * B * (
transpose(Ysub) * Gpow3 -
diag(Beta)) * D;
607 B = (X * hypTan) / numSamples -
elem_mult(
reshape(repeat(
sumcol(1 -
pow(hypTan, 2)), B.rows()), B.rows(), B.cols()), B) / numSamples * a1;
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;
631 B = B + myy * B * (
transpose(Ysub) * hypTan -
diag(Beta)) * D;
639 mat ex =
exp(-a2 * Usquared / 2);
642 B = (X * gauss) / numSamples -
elem_mult(
reshape(repeat(
sumcol(dGauss), B.rows()), B.rows(), B.cols()), B) / numSamples;
659 mat ex =
exp(-a2 * Usquared / 2);
662 B = (Xsub * gauss) / Xsub.cols() -
elem_mult(
reshape(repeat(
sumcol(dGauss), B.rows()), B.rows(), B.cols()), B) / Xsub.cols();
672 B = B + myy * B * (
transpose(Ysub) * gauss -
diag(Beta)) * D;
691 B = (Xsub * (
pow(
transpose(Xsub) * B, 2))) / Xsub.cols();
699 B = B + myy * B * (
transpose(Ysub) * Gskew -
diag(Beta)) * D;
717 A =
zeros(whiteningMatrix.cols(), numOfIC);
724 while (
round <= numOfIC) {
728 usedNlinearity = gOrig;
733 int endFinetuning = 0;
737 if (initialStateMode == 0)
739 w =
randu(vectorSize) - 0.5;
741 else w = whiteningMatrix * guess.get_col(
round);
753 while (i <= maxNumIterations + gabba) {
761 if (i == maxNumIterations + 1) {
767 if (numFailures > failureLimit) {
771 A = dewhiteningMatrix * B;
786 else if (i >= endFinetuning) wOld = w;
788 if (
norm(w - wOld) < epsilon ||
norm(w + wOld) < epsilon) {
790 if (fineTuningEnabled && notFine) {
794 wOld =
zeros(vectorSize);
795 wOld2 =
zeros(vectorSize);
796 usedNlinearity = gFine;
797 myy = myyK * myyOrig;
798 endFinetuning = maxFinetune + i;
806 B.set_col(
round - 1, w);
808 A.set_col(
round - 1, dewhiteningMatrix*w);
818 else if (stabilizationEnabled) {
820 if (stroke == 0.0 && (
norm(w - wOld2) < epsilon ||
norm(w + wOld2) < epsilon)) {
825 if (
mod(usedNlinearity, 2) == 0) {
833 else if (stroke != 0.0) {
838 if (myy == 1 &&
mod(usedNlinearity, 2) != 0) {
844 else if (notFine && !loong && i > maxNumIterations / 2) {
849 if (
mod(usedNlinearity, 2) == 0) {
863 switch (usedNlinearity) {
867 w = (X *
pow(
transpose(X) * w, 3)) / numSamples - 3 * w;
872 vec Gpow3 = X *
pow(Y, 3) / numSamples;
873 double Beta =
dot(w, Gpow3);
874 w = w - myy * (Gpow3 - Beta * w) / (3 - Beta);
879 w = (Xsub *
pow(
transpose(Xsub) * w, 3)) / Xsub.cols() - 3 * w;
885 double Beta =
dot(w, Gpow3);
886 w = w - myy * (Gpow3 - Beta * w) / (3 - Beta);
893 w = (X * hypTan - a1 *
sum(1 -
pow(hypTan, 2)) * w) / numSamples;
899 double Beta =
dot(w, X * hypTan);
900 w = w - myy * ((X * hypTan - Beta * w) / (a1 *
sum(1 -
pow(hypTan, 2)) - Beta));
906 w = (Xsub * hypTan - a1 *
sum(1 -
pow(hypTan, 2)) * w) / Xsub.cols();
912 double Beta =
dot(w, Xsub * hypTan);
913 w = w - myy * ((Xsub * hypTan - Beta * w) / (a1 *
sum(1 -
pow(hypTan, 2)) - Beta));
921 vec ex =
exp(-a2 * Usquared / 2);
924 w = (X * gauss -
sum(dGauss) * w) / numSamples;
931 vec ex =
exp(-a2 * Usquared / 2);
934 double Beta =
dot(w, X * gauss);
935 w = w - myy * ((X * gauss - Beta * w) / (
sum(dGauss) - Beta));
942 vec ex =
exp(-a2 * Usquared / 2);
945 w = (Xsub * gauss -
sum(dGauss) * w) / Xsub.cols();
952 vec ex =
exp(-a2 * Usquared / 2);
955 double Beta =
dot(w, Xsub * gauss);
956 w = w - myy * ((Xsub * gauss - Beta * w) / (
sum(dGauss) - Beta));
967 vec Gskew = X *
pow(Y, 2) / numSamples;
968 double Beta =
dot(w, Gskew);
969 w = w - myy * (Gskew - Beta * w / (-Beta));
974 w = (Xsub * (
pow(
transpose(Xsub) * w, 2))) / Xsub.cols();
980 double Beta =
dot(w, Gskew);
981 w = w - myy * (Gskew - Beta * w) / (-Beta);