main.cpp

4维简单胞映射程序: 对动力系统做全局分析

				ioscm.read((char *)(&gg),sizeof(long));
				ioscm.read((char *)(&cel),sizeof(long));
				if(gg==i)
				{
					temp[n]=cel;					  
					iogcm1.seekg(cel*sizeof(gcmdata));
					iogcm1.read((char *)(&gcm),sizeof(gcmdata));
					gcm.gr=-3;
					iogcm1.seekp(cel*sizeof(gcmdata));
					iogcm1.write((char *)(&gcm),sizeof(gcmdata));
					n++;
				}				  
			}			
			/////////////////// to determine whether each cell's image cells are in the temp[n]
			long prenum,afnum;
			prenum=0;
			afnum=1;          
			while(afnum>prenum)  ///////// to locate all candidate cells
			{
				prenum=n;
				for(long k=0;k<prenum;k++)
				{
					long te=temp[k];					  
					iogcm1.seekg(te*sizeof(gcmdata));
					long Iz;
					iogcm1.read((char *)(&gcm),sizeof(gcmdata));
					Iz=gcm.Iz;
					long Loc2=gcm.CpLoc;
					for(long kk=0;kk<Iz;kk++)
					{
						long Loc=Loc2+kk;
						iocp1.seekg(Loc*sizeof(cpdata));
						long image;
						iocp1.read((char *)(&image),sizeof(long));						
						long ggimage;
						iogcm1.clear();
						iogcm1.seekg((image)*sizeof(gcmdata));
						iogcm1.read((char *)(&gcm),sizeof(gcmdata));
						ggimage=gcm.gr;						
						if(ggimage==0)
						{
							iogcm1.seekg(image*sizeof(gcmdata));
							iogcm1.read((char *)(&gcm),sizeof(gcmdata));
							gcm.gr=-3;
							iogcm1.seekp(image*sizeof(gcmdata));
							iogcm1.write((char *)(&gcm),sizeof(gcmdata));
							temp[n]=image;
							n++;
						}
                        //*
						iosgr.seekg(image*sizeof(long));
						long sgg;
						iosgr.read((char *)(&sgg),sizeof(long));
						if(sgg>0&&sgg!=i)
						{
							flag=0;
							kk=Iz;
							k=prenum;
							prenum=n;
						}
						
					}
				}
				afnum=n;			
			}
			///////////////////////////////////////     determine 
			if(flag!=0)
			{
				for(long k=0;k<n;k++)
				{
					long tee=temp[k];
					iogcm1.seekg(tee*sizeof(gcmdata));
					long Iz;
					iogcm1.read((char *)(&gcm),sizeof(gcmdata));
					Iz=gcm.Iz;
					long Loc2;
					Loc2=gcm.CpLoc;
					for(long as=0;as<Iz;as++)
					{
						long Loc=Loc2+as;
						iocp1.seekg(Loc*sizeof(cpdata));
						long iima;
						iocp1.read((char *)(&iima),sizeof(long));
						iogcm1.seekg(iima*sizeof(gcmdata)+2*sizeof(long));
						long ggg;
						iogcm1.read((char *)(&ggg),sizeof(long));
						if(ggg!=-3)
							flag=0;
					}				
				}
			}
			
			if(flag!=0)         //  发现 persistent group  //
			{
				gn++;
				for(long h=0;h<n;h++)
				{
					iogcm1.seekp(temp[h]*sizeof(gcmdata)+2*sizeof(long));
					iogcm1.write((char *)(&gn),sizeof(long));
					long celp=temp[h];
					iopersistent.write((char *)(&celp),sizeof(long));
					iopersistent.write((char *)(&gn),sizeof(long));					   
				}				
			}	
			else
			{
				for(long ui=0;ui<n;ui++)
				{					
					iogcm1.seekg(temp[ui]*sizeof(gcmdata));
					iogcm1.read((char *)(&gcm),sizeof(gcmdata));
					gcm.gr=0;
					iogcm1.seekp(temp[ui]*sizeof(gcmdata));
					iogcm1.write((char *)(&gcm),sizeof(gcmdata));					  
				}
			}
			/////////  每次完后，gcm中gr应该清零，否则影响下一个group的判断  ////////			
			
		}		
	}
	//*/
    ///////////////////////////////////////////////////////
	
	
	
}


void determine()
{	
	void sort(long a[],long b[],long n);
	CMatrix mat;
	CMatrix cst;	
	long *flag;
	flag=new long [max];
	for(long i=0;i<max;i++)
	{
		if(i%2000==0)
			cout<<i<<endl;		
		flag[i]=-1;
	}	
	long tem[500000];/////////保存persistent cells     =================================????????????
	long n=0;
	fstream iopersistent;
	iopersistent.open("persistent.dat",ios::in|ios::out|ios::binary);  //////////	
	fstream iogcm,iocp;
	iogcm.open("gcm.dat",ios::in|ios::out|ios::binary);
	iocp.open("cp.dat",ios::in|ios::out|ios::binary);
    gcmdata gcm;
	cpdata cp;	
	fstream iopre;
	iopre.open("pre.dat",ios::out|ios::binary|ios::in);
	long preloc;
	iogcm.seekg((max-1)*sizeof(gcmdata));
	long jzmax;
	iogcm.read((char *)(&gcm),sizeof(gcmdata));
	preloc=gcm.PreLoc;
	jzmax=gcm.Jz;
	long totalpre=jzmax+preloc+1;
	long *fzpre;
	fzpre=new long [totalpre];
	iopre.seekg(0);
	for(long ttr=0;ttr<totalpre;ttr++)
	{
		iopre.read((char *)(&fzpre[ttr]),sizeof(long));
	}
	//////////////////////////////////////////////	
	iogcm.seekg((max-1)*sizeof(gcmdata));	
	iogcm.read((char *)(&gcm),sizeof(gcmdata));
	long Izmax;
	Izmax=gcm.Iz;
	long cpclo;
	cpclo=gcm.CpLoc;
	long totalcl;
	totalcl=Izmax+cpclo+1;
	//////
	cpdata *fcpl;
	fcpl=new cpdata [totalcl];
	iocp.seekg(0);
	for(ttr=0;ttr<totalcl;ttr++)
	{
		iocp.read((char *)(&fcpl[ttr]),sizeof(cpdata));
	}
	////////////////////////////////////////////
	gcmdata *fgcm;
	fgcm=new gcmdata [max];
	iogcm.seekg(0);
	for(ttr=0;ttr<max;ttr++)
	{
		iogcm.read((char *)(&fgcm[ttr]),sizeof(gcmdata));
	}
	
	//////////////////////// 输出 结果 的文件 //////////////////////////////////////////////////
	
	
	//////////////初始化result
	//long res=(gn+1)*sizeof(long)+gn*sizeof(float);
	float fres[max][3];
	for(i=0;i<max;i++)
	{
		if(i%2000==0)
			cout<<i<<endl;
		long tttt=0;
		float tff=0.0;
		///ioresult.write((char *)(&tttt),sizeof(long));		
		for(long j=0;j<3;j++)
		{
			fres[i][j]=0.0;
		}		
	}
	fstream iodm;             ////// 保存 cell 的, domicle 的数量,做分组用  
    iodm.open("todomiledata.dat",ios::in|ios::out|ios::binary);
	dom dm;
	dm.n=0;
	dm.total=0;
	//fstream ioder;
	//ioder.open("derminedornot.dat",ios::in|ios::out|ios::binary);
	long *fder;
	fder=new long [max];
	///fstream iocw;  //////////成尾数   计算是用
	//iocw.open("chengwei.dat",ios::in|ios::out|ios::binary);
	for(i=0;i<max;i++)
	{   /////////////////所有的cell 数量, total number  persistent group 也当成0
		if(i%2000==0)
			cout<<i<<endl;
		iodm.write((char *)(&dm),sizeof(dom)); /////////先保存 
		long ffff=0;
		fder[i]=ffff;		
	}	
	/////////////////////////////////////////////////////////////////////////////////
	
	
	//////////////////////////////////////////////////////////////////////////////
	long *ftemp;

	ftemp=new long [max];
	long ftn=0;

	for(i=2;i<gn+1;i++)      ////////////1为sink cell 只要取得其吸引域就可以了，不要算vj pj
	{
		ftn=0;
		///////// 确定 此组的cell
		//iotemp.close();
		long pend;
		//cur=iopersistent.tellg();
		iopersistent.seekg(0,ios::end);
		pend=iopersistent.tellg();
		iopersistent.seekg(0);
		//long per=0;
		//////////保留persistent cells
		n=0;
		long dfj=-1;///////不同批次preimage 的分界线
		////////////每次迭代计算preimage开始的地方
		//long nc=0;  //////////transcient cells 分成 nc
		while(pend!=iopersistent.tellg())
		{
			long cee,grg;
			iopersistent.read((char *)(&cee),sizeof(long));
			iopersistent.read((char *)(&grg),sizeof(long));
			if(grg==i)
			{
				/////改变 flag 	
				long ff=i;
				float pp=1.0;
				flag[cee]=ff;
				tem[n]=cee;				
				n++;	
				////////////////////////////////////  persistent cell为已知胞  ////////////////////////
				fder[cee]=ff;	
                                fres[cee][i-1]=pp;
				////////domicle
				iodm.seekg(cee*sizeof(dom));
				iodm.read((char *)(&dm),sizeof(dom));
				dm.n=0; ///////     用零来区分是否是 persistent cells 
				dm.total=dm.total+i;
				iodm.seekp(cee*sizeof(dom));
				iodm.write((char *)(&dm),sizeof(dom));				
				/////////////////////////////////////////////////////////////////////////////////////
			}
		}
		///////  注意,保存到temp中的transcient cell 的文件必需 每次循环时清空
		//iotemp.close();
		//iotemp.open("temptemp.dat",ios::trunc|ios::in|ios::out|ios::binary);  //////////
		//// 
		
		long numb=0;
		long handnum=0; 
		cout<<"first level preimage"<<endl;
		for(long k=0;k<n;k++)
		{
			long cee;
			cee=tem[k];   ///////// persistent cells
			long Jz;
            gcm=fgcm[cee];			
			Jz=gcm.Jz;
			long loc2=gcm.PreLoc;
			for(long p=0;p<Jz;p++)
			{
				long loc=loc2+p;
				long ffg,iima;
				iima=fzpre[loc];
				long dis;
				if(iima<0)///////////////////////////不取进来
				{
					ffg=i;
				}
				else
				{
				dis=fder[iima];
				ffg=flag[iima];
				}
				////////////////////////////////////////
				//long ddd;
				
				///////////////////////////////////////
				if(ffg!=i&&dis!=i)
				{					
					//iotemp.write((char *)(&iima),sizeof(long));
					ftemp[ftn]=iima;
					ftn++;
					flag[iima]=i;					
					iodm.seekg(iima*sizeof(dom));
					iodm.read((char *)(&dm),sizeof(dom));
					dm.n++;
					dm.total=dm.total+i;
					iodm.seekp(iima*sizeof(dom));
					iodm.write((char *)(&dm),sizeof(dom));
					numb++;
				}				
			}
		}		  
		//iotemp.write((char *)(&dfj),sizeof(long));/////第一批与后面的分界
		ftemp[ftn]=-1;
		ftn++;
		///////////////////////////////迭代
		long afn,bfn;
		afn=0;
		bfn=-1;
		long sta=0;
		while(afn>bfn)
		{
			cout<<"preimage of persistent  an"<<afn<<"   bn "<<bfn<<endl;
			bfn=afn;			
			//iotemp.seekg(dbe*sizeof(long));
			long cee;			
			cee=ftemp[sta];
			//iotemp.read((char *)(&cee),sizeof(long));			
			sta++;		
			while(cee!=-1)
			{				
				long Jz;
                gcm=fgcm[cee];				
				Jz=gcm.Jz;
				long loc2=gcm.PreLoc;
				for(long p=0;p<Jz;p++)
				{
					long loc=loc2+p;
					long fg,ima;
					ima=fzpre[loc];
					fg=flag[ima];
					if(fg!=i)
					{		
						numb++;
						afn++;
						ftemp[ftn]=ima;
						ftn++;						
						flag[ima]=i;
						iodm.seekg(ima*sizeof(dom));
						iodm.read((char *)(&dm),sizeof(dom));
						dm.n++;
						dm.total=dm.total+i;
						iodm.seekp(ima*sizeof(dom));
						iodm.write((char *)(&dm),sizeof(dom));
					}					
				}				
				cee=ftemp[sta];
				sta++;				
			}
			///////// 次批完
            ftemp[ftn]=-1;
			ftn++;			
		}
		
		////////上面获得了 此persistent group的 所有pre象,并已经按批次存放,用-1分隔开来  iotemp中
		cout<<numb<<endl;
		long dtem[5000000];     ///////////      保存临时像/ 第i批数据/////////////?????????????????
		long nu=0;		
		long cet;
		sta=0;
		cet=ftemp[sta];	
		sta++;
		while(sta<ftn)
		{
			cout<<" get i'th cell of i'th preimage"<<endl;
			///////       获取第i批数据及像       ////////  
			nu=0;
			long dtm3;
			if(cet!=-1)       ////////文件指针不能移出去
			{
				dtm3=fder[cet];
			}
			if(cet!=-1&&dtm3!=i)   
			{
				dtem[nu]=cet;
				nu++;
			}			
			cet=ftemp[sta];
			sta++;
			while(cet!=-1)  //////////            取得此大批preimage
			{				
				long dtm;
				dtm=fder[cet];				
				if(dtm!=i)        /////////好像肯定不等于的??????????????
				{
					dtem[nu]=cet;			/////////已经处理也要放进取,避免它的下几层胞丢失	
					nu++;
				}				
				cet=ftemp[sta];
				sta++;
			}
            ///////////////////////////////////////////////////////////////
			////////////////  取像  只取此pg 的吸引域里面的///////////
            if(nu>0)
			{				
				long an,bn;
				an=0;
				bn=-1;
				long psta=0;
				while(an>bn)
				{
					cout<<" i'th preimage an "<<an<<"   bn "<<bn<<endl;
					bn=an;
					for(long u=psta;u<nu;u++)   //////////获得每个像
					{
						if(nu>10000)
						{
							if(u%500==0)
								cout<<u<<"           "<<nu<<endl;
						}
						long cre=dtem[u];
						psta++;
						long dt; ////////////// 考察是否已经确定
                        gcm=fgcm[cre];						
						long Iz;
						Iz=gcm.Iz;
						long loc2=gcm.CpLoc;
						for(long y=0;y<Iz;y++)
						{
							long loc=loc2+y;
							cp=fcpl[loc];							
							long cim;
							cim=cp.Cz;
							//////// whether determined
							dt=fder[cim];							
							long tfgg;
							tfgg=flag[cim];
							if(dt!=i&&tfgg==i) //////  防止边界上,其它域里面的胞也保存进去
							{
								//////看是否已经保存了
								long t=0;
								while(dtem[t]!=cim&&t<nu)
								{
									t++;
								}
								if(t==nu) /////////没有保存
								{
									dtem[nu]=cim;
									nu++;
									an++;
								}
							}
						}
					}					
				}				
			}
			/////////////////////                                  计算
			/////////此大批共有nu个,都保存到dtem[]中了 ////////////////
			if(nu>0)
			{
				//////////         取得第i小批      //////////////////
				long nt=0;         /////已经处理cell数量
				while(nt<nu)     ////////////////表示还没有处理完此大批的cell
				{
					long cel;
					long tt=nu-1;
					if((nu-nt)<20000)
					{
						tt=0;
						while(dtem[tt]==-1)    /////////用-1 表示已经处理过的cell
						{
							tt++;
						}
					}
					else
					{
						tt=nu-1;
						while(dtem[tt]==-1)    /////////用-1 表示已经处理过的cell
						{
							tt--;
						}
					}
					cel=dtem[tt]; ///////// 碰到一个没有处理的cell
					////获取围着它的小批,注意只保存此组里面的, 边界 cell 
					long stem[5000000]; //////////////??????????????????????????????????????????????/       
					long snu=0;
                    long afn=0;
					long bfn=-1;
					stem[0]=cel;
					snu++;
					long psta=0;
					while(afn>bfn)  ///////    preimage   与    image of cel
					{
						cout<<"small  i'th afn "<<afn<<endl;
						bfn=afn;
						///               long tbg=0;						
						for(long k=psta;k<snu;k++)  //////////
						{		
							if(snu>10000)
							{
								if(k%500==0)
									cout<<k<<"           "<<snu<<endl;
							}
							long tce;
							tce=stem[k];
							psta++;
							////////   考察preimage
                            gcm=fgcm[tce]; 
							long Jz,Iz;							
							Iz=gcm.Iz;							
							long lociz=gcm.CpLoc;
							Jz=gcm.Jz;
							long locjz=gcm.PreLoc;
							
							for(long t=0;t<Iz;t++)
							{
								long loc=lociz+t;
								long icel;
								cp=fcpl[loc];								
								icel=cp.Cz;
								////////////查看是否是其他吸引域里面的,因为边界容易出现这种情况
                                long tt=0;