algorithm.cpp

mfold

  }
  else if (i>(number)) {
    ct->basepr[rep][i-(number)] = j - (number);
    ct->basepr[rep][j-(number)] = i - (number);
    i = i - (number);
    j = j - (number);
  }
  else {
    ct->basepr[rep][j-(number)] = i;
    ct->basepr[rep][i] = j - (number);
  }
  open = 0;
  //does i, j stack over i+1, j-1
  if ((i!=(number))&&(j!=((number)+1))) {
    if (en==(erg1(i, j, i+1, j-1, ct, data) + v->f(i+1, j-1))) {
      i++;
      j--;
      en = v->f(i, j);

      //output the stack value to dump file:\
	//if (filename!=0) {
	// 	out << "pair stack of "<<i<" and "<<j<<" onto "<<(i-1)<<
	//   	" and "<<(j+1)<<" dG = "<<(v->f(i-1, j+1)-v->f(i, j))
	//      <<"\n";

	//}

	goto sub300;
    }
  }
  //does i, j close a hairpin loop?
  if (en==erg3(i, j, ct, data, fce[i][j])) {

    //output the hairpin loop energy to dump file:
    //if (filename!=0) {
    //	out << "hairpin loop between nucleotides "<<i<<" and "<<j<<
    //		" dG = "<<erg3(i, j, ct, data, fce[i][j])<<"\n";
    //}
    goto sub100;
  }


  //ep = e corrected for forced base pairs

  ep = en;


  if ((i+2)>(j-3)) {
    //tidy up loose ends
    if ((ep==0)||(ep==penalty(i, j, ct, data))||((i!=(number))&&
					      (ep==erg4(i, j, i+1, 1, ct, data, lfce[i+1])+penalty(i, j, ct, data)))) goto sub100;
    else if ((j!=((number)+1))&&
	     (ep==(erg4(i, j, j-1, 2, ct, data, lfce[j-1])+penalty(i, j, ct, data))))
      goto sub100;
    else if ((i!=(number))&&(j!=((number)+1))&&
	     (ep==(erg4(i, j, i+1, 1, ct, data, lfce[i+1])+
		   erg4(i, j, j-1, 2, ct, data, lfce[j-1])+penalty(i, j, ct, data))))
      goto sub100;

    else {
      errmsg (100, 2);
      cin >> k;
    }
  }

  else if (i>=((number)-1)) {
    //up to one base hanging on i
    if (ep==w5[j-(number)-1]+penalty(i, j, ct, data)) {
      push (&stack, 1, j-(number)-1,
	    w5[j-(number)-1], 1);
      goto sub100;
    }
    else if (i!=(number)) {
      if (ep==(erg4(i, j, i+1, 1, ct, data, lfce[i+1])+penalty(i, j, ct, data)+
	       w5[j-(number)-1])){
	push (&stack, 1, j-(number)-1,
	      w5[j-(number)-1], 1);
	goto sub100;
      }
      else if (ep==(erg4(i, j, i+1, 1, ct, data, lfce[i+1]) +
		    erg4(i, j, j-1, 2, ct, data, lfce[j-1]) + w5[j-(number)-2]
		    +penalty(i, j, ct, data))) {
	push (&stack, 1, j-(number)-2,
	      w5[j-(number)-2], 1);
	goto sub100;
      }
    }
    else if (ep==(erg4(i, j, j-1, 2, ct, data, lfce[j-1])+
		  w5[j-(number)-2]+penalty(i, j, ct, data))) {
      push (&stack, 1, j-(number)-2, w5[j-(number)-2],
	    1);
      goto sub100;
    }
  }
  else if ((j==((number)+1))||(j==(number)+2)) {
    //up to one base hanging on j
    if (ep==w3[i+1]+penalty(i, j, ct, data)) {
      push (&stack, i+1, number, w3[i+1], 1);
      goto sub100;
    }
    else if (ep == (erg4(i, j, i+1, 1, ct, data, lfce[i+1])+penalty(i, j, ct, data)+ w3[i+2])) {
      push (&stack, i+2, number, w3[i+2], 1);
      goto sub100;
    }
    else if (j!=((number)+1)) {
      if (ep==(erg4(i, j, j-1, 2, ct, data, lfce[j-1])+penalty(i, j, ct, data)+w3[i+1])) {
	push(&stack, i+1, number, w3[i+1], 1);
	goto sub100;
      }
      else if (ep==(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] +
		    penalty(i, j, ct, data))) {
	push(&stack, i+2, number, w3[i+2], 1);
	goto sub100;
      }
    }
  }
  else {
    //perhaps i, j closes a multi-branch loop
    k = i +2;
    while (k<=(j-3)) {
      if (k!=number) {
	if (ep==(w->f(i+1, k)+w->f(k+1, j-1)+data->eparam[10]+
		 data->eparam[5] + penalty(i, j, ct, data))) {
	  //multi loop, no dangling ends on i, j
	  push (&stack, i+1, k, w->f(i+1, k), 0);
	  push (&stack, k+1, j-1, w->f(k+1, j-1), 0);
	  goto sub100;
	}
	else if (ep==(erg4(i, j, i+1, 1, ct, data, lfce[i+1])+penalty(i, j, ct, data)+w->f(i+2, k)+
		      w->f(k+1, j-1) +data->eparam[10]+data->eparam[6]+
		      data->eparam[5])) {
	  //multi loop i and j basepaired with i+1 dangling
	  push (&stack, i+2, k, w->f(i+2, k), 0);
	  push (&stack, k+1, j-1, w->f(k+1, j-1), 0);
	  goto sub100;
	}
	else if (ep==(erg4(i, j, j-1, 2, ct, data, lfce[j-1])+penalty(i, j, ct, data)+
		      w->f(i+1, k)+w->f(k+1, j-2)+data->eparam[10] +
		      data->eparam[6] + data->eparam[5])) {
	  //multi loop with j-1 dangling over a basepair of i, j
	  push (&stack, i+1, k, w->f(i+1, k), 0);
	  push (&stack, k+1, j-2, w->f(k+1, j-2), 0);
	  goto sub100;
	}
	else if (ep==(erg4(i, j, i+1, 1, ct, data, lfce[i+1])+
		      erg4(i, j, j-1, 2, ct, data, lfce[j-1]) + w->f(i+2, k) +
		      w->f(k+1, j-2) + data->eparam[10] + 2*(data->eparam[6]) +
		      data->eparam[5]+penalty(i, j, ct, data))) {
	  //multi loop with i+1 and j-1 dangling over pair i, j
	  push(&stack, i+2, k, w->f(i+2, k), 0);
	  push(&stack, k+1, j-2, w->f(k+1, j-2), 0);
	  goto sub100;
	}
      }
      else {
	if (ep==(w3[i+1]+w5[j-(number)-1]+penalty(i, j, ct, data))) {
	  //exterior loop, no dangling ends over i, j pair
	  push(&stack, i+1, number, w3[i+1], 1);
	  push(&stack, 1, (j-(number)-1),
	       w5[j-(number)-1], 1);
	  goto sub100;
	}
	else if (ep==(erg4(i, j, i+1, 1, ct, data, lfce[i+1])+penalty(i, j, ct, data)+w3[i+2]+
		      w5[j-(number)-1])) {
	  //exterior loop i+1 dangles over i, j pair
	  push(&stack, i+2, number, w3[i+2], 1);
	  push(&stack, 1, j-(number)-1,
	       w5[j-(number)-1], 1);
	  goto sub100;
	}
	else if (ep==(erg4(i, j, j-1, 2, ct, data, lfce[j-1])+penalty(i, j, ct, data)+w3[i+1]+
		      w5[j-(number)-2])) {
	  //exterior loop with j-1 dangling over i, j pair
	  push(&stack, i+1, number, w3[i+1], 1);
	  push(&stack, 1, j-(number)-2,
	       w5[j-(number)-2], 1);
	  goto sub100;
	}
	else if (ep==(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))) {
	  //exterior loop with j-1 and i+1 dangling over i, j pair
	  push(&stack, i+2, number, w3[i+2], 1);
	  push(&stack, 1, j-(number)-2,
	       w5[j-(number)-2], 1);
	  goto sub100;
	}
      }
      k++;
    }
    /////
    if (ct->intermolecular) {
      k = i +2;
      while (k<=(j-3)) {
	if (k!=number) {
	  if (ep==(w2->f(i+1, k)+w2->f(k+1, j-1)+
		   + penalty(i, j, ct, data))) {
	    //multi loop, no dangling ends on i, j
	    push (&stack, i+1, k, w2->f(i+1, k), 0);
	    push (&stack, k+1, j-1, w2->f(k+1, j-1), 0);
	    goto sub100;
	  }
	  else if (ep==(erg4(i, j, i+1, 1, ct, data, lfce[i+1])+penalty(i, j, ct, data)+w2->f(i+2, k)+
			w2->f(k+1, j-1))) {
	    //multi loop i and j basepaired with i+1 dangling
	    push (&stack, i+2, k, w2->f(i+2, k), 0);
	    push (&stack, k+1, j-1, w2->f(k+1, j-1), 0);
	    goto sub100;
	  }
	  else if (ep==(erg4(i, j, j-1, 2, ct, data, lfce[j-1])+penalty(i, j, ct, data)+
			w2->f(i+1, k)+w2->f(k+1, j-2))) {
	    //multi loop with j-1 dangling over a basepair of i, j
	    push (&stack, i+1, k, w2->f(i+1, k), 0);
	    push (&stack, k+1, j-2, w2->f(k+1, j-2), 0);
	    goto sub100;
	  }
	  else if (ep==(erg4(i, j, i+1, 1, ct, data, lfce[i+1])+
			erg4(i, j, j-1, 2, ct, data, lfce[j-1]) + w2->f(i+2, k) +
                     	w2->f(k+1, j-2) +penalty(i, j, ct, data))) {
	    //multi loop with i+1 and j-1 dangling over pair i, j
	    push(&stack, i+2, k, w2->f(i+2, k), 0);
	    push(&stack, k+1, j-2, w2->f(k+1, j-2), 0);
	    goto sub100;
	  }
	}
	k++;
      }
    }
  }
  //i, j must close a bulge or interior loop
  //         sub500:
  for (d=(j-i-3); d>=1; d--) {
    for (ip=(i+1); ip<=(j-1-d); ip++) {
      jp=d+ip;
      if (abs(ip-i+jp-j)<=(data->eparam[8])) {
	//if (en==(erg2(i, j, ip, jp, ct, data, 0,
	//			0)+v[ip][jp])) {
	if (en==(erg2(i, j, ip, jp, ct, data, fce[i][ip],
		      fce[jp][j])+v->f(ip, jp))) {

	  //output the internal/bulge loop energy to the dump file:
	  //if (filename != 0) {
	  // 	out << "bulge/internal loop starting with pair "<<i<<" and "<<j<<
	  //   	" ending with pair "<<ip<<" and "<<jp<<
	  //      " dG = "<<(v->f(i, j)-v->f(ip, jp))<<"\n";

	  //}

	  i=ip;
	  j=jp;
	  en=v->f(i, j);




	  goto sub300;
	}
      }
    }
  }




}



#define maxsort 90000 //The maximum number of basepairs within %cntrl8
//of the minimum free energy

void traceback(structure *ct, datatable *data, arrayclass *v, arrayclass *w, arrayclass *w2, int *w3, int *w5, int **fce,
	       bool *lfce, int vmin, int cntrl6, int cntrl8, int cntrl9) {

  bool flag, **mark;
  register int number;
  int ii = 0, sort;
  int rep, i;
  int iret, jret, numbp, count, count2, k1, k2, num, crit;
  //int heapi[maxsort+1], heapj[maxsort+1];
  int *heapi, *heapj;
  int cur, c, k, j, cntr;
  //stackstruct stack;
  int ji = 0; //(added during debugging)
  //int energy[maxsort+1];
  int *energy;
  //ofstream out("picked.out");






  ////////////////
  //cout << "traceback\n"<<flush;







  //mark keeps track of which pairs have been formed by the suboptimal routine


  sort = maxsort;
  number = ct->numofbases;





  //Construct heapi and heapj composed of pairs ij such that a structure
  //	containing ij has energy less than a given percent (ie:cntrl8) of
  //	the minimum folding energy

  if (vmin>=0) { //no viable structure found
    ct->numofstructures=1;
    ct->energy[1]=0;
    for (i=1; i<=ct->numofbases; i++) {
      ct->basepr[1][i]=0;
    }
    return;
  }

  //dynamically allocate space for mark:
  mark = new bool *[number + 1];
  for (i=0; i<=(number); i++)
    mark[i] = new bool [number + 1];


  //This is the traceback portion of the dynamic algorithm

  flag = true;
  rep = 1;

  for (count=1; count<=(number); count++) {
    for (count2=1; count2<=(number); count2++) {
      mark[count][count2]=false;
    }
  }


  crit = (int) (abs(vmin) * (double) cntrl8 / 100);

  //add limits to crit:
  if (crit>120) crit = 120;
  else if (crit<10) crit = 10;

  crit = crit + vmin;



  energy = new int [sort+1];
  heapi = new int [sort+1];
  heapj = new int [sort+1];

  num = 0;
  i = 1;
  j = 2;
  while (i<(number)) {
    if (num==sort) {
      //allocate more space for the heap
      delete[] heapi;
      delete[] heapj;
      delete[] energy;
      sort = 10*sort;
      heapi = new int [sort+1];
      heapj = new int [sort+1];
      energy = new int [sort+1];
      i = 1;
      j = 2;
      num = 0;

    }


    /////debugging:
    //if (i==1&&j==26) {
    // 	heapj[0] = (v->f(i, j)+v->f(j, i+number));
    //}


    //if ((v->f(i, j)+v->f(j, i+number))==crit) {
    // 	heapi[0] = crit;
    //}


    if ((v->f(i, j)+v->f(j, i+number))<=crit) {

      num++;
      heapi[num]=i;
      heapj[num]=j;
      energy[num] = (v->f(i, j)+v->f(j, i+number));
      j = j+cntrl9+1;
      if (j>number) {
	i++;
      	j=i+1;
      }
    }
    else {
      j++;
      if (j>number) {
      	i++;
	j=i+1;
      }
    }
  }

  ///debugging:
  //ofstream out2("heap.out");
  //out2 <<crit<<"\n";
  //for (i=1; i<=num; i++) {

  //	out2 << heapi[i]<<" "<<heapj[i]<<" "<<(v->f(heapi[i], heapj[i])+v->f(heapj[i], heapi[i]+number))<<"\n";

  //}





  //sort the base pair list:


  ///////////////////////////////////

  //make a heap:

  int q, up, ir;
  for (q=2; q<=num; q++) {
    cur = q;
    up = cur/2;
    while ((energy[cur]<energy[up])&&up>=1) {
      swap(&heapi[cur], &heapi[up]);
      swap(&heapj[cur], &heapj[up]);
      swap(&energy[cur], &energy[up]);
      cur = cur/2;
      up = up/2;
    }
  }

  //sort the heap:

  for (ir=num-1; ir>=1; ir--) {
    swap(&heapi[ir+1], &heapi[1]);
    swap(&heapj[ir+1], &heapj[1]);
    swap(&energy[ir+1], &energy[1]);

    up =1 ;
    c = 2;
    while (c<=ir) {
      if (c!=ir) {
	if (energy[c+1]<energy[c]) c++;
      }
      if (energy[c]<energy[up]) {
	swap(&heapi[c], &heapi[up]);
	swap(&heapj[c], &heapj[up]);
	swap(&energy[c], &energy[up]);
	up=c;
	c=2*c;
      }
      else c = ir+1;
    }
  }



  //out2 <<"\n";
  //for (i=1; i<=num; i++) {

  //	out2 << heapi[i]<<" "<<heapj[i]<<" "<<(v->f(heapi[i], heapj[i])+v->f(heapj[i], heapi[i]+number))<<"\n";

  //}

  cntr = num;


  //cout <<"2.5"<<flush;

  while (flag) {
    //This is routine to select the region of the structure to be
    //	folded, it allows for sub-optimal structure predictions
    //err=0;
    //Select the next valid unmarked basepair
    while (mark[heapi[cntr]][heapj[cntr]]) {
      if (cntr==1) {
      	flag=false;
	goto sub900;
      }
      cntr--;
    }
    iret = heapi[cntr];
    jret = heapj[cntr];
    rep++;

    //Traceback to find best structure on included fragment (ie:iret to jret)
    if (flag) {
      for (count=1; count<=2; count++) {
      	if (count==1) {
	  ii=iret;