算法及编程语言 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断

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

求助“胞映射”的程序

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

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

回复楼主的帖子

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

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

回复 2楼的帖子

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

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

我做的是简单胞映射

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

回复 4楼的帖子

我做的广义胞映射，二维的，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;////////////记录
//////////////////////
for(long i=0;i<max;i++)
{
if(i%20==0)
cout<<"get IzCpLoc cp.dat"<<i<<endl;
for(long k=0;k<snm;k++) /////////////每次都要置0-1
{
for(long e=0;e<2;e++)
{
if(e==0)
{
tempz=-1;
}
else
{
tempz=0;
}
}
}

long r=0;
for(long j=0;j<snm;j++)
{
tz=map(i,j);
r=0;
while(tempz!=tz && tempz!=-1)////////搜索到tz 或 tz为新值
{
r++;
}
tempz=tz;
tempz=tempz+1;
}

   r=0;
while(r<snm && tempz!=-1)
{
r++;
}

//Iz=r;
gcm.Iz=r;
long preIz;
if(i==0)
{
//CpLoc=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;
cp.Pz=(float)(tempz)/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 ;
for(b=0;b<totalpre;b++)
{
if(b%10000==0)
cout<<"write pre.dat "<<b<<endl;
pret=-1;

}
cout<<" get preimage"<<endl;
//////////////////////////////////////////read prel data
long *ptt;
ptt=new long ;
iogcm.seekg(0);
for(long tt=0;tt<max;tt++)
{
iogcm.read((char *)(&gcm),sizeof(gcmdata));
ptt=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;
long t=0;
long pce;
pce=pret;
while(pce!=-1)
{
t=t+1;
pce=pret;
}
pret=b;
}
}
////////////write data
for(tt=0;tt<totalpre;tt++)
{
if(tt%2000==0)
cout<<"save predata "<<tt<<endl;
iopre.write((char *)(&pret),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;       //==========================????????????????高维时注意
for(long i=0;i<max;i++)
{
sgr=0;
iosgr.write((char *)(&sgr),sizeof(long));
//sgr=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=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];//////select the first image cell
temz=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!=b)
{
u++;
}
for(long q=0;q<m;q++)
{
if(q<u)
{
iosgr.seekg(temz*sizeof(long));
long c=-2;
iosgr.write((char *)(&c),sizeof(long));
//sgr]=-2;
}
else
{
iosgr.seekp(temz*sizeof(long));
                        iosgr.write((char *)(&g),sizeof(long));
long tt=temz;
ioscm.write((char *)(&g),sizeof(long));/////////先写group
ioscm.write((char *)(&tt),sizeof(long));
}
}

}
else       //////吸引到已知一点上
{
for(long k=0;k<m;k++)
{
iosgr.seekg(temz*sizeof(long));
long c=-2;
                     iosgr.write((char *)(&c),sizeof(long));
}
}
}
}
}
cout<<"search for persistent groups"<<endl;
////////////////////////////////////////////////////////////////////////////////////////
{
long temp;       ////======================??????????????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=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=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
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;
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=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;
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*sizeof(gcmdata)+2*sizeof(long));
iogcm1.write((char *)(&gn),sizeof(long));
long celp=temp;
iopersistent.write((char *)(&celp),sizeof(long));
iopersistent.write((char *)(&gn),sizeof(long));
}
}
else
{
for(long ui=0;ui<n;ui++)
{
iogcm1.seekg(temp*sizeof(gcmdata));
iogcm1.read((char *)(&gcm),sizeof(gcmdata));
gcm.gr=0;
iogcm1.seekp(temp*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 ;
for(long i=0;i<max;i++)
{
   if(i%2000==0)
         cout<<i<<endl;
   flag=-1;
}
long tem;/////////保存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 ;
iopre.seekg(0);
for(long ttr=0;ttr<totalpre;ttr++)
{
   iopre.read((char *)(&fzpre),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 ;
iocp.seekg(0);
for(ttr=0;ttr<totalcl;ttr++)
{
   iocp.read((char *)(&fcpl),sizeof(cpdata));
}
////////////////////////////////////////////
gcmdata *fgcm;
fgcm=new gcmdata ;
iogcm.seekg(0);
for(ttr=0;ttr<max;ttr++)
{
   iogcm.read((char *)(&fgcm),sizeof(gcmdata));
}

//////////////////////// 输出结果的文件 //////////////////////////////////////////////////

//////////////初始化result
//long res=(gn+1)*sizeof(long)+gn*sizeof(float);
float fres;
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=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 ;
///fstream iocw;//////////成尾数计算是用
//iocw.open("chengwei.dat",ios::in|ios::out|ios::binary);
for(i=0;i<max;i++)
{ /////////////////所有的cell 数量, total numberpersistent group 也当成0
   if(i%2000==0)
         cout<<i<<endl;
   iodm.write((char *)(&dm),sizeof(dom)); /////////先保存
   long ffff=0;
   fder=ffff;
}
/////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////
long *ftemp;

ftemp=new long ;
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=ff;
            tem=cee;
            n++;
            ////////////////////////////////////persistent cell为已知胞////////////////////////
            fder=ff;
                           fres=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; ///////// persistent cells
         long Jz;
         gcm=fgcm;
         Jz=gcm.Jz;
         long loc2=gcm.PreLoc;
         for(long p=0;p<Jz;p++)
         {
            long loc=loc2+p;
            long ffg,iima;
            iima=fzpre;
            long dis;
            if(iima<0)///////////////////////////不取进来
            {
               ffg=i;
            }
            else
            {
            dis=fder;
            ffg=flag;
            }
            ////////////////////////////////////////
            //long ddd;

            ///////////////////////////////////////
            if(ffg!=i&&dis!=i)
            {
               //iotemp.write((char *)(&iima),sizeof(long));
               ftemp=iima;
               ftn++;
               flag=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=-1;
   ftn++;
   ///////////////////////////////迭代
   long afn,bfn;
   afn=0;
   bfn=-1;
   long sta=0;
   while(afn>bfn)
   {
         cout<<"preimage of persistentan"<<afn<<" bn "<<bfn<<endl;
         bfn=afn;
         //iotemp.seekg(dbe*sizeof(long));
         long cee;
         cee=ftemp;
         //iotemp.read((char *)(&cee),sizeof(long));
         sta++;
         while(cee!=-1)
         {
            long Jz;
            gcm=fgcm;
            Jz=gcm.Jz;
            long loc2=gcm.PreLoc;
            for(long p=0;p<Jz;p++)
            {
               long loc=loc2+p;
               long fg,ima;
               ima=fzpre;
               fg=flag;
               if(fg!=i)
               {
                     numb++;
                     afn++;
                     ftemp=ima;
                     ftn++;
                     flag=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++;
         }
         ///////// 次批完
         ftemp=-1;
         ftn++;
   }

   ////////上面获得了此persistent group的所有pre象,并已经按批次存放,用-1分隔开来iotemp中
   cout<<numb<<endl;
   long dtem; ///////////    保存临时像/ 第i批数据/////////////?????????????????
   long nu=0;
   long cet;
   sta=0;
   cet=ftemp;
   sta++;
   while(sta<ftn)
   {
         cout<<" get i'th cell of i'th preimage"<<endl;
         ///////    获取第i批数据及像    ////////
         nu=0;
         long dtm3;
         if(cet!=-1)    ////////文件指针不能移出去
         {
            dtm3=fder;
         }
         if(cet!=-1&&dtm3!=i)
         {
            dtem=cet;
            nu++;
         }
         cet=ftemp;
         sta++;
         while(cet!=-1)//////////          取得此大批preimage
         {
            long dtm;
            dtm=fder;
            if(dtm!=i)    /////////好像肯定不等于的??????????????
            {
               dtem=cet;          /////////已经处理也要放进取,避免它的下几层胞丢失
               nu++;
            }
            cet=ftemp;
            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;
                     psta++;
                     long dt; ////////////// 考察是否已经确定
                     gcm=fgcm;
                     long Iz;
                     Iz=gcm.Iz;
                     long loc2=gcm.CpLoc;
                     for(long y=0;y<Iz;y++)
                     {
                        long loc=loc2+y;
                        cp=fcpl;
                        long cim;
                        cim=cp.Cz;
                        //////// whether determined
                        dt=fder;
                        long tfgg;
                        tfgg=flag;
                        if(dt!=i&&tfgg==i) //////防止边界上,其它域里面的胞也保存进去
                        {
                           //////看是否已经保存了
                           long t=0;
                           while(dtem!=cim&&t<nu)
                           {
                                 t++;
                           }
                           if(t==nu) /////////没有保存
                           {
                                 dtem=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==-1) /////////用-1 表示已经处理过的cell
                     {
                        tt++;
                     }
               }
               else
               {
                     tt=nu-1;
                     while(dtem==-1) /////////用-1 表示已经处理过的cell
                     {
                        tt--;
                     }
               }
               cel=dtem; ///////// 碰到一个没有处理的cell
               ////获取围着它的小批,注意只保存此组里面的, 边界 cell
               long stem; //////////////??????????????????????????????????????????????/
               long snu=0;
               long afn=0;
               long bfn=-1;
               stem=cel;
               snu++;
               long psta=0;
               while(afn>bfn)/////// preimage 与 image of cel
               {
                     cout<<"smalli'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;
                        psta++;
                        //////// 考察preimage
                        gcm=fgcm;
                        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;
                           icel=cp.Cz;
                           ////////////查看是否是其他吸引域里面的,因为边界容易出现这种情况
                           long tt=0;
                           while(stem!=icel&&tt<snu)
                           {
                                 tt++;
                           }
                           if(tt==snu)/////查看是否是其他吸引域里面
                           {
                                 long iflgf;
                                 iflgf=flag;
                                 //////////已经确定的cell也不能加进去
                                 long ddi;
                                 ddi=fder;


                                 if(iflgf==i&&ddi!=i) ////
                                 {
                                    stem=icel;
                                    snu++;
                                    afn++;
                                 }
                           }
                        }
                     }
               }
               ///////////////// 此小批已经获得,放在stem[],有snu个 /////////////////////////////
               ///////// snu 个stem排序
               long *sstem;
               sstem=new long ;

               for(long uu=0;uu<snu;uu++)
               {
                     sstem=stem;
               }

               /////已经胞排好了,放在sstem中
               ////开始计算
               long threm;
               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==-1)///////////////从最后一个开始算
                        {
                           tt++;
                        }
                     }
                     else
                     {
                        tt=snu-1;
                        while(sstem==-1)///////////////从最后一个开始算
                        {
                           tt--;
                        }
                     }
                     cre=sstem;
                     threm=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;
                           psta++;
                           long Iz,loc,loc1;
                           gcm=fgcm;
                           Iz=gcm.Iz;
                           loc=gcm.CpLoc;
                           for(long t=0;t<Iz;t++)
                           {
                                 loc1=loc+t;
                                 cp=fcpl;
                                 long cie=cp.Cz;
                                 long det;
                                 det=fder;
                                 long flt;
                                 flt=flag;
                                 if(det!=i&&flt==i)
                                 {
                                    //////保存了吗
                                    long ttt=0;
                                    while(threm!=cie&&ttt<thrn)
                                    {
                                       ttt++;
                                    }
                                    if(ttt==thrn)
                                    {
                                       an++;
                                       threm=cie;
                                       thrn++;
                                    }
                                 }
                           }
                        }
                     }
                     //////////////不稳定组, //////////

                     //////计算,已经放到threm中
                     ////////pj////////////////

                     handnum=handnum+thrn;
                     cout<<"total: "<<numb<<"       handled : "<<handnum<<"    equation : "<<thrn<<endl;
                     float *xr,*xr2;
                     xr=new float ;/////////迭代用
                     xr2=new float ;
                     float err=1.0;
                     for(k=0;k<thrn;k++)
                     {
                        xr=0.0;
                     }
                     xr=0.65; ////////////////初值
                     if(thrn<1000)
                     {
                        mat.Realloc(thrn,thrn);
                        mat=0.0;
                     }
                     xlink **head;
                     xlink *p;
                     head=new xlink *;
                     if(thrn>=1000)
                     {
                        for(int j=0;j<thrn;j++)
                        {
                           head=new xlink ;
                           if(j>=1)
                                 head->num=0;

                        }
                     }

                     cst.Realloc(thrn,1);

                     cst=0.0;
                     float *cs;
                     cs=new float ;
                     if(thrn<1000)
                     {
                        for(k=0;k<thrn;k++)
                        {
                           if(thrn>1000&&thrn%500==0)
                                 cout<<"get xishu"<<k<<"    "<<thrn<<endl;
                           cs=0.0;
                           long ce=threm;
                           long Iz;
                           gcm=fgcm;
                           Iz=gcm.Iz;
                           long loc2=gcm.CpLoc;
                           mat.m_adValue=1.0;///////// 别忘记了
                           for(long t=0;t<Iz;t++)
                           {
                                 long loc=loc2+t;
                                 cp=fcpl;

                                 long ceim=cp.Cz;////////
                                 long tt=0;
                                 while(ceim!=threm&&tt<thrn) ///////????????tt<tnb
                                 {
                                    tt++;
                                 }
                                 if(tt==thrn) /////////象为已知胞, 也可能为未知胞
                                 { ////看到
                                    long ttfg;
                                    ttfg=flag;
                                    long dfdf;
                                    dfdf=fder;
                                    if(ttfg==i&&dfdf==i)
                                    {
                                       float pj;
                                                                                    pj=fres;
                                       cs=cs+cp.Pz*pj;////常数项
                                    }
                                 }
                                 else
                                 {
                                    if(tt==k)////////与ce相同
                                    {
                                       ////ce 到 ce的概率
                                       mat.m_adValue=mat.m_adValue-cp.Pz;
                                    }
                                    else
                                    {
                                       mat.m_adValue=-cp.Pz;
                                    }
                                 }
                           }
                        }
                     }
                     else
                     {
                        for(k=0;k<thrn;k++)
                        {
                           if(thrn>1000&&thrn%500==0)
                                 cout<<"get xishu"<<k<<"    "<<thrn<<endl;
                           cs=0.0;
                           long ce=threm;
                           long Iz;
                           gcm=fgcm;
                           Iz=gcm.Iz;
                           long loc2=gcm.CpLoc;
                           p=head;
                           p->num=k;
                           p->x=1.0;
                           //mat.m_adValue=1.0;///////// 别忘记了
                           for(long t=0;t<Iz;t++)
                           {
                                 long loc=loc2+t;
                                 cp=fcpl;
                                 long ceim=cp.Cz;////////
                                 long tt=0;
                                 while(ceim!=threm&&tt<thrn) ///////????????tt<tnb
                                 {
                                    tt++;
                                 }
                                 if(tt==thrn) /////////象为已知胞, 也可能为未知胞
                                 { ////看到
                                    long ttfg;
                                    ttfg=flag;
                                    long dfdf;
                                    dfdf=fder;
                                    if(ttfg==i&&dfdf==i)
                                    {
                                       float pj;
                                                                                    pj=fres;
                                       cs=cs+cp.Pz*pj;////常数项
                                    }
                                 }
                                 else
                                 {
                                    if(tt==k)////////与ce相同
                                    {
                                       ////ce 到 ce的概率
                                       //mat.m_adValue=mat.m_adValue-cp.Pz;
                                       p=head;
                                       p->x=p->x-cp.Pz;
                                    }
                                    else
                                    {
                                       //mat.m_adValue=-cp.Pz;
                                       p=head;
                                       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=cs;
                     }
                     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;
                                 p=p->next;/////第二个胞了，非tr
                                 while(p!=NULL)
                                 {
                                    tot=tot+(p->x)*xr;
                                    p=p->next;
                                 }
                                 p=head;
                                 float njt=p->x;
                                 xr2=(cs-tot)/njt;
                           }
                           err=0.0;
                           for(long jr=0;jr<thrn;jr++)
                           {
                                 err=err+sqrt((xr2-xr)*(xr2-xr));
                           }
                           for(jr=0;jr<thrn;jr++)
                           {
                                 xr=xr2;
                           }
                           cout<<err<<"          "<<err<<endl;
                        }
                        ////////////结果放在cst
                        for(long te=0;te<thrn;te++)
                        {
                           cst.m_adValue=xr;
                        }

                        for(long mt=0;mt<thrn;mt++)
                        {
                           xlink *q1,*q2;
                           q1=head->next;

                           while(q1!=NULL)
                           {
                                 q2=q1->next;
                                 delete q1;
                                 q1=q2;
                           }
                        }

                     }

                     ////存入数据
                     for(long tr=0;tr<thrn;tr++)
                     {
                        long tre=threm;
                        float pjr=cst.m_adValue;
                        fres=pjr;

                     }
                     /////////////vj////////////////
                     delete [] cs;
                     delete [] xr;
                     delete [] xr2;
                     ////////结果,改变sstem-1以及ioder
                     for(long u=0;u<thrn;u++)
                     {
                        long cme=threm;
                        long det=i;
                        fder=det;
                        long tr=0;
                        while(sstem!=cme&&tr<snu)
                        {
                           tr++;
                        }
                        sstem=-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;
                     long yu=i;
                     fder=yu;
                     /////////////////dtem 中 -1
                     long tt=0;
                     while(dtem!=cew&&tt<nu)
                     {
                        tt++;
                     }
                     dtem=-1;/////////////////////
               }
               nt=nt+snu;
            }
         }
         cet=ftemp;
         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),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<a)
            small=j;
   }
   if(small!=i)
   {
         temp1=a;
         temp2=b;
         a=a;
         a=temp1;
         b=b;
         b=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 static/image/common/back.gif
我做的针对四维状态空间，胞空间划分10*10*10*10;但这种划分太不细了。可是这样速度还是慢阿。
我的QQ：80 ...

师姐为什么加不了你呢？

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

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

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

页: [1]

声振论坛's Archiver

求助“胞映射”的 程序

回复 楼主 的帖子

回复 2楼 的帖子

回复 4楼 的帖子

求助“胞映射”的程序

回复楼主的帖子

回复 2楼的帖子

回复 4楼的帖子