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);
void set_sample_size(double fl_sampleSize)
Set sample size.
Definition of FastICA (Independent Component Analysis) for IT++.
Various functions on vectors and matrices - header file.
#define FICA_APPROACH_SYMM
Use symmetric approach : compute all ICs at a time.
void set_a1(double fl_a1)
Set parameter.
Mat< double > mat
Default Matrix Type.
Mat< Num_T > concat_horizontal(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Horizontal concatenation of two matrices.
double randu(void)
Generates a random uniform (0,1) number.
void set_mu(double fl_mu)
Set parameter.
static void remmean(mat inVectors, mat &outVectors, vec &meanValue)
Local functions for FastICA.
Definitions of eigenvalue decomposition functions.
bool separate(void)
Explicit launch of main FastICA function.
#define FICA_NONLIN_SKEW
Use skew non-linearity.
mat get_whitening_matrix()
Get the whitening matrix.
void set_init_guess(mat ma_initGuess)
Set initial guess matrix instead of random (default)
Vec< std::complex< double > > cvec
Definition of complex<double> vector type.
double norm(const cvec &v)
Calculate the 2-norm: norm(v)=sqrt(sum(abs(v).^2))
Sparse_Vec< double > sparse_vec
Type definition of a double sparse vector.
Mat< Num_T > elem_mult(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Element wise multiplication of two matrices.
static int pcamat(const mat vectors, const int numOfIC, int firstEig, int lastEig, mat &Es, vec &Ds)
Local functions for FastICA.
int mod(int k, int n)
Calculates the modulus, i.e. the signed reminder after division.
void set_max_num_iterations(int in_maxNumIterations)
Set maximum number of iterations.
T sum(const Vec< T > &v)
Sum of all elements in the vector.
void set_epsilon(double fl_epsilon)
Set convergence parameter .
#define FICA_NONLIN_TANH
Use tanh(x) non-linearity.
void set_nrof_independent_components(int in_nrIC)
Set number of independent components to separate.
static ivec getSamples(const int max, const double percentage)
Local functions for FastICA.
#define FICA_NONLIN_GAUSS
Use Gaussian non-linearity.
Minimum and maximum functions on vectors and matrices.
Definitions of signal processing functions.
double mean(const vec &v)
The mean value.
ITPP_EXPORT double round(double x)
Round to nearest integer, return result in double.
void sort(Vec< T > &data, SORTING_METHOD method=INTROSORT)
Sort the data vector in increasing order.
#define FICA_APPROACH_DEFL
Use deflation approach : compute IC one-by-one in a Gram-Schmidt-like fashion.
int get_nrof_independent_components()
Get number of independent components.
ITPP_EXPORT vec zeros(int size)
A Double vector of zeros.
const double eps
Constant eps.
T min(const Vec< T > &in)
Minimum value of vector.
void set_last_eig(int in_lastEig)
Set last eigenvalue index to take into account.
int max_index(const Vec< T > &in)
Return the postion of the maximum element in the vector.
void set_pca_only(bool in_PCAonly)
If true, only perform Principal Component Analysis (default = false)
bool svd(const mat &A, vec &S)
Get singular values s of a real matrix A using SVD.
vec exp(const vec &x)
Exp of the elements of a vector x.
Mat< std::complex< double > > cmat
Default Complex Matrix Type.
T max(const Vec< T > &v)
Maximum value of vector.
Trigonometric and hyperbolic functions - header file.
vec pow(const double x, const vec &y)
Calculates x to the power of y (x^y)
mat cov(const mat &X, bool is_zero_mean)
Covariance matrix calculation.
bin abs(const bin &inbin)
absolute value of bin
#define FICA_NONLIN_POW3
Use x^3 non-linearity.
int size(const Vec< T > &v)
Length of vector.
mat get_separating_matrix()
Get separating matrix.
static vec sumcol(const mat A)
Local functions for FastICA.
Mat< T > full(const Sparse_Mat< T > &s)
Convert a sparse matrix s into its dense representation.
void set_stabilization(bool in_stabilization)
Set stabilization mode true or off.
static mat mpower(const mat A, const double y)
Local functions for FastICA.
void set_max_fine_tune(int in_maxFineTune)
Set maximum number of iterations for fine tuning.
void set_non_linearity(int in_g)
Set non-linearity.
Resampling functions - header file.
Miscellaneous statistics functions and classes - header file.
ivec to_ivec(const Vec< T > &v)
Converts a Vec<T> to ivec.
#define FICA_INIT_GUESS
Set predefined start for Fast_ICA.
Definitions of special vectors and matrices.
#define it_warning(s)
Display a warning message.
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.
mat get_principal_eigenvectors()
Get nrIC first columns of the de-whitening matrix.
mat get_independent_components()
Get separated signals.
bool eig_sym(const mat &A, vec &d, mat &V)
Calculates the eigenvalues and eigenvectors of a symmetric real matrix.
static void whitenv(const mat vectors, const mat E, const mat D, mat &newVectors, mat &whiteningMatrix, mat &dewhiteningMatrix)
Local functions for FastICA.
void transpose(const Mat< T > &m, Mat< T > &out)
Transposition of the matrix m returning the transposed matrix in out.
#define FICA_INIT_RAND
Set random start for Fast_ICA.
Mat< T > reshape(const Mat< T > &m, int rows, int cols)
Reshape the matrix into an rows*cols matrix.
Vec< int > ivec
Definition of integer vector type.
mat get_dewhitening_matrix()
Get the de-whitening matrix.
void set_a2(double fl_a2)
Set parameter.
vec tanh(const vec &x)
Tan hyperbolic function.
vec sqrt(const vec &x)
Square root of the elements.
Definitions of Singular Value Decompositions.
Mat< T > diag(const Vec< T > &v, const int K=0)
Create a diagonal matrix using vector v as its diagonal.
static mat orth(const mat A)
Local functions for FastICA.
mat get_white_sig()
Get whitened signals.
Definition of classes for random number generators.
Fast_ICA(mat ma_mixed_sig)
Constructor.
void set_first_eig(int in_firstEig)
Set first eigenvalue index to take into account.
Sparse Vector Class definitions.
Num_T dot(const Vec< Num_T > &v1, const Vec< Num_T > &v2)
Inner (dot) product of two vectors v1 and v2.
#define FICA_TOL
Eigenvalues of the covariance matrix lower than FICA_TOL are discarded for analysis.
mat get_mixing_matrix()
Get mixing matrix.
void set_fine_tune(bool in_finetune)
Set fine tuning.
void set_approach(int in_approach)
Set approach : FICA_APPROACH_DEFL or FICA_APPROACH_SYMM (default)
Vec< double > vec
Definition of double vector type.
static void selcol(const mat oldMatrix, const vec maskVector, mat &newMatrix)
Local functions for FastICA.