dynamic.txt

mfold

	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++) {
      	for (k=0; k<6; k++) {
	  for (l=0; l<3; l++) {
	    savefile(data->dangle[i][j][k][l], &sav);
                
	  }
	}
      }
    }
    for (i=0; i<31; i++) {
      savefile(data->inter[i], &sav);
      savefile(data->bulge[i], &sav);
      savefile(data->hairpin[i], &sav);
    }
    for (i=0; i<6; i++) {
      for (j=0; j<6; j++) {
      	for (k=0; k<6; k++) {
	  for (l=0; l<6; l++) {
	    savefile(data->stack[i][j][k][l], &sav);
	    savefile(data->tstkh[i][j][k][l], &sav);
	    savefile(data->tstki[i][j][k][l], &sav);
	    savefile(data->coax[i][j][k][l], &sav);
	    savefile(data->tstackcoax[i][j][k][l], &sav);
	    savefile(data->coaxstack[i][j][k][l], &sav);
	    savefile(data->tstack[i][j][k][l], &sav);
	    savefile(data->tstkm[i][j][k][l], &sav);

	  }
	}
      }
    }
    for (i=0; i<maxtloop+1; i++) {
      for (j=0; j<2; j++) {
      	savefile(data->tloop[i][j], &sav);
      }
    }
    savefile(data->numoftloops, &sav);

    for (i=0; i<6; i++) {
      for (j=0; j<6; j++) {
      	for (k=0; k<6; k++) {
	  for (l=0; l<6; l++) {
	    for (m=0; m<6; m++) {
	      for (n=0; n<6; n++) {
		savefile(data->iloop11[i][j][k][l][m][n], &sav);
		for (o=0; o<6; o++) {
		  savefile(data->iloop21[i][j][k][l][m][n][o], &sav);
		  for (p=0; p<6; p++) {
		    savefile(data->iloop22[i][j][k][l][m][n][o][p], &sav);
		  }
		}
	      }
	    }
	  }
	}
      }
    }
    savefile(data->auend, &sav);
    savefile(data->gubonus, &sav);
    savefile(data->cint, &sav);
    savefile(data->cslope, &sav);
    savefile(data->c3, &sav);
    savefile(data->efn2a, &sav);
    savefile(data->efn2b, &sav);
    savefile(data->efn2c, &sav);
    savefile(data->numoftriloops, &sav);


    for (i=0; i<maxtloop+1; i++) {
      for (j=0; j<2; j++) {
      	savefile(data->triloop[i][j], &sav);
      }
    }
    convert.f = data->prelog;
    sav.write(convert.ch, 4);
    savefile(data->init, &sav);

    sav.close();
  }



  ////////////
  //cout <<"pretrace\n";

  for (i=0; i<2*number+3; i++) {
    delete[] work[i];
    delete[] wmb[i];
  }

  delete[] wmb;
  delete[] work;

  if (ct->intermolecular) {

    for (i=0; i<2*number+3; i++) {
      delete[] work2[i];
      delete[] wmb2[i];
    }

    delete[] wmb2;
    delete[] work2;

  }

  traceback(ct, data, &v, &w, w2, w3, w5, fce, lfce, vmin, cntrl6, cntrl8, cntrl9);





  delete[] lfce;
  de_allocate (fce, 2*number+1);







  delete[] w5;
  delete[] w3;

  if (ct->intermolecular) delete w2;



  return;
}