IT++ Logo
fastica.cpp
Go to the documentation of this file.
00001 
00062 #include <itpp/signal/fastica.h>
00063 #include <itpp/signal/sigfun.h>
00064 #include <itpp/signal/resampling.h>
00065 #include <itpp/base/algebra/eigen.h>
00066 #include <itpp/base/algebra/svd.h>
00067 #include <itpp/base/math/trig_hyp.h>
00068 #include <itpp/base/matfunc.h>
00069 #include <itpp/base/random.h>
00070 #include <itpp/base/sort.h>
00071 #include <itpp/base/specmat.h>
00072 #include <itpp/base/svec.h>
00073 #include <itpp/base/math/min_max.h>
00074 #include <itpp/stat/misc_stat.h>
00075 
00076 
00077 using namespace itpp;
00078 
00079 
00084 static void selcol(const mat oldMatrix, const vec maskVector, mat & newMatrix);
00085 static void pcamat(const mat vectors, const int numOfIC, int firstEig, int lastEig, mat & Es, vec & Ds);
00086 static void remmean(mat inVectors, mat & outVectors, vec & meanValue);
00087 static void whitenv(const mat vectors, const mat E, const mat D, mat & newVectors, mat & whiteningMatrix, mat & dewhiteningMatrix);
00088 static mat orth(const mat A);
00089 static mat mpower(const mat A, const double y);
00090 static ivec getSamples(const int max, const double percentage);
00091 static vec sumcol(const mat A);
00092 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);
00095 namespace itpp
00096 {
00097 
00098 // Constructor, init default values
00099 Fast_ICA::Fast_ICA(mat ma_mixedSig)
00100 {
00101 
00102   // Init default values
00103   approach = FICA_APPROACH_SYMM;
00104   g = FICA_NONLIN_POW3;
00105   finetune = true;
00106   a1 = 1.0;
00107   a2 = 1.0;
00108   mu = 1.0;
00109   epsilon = 0.0001;
00110   sampleSize = 1.0;
00111   stabilization = false;
00112   maxNumIterations = 100000;
00113   maxFineTune = 100;
00114   firstEig = 1;
00115 
00116   mixedSig = ma_mixedSig;
00117 
00118   lastEig = mixedSig.rows();
00119   numOfIC = mixedSig.rows();
00120   PCAonly = false;
00121   initState = FICA_INIT_RAND;
00122 
00123 }
00124 
00125 // Call main function
00126 void Fast_ICA::separate(void)
00127 {
00128 
00129   int Dim = numOfIC;
00130 
00131   mat mixedSigC;
00132   vec mixedMean;
00133 
00134   mat guess;
00135   if (initState == FICA_INIT_RAND)
00136     guess = zeros(Dim, Dim);
00137   else
00138     guess = mat(initGuess);
00139 
00140   VecPr = zeros(mixedSig.rows(), numOfIC);
00141 
00142   icasig = zeros(numOfIC, mixedSig.cols());
00143 
00144   remmean(mixedSig, mixedSigC, mixedMean);
00145 
00146   pcamat(mixedSigC, numOfIC, firstEig, lastEig, E, D);
00147 
00148   whitenv(mixedSigC, E, diag(D), whitesig, whiteningMatrix, dewhiteningMatrix);
00149 
00150 
00151   ivec NcFirst = to_ivec(zeros(numOfIC));
00152   vec NcVp = D;
00153   for (int i = 0; i < NcFirst.size(); i++) {
00154 
00155     NcFirst(i) = max_index(NcVp);
00156     NcVp(NcFirst(i)) = 0.0;
00157     VecPr.set_col(i, dewhiteningMatrix.get_col(i));
00158 
00159   }
00160 
00161   if (PCAonly == false) {
00162 
00163     Dim = whitesig.rows();
00164 
00165     if (numOfIC > Dim) numOfIC = Dim;
00166 
00167     fpica(whitesig, whiteningMatrix, dewhiteningMatrix, approach, numOfIC, g, finetune, a1, a2, mu, stabilization, epsilon, maxNumIterations, maxFineTune, initState, guess, sampleSize, A, W);
00168 
00169     icasig = W * mixedSig;
00170 
00171   }
00172 
00173   else { // PCA only : returns E as IcaSig
00174     icasig = VecPr;
00175   }
00176 }
00177 
00178 void Fast_ICA::set_approach(int in_approach) { approach = in_approach; if (approach == FICA_APPROACH_DEFL) finetune = true; }
00179 
00180 void Fast_ICA::set_nrof_independent_components(int in_nrIC) { numOfIC = in_nrIC; }
00181 
00182 void Fast_ICA::set_non_linearity(int in_g) { g = in_g; }
00183 
00184 void Fast_ICA::set_fine_tune(bool in_finetune) { finetune = in_finetune; }
00185 
00186 void Fast_ICA::set_a1(double fl_a1) { a1 = fl_a1; }
00187 
00188 void Fast_ICA::set_a2(double fl_a2) { a2 = fl_a2; }
00189 
00190 void Fast_ICA::set_mu(double fl_mu) { mu = fl_mu; }
00191 
00192 void Fast_ICA::set_epsilon(double fl_epsilon) { epsilon = fl_epsilon; }
00193 
00194 void Fast_ICA::set_sample_size(double fl_sampleSize) { sampleSize = fl_sampleSize; }
00195 
00196 void Fast_ICA::set_stabilization(bool in_stabilization) { stabilization = in_stabilization; }
00197 
00198 void Fast_ICA::set_max_num_iterations(int in_maxNumIterations) { maxNumIterations = in_maxNumIterations; }
00199 
00200 void Fast_ICA::set_max_fine_tune(int in_maxFineTune) { maxFineTune = in_maxFineTune; }
00201 
00202 void Fast_ICA::set_first_eig(int in_firstEig) { firstEig = in_firstEig; }
00203 
00204 void Fast_ICA::set_last_eig(int in_lastEig) { lastEig = in_lastEig; }
00205 
00206 void Fast_ICA::set_pca_only(bool in_PCAonly) { PCAonly = in_PCAonly; }
00207 
00208 void Fast_ICA::set_init_guess(mat ma_initGuess)
00209 {
00210   initGuess = ma_initGuess;
00211   initState = FICA_INIT_GUESS;
00212 }
00213 
00214 mat Fast_ICA::get_mixing_matrix() { if (PCAonly) { it_warning("No ICA performed."); return (zeros(1, 1));} else return A; }
00215 
00216 mat Fast_ICA::get_separating_matrix() { if (PCAonly) { it_warning("No ICA performed."); return(zeros(1, 1)); } else return W; }
00217 
00218 mat Fast_ICA::get_independent_components() { if (PCAonly) { it_warning("No ICA performed."); return(zeros(1, 1)); } else return icasig; }
00219 
00220 int Fast_ICA::get_nrof_independent_components() { return numOfIC; }
00221 
00222 mat Fast_ICA::get_principal_eigenvectors() { return VecPr; }
00223 
00224 mat Fast_ICA::get_whitening_matrix() { return whiteningMatrix; }
00225 
00226 mat Fast_ICA::get_dewhitening_matrix() { return dewhiteningMatrix; }
00227 
00228 mat Fast_ICA::get_white_sig() { return whitesig; }
00229 
00230 } // namespace itpp
00231 
00232 
00233 static void selcol(const mat oldMatrix, const vec maskVector, mat & newMatrix)
00234 {
00235 
00236   int numTaken = 0;
00237 
00238   for (int i = 0; i < size(maskVector); i++) if (maskVector(i) == 1) numTaken++;
00239 
00240   newMatrix = zeros(oldMatrix.rows(), numTaken);
00241 
00242   numTaken = 0;
00243 
00244   for (int i = 0; i < size(maskVector); i++) {
00245 
00246     if (maskVector(i) == 1) {
00247 
00248       newMatrix.set_col(numTaken, oldMatrix.get_col(i));
00249       numTaken++;
00250 
00251     }
00252   }
00253 
00254   return;
00255 
00256 }
00257 
00258 static void pcamat(const mat vectors, const int numOfIC, int firstEig, int lastEig, mat & Es, vec & Ds)
00259 {
00260 
00261   mat Et;
00262   vec Dt;
00263   cmat Ec;
00264   cvec Dc;
00265   double lowerLimitValue = 0.0,
00266                            higherLimitValue = 0.0;
00267 
00268   int oldDimension = vectors.rows();
00269 
00270   mat covarianceMatrix = cov(transpose(vectors), 0);
00271 
00272   eig_sym(covarianceMatrix, Dt, Et);
00273 
00274   int maxLastEig = 0;
00275 
00276   // Compute rank
00277   for (int i = 0; i < Dt.length(); i++) if (Dt(i) > FICA_TOL) maxLastEig++;
00278 
00279   // Force numOfIC components
00280   if (maxLastEig > numOfIC) maxLastEig = numOfIC;
00281 
00282   vec eigenvalues = zeros(size(Dt));
00283   vec eigenvalues2 = zeros(size(Dt));
00284 
00285   eigenvalues2 = Dt;
00286 
00287   sort(eigenvalues2);
00288 
00289   vec lowerColumns = zeros(size(Dt));
00290 
00291   for (int i = 0; i < size(Dt); i++) eigenvalues(i) = eigenvalues2(size(Dt) - i - 1);
00292 
00293   if (lastEig > maxLastEig) lastEig = maxLastEig;
00294 
00295   if (lastEig < oldDimension) lowerLimitValue = (eigenvalues(lastEig - 1) + eigenvalues(lastEig)) / 2;
00296   else lowerLimitValue = eigenvalues(oldDimension - 1) - 1;
00297 
00298   for (int i = 0; i < size(Dt); i++) if (Dt(i) > lowerLimitValue) lowerColumns(i) = 1;
00299 
00300   if (firstEig > 1) higherLimitValue = (eigenvalues(firstEig - 2) + eigenvalues(firstEig - 1)) / 2;
00301   else higherLimitValue = eigenvalues(0) + 1;
00302 
00303   vec higherColumns = zeros(size(Dt));
00304   for (int i = 0; i < size(Dt); i++) if (Dt(i) < higherLimitValue) higherColumns(i) = 1;
00305 
00306   vec selectedColumns = zeros(size(Dt));
00307   for (int i = 0; i < size(Dt); i++) selectedColumns(i) = (lowerColumns(i) == 1 && higherColumns(i) == 1) ? 1 : 0;
00308 
00309   selcol(Et, selectedColumns, Es);
00310 
00311   int numTaken = 0;
00312 
00313   for (int i = 0; i < size(selectedColumns); i++) if (selectedColumns(i) == 1) numTaken++;
00314 
00315   Ds = zeros(numTaken);
00316 
00317   numTaken = 0;
00318 
00319   for (int i = 0; i < size(Dt); i++)
00320     if (selectedColumns(i) == 1) {
00321       Ds(numTaken) = Dt(i);
00322       numTaken++;
00323     }
00324 
00325   return;
00326 
00327 }
00328 
00329 
00330 static void remmean(mat inVectors, mat & outVectors, vec & meanValue)
00331 {
00332 
00333   outVectors = zeros(inVectors.rows(), inVectors.cols());
00334   meanValue = zeros(inVectors.rows());
00335 
00336   for (int i = 0; i < inVectors.rows(); i++) {
00337 
00338     meanValue(i) = mean(inVectors.get_row(i));
00339 
00340     for (int j = 0; j < inVectors.cols(); j++) outVectors(i, j) = inVectors(i, j) - meanValue(i);
00341 
00342   }
00343 
00344 }
00345 
00346 static void whitenv(const mat vectors, const mat E, const mat D, mat & newVectors, mat & whiteningMatrix, mat & dewhiteningMatrix)
00347 {
00348 
00349   whiteningMatrix = zeros(E.cols(), E.rows());
00350   dewhiteningMatrix = zeros(E.rows(), E.cols());
00351 
00352   for (int i = 0; i < D.cols(); i++) {
00353     whiteningMatrix.set_row(i, std::pow(std::sqrt(D(i, i)), -1)*E.get_col(i));
00354     dewhiteningMatrix.set_col(i, std::sqrt(D(i, i))*E.get_col(i));
00355   }
00356 
00357   newVectors = whiteningMatrix * vectors;
00358 
00359   return;
00360 
00361 }
00362 
00363 static mat orth(const mat A)
00364 {
00365 
00366   mat Q;
00367   mat U, V;
00368   vec S;
00369   double eps = 2.2e-16;
00370   double tol = 0.0;
00371   int mmax = 0;
00372   int r = 0;
00373 
00374   svd(A, U, S, V);
00375   if (A.rows() > A.cols()) {
00376 
00377     U = U(0, U.rows() - 1, 0, A.cols() - 1);
00378     S = S(0, A.cols() - 1);
00379   }
00380 
00381   mmax = (A.rows() > A.cols()) ? A.rows() : A.cols();
00382 
00383   tol = mmax * eps * max(S);
00384 
00385   for (int i = 0; i < size(S); i++) if (S(i) > tol) r++;
00386 
00387   Q = U(0, U.rows() - 1, 0, r - 1);
00388 
00389   return (Q);
00390 }
00391 
00392 static mat mpower(const mat A, const double y)
00393 {
00394 
00395   mat T = zeros(A.rows(), A.cols());
00396   mat dd = zeros(A.rows(), A.cols());
00397   vec d = zeros(A.rows());
00398   vec dOut = zeros(A.rows());
00399 
00400   eig_sym(A, d, T);
00401 
00402   dOut = pow(d, y);
00403 
00404   diag(dOut, dd);
00405 
00406   for (int i = 0; i < T.cols(); i++) T.set_col(i, T.get_col(i) / norm(T.get_col(i)));
00407 
00408   return (T*dd*transpose(T));
00409 
00410 }
00411 
00412 static ivec getSamples(const int max, const double percentage)
00413 {
00414 
00415   vec rd = randu(max);
00416   sparse_vec sV;
00417   ivec out;
00418   int sZ = 0;
00419 
00420   for (int i = 0; i < max; i++) if (rd(i) < percentage) { sV.add_elem(sZ, i); sZ++; }
00421 
00422   out = to_ivec(full(sV));
00423 
00424   return (out);
00425 
00426 }
00427 
00428 static vec sumcol(const mat A)
00429 {
00430 
00431   vec out = zeros(A.cols());
00432 
00433   for (int i = 0; i < A.cols(); i++) { out(i) = sum(A.get_col(i)); }
00434 
00435   return (out);
00436 
00437 }
00438 
00439 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)
00440 {
00441 
00442   int vectorSize = X.rows();
00443   int numSamples = X.cols();
00444   int gOrig = g;
00445   int gFine = finetune + 1;
00446   double myyOrig = myy;
00447   double myyK = 0.01;
00448   int failureLimit = 5;
00449   int usedNlinearity = 0;
00450   double stroke = 0.0;
00451   int notFine = 1;
00452   int loong = 0;
00453   int initialStateMode = initState;
00454   double minAbsCos = 0.0, minAbsCos2 = 0.0;
00455 
00456   if (sampleSize * numSamples < 1000) sampleSize = (1000 / (double)numSamples < 1.0) ? 1000 / (double)numSamples : 1.0;
00457 
00458   if (sampleSize != 1.0) gOrig += 2;
00459   if (myy != 1.0) gOrig += 1;
00460 
00461   int fineTuningEnabled = 1;
00462 
00463   if (!finetune) {
00464     if (myy != 1.0) gFine = gOrig;
00465     else gFine = gOrig + 1;
00466     fineTuningEnabled = 0;
00467   }
00468 
00469   int stabilizationEnabled = stabilization;
00470 
00471   if (!stabilization && myy != 1.0) stabilizationEnabled = true;
00472 
00473   usedNlinearity = gOrig;
00474 
00475   if (initState == FICA_INIT_GUESS && guess.rows() != whiteningMatrix.cols()) {
00476     initialStateMode = 0;
00477 
00478   }
00479   else if (guess.cols() < numOfIC) {
00480 
00481     mat guess2 = randu(guess.rows(), numOfIC - guess.cols()) - 0.5;
00482     guess = concat_horizontal(guess, guess2);
00483   }
00484   else if (guess.cols() > numOfIC) guess = guess(0, guess.rows() - 1, 0, numOfIC - 1);
00485 
00486   if (approach == FICA_APPROACH_SYMM) {
00487 
00488     usedNlinearity = gOrig;
00489     stroke = 0;
00490     notFine = 1;
00491     loong = 0;
00492 
00493     A = zeros(vectorSize, numOfIC);
00494     mat B = zeros(vectorSize, numOfIC);
00495 
00496     if (initialStateMode == 0) B = orth(randu(vectorSize, numOfIC) - 0.5);
00497     else B = whiteningMatrix * guess;
00498 
00499     mat BOld = zeros(B.rows(), B.cols());
00500     mat BOld2 = zeros(B.rows(), B.cols());
00501 
00502     for (int round = 0; round < maxNumIterations; round++) {
00503 
00504       if (round == maxNumIterations - 1) {
00505 
00506         // If there is a convergence problem,
00507         // we still want ot return something.
00508         // This is difference with original
00509         // Matlab implementation.
00510         A = dewhiteningMatrix * B;
00511         W = transpose(B) * whiteningMatrix;
00512 
00513         return;
00514       }
00515 
00516       B = B * mpower(transpose(B) * B , -0.5);
00517 
00518       minAbsCos = min(abs(diag(transpose(B) * BOld)));
00519       minAbsCos2 = min(abs(diag(transpose(B) * BOld2)));
00520 
00521       if (1 - minAbsCos < epsilon) {
00522 
00523         if (fineTuningEnabled && notFine) {
00524 
00525           notFine = 0;
00526           usedNlinearity = gFine;
00527           myy = myyK * myyOrig;
00528           BOld = zeros(B.rows(), B.cols());
00529           BOld2 = zeros(B.rows(), B.cols());
00530 
00531         }
00532 
00533         else {
00534 
00535           A = dewhiteningMatrix * B;
00536           break;
00537 
00538         }
00539 
00540       } // IF epsilon
00541 
00542       else if (stabilizationEnabled) {
00543         if (!stroke && (1 - minAbsCos2 < epsilon)) {
00544 
00545           stroke = myy;
00546           myy /= 2;
00547           if (mod(usedNlinearity, 2) == 0) usedNlinearity += 1 ;
00548 
00549         }
00550         else if (stroke) {
00551 
00552           myy = stroke;
00553           stroke = 0;
00554           if (myy == 1 && mod(usedNlinearity, 2) != 0) usedNlinearity -= 1;
00555 
00556         }
00557         else if (!loong && (round > maxNumIterations / 2)) {
00558 
00559           loong = 1;
00560           myy /= 2;
00561           if (mod(usedNlinearity, 2) == 0) usedNlinearity += 1;
00562 
00563         }
00564 
00565       } // stabilizationEnabled
00566 
00567       BOld2 = BOld;
00568       BOld = B;
00569 
00570       switch (usedNlinearity) {
00571 
00572         // pow3
00573       case FICA_NONLIN_POW3 : {
00574         B = (X * pow(transpose(X) * B, 3)) / numSamples - 3 * B;
00575         break;
00576       }
00577       case(FICA_NONLIN_POW3+1) : {
00578         mat Y = transpose(X) * B;
00579         mat Gpow3 = pow(Y, 3);
00580         vec Beta = sumcol(pow(Y, 4));
00581         mat D = diag(pow(Beta - 3 * numSamples , -1));
00582         B = B + myy * B * (transpose(Y) * Gpow3 - diag(Beta)) * D;
00583         break;
00584       }
00585       case(FICA_NONLIN_POW3+2) : {
00586         mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
00587         B = (Xsub * pow(transpose(Xsub) * B, 3)) / Xsub.cols() - 3 * B;
00588         break;
00589       }
00590       case(FICA_NONLIN_POW3+3) : {
00591         mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B;
00592         mat Gpow3 = pow(Ysub, 3);
00593         vec Beta = sumcol(pow(Ysub, 4));
00594         mat D = diag(pow(Beta - 3 * Ysub.rows() , -1));
00595         B = B + myy * B * (transpose(Ysub) * Gpow3 - diag(Beta)) * D;
00596         break;
00597       }
00598 
00599       // TANH
00600       case FICA_NONLIN_TANH : {
00601         mat hypTan = tanh(a1 * transpose(X) * B);
00602         B = (X * hypTan) / numSamples - elem_mult(reshape(repeat(sumcol(1 - pow(hypTan, 2)), B.rows()), B.rows(), B.cols()), B) / numSamples * a1;
00603         break;
00604       }
00605       case(FICA_NONLIN_TANH+1) : {
00606         mat Y = transpose(X) * B;
00607         mat hypTan = tanh(a1 * Y);
00608         vec Beta = sumcol(elem_mult(Y, hypTan));
00609         vec Beta2 = sumcol(1 - pow(hypTan, 2));
00610         mat D = diag(pow(Beta - a1 * Beta2 , -1));
00611         B = B + myy * B * (transpose(Y) * hypTan - diag(Beta)) * D;
00612         break;
00613       }
00614       case(FICA_NONLIN_TANH+2) : {
00615         mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
00616         mat hypTan = tanh(a1 * transpose(Xsub) * B);
00617         B = (Xsub * hypTan) / Xsub.cols() -  elem_mult(reshape(repeat(sumcol(1 - pow(hypTan, 2)), B.rows()), B.rows(), B.cols()), B) / Xsub.cols() * a1;
00618         break;
00619       }
00620       case(FICA_NONLIN_TANH+3) : {
00621         mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B;
00622         mat hypTan = tanh(a1 * Ysub);
00623         vec Beta = sumcol(elem_mult(Ysub, hypTan));
00624         vec Beta2 = sumcol(1 - pow(hypTan, 2));
00625         mat D = diag(pow(Beta - a1 * Beta2 , -1));
00626         B = B + myy * B * (transpose(Ysub) * hypTan - diag(Beta)) * D;
00627         break;
00628       }
00629 
00630       // GAUSS
00631       case FICA_NONLIN_GAUSS : {
00632         mat U = transpose(X) * B;
00633         mat Usquared = pow(U, 2);
00634         mat ex = exp(-a2 * Usquared / 2);
00635         mat gauss = elem_mult(U, ex);
00636         mat dGauss = elem_mult(1 - a2 * Usquared, ex);
00637         B = (X * gauss) / numSamples - elem_mult(reshape(repeat(sumcol(dGauss), B.rows()), B.rows(), B.cols()), B) / numSamples;
00638         break;
00639       }
00640       case(FICA_NONLIN_GAUSS+1) : {
00641         mat Y = transpose(X) * B;
00642         mat ex = exp(-a2 * pow(Y, 2) / 2);
00643         mat gauss = elem_mult(Y, ex);
00644         vec Beta = sumcol(elem_mult(Y, gauss));
00645         vec Beta2 = sumcol(elem_mult(1 - a2 * pow(Y, 2), ex));
00646         mat D = diag(pow(Beta - Beta2 , -1));
00647         B = B + myy * B * (transpose(Y) * gauss - diag(Beta)) * D;
00648         break;
00649       }
00650       case(FICA_NONLIN_GAUSS+2) : {
00651         mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
00652         mat U = transpose(Xsub) * B;
00653         mat Usquared = pow(U, 2);
00654         mat ex = exp(-a2 * Usquared / 2);
00655         mat gauss = elem_mult(U, ex);
00656         mat dGauss = elem_mult(1 - a2 * Usquared, ex);
00657         B = (Xsub * gauss) / Xsub.cols() - elem_mult(reshape(repeat(sumcol(dGauss), B.rows()), B.rows(), B.cols()), B) / Xsub.cols();
00658         break;
00659       }
00660       case(FICA_NONLIN_GAUSS+3) : {
00661         mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B;
00662         mat ex = exp(-a2 * pow(Ysub, 2) / 2);
00663         mat gauss = elem_mult(Ysub, ex);
00664         vec Beta = sumcol(elem_mult(Ysub, gauss));
00665         vec Beta2 = sumcol(elem_mult(1 - a2 * pow(Ysub, 2), ex));
00666         mat D = diag(pow(Beta - Beta2 , -1));
00667         B = B + myy * B * (transpose(Ysub) * gauss - diag(Beta)) * D;
00668         break;
00669       }
00670 
00671       // SKEW
00672       case FICA_NONLIN_SKEW : {
00673         B = (X * (pow(transpose(X) * B, 2))) / numSamples;
00674         break;
00675       }
00676       case(FICA_NONLIN_SKEW+1) : {
00677         mat Y = transpose(X) * B;
00678         mat Gskew = pow(Y, 2);
00679         vec Beta = sumcol(elem_mult(Y, Gskew));
00680         mat D = diag(pow(Beta , -1));
00681         B = B + myy * B * (transpose(Y) * Gskew - diag(Beta)) * D;
00682         break;
00683       }
00684       case(FICA_NONLIN_SKEW+2) : {
00685         mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
00686         B = (Xsub * (pow(transpose(Xsub) * B, 2))) / Xsub.cols();
00687         break;
00688       }
00689       case(FICA_NONLIN_SKEW+3) : {
00690         mat Ysub = transpose(X.get_cols(getSamples(numSamples, sampleSize))) * B;
00691         mat Gskew = pow(Ysub, 2);
00692         vec Beta = sumcol(elem_mult(Ysub, Gskew));
00693         mat D = diag(pow(Beta , -1));
00694         B = B + myy * B * (transpose(Ysub) * Gskew - diag(Beta)) * D;
00695         break;
00696       }
00697 
00698 
00699       } // SWITCH usedNlinearity
00700 
00701     } // FOR maxIterations
00702 
00703     W = transpose(B) * whiteningMatrix;
00704 
00705 
00706   } // IF FICA_APPROACH_SYMM APPROACH
00707 
00708   // DEFLATION
00709   else {
00710 
00711     // FC 01/12/05
00712     A = zeros(whiteningMatrix.cols(), numOfIC);
00713     //    A = zeros( vectorSize, numOfIC );
00714     mat B = zeros(vectorSize, numOfIC);
00715     W = transpose(B) * whiteningMatrix;
00716     int round = 1;
00717     int numFailures = 0;
00718 
00719     while (round <= numOfIC) {
00720 
00721       myy = myyOrig;
00722 
00723       usedNlinearity = gOrig;
00724       stroke = 0;
00725 
00726       notFine = 1;
00727       loong = 0;
00728       int endFinetuning = 0;
00729 
00730       vec w = zeros(vectorSize);
00731 
00732       if (initialStateMode == 0)
00733 
00734         w = randu(vectorSize) - 0.5;
00735 
00736       else w = whiteningMatrix * guess.get_col(round);
00737 
00738       w = w - B * transpose(B) * w;
00739 
00740       w /= norm(w);
00741 
00742       vec wOld = zeros(vectorSize);
00743       vec wOld2 = zeros(vectorSize);
00744 
00745       int i = 1;
00746       int gabba = 1;
00747 
00748       while (i <= maxNumIterations + gabba) {
00749 
00750         w = w - B * transpose(B) * w;
00751 
00752         w /= norm(w);
00753 
00754         if (notFine) {
00755 
00756           if (i == maxNumIterations + 1) {
00757 
00758             round--;
00759 
00760             numFailures++;
00761 
00762             if (numFailures > failureLimit) {
00763 
00764               if (round == 0) {
00765 
00766                 A = dewhiteningMatrix * B;
00767                 W = transpose(B) * whiteningMatrix;
00768 
00769               } // IF round
00770 
00771               break;
00772 
00773             } // IF numFailures > failureLimit
00774 
00775             break;
00776 
00777           } // IF i == maxNumIterations+1
00778 
00779         } // IF NOTFINE
00780 
00781         else if (i >= endFinetuning) wOld = w;
00782 
00783         if (norm(w - wOld) < epsilon || norm(w + wOld) < epsilon) {
00784 
00785           if (fineTuningEnabled && notFine) {
00786 
00787             notFine = 0;
00788             gabba = maxFinetune;
00789             wOld = zeros(vectorSize);
00790             wOld2 = zeros(vectorSize);
00791             usedNlinearity = gFine;
00792             myy = myyK * myyOrig;
00793             endFinetuning = maxFinetune + i;
00794 
00795           } // IF finetuning
00796 
00797           else {
00798 
00799             numFailures = 0;
00800 
00801             B.set_col(round - 1, w);
00802 
00803             A.set_col(round - 1, dewhiteningMatrix*w);
00804 
00805             W.set_row(round - 1, transpose(whiteningMatrix)*w);
00806 
00807             break;
00808 
00809           } // ELSE finetuning
00810 
00811         } // IF epsilon
00812 
00813         else if (stabilizationEnabled) {
00814 
00815           if (stroke == 0.0 && (norm(w - wOld2) < epsilon || norm(w + wOld2) < epsilon)) {
00816 
00817             stroke = myy;
00818             myy /= 2.0 ;
00819 
00820             if (mod(usedNlinearity, 2) == 0) {
00821 
00822               usedNlinearity++;
00823 
00824             } // IF MOD
00825 
00826           }// IF !stroke
00827 
00828           else if (stroke != 0.0) {
00829 
00830             myy = stroke;
00831             stroke = 0.0;
00832 
00833             if (myy == 1 && mod(usedNlinearity, 2) != 0) {
00834               usedNlinearity--;
00835             }
00836 
00837           } // IF Stroke
00838 
00839           else if (notFine && !loong && i > maxNumIterations / 2) {
00840 
00841             loong = 1;
00842             myy /= 2.0;
00843 
00844             if (mod(usedNlinearity, 2) == 0) {
00845 
00846               usedNlinearity++;
00847 
00848             } // IF MOD
00849 
00850           } // IF notFine
00851 
00852         } // IF stabilization
00853 
00854 
00855         wOld2 = wOld;
00856         wOld = w;
00857 
00858         switch (usedNlinearity) {
00859 
00860           // pow3
00861         case FICA_NONLIN_POW3 : {
00862           w = (X * pow(transpose(X) * w, 3)) / numSamples - 3 * w;
00863           break;
00864         }
00865         case(FICA_NONLIN_POW3+1) : {
00866           vec Y = transpose(X) * w;
00867           vec Gpow3 = X * pow(Y, 3) / numSamples;
00868           double Beta = dot(w, Gpow3);
00869           w = w - myy * (Gpow3 - Beta * w) / (3 - Beta);
00870           break;
00871         }
00872         case(FICA_NONLIN_POW3+2) : {
00873           mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
00874           w = (Xsub * pow(transpose(Xsub) * w, 3)) / Xsub.cols() - 3 * w;
00875           break;
00876         }
00877         case(FICA_NONLIN_POW3+3) : {
00878           mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
00879           vec Gpow3 = Xsub * pow(transpose(Xsub) * w, 3) / (Xsub.cols());
00880           double Beta = dot(w, Gpow3);
00881           w = w - myy * (Gpow3 - Beta * w) / (3 - Beta);
00882           break;
00883         }
00884 
00885         // TANH
00886         case FICA_NONLIN_TANH : {
00887           vec hypTan = tanh(a1 * transpose(X) * w);
00888           w = (X * hypTan - a1 * sum(1 - pow(hypTan, 2)) * w) / numSamples;
00889           break;
00890         }
00891         case(FICA_NONLIN_TANH+1) : {
00892           vec Y = transpose(X) * w;
00893           vec hypTan = tanh(a1 * Y);
00894           double Beta = dot(w, X * hypTan);
00895           w = w - myy * ((X * hypTan - Beta * w) / (a1 * sum(1 - pow(hypTan, 2)) - Beta));
00896           break;
00897         }
00898         case(FICA_NONLIN_TANH+2) : {
00899           mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
00900           vec hypTan = tanh(a1 * transpose(Xsub) * w);
00901           w = (Xsub * hypTan - a1 * sum(1 - pow(hypTan, 2)) * w) / Xsub.cols();
00902           break;
00903         }
00904         case(FICA_NONLIN_TANH+3) : {
00905           mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
00906           vec hypTan = tanh(a1 * transpose(Xsub) * w);
00907           double Beta = dot(w, Xsub * hypTan);
00908           w = w - myy * ((Xsub * hypTan - Beta * w) / (a1 * sum(1 - pow(hypTan, 2)) - Beta));
00909           break;
00910         }
00911 
00912         // GAUSS
00913         case FICA_NONLIN_GAUSS : {
00914           vec u = transpose(X) * w;
00915           vec Usquared = pow(u, 2);
00916           vec ex = exp(-a2 * Usquared / 2);
00917           vec gauss = elem_mult(u, ex);
00918           vec dGauss = elem_mult(1 - a2 * Usquared, ex);
00919           w = (X * gauss - sum(dGauss) * w) / numSamples;
00920           break;
00921         }
00922         case(FICA_NONLIN_GAUSS+1) : {
00923           vec u = transpose(X) * w;
00924           vec Usquared = pow(u, 2);
00925 
00926           vec ex = exp(-a2 * Usquared / 2);
00927           vec gauss = elem_mult(u, ex);
00928           vec dGauss = elem_mult(1 - a2 * Usquared, ex);
00929           double Beta = dot(w, X * gauss);
00930           w = w - myy * ((X * gauss - Beta * w) / (sum(dGauss) - Beta));
00931           break;
00932         }
00933         case(FICA_NONLIN_GAUSS+2) : {
00934           mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
00935           vec u = transpose(Xsub) * w;
00936           vec Usquared = pow(u, 2);
00937           vec ex = exp(-a2 * Usquared / 2);
00938           vec gauss = elem_mult(u, ex);
00939           vec dGauss = elem_mult(1 - a2 * Usquared, ex);
00940           w = (Xsub * gauss - sum(dGauss) * w) / Xsub.cols();
00941           break;
00942         }
00943         case(FICA_NONLIN_GAUSS+3) : {
00944           mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
00945           vec u = transpose(Xsub) * w;
00946           vec Usquared = pow(u, 2);
00947           vec ex = exp(-a2 * Usquared / 2);
00948           vec gauss = elem_mult(u, ex);
00949           vec dGauss = elem_mult(1 - a2 * Usquared, ex);
00950           double Beta = dot(w, Xsub * gauss);
00951           w = w - myy * ((Xsub * gauss - Beta * w) / (sum(dGauss) - Beta));
00952           break;
00953         }
00954 
00955         // SKEW
00956         case FICA_NONLIN_SKEW : {
00957           w = (X * (pow(transpose(X) * w, 2))) / numSamples;
00958           break;
00959         }
00960         case(FICA_NONLIN_SKEW+1) : {
00961           vec Y = transpose(X) * w;
00962           vec Gskew = X * pow(Y, 2) / numSamples;
00963           double Beta = dot(w, Gskew);
00964           w = w - myy * (Gskew - Beta * w / (-Beta));
00965           break;
00966         }
00967         case(FICA_NONLIN_SKEW+2) : {
00968           mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
00969           w = (Xsub * (pow(transpose(Xsub) * w, 2))) / Xsub.cols();
00970           break;
00971         }
00972         case(FICA_NONLIN_SKEW+3) : {
00973           mat Xsub = X.get_cols(getSamples(numSamples, sampleSize));
00974           vec Gskew = Xsub * pow(transpose(Xsub) * w, 2) / Xsub.cols();
00975           double Beta = dot(w, Gskew);
00976           w = w - myy * (Gskew - Beta * w) / (-Beta);
00977           break;
00978         }
00979 
00980         } // SWITCH nonLinearity
00981 
00982         w /= norm(w);
00983         i++;
00984 
00985       } // WHILE i<= maxNumIterations+gabba
00986 
00987       round++;
00988 
00989     } // While round <= numOfIC
00990 
00991   } // ELSE Deflation
00992 
00993 } // FPICA
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Sat Jul 9 2011 15:21:32 for IT++ by Doxygen 1.7.4