杭州锐达数字技术有限公司
查看: 3012|回复: 9

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

[复制链接]
发表于 2008-3-18 20:04 | 显示全部楼层 |阅读模式

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

您需要 登录 才可以下载或查看,没有帐号?我要加入

x
:handshake 本人编了 一个程序,但是程序内有多个循环体套在一起,运行速度慢极了,哪位好心人能给出更好的程序?急
回复
分享到:

使用道具 举报

发表于 2008-3-18 21:16 | 显示全部楼层

回复 楼主 的帖子

能否将你的程序附上来?我也编了一个,Matlab编的,速度也一般。不知道你是做的哪一种胞映射,胞空间划分为多大?我们可以探讨。
 楼主| 发表于 2008-3-19 08:51 | 显示全部楼层

回复 2楼 的帖子

我做的针对四维状态空间,胞空间划分10*10*10*10;但这种划分太不细了。可是这样速度还是慢阿。
我的QQ:806885738,咱们交流一下
 楼主| 发表于 2008-3-19 17:15 | 显示全部楼层
我做的是简单胞映射
发表于 2008-3-19 23:06 | 显示全部楼层

回复 4楼 的帖子

我做的广义胞映射,二维的,duffing方程
 楼主| 发表于 2008-3-20 09:39 | 显示全部楼层
那你对简单胞映射应该很熟啊,给指点一下吧。主要是胞矢量坐标的选取原则(方法)。我是以区间段的个数定义胞矢量的坐标的,可以吗
发表于 2008-4-4 18:28 | 显示全部楼层
4维简单胞映射程序: 对动力系统做全局分析,程序未经验证
  1. #include <iostream.h>
  2. #include <fstream.h>
  3. #include <stdlib.h>
  4. #include <math.h>
  5. #include "mymatrix.h"

  6. const long step1=75;
  7. const long step2=90;////////////
  8. const long step3=75;
  9. const long step4=90;

  10. const long max=step1*step2*step3*step4+1;/////
  11. const float h1=5.05/step1;///////
  12. const float h2=6.06/step2;///////
  13. const float h3=5.05/step3;
  14. const float h4=6.06/step4;

  15. const float PI=3.14159265358979;
  16. float x01=-2.525;
  17. float x02=-3.03;
  18. float x03=-2.525;
  19. float x04=-3.03;


  20. long gn;     //     总组号
  21. long ugn;    ///总的不稳定组号
  22. long g_attr;
  23. struct dom
  24. {
  25.         long total;
  26.         long n;
  27. };

  28. struct gcmdata
  29. {
  30.         long Iz;
  31.         long Jz;
  32.         long gr;
  33.         long CpLoc;
  34.         long PreLoc;
  35.        
  36. };

  37. struct cpdata
  38. {
  39.         long Cz;
  40.         float Pz;
  41. };

  42. class lin
  43. {
  44. public:
  45.         lin *next;
  46.         long data;
  47.         lin();
  48. };

  49. lin::lin()
  50. {
  51.         next=NULL;
  52. }

  53. class xlink
  54. {
  55. public:
  56.         float x;
  57.         int num;
  58.         xlink *next;
  59.         xlink();
  60. };
  61. xlink::xlink()
  62. {
  63.         next=NULL;
  64. }

  65. void main()
  66. {
  67.         gn=0;
  68.         ugn=0;
  69.         long map(long z,long j,long n);   
  70.     void getbasedata();
  71.         void ab_per();
  72.         void determine();
  73.         void outputdata();
  74.         void fenzu();       
  75.     //////////////////////////////////////////////////////////////////       
  76.         ///*
  77.         getbasedata();       
  78.         ab_per();
  79.         /////////////////////////////////////////////////////////////////////////////////////
  80.        
  81.         //////////////////////////////////////////////////////////////////////////////////////
  82.         //gn=3;
  83.         //*
  84.         cout<<" start calculu"<<endl;
  85.         g_attr=gn;
  86.        
  87.         ofstream of("per_groupn.txt");
  88.         of<<"persistent group numbers: "<<gn<<endl; /////////读取数据时要参考
  89.         ////////////////////////////
  90.         if(gn<2)
  91.         {
  92.                 cout<<"there 's no stable zone"<<endl;
  93.                 exit(0);
  94.         }
  95.         ////////////////////////////
  96.         //*/
  97.     outputdata();
  98.         //determine();
  99.         //fenzu();
  100.        
  101.        
  102.         //*/
  103.        
  104. }


  105. long map(long m,long j)
  106. {
  107.         //*
  108.         float x1,x2,x3,x4;       
  109.         long z1,z2,z3,z4;
  110.         long j1,j2,j3,j4;
  111.         long number=2;
  112.         j4=(long)((j-1)/(number*number*number))+1;
  113.         j=j-(j4-1)*number*number*number;
  114.         j3=(long)((j-1)/(number*number))+1;
  115.         j=j-(j3-1)*(number*number);
  116.         j2=(long) ((j-1)/number)+1;
  117.     j1=j-(j2-1)*number;

  118.         z4=(long)((m-1)/(step1*step2*step3))+1;
  119.         m=m-(z4-1)*step1*step2*step3;
  120.         z3=(long)((m-1)/(step1*step2))+1;
  121.         m=m-(z3-1)*step1*step2;
  122.         z2=(long) ((m-1)/step1)+1;
  123.         z1=m-(z2-1)*step1;       
  124.         x1=z1*h1-h1+x01+j1*h1/number-0.5*h1/number;
  125.         x2=z2*h2-h2+x02+j2*h2/number-0.5*h2/number;
  126.         x3=z3*h3-h3+x03+j3*h3/number-0.5*h3/number;
  127.         x4=z4*h4-h4+x04+j4*h4/number-0.5*h4/number;
  128.        
  129.         //float x1,x2,x3,x4;
  130.         float h=0.08;
  131.         float mu=0.25;
  132.         float yt=0.04;
  133.         float v=0.1;
  134.         float t;
  135.         float k1,k2,k3,k4,l1,l2,l3,l4,s1,s2,s3,s4,d1,d2,d3,d4;       

  136.         ////////////////////////////
  137.         for(long i=0;i<30;i++)
  138.         {
  139.                 k1=x2;
  140.                 l1=mu*(1-x1*x1)*x2-(1+v)*x1+v*x3;
  141.                 s1=x4;
  142.                 d1=mu*(1-x3*x3)*x4+v*x1-(1+yt+v)*x3;

  143.                 k2=x2+h*l1/2;
  144.                 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);
  145.                 s2=x4+d1*h/2;
  146.                 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);


  147.                 k3=x2+h*l2/2;
  148.                 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);
  149.                 s3=x4+d2*h/2;
  150.                 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);

  151.                 k4=x2+h*l3/2;
  152.                 l4=mu*(1-(x1+h*k3)*(x1+h*k3))*(x2+h*l3)-(1+v)*(x1+h*k3)+v*(x3+h*s3);
  153.                 s4=x4+d3*h/2;
  154.                 d4=mu*(1-(x3+h*s3)*(x3+h*s3))*(x4+d3*h)+v*(x1+h*k3)-(1+yt+v)*(x3+h*s3);
  155.         
  156.                 x1=x1+h*(k1+2*k2+2*k3+k4)/6.0;
  157.                 x2=x2+h*(l1+2*l2+2*l3+l4)/6.0;
  158.                 x3=x3+h*(s1+2*s2+2*s3+s4)/6.0;
  159.                 x4=x4+h*(d1+2*d2+2*d3+d4)/6.0;
  160.             
  161.         }       
  162.    
  163.         z1=(long)((x1-x01)/h1)+1;
  164.         z2=(long)((x2-x02)/h2)+1;
  165.         z3=(long)((x3-x03)/h3)+1;
  166.         z4=(long)((x4-x04)/h4)+1;
  167.         if(z1>=step1 || z1<=0 || z2 >=step2 || z2<=0 || z3>=step3 || z3<=0 || m==0 || z4>=step4 || z4<=0)
  168.                 return 0;
  169.         else
  170.         {       
  171.                 return (z4-1)*step1*step2*step3+(z3-1)*step1*step2+(z2-1)*step1+z1;       
  172.        
  173.         }
  174.        
  175.         //*/
  176.        
  177. }

  178. void getbasedata()
  179. {       
  180.         fstream iogcm,iocp,iopre;
  181.         iocp.open("cp.dat",ios::out|ios::binary|ios::in);
  182.         iogcm.open("gcm.dat",ios::in|ios::out|ios::binary);
  183.         //ifstream igcm("gcm.dat",ios::in);
  184.         iopre.open("pre.dat",ios::out|ios::binary|ios::in);
  185.         gcmdata gcm;
  186.         cpdata cp;
  187.         ///*
  188.         long tz=0;
  189.         long snm=16;
  190.     long tempz[16][2];  ////////////记录       
  191.     //////////////////////       
  192.     for(long i=0;i<max;i++)
  193.         {
  194.                 if(i%20==0)
  195.                         cout<<"get Iz  CpLoc cp.dat"<<i<<endl;               
  196.                 for(long k=0;k<snm;k++)    /////////////每次都要置0  -1
  197.                 {
  198.                         for(long e=0;e<2;e++)
  199.                         {
  200.                                 if(e==0)
  201.                                 {
  202.                                         tempz[k][e]=-1;
  203.                                 }
  204.                                 else
  205.                                 {
  206.                                         tempz[k][e]=0;
  207.                                 }
  208.                         }
  209.                 }
  210.                
  211.                 long r=0;
  212.                 for(long j=0;j<snm;j++)
  213.                 {
  214.                         tz=map(i,j);
  215.                         r=0;
  216.                         while(tempz[r][0]!=tz && tempz[r][0]!=-1)  ////////搜索到tz 或 tz为新值
  217.                         {
  218.                                 r++;
  219.                         }
  220.                         tempz[r][0]=tz;
  221.                         tempz[r][1]=tempz[r][1]+1;
  222.                 }
  223.                
  224.         r=0;
  225.                 while(r<snm && tempz[r][0]!=-1)
  226.                 {
  227.                         r++;
  228.                 }
  229.                
  230.                 //Iz[i]=r;               
  231.                 gcm.Iz=r;
  232.                 long preIz;
  233.                 if(i==0)
  234.                 {
  235.                         //CpLoc[0]=0;                       
  236.                         gcm.CpLoc=0;
  237.                 }
  238.                 else
  239.                 {
  240.                        
  241.                         gcm.CpLoc=gcm.CpLoc+preIz;                               
  242.                 }                 
  243.                 ///
  244.                 preIz=r;
  245.                 gcm.Jz=0;
  246.                 gcm.PreLoc=0;
  247.                 gcm.gr=0;
  248.         //outgcm.seekp(i*sizeof(gcmdata));
  249.                 iogcm.write((char *)(&gcm),sizeof(gcmdata));               
  250.                 for(long t=0;t<r;t++)
  251.                 {
  252.                         cp.Cz=tempz[t][0];               
  253.                         cp.Pz=(float)(tempz[t][1])/snm;
  254.             ///////////////////
  255.                         iocp.write((char *)(&cp),sizeof(cpdata));
  256.                 }               
  257.         }
  258.        
  259.        
  260.         cout<<"The following part will get preimage data,please wait,first Jz"<<endl;
  261.         ////////////////////// preimage ////////////////////////////////////
  262.         ///iogcm.seekg(0+sizeof(long));       
  263.            for(long b=0;b<max;b++)
  264.         {
  265.                 if(b%2000==0)
  266.                         cout<<"get Jz  "<<b<<endl;
  267.                 iogcm.seekg(b*sizeof(gcmdata));
  268.                 iogcm.read((char *)(&gcm),sizeof(gcmdata));
  269.                 long I=gcm.Iz;
  270.                 for(long n=0;n<I;n++)
  271.                 {   
  272.                         long jz;
  273.                         long loc=gcm.CpLoc+n;
  274.                         iocp.seekg(loc*sizeof(cpdata));
  275.                         long cel;
  276.                         iocp.read((char *)(&cel),sizeof(long));
  277.                         iogcm.seekg(cel*sizeof(gcmdata)+sizeof(long),ios::beg);
  278.                         iogcm.read((char *)(&jz),sizeof(long));///////先读取jz                       
  279.                         jz++;
  280.                         iogcm.seekp(cel*sizeof(gcmdata)+sizeof(long),ios::beg);     
  281.                         iogcm.write((char *)(&jz),sizeof(long));                                               
  282.                 }               
  283.         }
  284.         ////////       计算   preloc    /////
  285.         cout<<endl<<"get preLoc   "<<endl;
  286.         long preloc;
  287.         for(b=0;b<max;b++)
  288.         {
  289.                 if(b%2000==0)
  290.                         cout<<"get preLoc   "<<b<<endl;
  291.                 if(b==0)
  292.                 {
  293.                         preloc=0;
  294.                         iogcm.seekp(0+4*sizeof(long));////        preLoc的位置
  295.             iogcm.write((char *)(&b),sizeof(long));
  296.                 }
  297.                 else
  298.                 {
  299.                         iogcm.seekg((b-1)*sizeof(gcmdata)+sizeof(long));////////   prejz
  300.                         long prejz;
  301.                         iogcm.read((char *)(&prejz),sizeof(long));
  302.             preloc=preloc+prejz;
  303.                         iogcm.seekp(b*sizeof(gcmdata)+4*sizeof(long));
  304.             iogcm.write((char *)(&preloc),sizeof(long));                       
  305.                 }
  306.         }       
  307.        
  308.         /////////////////////////////////////         保存 preimage cell ///
  309.         cout<<" preimage cell         "<<endl;
  310.     iogcm.seekg((max-1)*sizeof(gcmdata)+sizeof(long));
  311.         long jzmax;
  312.         iogcm.read((char *)(&jzmax),sizeof(long));
  313.         long totalpre=jzmax+preloc+1;///////////////////////////////////
  314.         /////////////////////////////////////////////////////////////////算法很糟糕/////////////////
  315.         cout<<"write preimage"<<endl;
  316.         long *pret;
  317.         pret=new long [totalpre];
  318.         for(b=0;b<totalpre;b++)
  319.         {
  320.                 if(b%10000==0)
  321.                         cout<<"write pre.dat    "<<b<<endl;
  322.                 pret[b]=-1;
  323.                
  324.         }       
  325.     cout<<" get preimage  "<<endl;
  326.         //////////////////////////////////////////read prel data
  327.         long *ptt;
  328.         ptt=new long [max];
  329.         iogcm.seekg(0);
  330.     for(long tt=0;tt<max;tt++)
  331.         {
  332.                 iogcm.read((char *)(&gcm),sizeof(gcmdata));
  333.                 ptt[tt]=gcm.PreLoc;               
  334.         }
  335.        
  336.         iocp.seekg(0);
  337.         iogcm.seekg(0);
  338.         for(b=0;b<max;b++)
  339.         {               
  340.                 if(b%500==0)
  341.                         cout<<"get preimage    "<<b<<endl;
  342.                 //cout<<b<<"             "<<b<<endl;               
  343.                 iogcm.read((char *)(&gcm),sizeof(gcmdata));
  344.                 long I=gcm.Iz;
  345.                 for(long n=0;n<I;n++)
  346.                 {               
  347.                         long cel;
  348.                         iocp.read((char *)(&cp),sizeof(cpdata));
  349.                         cel=cp.Cz;
  350.                         long preL;
  351.                         preL=ptt[cel];
  352.                         long t=0;
  353.                         long pce;                       
  354.                         pce=pret[preL+t];
  355.                         while(pce!=-1)
  356.                         {
  357.                                 t=t+1;
  358.                                 pce=pret[preL+t];
  359.                         }
  360.                         pret[preL+t]=b;                       
  361.                 }
  362.         }
  363.         ////////////write data
  364.         for(tt=0;tt<totalpre;tt++)
  365.         {
  366.                 if(tt%2000==0)
  367.                         cout<<"save predata    "<<tt<<endl;
  368.                 iopre.write((char *)(&pret[tt]),sizeof(long));
  369.         }   
  370.        
  371.     delete [] ptt;
  372.         delete [] pret;
  373.         iogcm.close();
  374.         iocp.close();
  375.         iopre.close();       
  376. }

  377. void ab_per()
  378. {       
  379.     ////select simple map group
  380.        
  381.         fstream iosgr,iogcm1,iocp1,ioscm,iopersistent;
  382.         iosgr.open("sgr.dat",ios::in|ios::out|ios::binary);
  383.     iogcm1.open("gcm.dat",ios::in|ios::out|ios::binary);
  384.         iocp1.open("cp.dat",ios::in|ios::out|ios::binary);
  385.         ioscm.open("tempscm.dat",ios::in|ios::out|ios::binary);////保存scm的结果
  386.         iopersistent.open("persistent.dat",ios::in|ios::out|ios::binary);
  387.         gcmdata gcm;
  388.         cpdata cp;
  389.         ///*
  390.         cout<<"simple cell method"<<endl;
  391.         long sgr;
  392.         long sgn=0;       
  393.         {
  394.                 long temz[10000000];         //==========================????????????????高维时注意
  395.                 for(long i=0;i<max;i++)
  396.                 {
  397.                         sgr=0;
  398.                         iosgr.write((char *)(&sgr),sizeof(long));
  399.                         //sgr[i]=0;
  400.                 }               
  401.                 long m,b;
  402.                 long g=0;
  403.                 for(long j=0;j<max;j++)
  404.                 {
  405.                         if(j%1000==0)
  406.                                 cout<<j<<"       simple cell method, get sgn  "<<endl;
  407.                         m=0;
  408.                         iosgr.seekg(j*sizeof(long));
  409.             iosgr.read((char *)(&sgr),sizeof(long));
  410.                         if(sgr!=0)
  411.                         {
  412.                                 continue;  ////////////////不是virgin cell,重新选取
  413.                         }
  414.                         else
  415.                         {   
  416.                                 sgr=-1;
  417.                                 iosgr.seekp(j*sizeof(long));
  418.                                 iosgr.write((char *)(&sgr),sizeof(long));                               
  419.                                 b=j;
  420.                                 temz[m]=b;
  421.                                 do
  422.                                 {
  423.                                         sgr=-1;
  424.                                         iosgr.seekp(b*sizeof(long));
  425.                                         iosgr.write((char *)(&sgr),sizeof(long));
  426.                                         m=m+1;
  427.                                         ////////////////
  428.                                         iogcm1.seekg(b*sizeof(gcmdata));
  429.                                         iogcm1.read((char *)(&gcm),sizeof(gcmdata));
  430.                                         long loc=gcm.CpLoc;
  431.                                         iocp1.seekg(loc*sizeof(cpdata));
  432.                                         iocp1.read((char *)(&cp),sizeof(cpdata));
  433.                                         b=cp.Cz;       
  434.                                         //b=Cz[CpLoc[b]];//////select the first image cell                               
  435.                                         temz[m]=b;
  436.                                         iosgr.seekg(b*sizeof(long));
  437.                                         iosgr.read((char *)(&sgr),sizeof(long));
  438.                                 }while(sgr==0);  ///                         
  439.                                 if(sgr==-1)     ///////////////新周期解
  440.                                 {
  441.                                         ///////////寻找i/////////(刚好在周期轨道上)
  442.                                         g++;
  443.                                         sgn=g;
  444.                                         long u=0;
  445.                                         while(temz[u]!=b)
  446.                                         {
  447.                                                 u++;
  448.                                         }
  449.                                         for(long q=0;q<m;q++)
  450.                                         {
  451.                                                 if(q<u)
  452.                                                 {
  453.                                                         iosgr.seekg(temz[q]*sizeof(long));
  454.                                                         long c=-2;
  455.                                                         iosgr.write((char *)(&c),sizeof(long));
  456.                                                         //sgr[temz[q]]=-2;
  457.                                                 }
  458.                                                 else
  459.                                                 {
  460.                                                         iosgr.seekp(temz[q]*sizeof(long));
  461.                             iosgr.write((char *)(&g),sizeof(long));
  462.                                                         long tt=temz[q];       
  463.                                                         ioscm.write((char *)(&g),sizeof(long));/////////先写group   
  464.                                                         ioscm.write((char *)(&tt),sizeof(long));                                                       
  465.                                                 }
  466.                                         }
  467.                                        
  468.                                 }
  469.                                 else         //////吸引到已知一点上
  470.                                 {
  471.                                         for(long k=0;k<m;k++)
  472.                                         {
  473.                                                 iosgr.seekg(temz[k]*sizeof(long));
  474.                                                 long c=-2;
  475.                         iosgr.write((char *)(&c),sizeof(long));                                                                                 
  476.                                         }
  477.                                 }
  478.                         }       
  479.                 }               
  480.         }       
  481.         cout<<"search for persistent groups"<<endl;
  482.         ////////////////////////////////////////////////////////////////////////////////////////
  483.         {
  484.                 long temp[10000000];          ////======================??????????????save candidate cells
  485.                 long n=0;
  486.                 long flag=1;
  487.                 //ioscm.clear();               
  488.                 for(long i=1;i<=sgn;i++)
  489.                 {
  490.                         cout<<sgn<<"          "<<i<<endl;
  491.                         ///////////////////////////////先判断组号是否已经是absorbing group了
  492.                         n=0;
  493.                         ioscm.seekg(0);
  494.                         long cel,gg;
  495.                         ioscm.read((char *)(&gg),sizeof(long));
  496.                         ioscm.read((char *)(&cel),sizeof(long));
  497.                         long cur,pend;
  498.                         cur=ioscm.tellg();
  499.                         ioscm.seekg(0,ios::end);
  500.                         pend=ioscm.tellg();
  501.                         ioscm.seekg(cur);
  502.                         while(pend!=ioscm.tellg())
  503.                         {
  504.                                 if(gg==i)
  505.                                 {
  506.                                         temp[n]=cel;
  507.                                         n++;                               
  508.                                 }
  509.                                 ioscm.read((char *)(&gg),sizeof(long));
  510.                                 ioscm.read((char *)(&cel),sizeof(long));                               
  511.                         }
  512.                         ioscm.seekg(0);
  513.                         cur=ioscm.tellg();
  514.                         ioscm.seekg(0,ios::end);
  515.                         pend=ioscm.tellg();
  516.                         ioscm.seekg(cur);
  517.                         n=0;
  518.                         while(pend!=ioscm.tellg())
  519.                         {
  520.                                 ioscm.read((char *)(&gg),sizeof(long));
  521.                                 ioscm.read((char *)(&cel),sizeof(long));
  522.                                 if(gg==i)
  523.                                 {
  524.                                         temp[n]=cel;                                          
  525.                                         iogcm1.seekg(cel*sizeof(gcmdata));
  526.                                         iogcm1.read((char *)(&gcm),sizeof(gcmdata));
  527.                                         gcm.gr=-3;
  528.                                         iogcm1.seekp(cel*sizeof(gcmdata));
  529.                                         iogcm1.write((char *)(&gcm),sizeof(gcmdata));
  530.                                         n++;
  531.                                 }                                  
  532.                         }                       
  533.                         /////////////////// to determine whether each cell's image cells are in the temp[n]
  534.                         long prenum,afnum;
  535.                         prenum=0;
  536.                         afnum=1;         
  537.                         while(afnum>prenum)  ///////// to locate all candidate cells
  538.                         {
  539.                                 prenum=n;
  540.                                 for(long k=0;k<prenum;k++)
  541.                                 {
  542.                                         long te=temp[k];                                          
  543.                                         iogcm1.seekg(te*sizeof(gcmdata));
  544.                                         long Iz;
  545.                                         iogcm1.read((char *)(&gcm),sizeof(gcmdata));
  546.                                         Iz=gcm.Iz;
  547.                                         long Loc2=gcm.CpLoc;
  548.                                         for(long kk=0;kk<Iz;kk++)
  549.                                         {
  550.                                                 long Loc=Loc2+kk;
  551.                                                 iocp1.seekg(Loc*sizeof(cpdata));
  552.                                                 long image;
  553.                                                 iocp1.read((char *)(&image),sizeof(long));                                               
  554.                                                 long ggimage;
  555.                                                 iogcm1.clear();
  556.                                                 iogcm1.seekg((image)*sizeof(gcmdata));
  557.                                                 iogcm1.read((char *)(&gcm),sizeof(gcmdata));
  558.                                                 ggimage=gcm.gr;                                               
  559.                                                 if(ggimage==0)
  560.                                                 {
  561.                                                         iogcm1.seekg(image*sizeof(gcmdata));
  562.                                                         iogcm1.read((char *)(&gcm),sizeof(gcmdata));
  563.                                                         gcm.gr=-3;
  564.                                                         iogcm1.seekp(image*sizeof(gcmdata));
  565.                                                         iogcm1.write((char *)(&gcm),sizeof(gcmdata));
  566.                                                         temp[n]=image;
  567.                                                         n++;
  568.                                                 }
  569.                         //*
  570.                                                 iosgr.seekg(image*sizeof(long));
  571.                                                 long sgg;
  572.                                                 iosgr.read((char *)(&sgg),sizeof(long));
  573.                                                 if(sgg>0&&sgg!=i)
  574.                                                 {
  575.                                                         flag=0;
  576.                                                         kk=Iz;
  577.                                                         k=prenum;
  578.                                                         prenum=n;
  579.                                                 }
  580.                                                
  581.                                         }
  582.                                 }
  583.                                 afnum=n;                       
  584.                         }
  585.                         ///////////////////////////////////////     determine
  586.                         if(flag!=0)
  587.                         {
  588.                                 for(long k=0;k<n;k++)
  589.                                 {
  590.                                         long tee=temp[k];
  591.                                         iogcm1.seekg(tee*sizeof(gcmdata));
  592.                                         long Iz;
  593.                                         iogcm1.read((char *)(&gcm),sizeof(gcmdata));
  594.                                         Iz=gcm.Iz;
  595.                                         long Loc2;
  596.                                         Loc2=gcm.CpLoc;
  597.                                         for(long as=0;as<Iz;as++)
  598.                                         {
  599.                                                 long Loc=Loc2+as;
  600.                                                 iocp1.seekg(Loc*sizeof(cpdata));
  601.                                                 long iima;
  602.                                                 iocp1.read((char *)(&iima),sizeof(long));
  603.                                                 iogcm1.seekg(iima*sizeof(gcmdata)+2*sizeof(long));
  604.                                                 long ggg;
  605.                                                 iogcm1.read((char *)(&ggg),sizeof(long));
  606.                                                 if(ggg!=-3)
  607.                                                         flag=0;
  608.                                         }                               
  609.                                 }
  610.                         }
  611.                        
  612.                         if(flag!=0)         //  发现 persistent group  //
  613.                         {
  614.                                 gn++;
  615.                                 for(long h=0;h<n;h++)
  616.                                 {
  617.                                         iogcm1.seekp(temp[h]*sizeof(gcmdata)+2*sizeof(long));
  618.                                         iogcm1.write((char *)(&gn),sizeof(long));
  619.                                         long celp=temp[h];
  620.                                         iopersistent.write((char *)(&celp),sizeof(long));
  621.                                         iopersistent.write((char *)(&gn),sizeof(long));                                          
  622.                                 }                               
  623.                         }       
  624.                         else
  625.                         {
  626.                                 for(long ui=0;ui<n;ui++)
  627.                                 {                                       
  628.                                         iogcm1.seekg(temp[ui]*sizeof(gcmdata));
  629.                                         iogcm1.read((char *)(&gcm),sizeof(gcmdata));
  630.                                         gcm.gr=0;
  631.                                         iogcm1.seekp(temp[ui]*sizeof(gcmdata));
  632.                                         iogcm1.write((char *)(&gcm),sizeof(gcmdata));                                          
  633.                                 }
  634.                         }
  635.                         /////////  每次完后,gcm中gr应该清零,否则影响下一个group的判断  ////////                       
  636.                        
  637.                 }               
  638.         }
  639.         //*/
  640.     ///////////////////////////////////////////////////////
  641.        
  642.        
  643.        
  644. }
复制代码
发表于 2008-4-4 18:28 | 显示全部楼层
  1. void determine()
  2. {   
  3.     void sort(long a[],long b[],long n);
  4.     CMatrix mat;
  5.     CMatrix cst;   
  6.     long *flag;
  7.     flag=new long [max];
  8.     for(long i=0;i<max;i++)
  9.     {
  10.         if(i%2000==0)
  11.             cout<<i<<endl;        
  12.         flag[i]=-1;
  13.     }   
  14.     long tem[500000];/////////保存persistent cells     =================================????????????
  15.     long n=0;
  16.     fstream iopersistent;
  17.     iopersistent.open("persistent.dat",ios::in|ios::out|ios::binary);  //////////   
  18.     fstream iogcm,iocp;
  19.     iogcm.open("gcm.dat",ios::in|ios::out|ios::binary);
  20.     iocp.open("cp.dat",ios::in|ios::out|ios::binary);
  21.     gcmdata gcm;
  22.     cpdata cp;   
  23.     fstream iopre;
  24.     iopre.open("pre.dat",ios::out|ios::binary|ios::in);
  25.     long preloc;
  26.     iogcm.seekg((max-1)*sizeof(gcmdata));
  27.     long jzmax;
  28.     iogcm.read((char *)(&gcm),sizeof(gcmdata));
  29.     preloc=gcm.PreLoc;
  30.     jzmax=gcm.Jz;
  31.     long totalpre=jzmax+preloc+1;
  32.     long *fzpre;
  33.     fzpre=new long [totalpre];
  34.     iopre.seekg(0);
  35.     for(long ttr=0;ttr<totalpre;ttr++)
  36.     {
  37.         iopre.read((char *)(&fzpre[ttr]),sizeof(long));
  38.     }
  39.     //////////////////////////////////////////////   
  40.     iogcm.seekg((max-1)*sizeof(gcmdata));   
  41.     iogcm.read((char *)(&gcm),sizeof(gcmdata));
  42.     long Izmax;
  43.     Izmax=gcm.Iz;
  44.     long cpclo;
  45.     cpclo=gcm.CpLoc;
  46.     long totalcl;
  47.     totalcl=Izmax+cpclo+1;
  48.     //////
  49.     cpdata *fcpl;
  50.     fcpl=new cpdata [totalcl];
  51.     iocp.seekg(0);
  52.     for(ttr=0;ttr<totalcl;ttr++)
  53.     {
  54.         iocp.read((char *)(&fcpl[ttr]),sizeof(cpdata));
  55.     }
  56.     ////////////////////////////////////////////
  57.     gcmdata *fgcm;
  58.     fgcm=new gcmdata [max];
  59.     iogcm.seekg(0);
  60.     for(ttr=0;ttr<max;ttr++)
  61.     {
  62.         iogcm.read((char *)(&fgcm[ttr]),sizeof(gcmdata));
  63.     }
  64.    
  65.     //////////////////////// 输出 结果 的文件 //////////////////////////////////////////////////
  66.    
  67.    
  68.     //////////////初始化result
  69.     //long res=(gn+1)*sizeof(long)+gn*sizeof(float);
  70.     float fres[max][3];
  71.     for(i=0;i<max;i++)
  72.     {
  73.         if(i%2000==0)
  74.             cout<<i<<endl;
  75.         long tttt=0;
  76.         float tff=0.0;
  77.         ///ioresult.write((char *)(&tttt),sizeof(long));        
  78.         for(long j=0;j<3;j++)
  79.         {
  80.             fres[i][j]=0.0;
  81.         }        
  82.     }
  83.     fstream iodm;             ////// 保存 cell 的, domicle 的数量,做分组用  
  84.     iodm.open("todomiledata.dat",ios::in|ios::out|ios::binary);
  85.     dom dm;
  86.     dm.n=0;
  87.     dm.total=0;
  88.     //fstream ioder;
  89.     //ioder.open("derminedornot.dat",ios::in|ios::out|ios::binary);
  90.     long *fder;
  91.     fder=new long [max];
  92.     ///fstream iocw;  //////////成尾数   计算是用
  93.     //iocw.open("chengwei.dat",ios::in|ios::out|ios::binary);
  94.     for(i=0;i<max;i++)
  95.     {   /////////////////所有的cell 数量, total number  persistent group 也当成0
  96.         if(i%2000==0)
  97.             cout<<i<<endl;
  98.         iodm.write((char *)(&dm),sizeof(dom)); /////////先保存
  99.         long ffff=0;
  100.         fder[i]=ffff;        
  101.     }   
  102.     /////////////////////////////////////////////////////////////////////////////////
  103.    
  104.    
  105.     //////////////////////////////////////////////////////////////////////////////
  106.     long *ftemp;

  107.     ftemp=new long [max];
  108.     long ftn=0;

  109.     for(i=2;i<gn+1;i++)      ////////////1为sink cell 只要取得其吸引域就可以了,不要算vj pj
  110.     {
  111.         ftn=0;
  112.         ///////// 确定 此组的cell
  113.         //iotemp.close();
  114.         long pend;
  115.         //cur=iopersistent.tellg();
  116.         iopersistent.seekg(0,ios::end);
  117.         pend=iopersistent.tellg();
  118.         iopersistent.seekg(0);
  119.         //long per=0;
  120.         //////////保留persistent cells
  121.         n=0;
  122.         long dfj=-1;///////不同批次preimage 的分界线
  123.         ////////////每次迭代计算preimage开始的地方
  124.         //long nc=0;  //////////transcient cells 分成 nc
  125.         while(pend!=iopersistent.tellg())
  126.         {
  127.             long cee,grg;
  128.             iopersistent.read((char *)(&cee),sizeof(long));
  129.             iopersistent.read((char *)(&grg),sizeof(long));
  130.             if(grg==i)
  131.             {
  132.                 /////改变 flag     
  133.                 long ff=i;
  134.                 float pp=1.0;
  135.                 flag[cee]=ff;
  136.                 tem[n]=cee;               
  137.                 n++;   
  138.                 ////////////////////////////////////  persistent cell为已知胞  ////////////////////////
  139.                 fder[cee]=ff;   
  140.                                 fres[cee][i-1]=pp;
  141.                 ////////domicle
  142.                 iodm.seekg(cee*sizeof(dom));
  143.                 iodm.read((char *)(&dm),sizeof(dom));
  144.                 dm.n=0; ///////     用零来区分是否是 persistent cells
  145.                 dm.total=dm.total+i;
  146.                 iodm.seekp(cee*sizeof(dom));
  147.                 iodm.write((char *)(&dm),sizeof(dom));               
  148.                 /////////////////////////////////////////////////////////////////////////////////////
  149.             }
  150.         }
  151.         ///////  注意,保存到temp中的transcient cell 的文件必需 每次循环时清空
  152.         //iotemp.close();
  153.         //iotemp.open("temptemp.dat",ios::trunc|ios::in|ios::out|ios::binary);  //////////
  154.         ////
  155.         
  156.         long numb=0;
  157.         long handnum=0;
  158.         cout<<"first level preimage"<<endl;
  159.         for(long k=0;k<n;k++)
  160.         {
  161.             long cee;
  162.             cee=tem[k];   ///////// persistent cells
  163.             long Jz;
  164.             gcm=fgcm[cee];            
  165.             Jz=gcm.Jz;
  166.             long loc2=gcm.PreLoc;
  167.             for(long p=0;p<Jz;p++)
  168.             {
  169.                 long loc=loc2+p;
  170.                 long ffg,iima;
  171.                 iima=fzpre[loc];
  172.                 long dis;
  173.                 if(iima<0)///////////////////////////不取进来
  174.                 {
  175.                     ffg=i;
  176.                 }
  177.                 else
  178.                 {
  179.                 dis=fder[iima];
  180.                 ffg=flag[iima];
  181.                 }
  182.                 ////////////////////////////////////////
  183.                 //long ddd;
  184.                
  185.                 ///////////////////////////////////////
  186.                 if(ffg!=i&&dis!=i)
  187.                 {                    
  188.                     //iotemp.write((char *)(&iima),sizeof(long));
  189.                     ftemp[ftn]=iima;
  190.                     ftn++;
  191.                     flag[iima]=i;                    
  192.                     iodm.seekg(iima*sizeof(dom));
  193.                     iodm.read((char *)(&dm),sizeof(dom));
  194.                     dm.n++;
  195.                     dm.total=dm.total+i;
  196.                     iodm.seekp(iima*sizeof(dom));
  197.                     iodm.write((char *)(&dm),sizeof(dom));
  198.                     numb++;
  199.                 }               
  200.             }
  201.         }         
  202.         //iotemp.write((char *)(&dfj),sizeof(long));/////第一批与后面的分界
  203.         ftemp[ftn]=-1;
  204.         ftn++;
  205.         ///////////////////////////////迭代
  206.         long afn,bfn;
  207.         afn=0;
  208.         bfn=-1;
  209.         long sta=0;
  210.         while(afn>bfn)
  211.         {
  212.             cout<<"preimage of persistent  an"<<afn<<"   bn "<<bfn<<endl;
  213.             bfn=afn;            
  214.             //iotemp.seekg(dbe*sizeof(long));
  215.             long cee;            
  216.             cee=ftemp[sta];
  217.             //iotemp.read((char *)(&cee),sizeof(long));            
  218.             sta++;        
  219.             while(cee!=-1)
  220.             {               
  221.                 long Jz;
  222.                 gcm=fgcm[cee];               
  223.                 Jz=gcm.Jz;
  224.                 long loc2=gcm.PreLoc;
  225.                 for(long p=0;p<Jz;p++)
  226.                 {
  227.                     long loc=loc2+p;
  228.                     long fg,ima;
  229.                     ima=fzpre[loc];
  230.                     fg=flag[ima];
  231.                     if(fg!=i)
  232.                     {        
  233.                         numb++;
  234.                         afn++;
  235.                         ftemp[ftn]=ima;
  236.                         ftn++;                        
  237.                         flag[ima]=i;
  238.                         iodm.seekg(ima*sizeof(dom));
  239.                         iodm.read((char *)(&dm),sizeof(dom));
  240.                         dm.n++;
  241.                         dm.total=dm.total+i;
  242.                         iodm.seekp(ima*sizeof(dom));
  243.                         iodm.write((char *)(&dm),sizeof(dom));
  244.                     }                    
  245.                 }               
  246.                 cee=ftemp[sta];
  247.                 sta++;               
  248.             }
  249.             ///////// 次批完
  250.             ftemp[ftn]=-1;
  251.             ftn++;            
  252.         }
  253.         
  254.         ////////上面获得了 此persistent group的 所有pre象,并已经按批次存放,用-1分隔开来  iotemp中
  255.         cout<<numb<<endl;
  256.         long dtem[5000000];     ///////////      保存临时像/ 第i批数据/////////////?????????????????
  257.         long nu=0;        
  258.         long cet;
  259.         sta=0;
  260.         cet=ftemp[sta];   
  261.         sta++;
  262.         while(sta<ftn)
  263.         {
  264.             cout<<" get i'th cell of i'th preimage"<<endl;
  265.             ///////       获取第i批数据及像       ////////  
  266.             nu=0;
  267.             long dtm3;
  268.             if(cet!=-1)       ////////文件指针不能移出去
  269.             {
  270.                 dtm3=fder[cet];
  271.             }
  272.             if(cet!=-1&&dtm3!=i)   
  273.             {
  274.                 dtem[nu]=cet;
  275.                 nu++;
  276.             }            
  277.             cet=ftemp[sta];
  278.             sta++;
  279.             while(cet!=-1)  //////////            取得此大批preimage
  280.             {               
  281.                 long dtm;
  282.                 dtm=fder[cet];               
  283.                 if(dtm!=i)        /////////好像肯定不等于的??????????????
  284.                 {
  285.                     dtem[nu]=cet;            /////////已经处理也要放进取,避免它的下几层胞丢失   
  286.                     nu++;
  287.                 }               
  288.                 cet=ftemp[sta];
  289.                 sta++;
  290.             }
  291.             ///////////////////////////////////////////////////////////////
  292.             ////////////////  取像  只取此pg 的吸引域里面的///////////
  293.             if(nu>0)
  294.             {               
  295.                 long an,bn;
  296.                 an=0;
  297.                 bn=-1;
  298.                 long psta=0;
  299.                 while(an>bn)
  300.                 {
  301.                     cout<<" i'th preimage an "<<an<<"   bn "<<bn<<endl;
  302.                     bn=an;
  303.                     for(long u=psta;u<nu;u++)   //////////获得每个像
  304.                     {
  305.                         if(nu>10000)
  306.                         {
  307.                             if(u%500==0)
  308.                                 cout<<u<<"           "<<nu<<endl;
  309.                         }
  310.                         long cre=dtem[u];
  311.                         psta++;
  312.                         long dt; ////////////// 考察是否已经确定
  313.                         gcm=fgcm[cre];                        
  314.                         long Iz;
  315.                         Iz=gcm.Iz;
  316.                         long loc2=gcm.CpLoc;
  317.                         for(long y=0;y<Iz;y++)
  318.                         {
  319.                             long loc=loc2+y;
  320.                             cp=fcpl[loc];                           
  321.                             long cim;
  322.                             cim=cp.Cz;
  323.                             //////// whether determined
  324.                             dt=fder[cim];                           
  325.                             long tfgg;
  326.                             tfgg=flag[cim];
  327.                             if(dt!=i&&tfgg==i) //////  防止边界上,其它域里面的胞也保存进去
  328.                             {
  329.                                 //////看是否已经保存了
  330.                                 long t=0;
  331.                                 while(dtem[t]!=cim&&t<nu)
  332.                                 {
  333.                                     t++;
  334.                                 }
  335.                                 if(t==nu) /////////没有保存
  336.                                 {
  337.                                     dtem[nu]=cim;
  338.                                     nu++;
  339.                                     an++;
  340.                                 }
  341.                             }
  342.                         }
  343.                     }                    
  344.                 }               
  345.             }
  346.             /////////////////////                                  计算
  347.             /////////此大批共有nu个,都保存到dtem[]中了 ////////////////
  348.             if(nu>0)
  349.             {
  350.                 //////////         取得第i小批      //////////////////
  351.                 long nt=0;         /////已经处理cell数量
  352.                 while(nt<nu)     ////////////////表示还没有处理完此大批的cell
  353.                 {
  354.                     long cel;
  355.                     long tt=nu-1;
  356.                     if((nu-nt)<20000)
  357.                     {
  358.                         tt=0;
  359.                         while(dtem[tt]==-1)    /////////用-1 表示已经处理过的cell
  360.                         {
  361.                             tt++;
  362.                         }
  363.                     }
  364.                     else
  365.                     {
  366.                         tt=nu-1;
  367.                         while(dtem[tt]==-1)    /////////用-1 表示已经处理过的cell
  368.                         {
  369.                             tt--;
  370.                         }
  371.                     }
  372.                     cel=dtem[tt]; ///////// 碰到一个没有处理的cell
  373.                     ////获取围着它的小批,注意只保存此组里面的, 边界 cell
  374.                     long stem[5000000]; //////////////??????????????????????????????????????????????/      
  375.                     long snu=0;
  376.                     long afn=0;
  377.                     long bfn=-1;
  378.                     stem[0]=cel;
  379.                     snu++;
  380.                     long psta=0;
  381.                     while(afn>bfn)  ///////    preimage   与    image of cel
  382.                     {
  383.                         cout<<"small  i'th afn "<<afn<<endl;
  384.                         bfn=afn;
  385.                         ///               long tbg=0;                        
  386.                         for(long k=psta;k<snu;k++)  //////////
  387.                         {        
  388.                             if(snu>10000)
  389.                             {
  390.                                 if(k%500==0)
  391.                                     cout<<k<<"           "<<snu<<endl;
  392.                             }
  393.                             long tce;
  394.                             tce=stem[k];
  395.                             psta++;
  396.                             ////////   考察preimage
  397.                             gcm=fgcm[tce];
  398.                             long Jz,Iz;                           
  399.                             Iz=gcm.Iz;                           
  400.                             long lociz=gcm.CpLoc;
  401.                             Jz=gcm.Jz;
  402.                             long locjz=gcm.PreLoc;
  403.                            
  404.                             for(long t=0;t<Iz;t++)
  405.                             {
  406.                                 long loc=lociz+t;
  407.                                 long icel;
  408.                                 cp=fcpl[loc];                                
  409.                                 icel=cp.Cz;
  410.                                 ////////////查看是否是其他吸引域里面的,因为边界容易出现这种情况
  411.                                 long tt=0;
  412.                                 while(stem[tt]!=icel&&tt<snu)
  413.                                 {
  414.                                     tt++;
  415.                                 }
  416.                                 if(tt==snu)/////查看是否是其他吸引域里面
  417.                                 {
  418.                                     long iflgf;
  419.                                     iflgf=flag[icel];
  420.                                     //////////已经确定的cell也不能加进去                                    
  421.                                     long ddi;                                    
  422.                                     ddi=fder[icel];
  423.                                     
  424.                                     
  425.                                     if(iflgf==i&&ddi!=i) ////
  426.                                     {
  427.                                         stem[snu]=icel;
  428.                                         snu++;
  429.                                         afn++;
  430.                                     }
  431.                                 }
  432.                             }                           
  433.                         }                        
  434.                     }
  435.                     /////////////////   此小批已经获得,放在stem[]  ,  有snu个 /////////////////////////////
  436.                     ///////// snu 个stem排序
  437.                     long *sstem;
  438.                     sstem=new long [snu];
  439.                     
  440.                     for(long uu=0;uu<snu;uu++)
  441.                     {
  442.                         sstem[uu]=stem[uu];
  443.                     }                    
  444.                     
  445.                     /////已经胞排好了,放在sstem[snu]中
  446.                     ////开始计算                    
  447.                     long threm[5000000];
  448.                     long thrn=0;
  449.                     long nhand=0;
  450.                     while(nhand<snu)
  451.                     {
  452.                         thrn=0;
  453.                         long cre;
  454.                         long tt=snu-1;
  455.                         if((snu-nhand)<200)
  456.                         {
  457.                             tt=0;                           
  458.                             while(sstem[tt]==-1)///////////////从最后一个开始算
  459.                             {
  460.                                 tt++;
  461.                             }
  462.                         }
  463.                         else
  464.                         {
  465.                             tt=snu-1;                           
  466.                             while(sstem[tt]==-1)///////////////从最后一个开始算
  467.                             {
  468.                                 tt--;
  469.                             }
  470.                         }
  471.                         cre=sstem[tt];
  472.                         threm[thrn]=cre;
  473.                         thrn++;
  474.                         long an=0;
  475.                         long bn=-1;
  476.                         long psta=0;
  477.                         while(an>bn)
  478.                         {
  479.                             bn=an;
  480.                             for(long u=psta;u<thrn;u++)
  481.                             {
  482.                                 long crm=threm[u];
  483.                                 psta++;
  484.                                 long Iz,loc,loc1;
  485.                                 gcm=fgcm[crm];
  486.                                 Iz=gcm.Iz;
  487.                                 loc=gcm.CpLoc;
  488.                                 for(long t=0;t<Iz;t++)
  489.                                 {
  490.                                     loc1=loc+t;
  491.                                     cp=fcpl[loc1];                                    
  492.                                     long cie=cp.Cz;
  493.                                     long det;                                    
  494.                                     det=fder[cie];
  495.                                     long flt;
  496.                                     flt=flag[cie];
  497.                                     if(det!=i&&flt==i)
  498.                                     {
  499.                                         //////保存了吗
  500.                                         long ttt=0;
  501.                                         while(threm[ttt]!=cie&&ttt<thrn)
  502.                                         {
  503.                                             ttt++;
  504.                                         }
  505.                                         if(ttt==thrn)
  506.                                         {
  507.                                             an++;
  508.                                             threm[thrn]=cie;
  509.                                             thrn++;
  510.                                         }
  511.                                     }
  512.                                 }
  513.                             }
  514.                         }
  515.                         //////////////不稳定组,   //////////
  516.                                                 
  517.                         //////计算,已经放到threm[thrn]中
  518.                         ////////pj////////////////
  519.                         
  520.                         handnum=handnum+thrn;
  521.                         cout<<"total  : "<<numb<<"         handled : "<<handnum<<"        equation : "<<thrn<<endl;                    
  522.                         float *xr,*xr2;
  523.                         xr=new float [thrn];/////////迭代用
  524.                         xr2=new float [thrn];
  525.                         float err=1.0;
  526.                         for(k=0;k<thrn;k++)
  527.                         {
  528.                             xr[k]=0.0;
  529.                         }
  530.                         xr[0]=0.65;   ////////////////初值
  531.                         if(thrn<1000)
  532.                         {
  533.                             mat.Realloc(thrn,thrn);
  534.                             mat=0.0;
  535.                         }
  536.                         xlink **head;
  537.                         xlink *p;
  538.                         head=new xlink *[thrn];
  539.                         if(thrn>=1000)
  540.                         {
  541.                             for(int j=0;j<thrn;j++)
  542.                             {
  543.                                 head[j]=new xlink [1];
  544.                                 if(j>=1)
  545.                                     head[j]->num=0;
  546.                                 
  547.                             }
  548.                         }
  549.                         
  550.                         cst.Realloc(thrn,1);
  551.                         
  552.                         cst=0.0;
  553.                         float *cs;
  554.                         cs=new float [thrn];
  555.                         if(thrn<1000)
  556.                         {
  557.                             for(k=0;k<thrn;k++)
  558.                             {
  559.                                 if(thrn>1000&&thrn%500==0)
  560.                                     cout<<"get xishu  "<<k<<"      "<<thrn<<endl;
  561.                                 cs[k]=0.0;
  562.                                 long ce=threm[k];
  563.                                 long Iz;
  564.                                 gcm=fgcm[ce];
  565.                                 Iz=gcm.Iz;
  566.                                 long loc2=gcm.CpLoc;
  567.                                 mat.m_adValue[k][k]=1.0;/////////   别忘记了
  568.                                 for(long t=0;t<Iz;t++)
  569.                                 {
  570.                                     long loc=loc2+t;
  571.                                     cp=fcpl[loc];                    
  572.                                     
  573.                                     long ceim=cp.Cz;////////
  574.                                     long tt=0;
  575.                                     while(ceim!=threm[tt]&&tt<thrn) ///////  ????????tt<tnb
  576.                                     {
  577.                                         tt++;
  578.                                     }                                
  579.                                     if(tt==thrn) /////////象为已知胞, 也可能为未知胞
  580.                                     {    ////看 到
  581.                                         long ttfg;
  582.                                         ttfg=flag[ceim];                                    
  583.                                         long dfdf;
  584.                                         dfdf=fder[ceim];                                    
  585.                                         if(ttfg==i&&dfdf==i)
  586.                                         {                                            
  587.                                             float pj;                                            
  588.                                                                                         pj=fres[ceim][i-1];
  589.                                             cs[k]=cs[k]+cp.Pz*pj;  ////常数项
  590.                                         }
  591.                                     }
  592.                                     else
  593.                                     {   
  594.                                         if(tt==k)////////与ce相同
  595.                                         {
  596.                                             ////ce 到 ce  的 概率
  597.                                             mat.m_adValue[k][k]=mat.m_adValue[k][k]-cp.Pz;
  598.                                         }
  599.                                         else
  600.                                         {
  601.                                             mat.m_adValue[k][tt]=-cp.Pz;
  602.                                         }
  603.                                     }
  604.                                 }                           
  605.                             }
  606.                         }
  607.                         else
  608.                         {
  609.                             for(k=0;k<thrn;k++)
  610.                             {
  611.                                 if(thrn>1000&&thrn%500==0)
  612.                                     cout<<"get xishu  "<<k<<"      "<<thrn<<endl;
  613.                                 cs[k]=0.0;
  614.                                 long ce=threm[k];
  615.                                 long Iz;
  616.                                 gcm=fgcm[ce];
  617.                                 Iz=gcm.Iz;
  618.                                 long loc2=gcm.CpLoc;                                
  619.                                 p=head[k];
  620.                                 p->num=k;
  621.                                 p->x=1.0;
  622.                                 //mat.m_adValue[k][k]=1.0;/////////   别忘记了
  623.                                 for(long t=0;t<Iz;t++)
  624.                                 {
  625.                                     long loc=loc2+t;
  626.                                     cp=fcpl[loc];                                    
  627.                                     long ceim=cp.Cz;////////
  628.                                     long tt=0;
  629.                                     while(ceim!=threm[tt]&&tt<thrn) ///////  ????????tt<tnb
  630.                                     {
  631.                                         tt++;
  632.                                     }                                
  633.                                     if(tt==thrn) /////////象为已知胞, 也可能为未知胞
  634.                                     {    ////看 到
  635.                                         long ttfg;
  636.                                         ttfg=flag[ceim];                                    
  637.                                         long dfdf;
  638.                                         dfdf=fder[ceim];                                    
  639.                                         if(ttfg==i&&dfdf==i)
  640.                                         {                                            
  641.                                             float pj;                                            
  642.                                                                                         pj=fres[ceim][i-1];
  643.                                             cs[k]=cs[k]+cp.Pz*pj;  ////常数项
  644.                                         }
  645.                                     }
  646.                                     else
  647.                                     {   
  648.                                         if(tt==k)////////与ce相同
  649.                                         {
  650.                                             ////ce 到 ce  的 概率
  651.                                             //mat.m_adValue[k][k]=mat.m_adValue[k][k]-cp.Pz;
  652.                                             p=head[k];
  653.                                             p->x=p->x-cp.Pz;
  654.                                         }
  655.                                         else
  656.                                         {
  657.                                             //mat.m_adValue[k][tt]=-cp.Pz;
  658.                                             p=head[k];                                            
  659.                                             while(p->next!=NULL)
  660.                                             {
  661.                                                 p=p->next;
  662.                                             }
  663.                                             p->next=new xlink();
  664.                                             p=p->next;
  665.                                             p->x=-cp.Pz;
  666.                                             p->num=tt;
  667.                                         }
  668.                                     }
  669.                                 }                           
  670.                             }
  671.                         }
  672.                         /////////上面获得系数                    
  673.                         for(long yu=0;yu<thrn;yu++)
  674.                         {
  675.                             cst.m_adValue[yu][0]=cs[yu];
  676.                         }
  677.                         if(thrn<1000)
  678.                         {
  679.                             mat.MatInv();                        
  680.                             cst=mat*cst;
  681.                             ///////////////////////                        
  682.                            
  683.                         }
  684.                         else
  685.                         {
  686.                             err=1.0;
  687.                             while(err>0.01)
  688.                             {
  689.                                 float tot=0.0;
  690.                                 for(long tr=0;tr<thrn;tr++)
  691.                                 {
  692.                                     if(tr%500==0)
  693.                                         cout<<tr<<"    "<<thrn<<endl;
  694.                                     tot=0.0;
  695.                                     p=head[tr];
  696.                                     p=p->next;/////第二个胞了,非tr
  697.                                     while(p!=NULL)
  698.                                     {
  699.                                         tot=tot+(p->x)*xr[p->num];
  700.                                         p=p->next;
  701.                                     }                                    
  702.                                     p=head[tr];                                    
  703.                                     float njt=p->x;
  704.                                     xr2[tr]=(cs[tr]-tot)/njt;
  705.                                 }
  706.                                 err=0.0;
  707.                                 for(long jr=0;jr<thrn;jr++)
  708.                                 {
  709.                                     err=err+sqrt((xr2[jr]-xr[jr])*(xr2[jr]-xr[jr]));
  710.                                 }
  711.                                 for(jr=0;jr<thrn;jr++)
  712.                                 {
  713.                                     xr[jr]=xr2[jr];
  714.                                 }
  715.                                 cout<<err<<"            "<<err<<endl;
  716.                             }
  717.                             ////////////结果放在cst
  718.                             for(long te=0;te<thrn;te++)
  719.                             {
  720.                                 cst.m_adValue[te][0]=xr[te];
  721.                             }
  722.                            
  723.                             for(long mt=0;mt<thrn;mt++)
  724.                             {
  725.                                 xlink *q1,*q2;
  726.                                 q1=head[mt]->next;                                
  727.                                 
  728.                                 while(q1!=NULL)
  729.                                 {                                    
  730.                                     q2=q1->next;
  731.                                     delete q1;
  732.                                     q1=q2;
  733.                                 }                           
  734.                             }


  735.                         }
  736.                         
  737.                         ////存入数据
  738.                         for(long tr=0;tr<thrn;tr++)
  739.                         {                           
  740.                             long tre=threm[tr];
  741.                             float pjr=cst.m_adValue[tr][0];
  742.                             fres[tre][i-1]=pjr;                           
  743.                            
  744.                         }
  745.                         /////////////vj////////////////                                                
  746.                         delete [] cs;
  747.                         delete [] xr;
  748.                         delete [] xr2;
  749.                         ////////结果,改变sstem  -1  以及ioder
  750.                         for(long u=0;u<thrn;u++)
  751.                         {
  752.                             long cme=threm[u];
  753.                             long det=i;
  754.                             fder[cme]=det;
  755.                             long tr=0;
  756.                             while(sstem[tr]!=cme&&tr<snu)
  757.                             {
  758.                                 tr++;
  759.                             }
  760.                             sstem[tr]=-1;
  761.                         }
  762.                         nhand=nhand+thrn;
  763.                         //for(long loop=0;loop<thrn;loop++)//////////////////清理
  764.                         delete [] head;
  765.                         //delete [] head;
  766.                     }
  767.                     //delete [] stt;
  768.                     //delete [] stn;
  769.                     delete [] sstem;
  770.                     ///////////////////////////////////////
  771.                     
  772.                     for(long u=0;u<snu;u++)
  773.                     {
  774.                         long cew;
  775.                         cew=stem[u];                                                                        
  776.                         long yu=i;
  777.                         fder[cew]=yu;                        
  778.                         /////////////////dtem 中 -1
  779.                         long tt=0;
  780.                         while(dtem[tt]!=cew&&tt<nu)
  781.                         {
  782.                             tt++;
  783.                         }
  784.                         dtem[tt]=-1;/////////////////////
  785.                     }                                       
  786.                     nt=nt+snu;
  787.                 }
  788.             }
  789.             cet=ftemp[sta];
  790.             sta++;            
  791.             ///////////////此批计算完//////////////////////
  792.         }        
  793.     }   
  794.     ////////////////////////

  795.     fstream ioresult;         ////// 保存结果 gr,Vj,pj
  796.     ioresult.open("result.dat",ios::in|ios::out|ios::binary);

  797.     ///long res=3*sizeof(float);/////////////////////////////////////////
  798.     for(long mt=0;mt<max;mt++)
  799.     {
  800.         for(long t=0;t<3;t++)
  801.         {
  802.             ioresult.write((char *)(&fres[mt][t]),sizeof(float));
  803.         }
  804.     }
  805.     delete [] fzpre;
  806.     delete [] flag;
  807.     delete [] fder;
  808.     delete [] fcpl;
  809.     delete [] fgcm;
  810.     delete [] ftemp;
  811. }
  812. ////////////排序算法
  813. void sort(long a[],long b[],long n)
  814. {
  815.     long i,j,small;
  816.     long temp1,temp2;
  817.     for(i=0;i<n-1;i++)
  818.     {
  819.         small=i;
  820.         for(j=i+1;j<n;j++)
  821.         {
  822.             if(a[j]<a[small])
  823.                 small=j;
  824.         }
  825.         if(small!=i)
  826.         {
  827.             temp1=a[i];
  828.             temp2=b[i];
  829.             a[i]=a[small];            
  830.             a[small]=temp1;
  831.             b[i]=b[small];
  832.             b[small]=temp2;
  833.         }
  834.     }
  835. }

  836. long pai(long m)
  837. {
  838.     long n=g_attr;//
  839.     long nt=1;
  840.     long mt=1;
  841.     if(m>n)
  842.     {
  843.         return -1;
  844.     }
  845.     else
  846.     {
  847.         for(long i=1;i<=m;i++)
  848.         {
  849.             nt=nt*n;
  850.             n=n-1;
  851.             
  852.         }
  853.         if(m==0)
  854.         {
  855.             mt=1;
  856.         }
  857.         else
  858.         {
  859.             for(i=1;i<=m;i++)
  860.             {
  861.                 mt=mt*i;
  862.             }
  863.         }
  864.         long re;
  865.         re=(long) (nt/mt);
  866.         return re;
  867.     }
  868. }

  869. void fenzu()
  870. {
  871.     ///// 根据  iodm 中的信息,改变ioresult 中的gr值   
  872.     fstream iodm;    ////// 保存 cell 的, domicle 的数量,做分组用  
  873.     iodm.open("todomiledata.dat",ios::in|ios::out|ios::binary);
  874.     dom dm;
  875.     fstream ioresult;         ////// 保存结果 gr,Vj,pj
  876.     ioresult.open("result.dat",ios::in|ios::out|ios::binary);
  877.     long res=sizeof(long)+2*g_attr*sizeof(float);//////////////////////////
  878.     long pai(long m);
  879.     for(long i=0;i<max;i++)
  880.     {
  881.         if(i%5000==0)
  882.             cout<<i<<endl;
  883.         long dn;
  884.         long total=0;
  885.         iodm.seekg(i*sizeof(dom));
  886.         iodm.read((char *)(&dm),sizeof(dom));
  887.         dn=dm.n;
  888.         total=dm.total;
  889.         long gpai=0;
  890.         if(dn==0)
  891.         {
  892.             ioresult.seekp(i*res);
  893.             ioresult.write((char *)(&total),sizeof(long));
  894.         }
  895.                
  896.     }
  897.     ///////////////////////////////////
  898.    
  899. }


  900. void outputdata()
  901. {
  902.     fstream iopersistent;   
  903.     iopersistent.open("persistent.dat",ios::in|ios::out|ios::binary);
  904.     //g_attr=2;
  905.    
  906.     long tr=0;
  907.     ofstream df("persistent.txt");
  908.     iopersistent.seekg(0,ios::beg);
  909.     long cur,pend;
  910.     cur=iopersistent.tellg();
  911.     iopersistent.seekg(0,ios::end);
  912.     pend=iopersistent.tellg();
  913.     iopersistent.seekg(cur);
  914.     while(pend!=iopersistent.tellg())
  915.     {
  916.         long cee,grg;
  917.         iopersistent.read((char *)(&cee),sizeof(long));
  918.         iopersistent.read((char *)(&grg),sizeof(long));
  919.         float x1,x2,x3,x4;   
  920.         long z1,z2,z3,z4,m;
  921.         m=cee;
  922.         z4=(long)((m-1)/(step1*step2*step3))+1;
  923.         m=m-(z4-1)*step1*step2*step3;
  924.         z3=(long)((m-1)/(step1*step2))+1;
  925.         m=m-(z3-1)*step1*step2;
  926.         z2=(long) ((m-1)/step1)+1;
  927.         z1=m-(z2-1)*step1;   
  928.         x1=z1*h1-h1+x01;
  929.         x2=z2*h2-h2+x02;
  930.         x3=z3*h3-h3+x03;
  931.         x4=z4*h4-h4+x04;
  932.         
  933.         df<<x1<<"           "<<x2<<"           "<<x3<<"           "<<x4<<endl;
  934.         
  935.         
  936.     }
  937.    
  938. }
复制代码
发表于 2013-11-21 17:05 | 显示全部楼层

师姐  为什么加不了你呢?
发表于 2018-4-25 11:23 | 显示全部楼层
cinly_zhu 发表于 2008-3-19 23:06
**** 作者被禁止或删除 内容自动屏蔽 ****

你好,我是刚开始研究胞映射相关课题的学生,求学长指点。能加下qq嘛?
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

快速回复 返回顶部 返回列表