求助“胞映射”的程序

梅梅 · 发表于 2008-3-18 20:04

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

:handshake 本人编了一个程序，但是程序内有多个循环体套在一起，运行速度慢极了，哪位好心人能给出更好的程序？急

cinly_zhu · 发表于 2008-3-18 21:16

能否将你的程序附上来？我也编了一个，Matlab编的，速度也一般。不知道你是做的哪一种胞映射，胞空间划分为多大？我们可以探讨。

梅梅 · 发表于 2008-3-19 08:51

我做的针对四维状态空间，胞空间划分10*10*10*10;但这种划分太不细了。可是这样速度还是慢阿。
我的QQ：806885738，咱们交流一下

梅梅 · 发表于 2008-3-19 17:15

我做的是简单胞映射

cinly_zhu · 发表于 2008-3-19 23:06

我做的广义胞映射，二维的，duffing方程

梅梅 · 发表于 2008-3-20 09:39

那你对简单胞映射应该很熟啊，给指点一下吧。主要是胞矢量坐标的选取原则（方法）。我是以区间段的个数定义胞矢量的坐标的，可以吗

风花雪月 · 发表于 2008-4-4 18:28

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

#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <math.h>
#include "mymatrix.h"
const long step1=75;
const long step2=90;////////////
const long step3=75;
const long step4=90;
const long max=step1*step2*step3*step4+1;/////
const float h1=5.05/step1;///////
const float h2=6.06/step2;///////
const float h3=5.05/step3;
const float h4=6.06/step4;
const float PI=3.14159265358979;
float x01=-2.525;
float x02=-3.03;
float x03=-2.525;
float x04=-3.03;
long gn; // 总组号
long ugn; ///总的不稳定组号
long g_attr;
struct dom
{
long total;
long n;
};
struct gcmdata
{
long Iz;
long Jz;
long gr;
long CpLoc;
long PreLoc;
};
struct cpdata
{
long Cz;
float Pz;
};
class lin
{
public:
lin *next;
long data;
lin();
};
lin::lin()
{
next=NULL;
}
class xlink
{
public:
float x;
int num;
xlink *next;
xlink();
};
xlink::xlink()
{
next=NULL;
}
void main()
{
gn=0;
ugn=0;
long map(long z,long j,long n);
void getbasedata();
void ab_per();
void determine();
void outputdata();
void fenzu();
//////////////////////////////////////////////////////////////////
///*
getbasedata();
ab_per();
/////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//gn=3;
//*
cout<<" start calculu"<<endl;
g_attr=gn;
ofstream of("per_groupn.txt");
of<<"persistent group numbers: "<<gn<<endl; /////////读取数据时要参考
////////////////////////////
if(gn<2)
{
cout<<"there 's no stable zone"<<endl;
exit(0);
}
////////////////////////////
//*/
outputdata();
//determine();
//fenzu();
//*/
}
long map(long m,long j)
{
//*
float x1,x2,x3,x4;
long z1,z2,z3,z4;
long j1,j2,j3,j4;
long number=2;
j4=(long)((j-1)/(number*number*number))+1;
j=j-(j4-1)*number*number*number;
j3=(long)((j-1)/(number*number))+1;
j=j-(j3-1)*(number*number);
j2=(long) ((j-1)/number)+1;
j1=j-(j2-1)*number;
z4=(long)((m-1)/(step1*step2*step3))+1;
m=m-(z4-1)*step1*step2*step3;
z3=(long)((m-1)/(step1*step2))+1;
m=m-(z3-1)*step1*step2;
z2=(long) ((m-1)/step1)+1;
z1=m-(z2-1)*step1;
x1=z1*h1-h1+x01+j1*h1/number-0.5*h1/number;
x2=z2*h2-h2+x02+j2*h2/number-0.5*h2/number;
x3=z3*h3-h3+x03+j3*h3/number-0.5*h3/number;
x4=z4*h4-h4+x04+j4*h4/number-0.5*h4/number;
//float x1,x2,x3,x4;
float h=0.08;
float mu=0.25;
float yt=0.04;
float v=0.1;
float t;
float k1,k2,k3,k4,l1,l2,l3,l4,s1,s2,s3,s4,d1,d2,d3,d4;
////////////////////////////
for(long i=0;i<30;i++)
{
k1=x2;
l1=mu*(1-x1*x1)*x2-(1+v)*x1+v*x3;
s1=x4;
d1=mu*(1-x3*x3)*x4+v*x1-(1+yt+v)*x3;
k2=x2+h*l1/2;
l2=mu*(1-(x1+h*k1/2)*(x1+h*k1/2))*(x2+h*l1/2)-(1+v)*(x1+h*k1/2)+v*(x3+h*s1/2);
s2=x4+d1*h/2;
d2=mu*(1-(x3+h*s1/2)*(x3+h*s1/2))*(x4+d1*h/2)+v*(x1+h*k1/2)-(1+yt+v)*(x3+h*s1/2);
k3=x2+h*l2/2;
l3=mu*(1-(x1+h*k2/2)*(x1+h*k2/2))*(x2+h*l2/2)-(1+v)*(x1+h*k2/2)+v*(x3+h*s2/2);
s3=x4+d2*h/2;
d3=mu*(1-(x3+h*s2/2)*(x3+h*s2/2))*(x4+d2*h/2)+v*(x1+h*k2/2)-(1+yt+v)*(x3+h*s2/2);
k4=x2+h*l3/2;
l4=mu*(1-(x1+h*k3)*(x1+h*k3))*(x2+h*l3)-(1+v)*(x1+h*k3)+v*(x3+h*s3);
s4=x4+d3*h/2;
d4=mu*(1-(x3+h*s3)*(x3+h*s3))*(x4+d3*h)+v*(x1+h*k3)-(1+yt+v)*(x3+h*s3);
x1=x1+h*(k1+2*k2+2*k3+k4)/6.0;
x2=x2+h*(l1+2*l2+2*l3+l4)/6.0;
x3=x3+h*(s1+2*s2+2*s3+s4)/6.0;
x4=x4+h*(d1+2*d2+2*d3+d4)/6.0;
}
z1=(long)((x1-x01)/h1)+1;
z2=(long)((x2-x02)/h2)+1;
z3=(long)((x3-x03)/h3)+1;
z4=(long)((x4-x04)/h4)+1;
if(z1>=step1 || z1<=0 || z2 >=step2 || z2<=0 || z3>=step3 || z3<=0 || m==0 || z4>=step4 || z4<=0)
return 0;
else
{
return (z4-1)*step1*step2*step3+(z3-1)*step1*step2+(z2-1)*step1+z1;
}
//*/
}
void getbasedata()
{
fstream iogcm,iocp,iopre;
iocp.open("cp.dat",ios::out|ios::binary|ios::in);
iogcm.open("gcm.dat",ios::in|ios::out|ios::binary);
//ifstream igcm("gcm.dat",ios::in);
iopre.open("pre.dat",ios::out|ios::binary|ios::in);
gcmdata gcm;
cpdata cp;
///*
long tz=0;
long snm=16;
long tempz[16][2]; ////////////记录
//////////////////////
for(long i=0;i<max;i++)
{
if(i%20==0)
cout<<"get Iz CpLoc cp.dat"<<i<<endl;
for(long k=0;k<snm;k++) /////////////每次都要置0 -1
{
for(long e=0;e<2;e++)
{
if(e==0)
{
tempz[k][e]=-1;
}
else
{
tempz[k][e]=0;
}
}
}
long r=0;
for(long j=0;j<snm;j++)
{
tz=map(i,j);
r=0;
while(tempz[r][0]!=tz && tempz[r][0]!=-1) ////////搜索到tz 或 tz为新值
{
r++;
}
tempz[r][0]=tz;
tempz[r][1]=tempz[r][1]+1;
}
r=0;
while(r<snm && tempz[r][0]!=-1)
{
r++;
}
//Iz[i]=r;
gcm.Iz=r;
long preIz;
if(i==0)
{
//CpLoc[0]=0;
gcm.CpLoc=0;
}
else
{
gcm.CpLoc=gcm.CpLoc+preIz;
}
///
preIz=r;
gcm.Jz=0;
gcm.PreLoc=0;
gcm.gr=0;
//outgcm.seekp(i*sizeof(gcmdata));
iogcm.write((char *)(&gcm),sizeof(gcmdata));
for(long t=0;t<r;t++)
{
cp.Cz=tempz[t][0];
cp.Pz=(float)(tempz[t][1])/snm;
///////////////////
iocp.write((char *)(&cp),sizeof(cpdata));
}
}
cout<<"The following part will get preimage data,please wait,first Jz"<<endl;
////////////////////// preimage ////////////////////////////////////
///iogcm.seekg(0+sizeof(long));
for(long b=0;b<max;b++)
{
if(b%2000==0)
cout<<"get Jz "<<b<<endl;
iogcm.seekg(b*sizeof(gcmdata));
iogcm.read((char *)(&gcm),sizeof(gcmdata));
long I=gcm.Iz;
for(long n=0;n<I;n++)
{
long jz;
long loc=gcm.CpLoc+n;
iocp.seekg(loc*sizeof(cpdata));
long cel;
iocp.read((char *)(&cel),sizeof(long));
iogcm.seekg(cel*sizeof(gcmdata)+sizeof(long),ios::beg);
iogcm.read((char *)(&jz),sizeof(long));///////先读取jz
jz++;
iogcm.seekp(cel*sizeof(gcmdata)+sizeof(long),ios::beg);
iogcm.write((char *)(&jz),sizeof(long));
}
}
//////// 计算 preloc /////
cout<<endl<<"get preLoc "<<endl;
long preloc;
for(b=0;b<max;b++)
{
if(b%2000==0)
cout<<"get preLoc "<<b<<endl;
if(b==0)
{
preloc=0;
iogcm.seekp(0+4*sizeof(long));//// preLoc的位置
iogcm.write((char *)(&b),sizeof(long));
}
else
{
iogcm.seekg((b-1)*sizeof(gcmdata)+sizeof(long));//////// prejz
long prejz;
iogcm.read((char *)(&prejz),sizeof(long));
preloc=preloc+prejz;
iogcm.seekp(b*sizeof(gcmdata)+4*sizeof(long));
iogcm.write((char *)(&preloc),sizeof(long));
}
}
///////////////////////////////////// 保存 preimage cell ///
cout<<" preimage cell "<<endl;
iogcm.seekg((max-1)*sizeof(gcmdata)+sizeof(long));
long jzmax;
iogcm.read((char *)(&jzmax),sizeof(long));
long totalpre=jzmax+preloc+1;///////////////////////////////////
/////////////////////////////////////////////////////////////////算法很糟糕/////////////////
cout<<"write preimage"<<endl;
long *pret;
pret=new long [totalpre];
for(b=0;b<totalpre;b++)
{
if(b%10000==0)
cout<<"write pre.dat "<<b<<endl;
pret[b]=-1;
}
cout<<" get preimage "<<endl;
//////////////////////////////////////////read prel data
long *ptt;
ptt=new long [max];
iogcm.seekg(0);
for(long tt=0;tt<max;tt++)
{
iogcm.read((char *)(&gcm),sizeof(gcmdata));
ptt[tt]=gcm.PreLoc;
}
iocp.seekg(0);
iogcm.seekg(0);
for(b=0;b<max;b++)
{
if(b%500==0)
cout<<"get preimage "<<b<<endl;
//cout<<b<<" "<<b<<endl;
iogcm.read((char *)(&gcm),sizeof(gcmdata));
long I=gcm.Iz;
for(long n=0;n<I;n++)
{
long cel;
iocp.read((char *)(&cp),sizeof(cpdata));
cel=cp.Cz;
long preL;
preL=ptt[cel];
long t=0;
long pce;
pce=pret[preL+t];
while(pce!=-1)
{
t=t+1;
pce=pret[preL+t];
}
pret[preL+t]=b;
}
}
////////////write data
for(tt=0;tt<totalpre;tt++)
{
if(tt%2000==0)
cout<<"save predata "<<tt<<endl;
iopre.write((char *)(&pret[tt]),sizeof(long));
}
delete [] ptt;
delete [] pret;
iogcm.close();
iocp.close();
iopre.close();
}
void ab_per()
{
////select simple map group
fstream iosgr,iogcm1,iocp1,ioscm,iopersistent;
iosgr.open("sgr.dat",ios::in|ios::out|ios::binary);
iogcm1.open("gcm.dat",ios::in|ios::out|ios::binary);
iocp1.open("cp.dat",ios::in|ios::out|ios::binary);
ioscm.open("tempscm.dat",ios::in|ios::out|ios::binary);////保存scm的结果
iopersistent.open("persistent.dat",ios::in|ios::out|ios::binary);
gcmdata gcm;
cpdata cp;
///*
cout<<"simple cell method"<<endl;
long sgr;
long sgn=0;
{
long temz[10000000]; //==========================????????????????高维时注意
for(long i=0;i<max;i++)
{
sgr=0;
iosgr.write((char *)(&sgr),sizeof(long));
//sgr[i]=0;
}
long m,b;
long g=0;
for(long j=0;j<max;j++)
{
if(j%1000==0)
cout<<j<<" simple cell method, get sgn "<<endl;
m=0;
iosgr.seekg(j*sizeof(long));
iosgr.read((char *)(&sgr),sizeof(long));
if(sgr!=0)
{
continue; ////////////////不是virgin cell,重新选取
}
else
{
sgr=-1;
iosgr.seekp(j*sizeof(long));
iosgr.write((char *)(&sgr),sizeof(long));
b=j;
temz[m]=b;
do
{
sgr=-1;
iosgr.seekp(b*sizeof(long));
iosgr.write((char *)(&sgr),sizeof(long));
m=m+1;
////////////////
iogcm1.seekg(b*sizeof(gcmdata));
iogcm1.read((char *)(&gcm),sizeof(gcmdata));
long loc=gcm.CpLoc;
iocp1.seekg(loc*sizeof(cpdata));
iocp1.read((char *)(&cp),sizeof(cpdata));
b=cp.Cz;
//b=Cz[CpLoc[b]];//////select the first image cell
temz[m]=b;
iosgr.seekg(b*sizeof(long));
iosgr.read((char *)(&sgr),sizeof(long));
}while(sgr==0); ///
if(sgr==-1) ///////////////新周期解
{
///////////寻找i/////////（刚好在周期轨道上）
g++;
sgn=g;
long u=0;
while(temz[u]!=b)
{
u++;
}
for(long q=0;q<m;q++)
{
if(q<u)
{
iosgr.seekg(temz[q]*sizeof(long));
long c=-2;
iosgr.write((char *)(&c),sizeof(long));
//sgr[temz[q]]=-2;
}
else
{
iosgr.seekp(temz[q]*sizeof(long));
iosgr.write((char *)(&g),sizeof(long));
long tt=temz[q];
ioscm.write((char *)(&g),sizeof(long));/////////先写group
ioscm.write((char *)(&tt),sizeof(long));
}
}
}
else //////吸引到已知一点上
{
for(long k=0;k<m;k++)
{
iosgr.seekg(temz[k]*sizeof(long));
long c=-2;
iosgr.write((char *)(&c),sizeof(long));
}
}
}
}
}
cout<<"search for persistent groups"<<endl;
////////////////////////////////////////////////////////////////////////////////////////
{
long temp[10000000]; ////======================??????????????save candidate cells
long n=0;
long flag=1;
//ioscm.clear();
for(long i=1;i<=sgn;i++)
{
cout<<sgn<<" "<<i<<endl;
///////////////////////////////先判断组号是否已经是absorbing group了
n=0;
ioscm.seekg(0);
long cel,gg;
ioscm.read((char *)(&gg),sizeof(long));
ioscm.read((char *)(&cel),sizeof(long));
long cur,pend;
cur=ioscm.tellg();
ioscm.seekg(0,ios::end);
pend=ioscm.tellg();
ioscm.seekg(cur);
while(pend!=ioscm.tellg())
{
if(gg==i)
{
temp[n]=cel;
n++;
}
ioscm.read((char *)(&gg),sizeof(long));
ioscm.read((char *)(&cel),sizeof(long));
}
ioscm.seekg(0);
cur=ioscm.tellg();
ioscm.seekg(0,ios::end);
pend=ioscm.tellg();
ioscm.seekg(cur);
n=0;
while(pend!=ioscm.tellg())
{
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的判断 ////////
}
}
//*/
///////////////////////////////////////////////////////
}

复制代码

风花雪月 · 发表于 2008-4-4 18:28

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;
while(stem[tt]!=icel&&tt<snu)
{
tt++;
}
if(tt==snu)/////查看是否是其他吸引域里面
{
long iflgf;
iflgf=flag[icel];
//////////已经确定的cell也不能加进去
long ddi;
ddi=fder[icel];
if(iflgf==i&&ddi!=i) ////
{
stem[snu]=icel;
snu++;
afn++;
}
}
}
}
}
///////////////// 此小批已经获得,放在stem[] , 有snu个 /////////////////////////////
///////// snu 个stem排序
long *sstem;
sstem=new long [snu];
for(long uu=0;uu<snu;uu++)
{
sstem[uu]=stem[uu];
}
/////已经胞排好了,放在sstem[snu]中
////开始计算
long threm[5000000];
long thrn=0;
long nhand=0;
while(nhand<snu)
{
thrn=0;
long cre;
long tt=snu-1;
if((snu-nhand)<200)
{
tt=0;
while(sstem[tt]==-1)///////////////从最后一个开始算
{
tt++;
}
}
else
{
tt=snu-1;
while(sstem[tt]==-1)///////////////从最后一个开始算
{
tt--;
}
}
cre=sstem[tt];
threm[thrn]=cre;
thrn++;
long an=0;
long bn=-1;
long psta=0;
while(an>bn)
{
bn=an;
for(long u=psta;u<thrn;u++)
{
long crm=threm[u];
psta++;
long Iz,loc,loc1;
gcm=fgcm[crm];
Iz=gcm.Iz;
loc=gcm.CpLoc;
for(long t=0;t<Iz;t++)
{
loc1=loc+t;
cp=fcpl[loc1];
long cie=cp.Cz;
long det;
det=fder[cie];
long flt;
flt=flag[cie];
if(det!=i&&flt==i)
{
//////保存了吗
long ttt=0;
while(threm[ttt]!=cie&&ttt<thrn)
{
ttt++;
}
if(ttt==thrn)
{
an++;
threm[thrn]=cie;
thrn++;
}
}
}
}
}
//////////////不稳定组, //////////
//////计算,已经放到threm[thrn]中
////////pj////////////////
handnum=handnum+thrn;
cout<<"total : "<<numb<<" handled : "<<handnum<<" equation : "<<thrn<<endl;
float *xr,*xr2;
xr=new float [thrn];/////////迭代用
xr2=new float [thrn];
float err=1.0;
for(k=0;k<thrn;k++)
{
xr[k]=0.0;
}
xr[0]=0.65; ////////////////初值
if(thrn<1000)
{
mat.Realloc(thrn,thrn);
mat=0.0;
}
xlink **head;
xlink *p;
head=new xlink *[thrn];
if(thrn>=1000)
{
for(int j=0;j<thrn;j++)
{
head[j]=new xlink [1];
if(j>=1)
head[j]->num=0;
}
}
cst.Realloc(thrn,1);
cst=0.0;
float *cs;
cs=new float [thrn];
if(thrn<1000)
{
for(k=0;k<thrn;k++)
{
if(thrn>1000&&thrn%500==0)
cout<<"get xishu "<<k<<" "<<thrn<<endl;
cs[k]=0.0;
long ce=threm[k];
long Iz;
gcm=fgcm[ce];
Iz=gcm.Iz;
long loc2=gcm.CpLoc;
mat.m_adValue[k][k]=1.0;///////// 别忘记了
for(long t=0;t<Iz;t++)
{
long loc=loc2+t;
cp=fcpl[loc];
long ceim=cp.Cz;////////
long tt=0;
while(ceim!=threm[tt]&&tt<thrn) /////// ????????tt<tnb
{
tt++;
}
if(tt==thrn) /////////象为已知胞, 也可能为未知胞
{ ////看到
long ttfg;
ttfg=flag[ceim];
long dfdf;
dfdf=fder[ceim];
if(ttfg==i&&dfdf==i)
{
float pj;
pj=fres[ceim][i-1];
cs[k]=cs[k]+cp.Pz*pj; ////常数项
}
}
else
{
if(tt==k)////////与ce相同
{
////ce 到 ce 的概率
mat.m_adValue[k][k]=mat.m_adValue[k][k]-cp.Pz;
}
else
{
mat.m_adValue[k][tt]=-cp.Pz;
}
}
}
}
}
else
{
for(k=0;k<thrn;k++)
{
if(thrn>1000&&thrn%500==0)
cout<<"get xishu "<<k<<" "<<thrn<<endl;
cs[k]=0.0;
long ce=threm[k];
long Iz;
gcm=fgcm[ce];
Iz=gcm.Iz;
long loc2=gcm.CpLoc;
p=head[k];
p->num=k;
p->x=1.0;
//mat.m_adValue[k][k]=1.0;///////// 别忘记了
for(long t=0;t<Iz;t++)
{
long loc=loc2+t;
cp=fcpl[loc];
long ceim=cp.Cz;////////
long tt=0;
while(ceim!=threm[tt]&&tt<thrn) /////// ????????tt<tnb
{
tt++;
}
if(tt==thrn) /////////象为已知胞, 也可能为未知胞
{ ////看到
long ttfg;
ttfg=flag[ceim];
long dfdf;
dfdf=fder[ceim];
if(ttfg==i&&dfdf==i)
{
float pj;
pj=fres[ceim][i-1];
cs[k]=cs[k]+cp.Pz*pj; ////常数项
}
}
else
{
if(tt==k)////////与ce相同
{
////ce 到 ce 的概率
//mat.m_adValue[k][k]=mat.m_adValue[k][k]-cp.Pz;
p=head[k];
p->x=p->x-cp.Pz;
}
else
{
//mat.m_adValue[k][tt]=-cp.Pz;
p=head[k];
while(p->next!=NULL)
{
p=p->next;
}
p->next=new xlink();
p=p->next;
p->x=-cp.Pz;
p->num=tt;
}
}
}
}
}
/////////上面获得系数
for(long yu=0;yu<thrn;yu++)
{
cst.m_adValue[yu][0]=cs[yu];
}
if(thrn<1000)
{
mat.MatInv();
cst=mat*cst;
///////////////////////
}
else
{
err=1.0;
while(err>0.01)
{
float tot=0.0;
for(long tr=0;tr<thrn;tr++)
{
if(tr%500==0)
cout<<tr<<" "<<thrn<<endl;
tot=0.0;
p=head[tr];
p=p->next;/////第二个胞了，非tr
while(p!=NULL)
{
tot=tot+(p->x)*xr[p->num];
p=p->next;
}
p=head[tr];
float njt=p->x;
xr2[tr]=(cs[tr]-tot)/njt;
}
err=0.0;
for(long jr=0;jr<thrn;jr++)
{
err=err+sqrt((xr2[jr]-xr[jr])*(xr2[jr]-xr[jr]));
}
for(jr=0;jr<thrn;jr++)
{
xr[jr]=xr2[jr];
}
cout<<err<<" "<<err<<endl;
}
////////////结果放在cst
for(long te=0;te<thrn;te++)
{
cst.m_adValue[te][0]=xr[te];
}
for(long mt=0;mt<thrn;mt++)
{
xlink *q1,*q2;
q1=head[mt]->next;
while(q1!=NULL)
{
q2=q1->next;
delete q1;
q1=q2;
}
}
}
////存入数据
for(long tr=0;tr<thrn;tr++)
{
long tre=threm[tr];
float pjr=cst.m_adValue[tr][0];
fres[tre][i-1]=pjr;
}
/////////////vj////////////////
delete [] cs;
delete [] xr;
delete [] xr2;
////////结果,改变sstem -1 以及ioder
for(long u=0;u<thrn;u++)
{
long cme=threm[u];
long det=i;
fder[cme]=det;
long tr=0;
while(sstem[tr]!=cme&&tr<snu)
{
tr++;
}
sstem[tr]=-1;
}
nhand=nhand+thrn;
//for(long loop=0;loop<thrn;loop++)//////////////////清理
delete [] head;
//delete [] head;
}
//delete [] stt;
//delete [] stn;
delete [] sstem;
///////////////////////////////////////
for(long u=0;u<snu;u++)
{
long cew;
cew=stem[u];
long yu=i;
fder[cew]=yu;
/////////////////dtem 中 -1
long tt=0;
while(dtem[tt]!=cew&&tt<nu)
{
tt++;
}
dtem[tt]=-1;/////////////////////
}
nt=nt+snu;
}
}
cet=ftemp[sta];
sta++;
///////////////此批计算完//////////////////////
}
}
////////////////////////
fstream ioresult; ////// 保存结果 gr,Vj,pj
ioresult.open("result.dat",ios::in|ios::out|ios::binary);
///long res=3*sizeof(float);/////////////////////////////////////////
for(long mt=0;mt<max;mt++)
{
for(long t=0;t<3;t++)
{
ioresult.write((char *)(&fres[mt][t]),sizeof(float));
}
}
delete [] fzpre;
delete [] flag;
delete [] fder;
delete [] fcpl;
delete [] fgcm;
delete [] ftemp;
}
////////////排序算法
void sort(long a[],long b[],long n)
{
long i,j,small;
long temp1,temp2;
for(i=0;i<n-1;i++)
{
small=i;
for(j=i+1;j<n;j++)
{
if(a[j]<a[small])
small=j;
}
if(small!=i)
{
temp1=a[i];
temp2=b[i];
a[i]=a[small];
a[small]=temp1;
b[i]=b[small];
b[small]=temp2;
}
}
}
long pai(long m)
{
long n=g_attr;//
long nt=1;
long mt=1;
if(m>n)
{
return -1;
}
else
{
for(long i=1;i<=m;i++)
{
nt=nt*n;
n=n-1;
}
if(m==0)
{
mt=1;
}
else
{
for(i=1;i<=m;i++)
{
mt=mt*i;
}
}
long re;
re=(long) (nt/mt);
return re;
}
}
void fenzu()
{
///// 根据 iodm 中的信息,改变ioresult 中的gr值
fstream iodm; ////// 保存 cell 的, domicle 的数量,做分组用
iodm.open("todomiledata.dat",ios::in|ios::out|ios::binary);
dom dm;
fstream ioresult; ////// 保存结果 gr,Vj,pj
ioresult.open("result.dat",ios::in|ios::out|ios::binary);
long res=sizeof(long)+2*g_attr*sizeof(float);//////////////////////////
long pai(long m);
for(long i=0;i<max;i++)
{
if(i%5000==0)
cout<<i<<endl;
long dn;
long total=0;
iodm.seekg(i*sizeof(dom));
iodm.read((char *)(&dm),sizeof(dom));
dn=dm.n;
total=dm.total;
long gpai=0;
if(dn==0)
{
ioresult.seekp(i*res);
ioresult.write((char *)(&total),sizeof(long));
}
}
///////////////////////////////////
}
void outputdata()
{
fstream iopersistent;
iopersistent.open("persistent.dat",ios::in|ios::out|ios::binary);
//g_attr=2;
long tr=0;
ofstream df("persistent.txt");
iopersistent.seekg(0,ios::beg);
long cur,pend;
cur=iopersistent.tellg();
iopersistent.seekg(0,ios::end);
pend=iopersistent.tellg();
iopersistent.seekg(cur);
while(pend!=iopersistent.tellg())
{
long cee,grg;
iopersistent.read((char *)(&cee),sizeof(long));
iopersistent.read((char *)(&grg),sizeof(long));
float x1,x2,x3,x4;
long z1,z2,z3,z4,m;
m=cee;
z4=(long)((m-1)/(step1*step2*step3))+1;
m=m-(z4-1)*step1*step2*step3;
z3=(long)((m-1)/(step1*step2))+1;
m=m-(z3-1)*step1*step2;
z2=(long) ((m-1)/step1)+1;
z1=m-(z2-1)*step1;
x1=z1*h1-h1+x01;
x2=z2*h2-h2+x02;
x3=z3*h3-h3+x03;
x4=z4*h4-h4+x04;
df<<x1<<" "<<x2<<" "<<x3<<" "<<x4<<endl;
}
}

复制代码

水风而动 · 发表于 2013-11-21 17:05

梅梅发表于 2008-3-19 08:51

登录/注册后可看大图

我做的针对四维状态空间，胞空间划分10*10*10*10;但这种划分太不细了。可是这样速度还是慢阿。
我的QQ：80 ...

师姐为什么加不了你呢？

YjF · 发表于 2018-4-25 11:23

cinly_zhu 发表于 2008-3-19 23:06
**** 作者被禁止或删除内容自动屏蔽 ****

你好，我是刚开始研究胞映射相关课题的学生，求学长指点。能加下qq嘛？

账号		自动登录	找回密码
密码			我要加入

[经典算法] 求助“胞映射”的程序

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

回复楼主的帖子

回复 2楼的帖子

回复 4楼的帖子

[经典算法] 求助“胞映射”的 程序

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

回复 楼主 的帖子

回复 2楼 的帖子

回复 4楼 的帖子

[经典算法] 求助“胞映射”的程序

回复楼主的帖子

回复 2楼的帖子

回复 4楼的帖子