algorithm.cpp

mfold

	  ji=jret;
	}
      	if (count==2) {
	  ii=jret;
	  ji=iret+(number);
	}

	trace(ct, data, ii, ji, v, w, w2, rep-1, lfce, fce, w3, w5);







	if (count==2) {





	  ct->energy[rep-1] = v->f(iret, jret) + v->f(jret, iret+(number));
	  //count the number of new base pairs not within window of existing
	  //base pairs
	  numbp = 0;
	  for (k=1; k<=number; k++) {
	    if (k<(ct->basepr[rep-1][k])) {
	      if (!(mark[k][ct->basepr[rep-1][k]])) numbp++;
	    }
	  }
	  for (k=1; k<=(number); k++) {
	    if (k<ct->basepr[rep-1][k]) {
	      //Mark "traced back" base pairs and also base pairs
	      //	which are within a window of cntrl9
	      mark[k][ct->basepr[rep-1][k]]=true;
	      if (cntrl9>0) {
		for (k1=-(cntrl9); k1<=cntrl9; k1++) {
		  for (k2=-cntrl9; k2<=cntrl9; k2++) {
		    if (((k+k1)>0)&&((k+k1)<(ct->basepr[rep-1][k]))&&
			((ct->basepr[rep-1][k]+k2)<=(number)))
		      mark[k+k1][(ct->basepr[rep-1][k])+k2]=true;
		  }
		}
	      }
	    }
	  }
	  if (numbp<=cntrl9) {
	    rep--;

	    goto sub900;
	  }
	  else {

	    //place the structure name (from ctlabel[1]) into each structure
	    strcpy(ct->ctlabel[rep-1], ct->ctlabel[1]);
            

            //debugging
            // 	out << (rep-1) << " picked pair: "<<iret<<" "<<jret<<"\n";

	  }
	  if (rep>cntrl6) flag=false;
	}
      }
    sub900:
      continue;
    }
  }

  (ct->numofstructures)=rep-1;


  //for debugging:
  //energydump (ct, v, 1, "dump.out");
  //energydump (ct, data, v, 1, "dump-au.out");








  de_allocate (mark, number+1);
  delete[] energy;
  delete[] heapi;
  delete[] heapj;

}





void force(structure *ct, arrayclass *v, int **fce, bool *lfce){
  int i;
  register int number;

  number = ct->numofbases;

  for (i=1; i<=ct->nnopair; i++) {
    forcesingle(ct->nopair[i], ct, v);
  }

  for (i=1; i<=ct->npair; i++) {
    forcepair(ct->pair[i][0], ct->pair[i][1], ct, v);
  }

  for (i=1; i<=ct->ndbl; i++) {
    forcedbl(ct->dbl[i], ct, fce, lfce);
  }


  //u's in gu pairs must be double stranded
  for (i=0; i<ct->ngu; i++) forcedbl(ct->gu[i], ct, fce, lfce);



  if (ct->intermolecular) {//this indicates an intermolecular folding
    //for (i=0; i<=3; i++) {
    // 	forcesingle(ct->inter[i])//don't allow the intermolecular indicators to pair
    //}
    for (i=0; i<3; i++) {
      forceinter(ct->inter[i], ct, fce);
    }


    //mark the fce array so that initiation is added to multibranch loops:
    //fce[ct->inter[0]][ct->inter[1]] = 2;
    //fce[ct->inter[1]][ct->inter[2]] = 2;
    //fce[ct->inter[0]][ct->inter[1]+number] = 2;
    //fce[ct->inter[1]][ct->inter[2]+number] = 2;

    fce[ct->inter[1]][ct->inter[1]] = 2;


  }

  //Double up the sequence
  for (i=1; i<=number; i++) {
    ct->numseq[(number)+i] = ct->numseq[i];
  }


}



#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]);