matrix.h

fastica C++实现,很简单的,很容易看的

#include "stdafx.h"

//// allocate memory for matrix with r row and c column
//
double** alloc(int r, int c)
{
	double** a=new double*[r];
	for(int i=0;i<r;i++) a[i]=new double[c];
	return a;
}

//// allocate memory for square matrix with n row and n column
//
double** alloc(int n)
{
	return alloc(n,n);
}


////  get sample data from filename with c variables and r samples,
////  return data for matrix with r rows (number of samples) 
////  and c columns (number of variables)
//
double** readsample(char* filename, int r, int c)
{
	ifstream fin(filename);
	double** a=alloc(r,c);
	for(int i=0;i<r;i++) for(int j=0;j<c;j++) fin>>a[i][j];
	fin.close();
	return a;
}

//// write matrix a with r rows and c columns to filename 
//
void write(char* filename, double** a, int r, int c)
{
	FILE *stream;
	stream=fopen(filename,"w");
	for(int i=0;i<r;i++) 
	{
		for(int j=0;j<c;j++) fprintf(stream,"%16.12f ",a[i][j]);//fout<<setw(10)<<a[i][j]<<"  ";
		fprintf(stream,"\r\n");
	}
	fclose(stream);
}

//// write square matrix a with r rows and r columns to filename 
//
void write(char* filename, double** a, int n)
{
	write(filename,a,n,n);
}

////  write array a with length n to filename 
//
void write(char* filename, double* a, int r)
{
	ofstream fout(filename);
	for(int i=0;i<r;i++) 
	{
		fout<<setw(6)<<a[i]<<"  ";
	}
	fout.close();
}

//// clear matrix a (rXc)
//
void clear(double** &a, int r, int c)
{
	for(int i=0;i<r;i++) delete[]a[i];
	delete[] a;
	a=NULL;
}

//// clear matrix a (nXn)
//
void clear(double** &a, int n)
{
	clear(a,n,n);
}

//// return matrix a add b (rXc)
// 
double** add(double** a,double** b, int r, int c)
{
	double** x=alloc(r,c);
	for(int i=0;i<r;i++)
		for(int j=0;j<c;j++)
			x[i][j]=a[i][j]+b[i][j];
	return x;
}

////  return the sum of a and b (nXn)
//
double** add(double** a,double** b, int n)
{
	return add(a,b,n,n);
}

////  a-b
double** sub(double** a,double** b, int r, int c)
{
	double** x=alloc(r,c);
	for(int i=0;i<r;i++)
		for(int j=0;j<c;j++)
			x[i][j]=a[i][j]-b[i][j];
	return x;
}

////  a-b
double** sub(double** a,double** b, int n)
{
	return sub(a,b,n,n);
}

////  a(rXc) * b(cXc2)
double** mult(double** a,double** b, int r, int c, int c2)
{
	double** x=alloc(r,c2);
	for(int i=0;i<r;i++)
	{
		for(int j=0;j<c2;j++)
		{
			x[i][j]=0;
			for(int k=0;k<c;k++)
			{
				x[i][j]+=a[i][k]*b[k][j];
			}
		}
	}
	return x;
}


////  a*b nXn
double** mult(double** a,double** b, int n)
{
	return mult(a,b,n,n,n);
}

////  a(const)*b rXc
double** mult(double a,double** b, int r, int c)
{
	double** x=alloc(r,c);
	for(int i=0;i<r;i++)
	{
		for(int j=0;j<c;j++)
		{
			x[i][j]=0;
			for(int k=0;k<c;k++)
			{
				x[i][j]+=a*b[k][j];
			}
		}
	}
	return x;
}

////  a(const)*b nXn
double** mult(double a, double** b,int n)
{
	return mult(a,b,n,n);
}


//// transpose a rXc
double** T(double** a,int r, int c)
{
	double** x=alloc(c,r);
	for(int i=0;i<c;i++)
		for(int j=0;j<r;j++)
			x[i][j]=a[j][i];
	return x;
}

double** T(double** a,int c)
{
	return T(a,c,c);
}

////  m: square matrix    n: the rank of matrix m
////  d: diagonal element   e: off diagonal element
//
void householder(double** m, int n, double *d, double *e)
{
	int l,k,j,i;
	double scale, hh, h, g, f;

	for(i=n-1;i>0;i--)
	{
		l=i-1;
		h=scale=0;
		if(l>0)
		{
			for(k=0;k<l+1;k++)
			{
				scale+=fabs(m[i][k]);
			}
			if(scale==0.0)
			{
				e[i]=m[i][l];
			}
			else
			{
				for(k=0;k<l+1;k++)
				{
					m[i][k]/=scale;					
					h+=m[i][k]*m[i][k];					
				}
				f=m[i][l];
				g=(f>=0.0?-sqrt(h):sqrt(h));
				e[i]=scale*g;
				h-=f*g;
				m[i][l]=f-g;
				f=0.0;

				for(j=0;j<l+1;j++)
				{
					m[j][i]=m[i][j]/h;
					g=0.0;
					for(k=0;k<j+1;k++)
					{
						g+=m[j][k]*m[i][k];
					}
					for(k=j+1;k<l+1;k++)
					{
						g+=m[k][j]*m[i][k];
					}
					e[j]=g/h;
					f+=e[j]*m[i][j];
				}
				hh=f/(h+h);
				for(j=0;j<l+1;j++)
				{
					f=m[i][j];
					e[j]=g=e[j]-hh*f;
					for(k=0;k<j+1;k++)
					{
						m[j][k]-=( f*e[k]+g*m[i][k] );
					}
				}
			}
		}
		else
		{
			e[i]=m[i][l];
		}
		d[i]=h;
	}

	d[0]=0.0;
	e[0]=0.0;

	for(i=0;i<n;i++)
	{
		l=i;
		if(d[i]!=0.0)
		{
			for(j=0;j<l;j++)
			{
				g=0.0;
				for(k=0;k<l;k++)
				{
					g+=m[i][k]*m[k][j];
				}
				for(k=0;k<l;k++)
				{
					m[k][j]-=g*m[k][i];
				}
			}
		}
		d[i]=m[i][i];
		m[i][i]=1.0;
		for(j=0;j<l;j++)
		{
			m[j][i]=m[i][j]=0.0;
		}
	}
}

////    z: eigenvector matrix    n:  the rank of matrix
////    d: diagonal element=> eigenvalues  e: off-diagonal element 
//
void ql(double** z, int n,  double* d, double* e)
{
	int m,l,iter,i,k;
	double s,r,p,g,f,dd,c,b;
	for(i=1;i<n;i++)	
		e[i-1]=e[i];	
	e[n-1]=0.0;
	for(l=0;l<n;l++)
	{
		iter=0;
		do
		{
			for(m=l;m<n-1;m++)
			{
				dd=fabs(d[m])+fabs(d[m+1]);
				if(fabs(e[m])+dd==dd)  break;
			}
			if(m!=l)
			{
				if(iter++==30) 	return;
				g=(d[l+1]-d[l])/(e[l]+e[l]);
				r=sqrt(g*g+1.0);
				g=d[m]-d[l]+e[l]/( g + g>0?(r>0?r:-r):(r>0?-r:r) );
				s=c=1.0;
				p=0.0;
				for(i=m-1;i>=l;i--)
				{
					f=s*e[i];
					b=c*e[i];
					e[i+1]=(r=sqrt(f*f+g*g));
					if(r==0.0)
					{
						d[i+1]-=p;
						e[m]=0.0;
						break;
					}
					s=f/r;  c=g/r;  g=d[i+1]-p;	r=(d[i]-g)*s+2.0*c*b;  d[i+1]=g+(p=s*r);  g=c*r-b;
					for(k=0;k<n;k++)
					{
						f=z[k][i+1];
						z[k][i+1]=s*z[k][i]+c*f;
						z[k][i]=c*z[k][i]-s*f;
					}
				}
				if(r==0.0 && i>=l)  continue;
				d[l]-=p;  e[l]=g;  e[m]=0.0;
			}
		}while(m!=l);
	}
}

////  make samples a(r:number of samples, c: number of variables) meanvalues be 0
void zeromean1(double** a, int r, int c)
{
	for(int j=0;j<c;j++)
	{
		double m=0;
		for(int i=0;i<r;i++)
		{
			m+=a[i][j];
		}
		m/=r;
		for(i=0;i<r;i++)
		{
			a[i][j]-=m;
		}
	}
}

////  return samples  with 0 meanvalues 
////  a: samples  r:  number of samples  c: number of variables
//
double** zeromean(double** a, int r, int c)
{
	double** x=alloc(r,c);
	for(int j=0;j<c;j++)
	{
		double m=0;
		for(int i=0;i<r;i++)
		{
			m+=a[i][j];
		}
		m/=r;
		for(i=0;i<r;i++)
		{
			x[i][j]=a[i][j]-m;
		}
	}
	return x;
}

////    make samples a(0 meanvalues) be normalized
////    r:number of samples, c: number of variables
//
void normalize1(double** a, int r, int c)
{
	for(int j=0;j<c;j++)
	{
		double v=0;
		for(int i=0;i<r;i++)
		{
			v+=a[i][j]*a[i][j];
		}
		v=sqrt(v/(r-1));
		for(i=0;i<r;i++)
		{
			a[i][j]/=v;
		}
	}
}
////    return normalized samples 
////    a: samples(0 meanvalues)  r: number of samples  c: variables
//
double** normalize(double** a, int r, int c)
{	
	double** x=alloc(r,c);
	for(int j=0;j<c;j++)
	{
		double v=0;
		for(int i=0;i<r;i++)
		{
			v+=a[i][j]*a[i][j];
		}
		v=sqrt(v/(r-1));
		for(i=0;i<r;i++)
		{
			x[i][j]=a[i][j]/v;
		}
	}
	return x;
}

////   return the covariance of samples a normalized
////    a: samples  r: number of samples  c: variables
//
double** cov(double** a, int r, int c)
{
	double** covar=alloc(c);
	for(int j=0;j<c;j++)
	{
		for(int k=0;k<c;k++)
		{
			covar[j][k]=0;
			for(int i=0;i<r;i++)
			{
				covar[j][k]+=a[i][j]*a[i][k];
			}
			covar[j][k]/=(r-1);
		}
	}
	return covar;
}

//// permute eigenvector matrix q and eigenvalues d by descend order
//
void permute(double** q, double* d, int c)
{
	for(int n=0;n<c-1;n++)
	{
		for(int m=n+1;m<c;m++)
		{
			if(d[n]<d[m] )
			{
				swap(d[n],d[m]);
				for(int k=0;k<c;k++) swap(q[k][n],q[k][m]);
			}
		}
	}
}

//// V=pow(D,-1/2)*T(E) 
//
double** vforwhiten(double** e, double* d, int n)
{
	double** v=alloc(n);
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			v[i][j]=e[j][i]/sqrt(d[i]);
		}
	}
	return v;
}

////  z=Vx  . But owing to the structure of the samples, the right is z=x*v'. v' denotes the transpose of v
////  v: cXc matrix   x: normlized samples ( r: samples  c: variables)
//  
double** whiten(double** v, double** x, int r, int c)
{
	double** vT=T(v,c,c);
	double** w=mult(x,vT,r,c,c);
	clear(vT,c);
	return w;
}

double** whitensignal(double** a,int r, int c)
{
	double** x=zeromean(a,r,c);
	normalize1(x,r,c);
	double** Rx=cov(x,r,c);  
	double* d=new double[c];
	double* e=new double[c];
	householder(Rx,c,d,e);
	ql(Rx,c,d,e);
	permute(Rx,d,c); 

	double** VV=vforwhiten(Rx,d,c);  // whitening matrix
	double** WW=whiten(VV,x,r,c);  // whiten signal

	clear(x,r,c);
	clear(Rx,c);
	delete[]d;
	delete[]e;
	clear(VV,c);
	return WW;
	///////////////////////////////////////////////////////////////////////////////////////////

}

double** pca(double** a,int r,int c,double** &E, double* &d)
{	
	double** x=zeromean(a,r,c);
	normalize1(x,r,c);
	E=cov(x,r,c);// Rx->E
	d=new double[c];
	double* e=new double[c];
	householder(E,c,d,e);
	ql(E,c,d,e);
	permute(E,d,c);
	double** YY=mult(x,E,r,c,c);  // PCA signal

	delete[]e;
	clear(x,r,c);

	return YY;
}


double** orthonormalize(double** a,int c)
{
	double** aT=T(a,c);
	aT=mult(a,aT,c);

	double* d=new double[c];
	double* e=new double[c];
	householder(aT,c,d,e);
	ql(aT,c,d,e);// aT=E
	permute(aT,d,c);

	double** ww_1_2=alloc(c);
	for(int i=0;i<c;i++)
	{
		for(int j=0;j<=i;j++)
		{
			ww_1_2[i][j]=0;
			for(int k=0;k<c;k++)
			{
				ww_1_2[i][j]+=aT[i][k]*aT[j][k]/sqrt(d[k]);
			}
			ww_1_2[j][i]=ww_1_2[i][j];
		}
	}

	double** w=mult(ww_1_2,a,c);
	
	clear(aT,c);
	clear(ww_1_2,c);
	delete[]d;
	delete[]e;

	return w;
}



double** orthonormalize(double** a,int r, int c)
{
	double** aT=T(a,r);
	double** E=mult(a,aT,r,c,r);

	double* d=new double[r];
	double* e=new double[r];
	householder(E,r,d,e);
	ql(E,r,d,e);// 
	permute(E,d,r);

	double** ww_1_2=alloc(r);
	for(int i=0;i<r;i++)
	{
		for(int j=0;j<r;j++)
		{
			ww_1_2[i][j]=0;
			for(int k=0;k<c;k++)
			{
				ww_1_2[i][j]+=E[i][k]*E[j][k]/sqrt(d[k]);
			}
		}
	}

	double** w=mult(ww_1_2,a,r,r,c);
	
	clear(aT,c,r);
	clear(ww_1_2,r);
	clear(E,r);
	delete[]d;
	delete[]e;

	return w;
}

double** createrand(int c)
{
	double** a=alloc(c);
	for(int i=0;i<c;i++)
	{
		for(int j=0;j<c;j++)
		{
			a[i][j]=(rand()%1000000)/1000000.0;
		}
	}
	return a;
}


double** createrand(int r,int c)
{
	double** a=alloc(r,c);
	for(int i=0;i<r;i++)
	{
		for(int j=0;j<c;j++)
		{
			a[i][j]=(rand()%1000000)/1000000.0;
		}
	}
	return a;
}


double** createidentical(int r,int c)
{
	double** a=alloc(r,c);
	for(int i=0;i<r;i++)
	{
		for(int j=0;j<c;j++)
		{
			 a[i][j]=0;
		}
	}
	for(i=0;i<r;i++) a[i][i]=1;

	return a;
}

double compare(double** newW,double** oldW, int r, int m)
{
	double d=0;
	double a=0;
	for(int i=0;i<r;i++)
	{
		for(int j=0;j<m;j++)
		{
			d+=(newW[i][j]-oldW[i][j])*(newW[i][j]-oldW[i][j]);
			a+=oldW[i][j]*oldW[i][j];
		}
	}
	return sqrt(d/a/r/m);
}

void setvalue(double** from,double** to, int c)
{
	for(int i=0;i<c;i++)
	{
		for(int j=0;j<c;j++)
		{
			to[i][j]=from[i][j];
		}
	}
}

void setvalue(double** from,double** to,int r, int c)
{
	for(int i=0;i<r;i++)
	{
		for(int j=0;j<c;j++)
		{
			to[i][j]=from[i][j];
		}
	}
}

void orthnor(double** a,int r, int c)
{
	double** b=orthonormalize(a,r,c);
	setvalue(b,a,r,c);
	clear(b,r,c);
}

double g(double x)
{
	//return x*x*x;
	return tanh(x);
}
double gg(double x)
{
	//return 3*x*x;
	double gp=tanh(x);
	return 1-gp*gp;
}

double** G(double** a, int r, int c)
{
	double** ga=alloc(r,c);
	for(int i=0;i<r;i++)
	{
		for(int j=0;j<c;j++)
		{
			ga[i][j]=g(a[i][j]);
		}
	}
	return ga;
}

double** GP(double** a, int r, int c)
{
	double** gga=alloc(r,c);
	for(int i=0;i<r;i++)
	{
		for(int j=0;j<c;j++)
		{
			gga[i][j]=gg(a[i][j]);
		}
	}
	return gga;
}

double** fastica(double** a, int r, int c, int m, double** &neww)
{
	double** z=whitensignal(a,r,c);//write("z.txt",z,r,c);
	double** oldw=createrand(c,m);
	//double** oldw=createidentical(c,m);
	neww=orthonormalize(oldw,c,m);
	//setvalue(neww,oldw,c,m);write("oldw.txt",oldw,c,m);

	int iter=0;double zelta=0;
	do
	{
		iter++;    //cout<<iter<<": ";
		setvalue(neww,oldw,c,m);

		double** A=mult(z,oldw,r,c,m);
		double** g=  G(A,r,m);
		double** gp=GP(A,r,m);

		//double* gpp=new double[m];
		
		for(int j=0;j<c;j++)//    wj = E{z*g(wjT*z)} -E{gp(wjT*z)}*wj
		{
			for(int l=0;l<m;l++)
			{
				neww[j][l]=0;
				for(int i=0;i<r;i++)
				{
					neww[j][l]+=(z[i][j]*g[i][l]-oldw[j][l]*gp[i][l]);
				}
				neww[j][l]/=r;  //cout<<neww[j][l]<<"|"<<oldw[j][l]<<"    ";
			}
		}

		clear(A,r,m);
		clear(g,r,m);
		clear(gp,r,m);

		//double** ww=orthonormalize(neww,c,m);
		//clear(neww,c,m);
		//neww=ww;
		orthnor(neww,c,m);

		zelta=compare(neww,oldw,c,m);   cout<<iter<<":"<<zelta<<endl;
		
	}while( zelta>1e-7 && iter<300 );// 

	double** res=mult(z,neww,r,c,m);

	clear(z,r,c);
	clear(oldw,c,m);

	return res;

}