replik.c

模拟RNA的复制

     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      temp_seq[k]=0;
     }
     length=sq_length[xd5][yd5];
     sq_length[xd4][yd4]=sq_length[xd5][yd5];
     t_n[xd4][yd4]=t_n[xd5][yd5]; r_n[xd4][yd4]=r_n[xd5][yd5]; q_n[xd4][yd4]=q_n[xd5][yd5];
     for(k=1; k<=length; k++)
     {
      sq_seq[xd4][yd4][k]=sq_seq[xd5][yd5][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd4][yd4][k]=0;
     }
     length=temp_l;
     sq_length[xd5][yd5]=temp_l;
     t_n[xd5][yd5]=temp_t; r_n[xd5][yd5]=temp_r; q_n[xd5][yd5]=temp_q;
     for(k=1; k<=length; k++)
     {
      sq_seq[xd5][yd5][k]=temp_seq[k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd5][yd5][k]=0;
     }
     length=sq_length[xd6][yd6];
     temp_l=length;
     temp_t=t_n[xd6][yd6]; temp_r=r_n[xd6][yd6]; temp_q=q_n[xd6][yd6];
     for(k=1; k<=length; k++)
     {
      temp_seq[k]=sq_seq[xd6][yd6][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      temp_seq[k]=0;
     }
     length=sq_length[xd1][yd1];
     sq_length[xd6][yd6]=sq_length[xd1][yd1];
     t_n[xd6][yd6]=t_n[xd1][yd1]; r_n[xd6][yd6]=r_n[xd1][yd1]; q_n[xd6][yd6]=q_n[xd1][yd1];
     for(k=1; k<=length; k++)
     {
      sq_seq[xd6][yd6][k]=sq_seq[xd1][yd1][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd6][yd6][k]=0;
     }
     length=sq_length[xd3][yd3];
     sq_length[xd1][yd1]=sq_length[xd3][yd3];
     t_n[xd1][yd1]=t_n[xd3][yd3]; r_n[xd1][yd1]=r_n[xd3][yd3]; q_n[xd1][yd1]=q_n[xd3][yd3];
     for(k=1; k<=length; k++)
     {
      sq_seq[xd1][yd1][k]=sq_seq[xd3][yd3][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd1][yd1][k]=0;
     }
     length=temp_l;
     sq_length[xd3][yd3]=temp_l;
     t_n[xd3][yd3]=temp_t; r_n[xd3][yd3]=temp_r; q_n[xd3][yd3]=temp_q;
     for(k=1; k<=length; k++)
     {
      sq_seq[xd3][yd3][k]=temp_seq[k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd3][yd3][k]=0;
     }
  }
 
 void mc_step(void)
  {
   mc_stepi();

   for (j=1; j<=SIZE; j++)
   {
    xy_choose();
    if(randd()<DIFFRATE)
    {
     toff_marg();
    }

    length=sq_length[x][y];

    if(!length) {continue;} /* if empty site -> continue */

    if(randd()<DEATH)        /* if death occurs -> overwrite, continue */
    {
     for(k=1; k<=length; k++)
     {
      sq_seq[x][y][k]=0;
     }
     t_n[x][y]=0;
     r_n[x][y]=0;
     q_n[x][y]=0;
     sq_length[x][y]=0;
     n_denz=(n_denz*SIZE-length)/SIZE;
     continue;
    }

    if (length<CATCH) {continue;}  /* if too short to be a template -> continue */

    neighbourhood();

    if (!(life*empty)) {continue;} /* if there is no at least one      */
				                   /* empty and one filled site        */
				                   /* among the neighbours -> continue */
    rd=0.3333*life;
    if(randd()>rd) {continue;}

    s1=randl(life)+1;
    s1=lifes[s1];
    s2=randl(empty)+1;
    s2=empties[s2];

    if (s1==1) { xs1=left[x];  ys1=y;       }
    if (s1==2) { xs1=right[x]; ys1=y;       }
    if (s1==3) { xs1=x;        ys1=up[y];   } /* coordinates of replicase */
    if (s1==4) { xs1=x;        ys1=down[y]; }

    if (s2==1) { xs2=left[x];  ys2=y;       }
    if (s2==2) { xs2=right[x]; ys2=y;       }
    if (s2==3) { xs2=x;        ys2=up[y];   } /* coordinates of new chain    */
    if (s2==4) { xs2=x;        ys2=down[y]; }


    /* longest continuous sequences of monomers in replicase */

    max_B=r_n[xs1][ys1];
    max_C=q_n[xs1][ys1];


    r=fgv_replicase[max_B];

    q=fgv_fidelity[max_C];


    /* longest continuous sequences of monomers in template */

    max_A=t_n[x][y];


    t=fgv_template[max_A];

    length;

    rd=randd();

    n_factor=(NUKLPAR-n_denz)/NUKLPAR;
    if(rd<n_factor*(t*r/length))                       /* replication */
    {
     sq_length[xs2][ys2]=length;
     for(k=1; k<=length; k++)
     {
      rd=randd();
      if(rd<q)
      {
       sq_seq[xs2][ys2][k]=sq_seq[x][y][k];
      } else
      {
       nukl=sq_seq[x][y][k];
       rl=randl(3)+1;
       switch (nukl)
       {
       case 1:
        if(rl==1) sq_seq[xs2][ys2][k]=2;
        if(rl==2) sq_seq[xs2][ys2][k]=3;
        if(rl==3) sq_seq[xs2][ys2][k]=4;
        break;

       case 2:
        if(rl==1) sq_seq[xs2][ys2][k]=1;
        if(rl==2) sq_seq[xs2][ys2][k]=3;
        if(rl==3) sq_seq[xs2][ys2][k]=4;
       break;

       case 3:
        if(rl==1) sq_seq[xs2][ys2][k]=1;
        if(rl==2) sq_seq[xs2][ys2][k]=2;
        if(rl==3) sq_seq[xs2][ys2][k]=4;
       break;

       case 4:
        if(rl==1) sq_seq[xs2][ys2][k]=1;
        if(rl==2) sq_seq[xs2][ys2][k]=2;
        if(rl==3) sq_seq[xs2][ys2][k]=3;
       break;
       }
      }
     }
     rd=randd();
     if(rd<b)
     {
      length_o=length;
      length=randl(length-1)+1;
      sq_length[xs2][ys2]=length;
      if(randd()<0.5)
      {
       for(k=(length+1); k<=length_o; k++)
       {
        sq_seq[xs2][ys2][k]=0;
       }
      } else
      {
       for(k=1; k<=length; k++)
       {
        sq_seq[xs2][ys2][k]=sq_seq[xs2][ys2][length_o-length+k];
       }
       for(k=(length+1); k<=length_o; k++)
       {
        sq_seq[xs2][ys2][k]=0;
       }
      }
     }
     if((rd<(a+b)) && (rd>b))
     {
      length++;
      if(randd()<0.5)
      {
      sq_length[xs2][ys2]++;
      sq_seq[xs2][ys2][length]=randl(4)+1;
      } else
      {
      for(k=length; k>=2; k--)
      sq_seq[xs2][ys2][k]=sq_seq[xs2][ys2][k-1];
      sq_seq[xs2][ys2][1]=randl(4)+1;
      sq_length[xs2][ys2]++;
      }
     }
     n_denz=(n_denz*SIZE+length)/SIZE;
     trq_update();
    }
   }
  }

 void save(void)
 {
  tf2=fopen("pos.dat","w+t");
  fprintf(tf2,"%d",position);
  fflush(tf2);
  fclose(tf2);
  position++;
 }

 void result1(void)
 {
  double elokszama=0;
  long int lepesszam;
  double a_lanchossz=0, S_lanchossz=0, max_lanchossz=0, min_lanchossz=1000;
  double A_atlag=0, B_atlag=0, C_atlag=0, D_atlag=0,
         As_atlag=0, Bs_atlag=0, Cs_atlag=0, Ds_atlag=0,
         t_atlag=0, r_atlag=0, q_atlag=0, t_r_l_atlag; 
  double A_szoras=0, B_szoras=0, C_szoras=0, D_szoras=0,
         As_szoras=0, Bs_szoras=0, Cs_szoras=0, Ds_szoras=0,
         t_szoras=0, r_szoras=0, q_szoras=0;
  double A_min=1000, B_min=1000, C_min=1000, D_min=1000,
         As_min=1000, Bs_min=1000, Cs_min=1000, Ds_min=1000,
         t_min=1000, r_min=1000, q_min=1000;  
  double A_max=0, B_max=0, C_max=0, D_max=0,
         As_max=0, Bs_max=0, Cs_max=0, Ds_max=0,
         t_max=0, r_max=0, q_max=0, t_r_l_max=0;
  double denzitas=0;
  long int dat, len, max, maxA, maxB, maxC, maxD, minA, minB, minC, minD;
  long int hossz;
  double ar;
  FILE * tf;
  FILE * tfd;

  lepesszam=i;

  for(j=1; j<=SIDE; j++)
  {
   for(k=1; k<=SIDE; k++)
   {
    hossz=sq_length[j][k];
    if(hossz) 
    { 
     denzitas+=1; elokszama+=1; 
     a_lanchossz+=hossz;
     maxA=0, maxB=0, maxC=0, maxD=0; 
     for(l=1; l<=hossz; l++)
     {
      if(sq_seq[j][k][l]==1) {A_atlag+=1; maxA++;}
      if(sq_seq[j][k][l]==2) {B_atlag+=1; maxB++;}
      if(sq_seq[j][k][l]==3) {C_atlag+=1; maxC++;}
      if(sq_seq[j][k][l]==4) {D_atlag+=1; maxD++;}
     }
     max=t_n[j][k];
     As_atlag+=max;
     if(hossz>=CATCH) t_atlag+=fgv_template[max];
     max=r_n[j][k];
     Bs_atlag+=max;
     if(hossz>=CATCH) r_atlag+=fgv_replicase[max];
     max=q_n[j][k];
     Cs_atlag+=max;
     q_atlag+=fgv_fidelity[max];

     len=0; max=0;
     for(l=1; l<=hossz; l++)
     {
      dat=sq_seq[j][k][l];
      if(dat==4)
      {
       len++;
      }
      if(dat!=4)
      {
       if(len>max) max=len;
       len=0;
      }
     }
     Ds_atlag+=max;
    }
   }
  }

  a_lanchossz=a_lanchossz/elokszama;
  A_atlag=A_atlag/elokszama;
  B_atlag=B_atlag/elokszama;
  C_atlag=C_atlag/elokszama;
  D_atlag=D_atlag/elokszama;

  As_atlag=As_atlag/elokszama;
  Bs_atlag=Bs_atlag/elokszama;
  Cs_atlag=Cs_atlag/elokszama;
  Ds_atlag=Ds_atlag/elokszama;

  t_atlag=t_atlag/elokszama;
  r_atlag=r_atlag/elokszama;
  q_atlag=q_atlag/elokszama;

  denzitas=denzitas/SIZE;


  tf=fopen("result1.dat","a+t"); 

  fprintf(tf,"%10d", lepesszam); 

  fprintf(tf,"%8.2f",a_lanchossz); 
  fprintf(tf,"%6.2f",A_atlag); 
  fprintf(tf,"%6.2f",B_atlag); 
  fprintf(tf,"%6.2f",C_atlag); 
  fprintf(tf,"%6.2f",D_atlag); 
  fprintf(tf,"%6.2f",As_atlag); 
  fprintf(tf,"%6.2f",Bs_atlag); 
  fprintf(tf,"%6.2f",Cs_atlag); 
  fprintf(tf,"%6.2f",Ds_atlag); 
  fprintf(tf,"%8.4f",t_atlag); 
  fprintf(tf,"%8.4f",r_atlag); 
  fprintf(tf,"%8.4f",q_atlag); 
  fprintf(tf,"%7.3f\n",denzitas); 

  fflush(tf);
  fclose(tf);
 }

 void result2(void)
 {
  double elokszama=0;
  long int lepesszam;
  long int dat, len, max;
  restomb1 arany;
  restomb2 composition;
  long int hossz;
  double ar;

  FILE * tf;
  FILE * tfd;

  lepesszam=i;

  for(j=0; j<=MAXL_N; j++) arany[j]=0;
  for(j=0; j<=MAXL_N; j++)
  { 
   for(k=0; k<=8; k++) composition[j][k]=0;
  }

  for(j=1; j<=SIDE; j++)
  {
   for(k=1; k<=SIDE; k++)
   {
    hossz=sq_length[j][k];
    if(hossz) 
    {
     arany[hossz]+=1;

     for(l=1; l<=hossz; l++)
     {
      if(sq_seq[j][k][l]==1) composition[hossz][1]+=1;
      if(sq_seq[j][k][l]==2) composition[hossz][2]+=1;
      if(sq_seq[j][k][l]==3) composition[hossz][3]+=1;
      if(sq_seq[j][k][l]==4) composition[hossz][4]+=1;
     }
     composition[hossz][5]+=t_n[j][k];
     composition[hossz][6]+=r_n[j][k];
     composition[hossz][7]+=q_n[j][k];
     len=0; max=0;
     for(l=1; l<=hossz; l++)
     {
      dat=sq_seq[j][k][l];
      if(dat==4)
      {
       len++;
      }
      if(dat!=4)
      {
       if(len>max) max=len;
       len=0;
      }
     }
     if(hossz==1) max=len;
     composition[hossz][8]+=max;
    }
    else
    {
     arany[0]+=1;
    }
   }
  }
  for(j=0; j<=MAXL_N; j++)
  {
   for(k=1; k<=8; k++)
   {
    if (arany[j]) composition[j][k]=composition[j][k]/arany[j];
   }
  }
  for(j=0; j<=MAXL_N; j++)
  {
   arany[j]=arany[j]/SIZE;
  }

  tf=fopen("result2.dat","a+t"); 

  for(j=1; j<=(MAXL_N-1); j++)
  {
   fprintf(tf,"%8.4f", composition[j][1]);
  }
  fprintf(tf,"%8.4f\n", composition[MAXL_N][1]);
  for(j=1; j<=(MAXL_N-1); j++)
  {
   fprintf(tf,"%8.4f", composition[j][2]);
  }
  fprintf(tf,"%8.4f\n", composition[MAXL_N][2]);
  for(j=1; j<=(MAXL_N-1); j++)
  {
   fprintf(tf,"%8.4f", composition[j][3]);
  }
  fprintf(tf,"%8.4f\n", composition[MAXL_N][3]);
  for(j=1; j<=(MAXL_N-1); j++)
  {
   fprintf(tf,"%8.4f", composition[j][4]);
  }
  fprintf(tf,"%8.4f\n", composition[MAXL_N][4]);

  for(j=1; j<=(MAXL_N-1); j++)
  {
   fprintf(tf,"%8.4f", composition[j][5]);
  }
  fprintf(tf,"%8.4f\n", composition[MAXL_N][5]);
  for(j=1; j<=(MAXL_N-1); j++)
  {
   fprintf(tf,"%8.4f", composition[j][6]);
  }
  fprintf(tf,"%8.4f\n", composition[MAXL_N][6]);
  for(j=1; j<=(MAXL_N-1); j++)
  {
   fprintf(tf,"%8.4f", composition[j][7]);
  }
  fprintf(tf,"%8.4f\n", composition[MAXL_N][7]);
  for(j=1; j<=(MAXL_N-1); j++)
  {
   fprintf(tf,"%8.4f", composition[j][8]);
  }
  fprintf(tf,"%8.4f\n", composition[MAXL_N][8]);

  fprintf(tf,"%8d", lepesszam); 

  for(j=0; j<=(MAXL_N-1); j++)
  {
   fprintf(tf,"%9.5f", arany[j]);
  }
  fprintf(tf,"%9.5f\n", arany[MAXL_N]);

  fflush(tf);
  fclose(tf);
 }

 void result3(void)
 {
  tf=fopen("result3.dat","a+t"); 
  for(j=1; j<=SIDE; j++)
  {
   for(k=1; k<=SIDE; k++)
   {
    for(l=1; l<=(sq_length[j][k]); l++)
    {
     fprintf(tf," %1d", sq_seq[j][k][l]);
    }
    fprintf(tf," 0\n");
   }
  }
  fflush(tf);
  fclose(tf);
 }
 
 main()
 {
  inic();                               /* initialisation */
  for(i=0; i<=STEPNUMBER; i++)          /* Monte-Carlo cycle */
  {
   if((i%SAVEL)==0) {save();}
   if((i%RESULTL1)==0) {result1();}
   if((i%RESULTL2)==0) {result2();}
   if((i%RESULTL3)==0) {result3();}
   mc_step();
  }
 }