print.cpp

fastica 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);
	double** oldw=createrand(c,m);
	neww=orthonormalize(oldw,c,m);

	int iter=0;double zelta=0;
	do
	{
		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);
		
		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;  
			}
		}

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

		orthnor(neww,c,m);

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

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

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

	return res;

}

////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////


void testwhitensignal()
{
	int r=100,c=4;
	double** a=readsample("data.txt",r,c);
	write("inputdata输入数据.txt",a,r,c);

	double** x=zeromean(a,r,c);
	write("0mean0均值数据.txt",x,r,c);

	double** x1=normalize(x,r,c);
	write("norm归一化数据.txt",x1,r,c);

	double** Rx=cov(x1,r,c);
	write("Rx归一化数据的关联矩阵.txt",Rx,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** EE=Rx;  // eigenvector matrix
	double* DD=d;    // eigenvalues
	write("EE关联矩阵的本征矢矩阵.txt",EE,c);
	write("DD关联矩阵的本征值.txt",DD,c);

	double** VV=vforwhiten(EE,DD,c);  // whitening matrix
	write("VV白化矩阵.txt",VV,c);

	double** WW=whiten(VV,x1,r,c);  // whiten signal
	write("WW白化矩阵白化的数据.txt",WW,r,c);


	double** YY=mult(x1,EE,r,c,c);  // PCA signal
	write("PCA数据.txt",YY,r,c);

	double** normYY=normalize(YY,r,c);  // normalized PCA signal == whiten signal
	write("PCA数据归一化.txt",normYY,r,c);
	//////////////////////////////////////////////////////////////////////////


	double** ws=whitensignal(a,r,c);
	write("whitesignal函数白化的数据.txt",ws,r,c);
	//////////////////////////////////////////////////////////////////////////

}

void testpca()
{
	int r=100,c=4;
	double** a=readsample("data.txt",r,c);
	double** E=alloc(c);
	double* d=new double[c];
	double** pcad=pca(a,r,c,E,d);
	write("pcadata.txt",pcad,r,c);
	write("pcaE.txt",E,c);
	write("pcad.txt",d,c);
}

void testorth()
{
	int r=10;
	int c=10;
	double** a=createrand(r,c);
	write("orth0.txt",a,r,c);
	write("orth.txt",orthonormalize(a,r,c),r,c);
}

void testfastica()
{	
	int r=100,c=4;
	double** a=readsample("data.txt",r,c);

	int m=4;
	double** neww=NULL;
	double** s=fastica(a,r,c,m,neww);
	write("ica_signall.txt",s,r,m);
	write("WWW.txt",neww,c,m);

}

int main(int argc, char* argv[])
{
	//testwhitensignal();
	//testpca();
	//testorth();
	testfastica();
	return 0;
}


*********************  data.txt   *************************
4.627238406926    3.752335619713    2.877432832500    -2.250438085748 
7.270536476463    9.558582591152    7.771576003051    -0.274960259558 
7.752143168529    9.978679112322    8.529400224987    -0.391417625482 
7.002878912580    6.802325546140    8.083828572249    -1.670058256068 
6.066774374083    3.159137667640    7.478184218484    -3.068620215399 
6.125636313564    2.591104461992    7.893293553817    -3.407747668953 
5.391177135975    1.731541797777    3.534565637893    -4.478836348608 
6.916235163349    5.728677423514    5.463465509035    -3.232781594479 
7.473377089811    6.896691891648    6.445771320731    -2.901337837596 
6.714479858882    4.184405680929    6.131414201699    -3.837492671446 
5.496903810874    0.156398111113    5.375611149057    -5.189244113704 
7.086125870699    2.598453336204    11.446580156735    -3.696912793347 
7.966155225095    4.959010300643    12.815739171367    -2.882667171027 
9.400158092841    9.015944946329    14.748716596833    -1.489835791646 
10.116324018096    10.944381967795    15.971003709309    -0.796974781002 
2.546648719997    5.426816006685    5.432583821384    2.058707621236 
-0.767029145663    -2.710974807355    -1.374097664215    -1.263436950498 
-1.497422923607    -5.105894887464    -1.593706824075    -2.005477065896 
-0.908433351689    -3.556981726227    -0.498023150712    -1.438360123031 
0.508716108666    0.452374918590    1.418952930592    -0.060251992966 
1.496813526808    3.141064915238    2.897302534961    0.864868008078 
3.216265203722    5.978644400175    9.099882705037    2.490882808796 
2.027207352691    2.041884103545    8.374224330016    1.171715765126 
1.166151238932    -0.970386833029    7.959563163193    0.138039136041 
1.583546775691    -0.214386552045    8.804186865746    0.334897204622 
3.059268286666    3.642623911654    10.686002902534    1.537247147465 
2.408953791063    3.043958709411    6.413837787457    0.556348446948 
2.603284660766    2.891288460112    6.966974957493    0.359087184818 
1.723558171665    -0.573407657295    6.420382256811    -0.976412715897 
0.895071220971    -3.977765642535    5.898348956521    -2.327373423246 
-5.726276180054    -7.460453782183    -3.923201854261    0.897989826239 
-4.162591719821    -3.881804506893    -2.108351449981    1.800981263924 
-0.418777603729    4.136549882831    5.863644755630    4.813992494714 
0.323807079973    5.055990773899    6.801421367771    4.755949116847 
-0.134569377199    2.276188109634    6.510705372549    3.428115931885 
-0.833277352261    -1.318591796258    5.952805662877    1.792820737635 
-0.549320875938    -2.056025497075    6.351671219139    1.075442031141 
-0.946465910226    -2.920621204212    2.039693025708    -0.384743557232 
1.115953398650    1.509698848573    4.169125773364    0.556601465930 
2.375275594880    3.455448898764    5.473953218472    0.641078901690 
2.353312381236    1.488800032725    5.477911628961    -0.604688134286 
1.798943874552    -2.136172671593    4.932005329508    -2.426509005559 
3.986408797273    0.408774592335    11.117835328696    -1.544389500769 
5.511094532452    2.923462562775    12.623125542296    -1.356800452202 
7.739015258823    7.512455386933    14.821536427130    -0.491262212340 
2.430681052502    6.862902247291    5.997026072735    3.257182386059 
2.773420613610    5.754053877407    6.298170735662    2.207326819049 
0.307381967313    -1.785065449489    -0.216370436559    -1.655925798285 
0.296363578608    -3.958210771851    -0.269847467170    -3.061698674422 
1.571472614503    -2.259655988268    0.966608021716    -3.171831388712 
3.756665989955    2.192208434342    3.119735455975    -2.355410748906 
5.614499761015    5.694554404670    4.954806246099    -1.843092312853 
8.203690438684    9.430173218649    11.538217510294    -0.569608277969 
7.749187750974    6.087720526013    11.087297759746    -2.303761505439 
7.438490284553    3.235393819729    10.796957608339    -3.852173566987 
8.298594199251    3.961814842891    11.696374148363    -4.182395971639 
10.233762939104    7.986730860257    13.691798518896    -3.385196037751 
12.112418026184    11.921616118208    15.653426436343    -2.587716757419 
10.821810238112    8.436152118815    10.464991194584    -4.898850835935 
10.315660866147    5.390976172290    10.091629559460    -6.361638748463 
2.686775418648    -1.533137021376    -0.857119559109    -4.442964198735 
3.022808769679    -1.862962692708    -0.334355012322    -4.928710740977 
4.555486645497    1.494579748491    1.412892498729    -4.148216659740 
8.318271797833    9.638052224014    9.423346876890    -1.067516551944 
9.084964624288    10.893619271661    10.460603554161    -0.913021981115 
8.528536894138    8.276371976066    10.202324625275    -2.012685590766 
7.520140041385    4.397486222329    9.519010547978    -3.496983387545 
7.316750218738    3.024651039680    9.666711923864    -4.111253096525 
6.382829331755    1.549614773136    5.103639877360    -5.394009538537 
7.964785100933    5.701487572796    7.084878618452    -4.102454222239 
8.788236369205    7.653769099403    8.329304924257    -3.515173242799 
8.261494706272    5.624733169545    8.243327252327    -4.228612810889 
7.029764932243    1.542627275445    7.470045574744    -5.602830226468 
8.376397781390    3.246975632549    13.295585590043    -4.360187571041 
9.016520563571    4.879672878885    14.422512317758    -3.791662807691 
3.482708662571    5.423111849518    5.906645752774    1.066607806251 
4.431340024068    8.044633405757    7.360166366076    1.988854414598 
4.078429823170    6.778834512967    7.517089485029    1.621919630662 
0.807904679767    -1.229806305607    0.753472672742    -1.657292341177 
-0.134754201991    -4.259813090849    0.322086582300    -2.610377127311 
0.184649760512    -3.515948792096    1.149244065036    -2.310197246890 
1.525874514331    0.271303257956    2.991915304046    -1.003964162758 
2.700557181215    3.527322287535    4.659013953358    0.113150577532 
4.691370031712    7.188385136915    11.135640966644    2.017240725899 
3.600601782370    3.557624220561    10.511452357657    0.804313694782 
2.568301970322    0.044370026664    9.929189731843    -0.391501034345 
2.700221143372    -0.041823639750    10.492404254610    -0.469951369668 
4.037916987368    3.416687154625    12.240645290594    0.605499571635 
5.522627973420    7.237329070506    14.113354524446    1.771598294520 
3.986716589395    3.914235345649    8.936267418847    -0.143210049666 
-3.701451505488    -2.560974597852    -1.892692577873    2.163861275760 
-4.650035524012    -6.306377538263    -2.529324017902    0.706550820811 
-4.598310715704    -7.146363747858    -2.192935131007    0.181327942920 
-3.225894428952    -4.121582549916    -0.563656773610    0.907291523621 
0.591992608461    4.138968532502    7.488120332960    4.008645624826 
1.596695178454    5.864578834440    8.693676474819    4.226874241820 
1.335155181728    3.694737117708    8.605354530384    3.109771493157 
0.575953187755    -0.062635760959    7.992357682933    1.427474482823 
0.591236433633    -1.587955113663    8.127731162877    0.454375431203 
-0.038360699740    -3.132712988198    3.588213478893    -1.225980734928