print.cpp

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

#include <stdio.h>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <fstream>
using namespace std;



//// 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]);
		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);