dynamic.txt

mfold

#define minloop 3 //The minimum substructure with a basepair possible


//This is the dynamic algorithm of Zuker:
void dynamic(structure* ct, datatable* data, int cntrl6, int cntrl8, int cntrl9,
	     TProgressDialog* update, char* save)
{

  //cout << "2!!!"<<flush;

  //int **v, **w,
  int vmin, d, ip, jp, ii, jj;
  int e[6], i;
  int k, j, l, m, n, o, p;
  int inc[6][6]={{0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 0}, {0, 0, 0, 1, 0, 0}, {0, 0, 1, 0, 1, 0},
		 {0, 1, 0, 1, 0, 0}, {0, 0, 0, 0, 0, 0}};
  bool *lfce; //[maxbases+1][maxbases+1];
  int before, after, key;
  int **fce;
  int **work = NULL, **work2 = NULL;
  int *w5, *w3, **wmb, **wmb2 = NULL;
  register int number, jmt;
  output convert;

  //array *v, *w;
  //inc is an array that saves time by showing which bases can pair before
  //	erg is called

  //number is the number of bases being folded
  //v[i][j] is the best energy for subsequence i to j when i and j are paired
  //w[i][j] is the best energy for subsequence i to j

  number = (ct->numofbases); //place the number of bases in a registered integer
  key = 0;
  //w = new array(number);
  //v = new array(number);


  //////////////////
  //cout << "preallocate\n";



  arrayclass w(number);
  arrayclass v(number);

  //add a second array for intermolecular folding:
  arrayclass *w2 = NULL;
  if (ct->intermolecular) {
    w2 = new arrayclass(number);
    work2 = new int *[2*number+3];
    wmb2 = new int *[2*number+3];
    for (i=0; i<2*number+3; i++) {
      work2[i] = new int [3];
      wmb2[i] = new int [3];
      for (j=0; j<=2; j++) wmb2[i][j] = infinity;
    }

  }


  fce = new int *[2*number+1];
  for (i=0; i<=2*number; i++)
    fce[i] = new int [2*number + 1];

  for (i=1; i<=2*number; i++) {
    for (j=1; j<=2*number; j++) {
      fce[i][j] = 0;
    }
  }


  lfce = new bool [2*number+1];

  for (i=0; i<=2*number; i++) lfce[i] = false;




  work = new int *[2*number+3];
  wmb = new int *[2*number+3];
  for (i=0; i<2*number+3; i++) {
    work[i] = new int [3];
    wmb[i] = new int [3];
    for (j=0; j<=2; j++) wmb[i][j] = infinity;
  }




  w5 = new int [number+1];
  w3 = new int [number+2];

  for (i=0; i<=number; i++) {
    w5[i] = 0;
    w3[i] = 0;
  }
  w3[number+1] = 0;



  force(ct, &v, fce, lfce);


  //This is the fill routine:


  vmin=infinity;
  for (j=1; j<=(2*(number)-1); j++) {

    if (((j%10)==0)&&update) update->update((100*j)/(2*ct->numofbases));

    for (i=min(j, number); i>=max(1, j-(number)+1); i--) {

      //debug
      //if(i==25&&j==26) {
      // 	v.f(0, 0) = 0;
      //}


      wmb[i][j%3] = infinity;

      if (ct->templated) {
   	if (i>ct->numofbases) ii = i - ct->numofbases;
	else ii = i;
	if (j>ct->numofbases) jj = j - ct->numofbases;
	else jj = j;
	if (jj<ii) {
	  p = jj;
	  jj = ii;
	  ii = p;
	}
   	if (!ct->tem[jj][ii]) goto sub2;
      }

      //Compute v[i][j], the minimum energy of the substructure from i to j,
      //inclusive, where i and j are base paired

      //if fce[i][j] = 1 or 3, these need to be set to infinity
      //	because 1 indicates one of the bases should be single-stranded
      //		and 3 indicates one of the bases needs to be paired to a different base

      if (v.f(i, j)==1) {
	v.f(i, j) = infinity + 50;
	goto sub2;
      }
      if (v.f(i, j)==3) {
   	v.f(i, j)= infinity+50;
	//w[i][j]= infinity;
	goto sub3;
      }
      if (v.f(i, j)==2) {//forcing a pair between i and j
   	key = 2;
      }

      //v[i][j] = infinity; //added 8/18/97


      if (j<=(number)) {
	if ((j-i)<=minloop) goto sub3;
      }
      else {
	if (i==(number)||j==((number)+1)) goto sub1;
      }

      if (inc[ct->numseq[i]][ct->numseq[j]]==0) goto sub2;
      //force u's into gu pairs
      for (ip=0; ip<ct->ngu; ip++) {
   	if (ct->gu[ip]==i) {
	  if (ct->numseq[j]!=3) {
	    v.f(i, j) = infinity;
	    goto sub2;
	  }
	}
	//else if ((ct->gu[ip]+number)==i) {
	// 	if (ct->numseq[j]!=3) {
	//   	v.f(i, j) = infinity;
	//      goto sub2;
	//   }
	//}
	else if (ct->gu[ip]==j) {
	  if (ct->numseq[i]!=3) {
	    v.f(i, j) = infinity;
	    goto sub2;
	  }
	}
	else if ((ct->gu[ip]+number)==j) {
	  if (ct->numseq[i]!=3) {
	    v.f(i, j) = infinity;
            goto sub2;
	  }
	}

      }

      //now check to make sure that this isn't an isolated pair:
      //	(consider a pair separated by a bulge as not! stacked)

      //before = 0 if a stacked pair cannot form 5' to i
      if (i>1&&j<(2*number)) {
	before = inc[ct->numseq[i-1]][ct->numseq[j+1]];
	//		What follows is to allow a bulge
	//		if (before == 0) {
	//			if (i>2) {if (incp[ct->numseq[i-2]][ct->numseq[j+1]]==1) before = 1; }
	//			else if (j-1<(2*number)) {if (incp[ct->numseq[i-1]][ct->numseq[j+2]]==1) before = 1; }
	//		}
      }
      else before = 0;

      //after = 0 if a stacked pair cannot form 3' to i
      if ((j-i)>5) {
	after = inc[ct->numseq[i+1]][ct->numseq[j-1]];
	//		What follows is to allow a bulge
	//		if (after == 0) {
	//			if ((incp[ct->numseq[i+2]][ct->numseq[j-1]]==1)&&((j-i)>6)) after = 1;
	//			else if ((incp[ct->numseq[i+1]][ct->numseq[j-2]]==1)&&((j-i)>6)) after = 1;
	//		}
      }
      else after = 0;

      //if there are no stackable pairs to i.j then don't allow a pair i, j
      if ((before==0)&&(after==0)) {
	goto sub2;
      }


      //Perhaps i and j close a hairpin:
      v.f(i, j)=min(v.f(i, j), erg3(i, j, ct, data, fce[i][j]));

      if ((j-i-1)>=(minloop+2)||j>(number))
      	//Perhaps i, j stacks over i+1, j-1
	v.f(i, j)=min(v.f(i, j), (erg1(i, j, i+1, j-1, ct, data)+v.f(i+1, j-1)));
      //Perhaps i, j closes an interior or bulge loop, search for the best
      //possibility
      if (((j-i-1)>=(minloop+3))||(j>(number))) {
      	for (d=(j-i-3); d>=1; d--) {
	  for (ip=(i+1); ip<=(j-1-d); ip++) {
	    jp = d+ip;
	    if ((j-i-2-d)>(data->eparam[7])) goto sub1;
	    if (abs(ip-i+jp-j)<=(data->eparam[8])) {
	      if (ip>(number)) {
		//if (jp<=number) {

		//	v.f(i, j)=min(v.f(i, j), (erg2(i, j, ip, jp, ct, data, fce[i][ip-number],
		//		fce[jp][j])+
		//		v.f(ip-(number), jp)));

		//}
		//else {
		v.f(i, j)=min(v.f(i, j), (erg2(i, j, ip, jp, ct, data, fce[i][ip-number],
					    fce[jp-number][j-number])+
				       v.f(ip-(number), jp-(number))));
		//}
	      }
	      else {
		if (jp<=number) {

		  v.f(i, j)=min(v.f(i, j), (erg2(i, j, ip, jp, ct, data, fce[i][ip],
					      fce[jp][j])+
					 v.f(ip, jp)));
		}
		else {

		  v.f(i, j)=min(v.f(i, j), (erg2(i, j, ip, jp, ct, data, fce[i][ip],
					      fce[jp-number][j-number])+
					 v.f(ip, jp)));

		}
	      }
	    }
	  }
	}
      }

      //Perhaps i, j closes a multi-loop, search for the best possibility

    sub1:

      if (((j-i-1)>=(2*minloop+4))||(j>(number))) {
	for (ii=1; ii<=4; ii++) e[ii]=infinity;

	if (j>number) {
	  e[1] = w3[i+1] + w5[j-number-1] + penalty(i, j, ct, data);
	  //if (e[1]==-1204) {
	  //    	v.f(0, 0) = 0;
	  //}


	  if (i!=number) e[2] = erg4(i, j, i+1, 1, ct, data, lfce[i+1]) + penalty(i, j, ct, data)+w3[i+2] + w5[j-number-1];
	  if (j!=(number+1)) e[3] = erg4(i, j, j-1, 2, ct, data, lfce[j-1]) +penalty(i, j, ct, data)+ w3[i+1] + w5[j-number-2];
	  if ((i!=number)&&(j!=(number+1))) {
	    e[4] = erg4(i, j, i+1, 1, ct, data, lfce[i+1]) +
	      erg4(i, j, j-1, 2, ct, data, lfce[j-1]) + w3[i+2] + w5[j-number-2]
	      +penalty(i, j, ct, data);

	  }
	}




	if ((j-i)>(2*minloop+4)) {
	  //no dangling ends on i-j pair:
	  e[1] = min(e[1], wmb[i+1][(j-1)%3]+data->eparam[5]+data->eparam[10]
		     + penalty(i, j, ct, data));


	  //i+1 dangles on i-j pair:
	  if (i!=number) e[2] = min(e[2], erg4(i, j, i+1, 1, ct, data, lfce[i+1]) + penalty(i, j, ct, data) +
				    wmb[i+2][(j-1)%3] + data->eparam[5] + data->eparam[6] + data->eparam[10]);
	  //j-1 dangles
	  if (j!=(number+1)) e[3] = min(e[3], erg4(i, j, j-1, 2, ct, data, lfce[j-1]) + penalty(i, j, ct, data) +
					wmb[i+1][(j-2)%3] + data->eparam[5] + data->eparam[6] + data->eparam[10]);
	  //both i+1 and j-1 dangle
	  if ((i!=number)&&(j!=(number+1))) {
	    e[4] = min(e[4],
		       erg4(i, j, i+1, 1, ct, data, lfce[i+1]) +  erg4(i, j, j-1, 2, ct, data, lfce[j-1]) +
		       wmb[i+2][(j-2)%3] + data->eparam[5] + 2*data->eparam[6] + data->eparam[10]
		       +penalty(i, j, ct, data));

	  }

	  if (ct->intermolecular) {

	    //intermolecular, so consider wmb2,
	    //don't add the multiloop penalties because this is a exterior loop

	    e[1] = min(e[1], wmb2[i+1][(j-1)%3]
		       + penalty(i, j, ct, data));


	    //i+1 dangles on i-j pair:
	    if (i!=number) e[2] = min(e[2], erg4(i, j, i+1, 1, ct, data, lfce[i+1]) + penalty(i, j, ct, data) +
				      wmb2[i+2][(j-1)%3]);
	    //j-1 dangles
	    if (j!=(number+1)) e[3] = min(e[3], erg4(i, j, j-1, 2, ct, data, lfce[j-1]) + penalty(i, j, ct, data) +
					  wmb2[i+1][(j-2)%3] );
	    //both i+1 and j-1 dangle
	    if ((i!=number)&&(j!=(number+1))) {
	      e[4] = min(e[4],
			 erg4(i, j, i+1, 1, ct, data, lfce[i+1]) +  erg4(i, j, j-1, 2, ct, data, lfce[j-1]) +
			 wmb2[i+2][(j-2)%3]
			 +penalty(i, j, ct, data));

	    }
	  }


	}



	v.f(i, j)=min(v.f(i, j), e[1]);
	v.f(i, j)=min(v.f(i, j), e[2]);
	v.f(i, j)=min(v.f(i, j), e[3]);
	v.f(i, j)=min(v.f(i, j), e[4]);
      }

      //Compute w[i][j]: best energy between i and j where i, j does not have
      //	to be a base pair
      //(an exterior loop when it contains n and 1 (ie:n+1)   )
    sub2:


      //w[i][j]=infinity;

      if (key==2) {//force a pair between i and j
	w.f(i, j) = v.f(i, j);
	key = 0;
	goto sub3;
      }


      for (ii=1; ii<=5; ii++) e[ii] = infinity;


      if (i!=number) {
      	//calculate the energy of i stacked onto the pair of i+1, j

	e[1] = v.f(i+1, j) + data->eparam[10] + data->eparam[6] +
	  erg4(j, i+1, i, 2, ct, data, lfce[i])+penalty(i+1, j, ct, data);


	if (!lfce[i]) {
	  if (fce[i][i]!=2)
	    e[4] = w.f(i+1, j) + data->eparam[6];
	  //this is for when i represents the center of an intermolecular linker:
	  else e[4] = w.f(i+1, j) + data->eparam[6] + infinity;
	}
      }
      if (j!=((number)+1)) {
      	//calculate the energy of j stacked onto the pair of i, j-1
	if (j!=1) {
	  e[2] = v.f(i, j-1) + data->eparam[10] + data->eparam[6] +
	    erg4(j-1, i, j, 1, ct, data, lfce[j])+penalty(i, j-1, ct, data);

	  if (!lfce[j]) {
	    if (fce[j][j]!=2) {
	      e[5] = w.f(i, j-1) + data->eparam[6];
	    }
	    else e[5] = w.f(i, j-1) + data->eparam[6] + infinity;

	  }
	}
      }
      if ((i!=(number))&&(j!=((number)+1))) {
      	//calculate i and j stacked onto the pair of i+1, j-1
	if (j!=1) {
	  e[3] = v.f(i+1, j-1) + data->eparam[10] + 2*(data->eparam[6]) +
	    erg4(j-1, i+1, i, 2, ct, data, lfce[i]) +
	    erg4(j-1, i+1, j, 1, ct, data, lfce[j])
	    +penalty(j-1, i+1, ct, data);
	}
      }
      w.f(i, j) = min(((data->eparam[10])+v.f(i, j)+penalty(j, i, ct, data)), e[1]);
      w.f(i, j) = min(w.f(i, j), e[2]);
      w.f(i, j) = min(w.f(i, j), e[3]);
      w.f(i, j) = min(w.f(i, j), e[4]);
      w.f(i, j) = min(w.f(i, j), e[5]);

      if (ct->intermolecular) {

      	wmb2[i][j%3] = infinity;
      	//keep track of w2:
	for (ii=1; ii<=5; ii++) e[ii] = 2*infinity;


	if (i!=number) {
	  //calculate the energy of i stacked onto the pair of i+1, j

	  e[1] = v.f(i+1, j) +
	    erg4(j, i+1, i, 2, ct, data, lfce[i])+penalty(i+1, j, ct, data)+infinity;


	  if (!lfce[i]) {
	    if (fce[i][i]!=2)
	      e[4] = w2->f(i+1, j);
	    //this is for when i represents the center of an intermolecular linker:
            else e[4] = w2->f(i+1, j) - infinity + data->init;
	  }
	}
      	if (j!=((number)+1)) {
	  //calculate the energy of j stacked onto the pair of i, j-1
	  if (j!=1) {
	    e[2] = v.f(i, j-1)  + infinity +
	      erg4(j-1, i, j, 1, ct, data, lfce[j])+penalty(i, j-1, ct, data);

	    if (!lfce[j]) {
	      if (fce[j][j]!=2) {
		e[5] = w2->f(i, j-1);
	      }
	      else e[5] = w2->f(i, j-1) - infinity + data->init;

	    }
	  }
      	}
      	if ((i!=(number))&&(j!=((number)+1))) {
	  //calculate i and j stacked onto the pair of i+1, j-1
	  if (j!=1) {
	    e[3] = v.f(i+1, j-1)  + infinity +
	      erg4(j-1, i+1, i, 2, ct, data, lfce[i]) +
	      erg4(j-1, i+1, j, 1, ct, data, lfce[j])
	      +penalty(j-1, i+1, ct, data);
	  }
      	}
      	w2->f(i, j) = min((infinity+v.f(i, j)+penalty(j, i, ct, data)), e[1]);
      	w2->f(i, j) = min(w2->f(i, j), e[2]);
      	w2->f(i, j) = min(w2->f(i, j), e[3]);
      	w2->f(i, j) = min(w2->f(i, j), e[4]);
      	w2->f(i, j) = min(w2->f(i, j), e[5]);




      }


      /*if (((j-i-1)>(2*minloop+2))||(j>(number))) {
      // search for an open bifuraction:
      for (k=i; k<=j-1; k++) {
      if (k==(number)) w.f(i, j)=min(w.f(i, j),
      w3[i]+w5[j-(number)]);
      else w.f(i, j) = min(w.f(i, j),
      w.f(i, k)+work[k+1][j%3]);
      }
      }  */
      ////fill wmb:

      if (((j-i-1)>(2*minloop+2))||j>number) {