algorithm.cpp

mfold

	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) {
	jmt = j%3;


	//search for an open bifurcation:
	for (k=i; k<=j-1; k++) {
	  if (k!=number) wmb[i][jmt] = min(wmb[i][jmt], w.f(i, k)+work[k+1][jmt]);
	  else wmb[i][jmt] = min(wmb[i][jmt], w3[i]+w5[j-number]);


	}

	if (ct->intermolecular) {
	  //intermolecular folding:

	  //search for an open bifurcation:
	  for (k=i; k<=j-1; k++) {
	    if (k!=number) {

	      wmb2[i][jmt] = min(wmb2[i][jmt], w2->f(i, k)+work[k+1][jmt]);
	    }
            	



	  }

	  w2->f(i, j) = min(w2->f(i, j), wmb2[i][j%3]);


	}


      }

      w.f(i, j) = min(w.f(i, j), wmb[i][j%3]);




      //Fill in work, the best energy for columns j, j-1, j-2
    sub3:

      work[i][j%3] = w.f(i, j);

      if (ct->intermolecular)
      	work2[i][j%3] = w2->f(i, j);





      //Calculate vmin, the best energy for the entire sequence
      if (j>(number)) {
      	vmin = min(vmin, v.f(i, j)+v.f(j-(number), i));

      }
    }

    //Compute w5[i], the energy of the best folding from 1->i, and
    //w3[i], the energy of the best folding from i-->numofbases
    if (j==(number)) {

      for (i=0; i<=(minloop+1); i++) {
      	if (lfce[i]) w5[i] = infinity;
	else if (i!=0) w5[i]=w5[i-1]; //condition added 8/18/97
	else w5[i] = 0; //added 8/18/97
	//w5[i]=0;
      }
      for (i=minloop+2; i<=(number); i++) {
      	if (lfce[i]) w5[i] = infinity;
	//else if (fce[i][i]==2) {
	//	w5[i] = w5[i-1] + data->init; //add the initiation
	//}															//for an intermolecular
	//interaction
	else w5[i] = w5[i-1];
	//w5[i]=w5[i-1];
      	for (k=1; k<=5; k++) e[k] = infinity; //e[k]=0;
      	for (k=0; k<=(i-4); k++) {
	  e[1] = min(e[1], (w5[k]+v.f(k+1, i)+penalty(i, k+1, ct, data)));

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

	  e[3] = min(e[3], (w5[k]+erg4(i-1, k+1, i, 1, ct, data, lfce[i])+v.f(k+1, i-1)+penalty(i-1, k+1, ct, data)));

	  e[4] = min(e[4], (w5[k]+erg4(i-1, k+2, k+1, 2, ct, data, lfce[k+1])+
			   erg4(i-1, k+2, i, 1, ct, data, lfce[i]) + v.f(k+2, i-1)+
			   penalty(i-1, k+2, ct, data)));

      	}
      	w5[i] = min(w5[i], e[1]);
      	w5[i] = min(w5[i], e[2]);
      	w5[i] = min(w5[i], e[3]);
      	w5[i] = min(w5[i], e[4]);
      }
      w3[0] = 0;
      w3[number+1] = 0;
      for (i=(number); i>=(number-minloop); i--)    //number+1 ... number-minloop
      	if (lfce[i]) w3[i] = infinity;
	else w3[i]=w3[i+1];
      //w3[i]=0;
      for (i=((number)-minloop-1); i>=1; i--) {
      	if (lfce[i]) w3[i] = infinity;
	//else if (fce[i][i]==2) {
	//	w3[i] = w3[i+1] + data->init; //add the initiation
	//}															//for an intermolecular
	//interaction
	else w3[i] = w3[i+1];
	//w3[i]=w3[i+1];
      	for (k=1; k<=5; k++) e[k] = infinity;
      	for (k=((number)+1); k>=(i+4); k--) {
	  e[1] = min(e[1], (v.f(i, k-1)+w3[k]+penalty(k-1, i, ct, data)));


	  e[2] = min(e[2], (v.f(i+1, k-1)+erg4(k-1, i+1, i, 2, ct, data, lfce[i])+penalty(k-1, i+1, ct, data) + w3[k]));

	  e[3] = min(e[3], (v.f(i, k-2)+erg4(k-2, i, k-1, 1, ct, data, lfce[k-1]) + penalty(k-2, i, ct, data) + w3[k]));

	  e[4] = min(e[4], (v.f(i+1, k-2)+erg4(k-2, i+1, i, 2, ct, data, lfce[i]) +
			   erg4(k-2, i+1, k-1, 1, ct, data, lfce[k-1])+w3[k]+
			   penalty(k-2, i+1, ct, data)));

      	}
      	w3[i] = min(w3[i], e[1]);
      	w3[i] = min(w3[i], e[2]);
      	w3[i] = min(w3[i], e[3]);
      	w3[i] = min(w3[i], e[4]);
      }
    }
    //fill in some work array values:
    if (j>=(number)) {
      for (k=j+1; k>=((number)+1); k--)
      	work[k][(j+1)%3] = w.f(k-(number), j+1-(number));

      if (ct->intermolecular) work2[k][(j+1)%3] = w2->f(k-(number), j+1-(number));

    }
  }

  //for (i=((number)+1); i<=(2*(number)-1); i++) {
  //	for (j=((number)+1); j<=(2*(number)); j++) {
  //   	v.f(i, j)=v.f(i-(number), j-(number));
  //      w.f(i, j)=w.f(i-(number), j-(number));
  //   }
  //}



  ////////////////
  //cout << "presave"<<"\n";



  ///debugging:
  //if (ct->intermolecular) {
  //ofstream dump("arrays.out");


  //output w and v along diagonals:
  //for (l=1; l<ct->numofbases; l++) {
  // 	for (i=1; i+l<=2*ct->numofbases; i++) {
  //    	j = i+l;
  //      dump << i <<" "<<j<<" "<<v.f(i, j)<<" "<<w.f(i, j)<<" "<<w2->f(i, j)<<"\n";
  //   }

  //}

  //for (j=1; j<=(2*(number)-1); j++) {
  //for (i=min(j, number); i>=max(1, j-(number)+1); i--) {
  //dump << i <<" "<<j<<" "<<v.f(i, j)<<" "<<w.f(i, j)<<" "<<w2->f(i, j)<<"\n";
  //}
  //}
  //dump.close();
  //}





  if (save!=0) {
    ofstream sav(save, ios::binary);

    //sav << ct->ctlabel[1];
    sav.write((ct->ctlabel[1]), ctheaderlength);

    convert.ch[0] = 0;
    convert.ch[1] = 0;
    convert.ch[2] = 0;
    convert.ch[3] = 0;

    convert.i = ct->numofbases;
    sav.write(convert.ch, 2);

    //sav << ct->numofbases << "\n";

    for (i=1; i<=ct->numofbases; i++) {
      convert.i = int (ct->numseq[i]);
      sav.write(convert.ch, 2);
      sav.write(ct->nucs+i, 1);
      convert.i = int (ct->hnumber[i]);
      sav.write(convert.ch, 2);

      //sav << ct->numseq[i]<<"\n";
    }

    convert.i = ct->npair;
    sav.write(convert.ch, 2);
    //sav << ct->npair << "\n";

    convert.i = ct->ndbl;
    sav.write(convert.ch, 2);
    //sav << ct->ndbl << "\n";

    convert.i = ct->nnopair;
    sav.write(convert.ch, 2);

    convert.i = ct->ngu;
    sav.write(convert.ch, 2);

    //sav << ct->nnopair << "\n";

    for (i=1; i<=ct->npair; i++) {
      convert.i = ct->pair[i][0];
      sav.write(convert.ch, 2);
      convert.i = ct->pair[i][1];
      sav.write(convert.ch, 2);

      //sav << ct->pair[i][0]<<"\n";
      //sav << ct->pair[i][1]<<"\n";
    }

    for (i=1; i<=ct->ndbl; i++) {
      convert.i = ct->dbl[i];
      sav.write(convert.ch, 2);
      //sav << ct->dbl[i] << "\n";
    }

    for (i=1; i<=ct->nnopair; i++) {
      convert.i = ct->nopair[i];
      sav.write(convert.ch, 2);
      //sav << ct->nopair[i] << "\n";
    }
    for (i=0; i<ct->ngu; i++) {
      convert.i = ct->gu[i];
      sav.write(convert.ch, 2);
    }

    if (ct->intermolecular) {
      //sav << 1 <<"\n";
      convert.i = int (1);
      sav.write(convert.ch, 1);
      for (i=0; i<3; i++) {
      	convert.i = ct->inter[i];
	sav.write(convert.ch, 2);
       	//sav << ct->inter[i] << "\n";
      }
    }
    else {
      convert.i = int (0);
      sav.write(convert.ch, 1);
      //sav << 0 <<"\n";
    }

    for (i=0; i<=ct->numofbases; i++) {
      //sav <<w5[i] << "\n";
      if (w5[i] < 0) {
      	convert.i = 0;
       	sav.write(convert.ch, 1);
	convert.i = -w5[i];
      }
      else {
      	convert.i = 1;
      	sav.write(convert.ch, 1);
      	convert.i = w5[i];
      }
      sav.write(convert.ch, 2);
    }
    for (i=0; i<=ct->numofbases+1; i++) {
      //sav <<w3[i] << "\n";

      if (w3[i] < 0) {
       	convert.i = 0;
      	sav.write(convert.ch, 1);
	convert.i = -w3[i];
      }
      else {
      	convert.i = 1;
      	sav.write(convert.ch, 1);
      	convert.i = w3[i];
      }
      sav.write(convert.ch, 2);
    }
    for (i=0; i<=ct->numofbases; i++) {
      for (j=0; j<=ct->numofbases; j++) {
      	//sav << v.dg[i][j] << "\n";
	//sav << w.dg[i][j] << "\n";
	if (v.dg[i][j] < 0) {
	  convert.i = 0;
	  sav.write(convert.ch, 1);
	  convert.i = -v.dg[i][j];
	}
	else {
	  convert.i = 1;
	  sav.write(convert.ch, 1);
	  convert.i = v.dg[i][j];
	}
	sav.write(convert.ch, 2);

	if (w.dg[i][j] < 0) {
	  convert.i = 0;
	  sav.write(convert.ch, 1);
	  convert.i = -w.dg[i][j];
	}
	else {
	  convert.i = 1;
	  sav.write(convert.ch, 1);
	  convert.i = w.dg[i][j];
	}
	sav.write(convert.ch, 2);

	if (ct->intermolecular) {
	  if (w2->dg[i][j] < 0) {
	    convert.i = 0;
	    sav.write(convert.ch, 1);
	    convert.i = -w2->dg[i][j];
	  }
	  else {
	    convert.i = 1;
	    sav.write(convert.ch, 1);
	    convert.i = w2->dg[i][j];
	  }
	  sav.write(convert.ch, 2);
	}

      }
    }
    //sav << vmin << "\n";
    if (vmin<0) {
      convert.i=0;
      sav.write(convert.ch, 1);
      convert.i = -vmin;
    }
    else {
      convert.i = 1;
      sav.write(convert.ch, 1);
      convert.i = vmin;

    }
    sav.write(convert.ch, 2);

    //also save the thermodynamic parameters
    //convert.i = ct->pair[i][0];
    //sav.write(convert.ch, 2);
    for (i=0; i<5; i++) {
      savefile(data->poppen[i], &sav);
    }
    savefile(data->maxpen, &sav);

    for (i=0; i<11; i++) {
      savefile(data->eparam[i], &sav);
    }
    for (i=0; i<6; i++) {
      for (j=0; j<6; j++) {