algorithm.cpp

mfold

				       [ct->numseq[helix[ip+1][0]]]);
		}
	      }

	      if (ip!=(sum-2)) {
		if ((helix[ip+2][1]-helix[ip+1][0])>1) {
		  coax[ip][ip+1] = min(coax[ip][ip+1],
				       data->tstackcoax[ct->numseq[helix[ip][0]+1]]
				       [ct->numseq[helix[ip+1][0]+1]]
				       [ct->numseq[helix[ip+1][1]]]
				       [ct->numseq[helix[ip+1][0]]] +
				       data->coaxstack[ct->numseq[helix[ip][0]]]
				       [ct->numseq[helix[ip][1]]]
				       [ct->numseq[helix[ip][0]+1]]
				       [ct->numseq[helix[ip+1][0]+1]]);
		}
	      }
	      else {
		//ip = sum - 2
		if (helix[sum-1][0]<ct->numofbases) {
		  coax[ip][ip+1] = min(coax[ip][ip+1],
				       data->tstackcoax[ct->numseq[helix[ip][0]+1]]
				       [ct->numseq[helix[ip+1][0]+1]]
				       [ct->numseq[helix[ip+1][1]]]
				       [ct->numseq[helix[ip+1][0]]] +
				       data->coaxstack[ct->numseq[helix[ip][0]]]
				       [ct->numseq[helix[ip][1]]]
				       [ct->numseq[helix[ip][0]+1]]
				       [ct->numseq[helix[ip+1][0]+1]]);
		}
	      }
	    }
	    else {//no possible stacks
	      coax[ip][ip+1] = coax[ip][ip] + coax[ip+1][ip+1];
	    }
	  }
	}
	else if (k>1) {
	  for (i=0; (i+k)<sum; i++) {
	    coax[i][i+k] = coax[i][i]+coax[i+1][i+k];
	    for (j=1; j<k; j++) {
	      coax[i][i+k] = min(coax[i][i+k],
				 coax[i][i+j]+coax[i+j+1][i+k]);
	    }
	  }


	}
	if (k==(sum-1)) {
	  ct->energy[count] = ct->energy[count] + coax[0][sum-1];
	  //cout << "count = "<<count<<"\n"<<flush;
	  //for (i=0; i<sum; i++) {
	  // 	for (j=i; j<sum; j++) {
	  // 		cout << "i "<<i<<"  j "<<j<<"  coax[i][j] = "<<coax[i][j]<<"\n";
	  // 	}
	  //}

	  //cin >> i;
	}







      }




      for (k=0; k<sum; k++) delete[] helix[k];
      delete[] helix;
      for (k=0; k<sum; k++) delete[] coax[k];
      delete[] coax;

      goto subroutine;
      /*
      // whittle away at 5' end:
      while (ct->basepr[count][i]==0&&ct->basepr[count][i+1]==0) {
      i++;
      if (i>=(j-1)) goto subroutine;
      }
      //whittle away at 3' end
      while (ct->basepr[count][j]==0&&ct->basepr[count][j-1]==0) {
      j--;
      if (i>=j-1) goto subroutine;
      }
      //i dangles on i+1
      if (ct->basepr[count][i]==0&&ct->basepr[count][i+1]>(i+1)) {
      ct->energy[count] = ct->energy[count]+
      min(0, erg4(ct->basepr[count][i+1], (i+1), i, 2, ct, data));
      i++;
      }
      //j dangles on j-1
      if ((ct->basepr[count][j]==0&&ct->basepr[count][j-1]!=0)&&
      ct->basepr[count][j-1]<(j-1)) {
      ct->energy[count]=ct->energy[count] +
      min(0, erg4(j-1, ct->basepr[count][j-1], j, 1, ct, data));
      j--;
      }
      //structure bifurcates
      if (ct->basepr[count][i]!=j) {
      k = ct->basepr[count][i];
      if (ct->basepr[count][k+1]!=0) {
      push (&stack, i, k, open, 0);
      push (&stack, k+1, j, open, 0);
      }
      else if(ct->basepr[count][k+2]==0) {
      push (&stack, i, k+1, open, 0);
      push (&stack, k+2, j, open, 0);
      }
      else if (erg4(k, i, k+1, 1, ct, data)<=
      erg4((ct->basepr[count][k+2]), k+2, k+1, 2, ct, data)) {
      push (&stack, i, k+1, open, 0);
      push (&stack, k+2, j, open, 0);
      }
      else {
      push (&stack, i, k, open, 0);
      push (&stack, k+1, j, open, 0);
      }
      goto subroutine;
      }
      push (&stack, i, j, open, null);
      goto subroutine;
      */
    }
  }

  for (i=0; i<=ct->numofbases; i++)
    delete[] fce[i];

  delete[] fce;


  return;
}



//void trace() {


//}

void trace(structure *ct, datatable *data, int ii, int ji,
	   arrayclass *v, arrayclass *w, arrayclass *w2, int rep, bool *lfce, int **fce, int *w3,
	   int *w5)
{
  int k, i, j, en, open, stz, ep, d, ip, jp;
  stackstruct stack;
  register int number;
  ofstream out;


  //if filename is not empty, open the file to add an energy
  //breakdown -- this is mostly important for debugging
  //if (filename != 0) {
  //   out.open(filename);
  //}


  number = ct->numofbases;
  //Here is the trace algorithm:
  //When count = 1, it provides the best structure on the included
  //	fragment (ie: iret to jret)
  //When count = 2, it provides the best structure on the excluded
  //	fragment


  //Set the stack counter to zero
  stack.sp = 0;

  //Zero the basepr array:
  if (ji<=(number)) {
    for (k=ii; k<=ji; k++) ct->basepr[rep][k] = 0;
  }
  else {
    for (k=1; k<=(ji-(number)); k++) ct->basepr[rep][k]=0;
    for (k=ii; k<=(number); k++) ct->basepr[rep][k]=0;
  }
  //Add forced pairs to the basepair array
  for (k=1; k<=ct->npair; k++) {
    ct->basepr[rep][ct->pair[k][0]]=ct->pair[k][1];
    ct->basepr[rep][ct->pair[k][1]]=ct->pair[k][0];
  }
  push (&stack, ii, ji, v->f(ii, ji), 0);
         
  j = 0;
 sub100:
  i = j;
  while(i==j) {
    //Take a fragment from the stack, ie: i to j
    //	with expected energy e
    //	open = 0 =>multi loop
    // open = 1 =>exterior loop

    pull (&stack, &i, &j, &en, &open, &stz);
    if (stz!=0) return;
  }


  //if (i==310&&j==330) {
  //	ct->energy[0] = 0;
  //}

  //If i and j are paired, goto sub300
  if (en==v->f(i, j)) goto sub300;

  if (open==0) {

    if (!ct->intermolecular) {
      while (en==(w->f(i+1, j)+data->eparam[6])) {
	//whittle away at 5' end
	i++;
	en=w->f(i, j);
	if(i>=j) goto sub100;
      }
      while (en==(w->f(i, j-1)+data->eparam[6])) {
	//whittle away at 3' end
	j--;
	en=w->f(i, j);
	if (i>=j) goto sub100;
      }
    }
    else {
      while (en==(w2->f(i+1, j))&&fce[i][i]!=2) {
	//whittle away at 5' end
	i++;
	en=w2->f(i, j);
	if(i>=j) goto sub100;
      }
      while (en==(w2->f(i, j-1))&&fce[j][j]!=2) {
	//whittle away at 3' end
	j--;
	en=w2->f(i, j);
	if (i>=j) goto sub100;
      }

      while (en==(w2->f(i+1, j)+infinity+data->init)&&fce[i][i]==2) {
	//whittle away at 5' end
	i++;
	en=w2->f(i, j);
	if(i>=j) goto sub100;
      }
      while (en==(w2->f(i, j-1)+infinity+data->init)&&fce[j][j]==2) {
	//whittle away at 3' end
	j--;
	en=w2->f(i, j);
	if (i>=j) goto sub100;
      }
      while (en==(w->f(i+1, j)+data->eparam[6])) {
	//whittle away at 5' end
	i++;
	en=w->f(i, j);
	if(i>=j) goto sub100;
      }
      while (en==(w->f(i, j-1)+data->eparam[6])) {
	//whittle away at 3' end
	j--;
	en=w->f(i, j);
	if (i>=j) goto sub100;
      }

    }

    if (en==(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))) {

      //i dangles over i+1, j
      i++;
      en=v->f(i, j);
    }
    else if (en==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)) {
      //j dangles over i, j-1
      j--;
      en=v->f(i, j);
    }
    else if (en==(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))) {
      //both i and j dangle over i+1, j-1
      i++;
      j--;
      en=v->f(i, j);
    }
    //check for stem closing a multi-loop
    if (en==(v->f(i, j)+data->eparam[10]+penalty(j, i, ct, data))) {
      en=v->f(i, j);
    }

    if (ct->intermolecular) {
      if (en==(v->f(i+1, j)+infinity+
	       erg4(j, i+1, i, 2, ct, data, lfce[i])+penalty(i+1, j, ct, data))) {

	//i dangles over i+1, j
	i++;
	en=v->f(i, j);
      }
      else if (en==v->f(i, j-1)+infinity+
	       erg4(j-1, i, j, 1, ct, data, lfce[j])+penalty(i, j-1, ct, data)) {
	//j dangles over i, j-1
	j--;
	en=v->f(i, j);
      }
      else if (en==(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))) {
	//both i and j dangle over i+1, j-1
	i++;
	j--;
	en=v->f(i, j);
      }
      //check for stem closing a multi-loop
      if (en==(v->f(i, j)+infinity+penalty(j, i, ct, data))) {
	en=v->f(i, j);
      }

    }


  }
  else {
    while ((j==number)&&(en==w3[i+1])){//&&(fce[i][i]!=2)) {
      //whittle away at 5' end
      i++;
      if (i>=j) goto sub100;
    }

    while ((i==1)&&(en==w5[j-1])){//&&(fce[j][j]!=2)) {
      //whittle away at 3' end
      j--;
      if (i>=j) goto sub100;
    }

    //remove the initiation penalty from exterior loops:
    /*while ((j==number)&&(en==(w3[i+1]+data->init))&&(fce[i][i]==2)) {
    //whittle away at 5' end
    en = w3[i+1];
    i++;
    if (i>=j) goto sub100;
    }

    while ((i==1)&&(en==(w5[j-1]+data->init))&&(fce[j][j]==2)) {
    //whittle away at 3' end
    en = w5[j-1];
    j--;
    if (i>=j) goto sub100;
    }
    while ((j==number)&&(en==w3[i+1])&&(fce[i][i]!=2)) {
    //whittle away at 5' end
    i++;
    if (i>=j) goto sub100;
    }

    while ((i==1)&&(en==w5[j-1])&&(fce[j][j]!=2)) {
    //whittle away at 3' end
    j--;
    if (i>=j) goto sub100;
    } */
    if (en==(v->f(i+1, j)+erg4(j, i+1, i, 2, ct, data, lfce[i])+
	     penalty(i+1, j, ct, data))) {
      //i dangles over i+1, j
      i++;
      en=v->f(i, j);
    }
    else if (en==(v->f(i, j-1)+erg4(j-1, i, j, 1, ct, data, lfce[j])+penalty(j-1, i, ct, data))) {
      //j dangles over i, j-1
      j--;
      en = v->f(i, j);
    }
    else if (en==(v->f(i+1, j-1)+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))) {
      //both j and i dangle over i+1, j-1
      i++;
      j--;
      en=v->f(i, j);
    }
    else if (en==(v->f(i, j)+penalty(i, j, ct, data))) {
      en = v->f(i, j);
    }
  }
  //At this point, the ends don't base pair and don't simply stack
  //on a helix, meaning the structure must bifurcate
  if (en!=v->f(i, j)) {
    k=i;
    while (k<(j-2)) {
      if ((open==0)&&(en==(w->f(i, k)+w->f(k+1, j)))) {
	//best structure is split into best structures on i, k and
	//k+1, j
	push (&stack, i, k, w->f(i, k), 0);
	push (&stack, k+1, j, w->f(k+1, j), 0);
	goto sub100;
      }
      //else if ((open==0)&&(en==(w->f(i, k)+w->f(k+1, j)+penalty(i, j, ct, data)))){
      //	//best structure is split into best structures on i, k and
      //   	//k+1, j
      //   push (&stack, i, k, w->f(i, k), 0);
      //   push (&stack, k+1, j, w->f(k+1, j), 0);
      //   goto sub100;
      //}
      //else if ((open==0)&&(en==(w->f(i, k)+w->f(k+1, j)-penalty(i, j, ct, data)))){
      //	//best structure is split into best structures on i, k and
      //   	//k+1, j
      //   push (&stack, i, k, w->f(i, k), 0);
      //   push (&stack, k+1, j, w->f(k+1, j), 0);
      //   goto sub100;
      //}
      else if ((open==1)&&(i==1)) {
	if (en==(w5[k]+v->f(k+1, j)+penalty(j, k+1, ct, data))) {
	  //best structure on 1, j splits into 1, k and pair k+1, j
	  push (&stack, 1, k, w5[k], 1);
	  push (&stack, k+1, j, v->f(k+1, j), 0);
	  goto sub100;
	}
	//else if (en==(w5[k]+v->f(k+1, j))) {
	//	//best structure on 1, j splits into 1, k and pair k+1, j
	//   push (&stack, 1, k, w5[k], 1);
	//   push (&stack, k+1, j, v->f(k+1, j), 0);
	//   goto sub100;
	//}
	else if (en==(w5[k]+erg4(j, k+2, k+1, 2, ct, data, lfce[k+1])
		      +v->f(k+2, j)+penalty(j, k+2, ct, data))) {
	  //best structure splits on 1, k splits into 1, k and
	  //	k+2 paired to j and k+1 stacks onto the pair
	  push (&stack, 1, k, w5[k], 1);
	  push (&stack, k+2, j, v->f(k+2, j), 0);
	  goto sub100;
	}
	else if (en==(w5[k]+erg4(j-1, k+1, j, 1, ct, data, lfce[j])+
		      v->f(k+1, j-1)+penalty(j-1, k+1, ct, data))) {
	  //best structure splits into best structure on 1, k
	  //	and the pair k+1, j-1 with j stacked onto the pair
	  push(&stack, 1, k, w5[k], 1);
	  push(&stack, k+1, j-1, v->f(k+1, j-1), 0);
	  goto sub100;

	}
	else if (en==(w5[k]+erg4(j-1, k+2, k+1, 2, ct, data, lfce[k+1]) +
		      erg4(j-1, k+2, j, 1, ct, data, lfce[j]) + v->f(k+2, j-1)
		      +penalty(j-1, k+2, ct, data))) {
	  //best structure on 1, j splits into 1, k and the pair
	  //	of k+2, j-1 and k+1 and j both stack
	  push(&stack, 1, k, w5[k], 1);
	  push(&stack, k+2, j-1, v->f(k+2, j-1), 0);
	  goto sub100;
	}
	// else if (en==(w5[k]+erg4(j-1, k+2, k+1, 2, ct, data, lfce[k+1]) +
	// 	erg4(j-1, k+2, j, 1, ct, data, lfce[j]) + v->f(k+2, j-1)
	//    )) {
	//    //best structure on 1, j splits into 1, k and the pair
	//	of k+2, j-1 and k+1 and j both stack
	//    push(&stack, 1, k, w5[k], 1);
	//    push(&stack, k+2, j-1, v->f(k+2, j-1), 0);
	//    goto sub100;
	// }
      }
      else if ((open==1)&&(j==number)) {
	if (en==(v->f(i, k+2)+w3[k+3]+penalty(k+2, i, ct, data))) {
	  //best structure on i, n splits into pair of i, k+2 and
	  //	k+3, number
	  push (&stack, i, k+2, v->f(i, k+2), 0);
	  push(&stack, k+3, number, w3[k+3], 1);
	  goto sub100;
	}
	else if (en==(v->f(i+1, k+2)+erg4(k+2, i+1, i, 2, ct, data, lfce[i]) +
		      w3[k+3]+penalty(i+1, k+2, ct, data))) {
	  //best structure on i, n splits into best structure of
	  //	pair i+1, k+2 with i stacked and k+3, number
	  push (&stack, i+1, k+2, v->f(i+1, k+2), 0);
	  push (&stack, k+3, number, w3[k+3], 1);
	  goto sub100;
	}
	else if (en==(v->f(i, k+1)+erg4(k+1, i, k+2, 1, ct, data, lfce[k+2])+
		      w3[k+3]+penalty(i, k+1, ct, data))) {
	  //best structure on i, n splits into pair i, k+1 with
	  //	k+2 stacked and k+3, number
	  push (&stack, i, k+1, v->f(i, k+1), 0);
	  push (&stack, k+3, number, w3[k+3], 1);
	  goto sub100;
	}
	else if (en==(v->f(i+1, k+1)+erg4(k+1, i+1, i, 2, ct, data, lfce[i])+
		      erg4(k+1, i+1, k+2, 1, ct, data, lfce[k+2])+w3[k+3]
		      +penalty(k+1, i+1, ct, data))) {
	  //best structure on i, number splits into pair of i+1, k+1
	  //	with both i and k+2 stacked and k+3, n
	  push (&stack, i+1, k+1, v->f(i+1, k+1), 0);
	  push (&stack, k+3, number, w3[k+3], 1);
	  goto sub100;
	}
      }
      k++;
    }

    errmsg(100, 1);
  }
  //Base pair found
 sub300:
  if (j<=(number)) {
    ct->basepr[rep][i]=j;
    ct->basepr[rep][j]=i;