/*
*  v 1.1 of conformation generation program
*
* Compile with cc -o cgofr cgofr.c -lm 
*   (on HP, use c89 rather than cc)  
*
* Tim Dellinger and Jennifer Gerbi 12-14-98 for 
* UIUC Matse 390AS group project
* www.students.uiuc.edu/~gerbi/project.html
*
* this version optimized to use a 6-member-chain. The
* array sizes can be reduced if the user wants to use
* a 5-chain or a 4-chain.
*
* Input files needed: (all in the format of a return after
*                      each number in the file.)
*   The program prompts for sequence size input from user.
* 
*   input.dat: the sequence chain. 1=H, 2=P. 
*
*   inputinst.dat: the list of possible conformations. 
*     This is the set of permutations of 0,1,2,3 (left, up,
*     right, down), for one less than the chain member number
*     (since the first one is fixed). note the first number in the
*     file is the number of conformations to be tested, and the second
*     number is the number of instructions in each conformation (again,
*     this is one less than the sequence size.)
* 
*   energies.dat: the energy table giving HH, HP, PP, etc. energy
*     interaction integer magnitudes.
*
* Output files produced:
*   
*  output.dat: the list of all valid conformations and their energies.
*
*/

/* global variables:
 *
 * buildconform: the array we read in from mathmatica that gives us
 * all the permutations of the possible conformation building
 * directions: left,up,right,and down.
 *
 * conform: the array of positions where the conformation is built.
 *
 * permlength: the number of instructions in each conformation building
 *             information set.
 * sequence: the array holding the monomer sequence, as listed in
 * definitions just above.
 *
 */

#include <stdio.h>
#include <stddef.h>
#include <math.h>
#define FAIL 0
#define MAXNUMTWO 2200       /* min is twice the sequence size*/
#define MAXNUM 1100          /* min for this is the sequence size*/
#define MAXPERMS 2200        /* max numbers allowed, change if needed*/
#define LEFT 0
#define UP 1
#define RIGHT 2
#define DOWN 3
#define EMPTY 0
#define H 1
#define P 2
#define HH 3
#define M 4

/* global variables */

static int buildconform[MAXNUM][MAXPERMS];
int energy[5][5];
static int conform[MAXNUMTWO][MAXNUMTWO], sequence[MAXNUM]; 
int permlength;

main()
{

int i,ii,jj,j,k,n,m,testnum;
int numperms;
int result,totalenergy;
int posy,posx,FLAG,check;
int sequencesize;

     /*****fcn to read in sequence array *****/
int fgetints(char *name, int max);

     /******fcn to read in the conformation building instructions****/
int getinstructions(char *name);
     
     /*****fcn to read in the energy table ******/
int gettable(char *name);

    /*******fcn to calculate the energy of each conformation*****/
int calcenergy(int irow, int sequencesize);

    /******fcn to write out the output file "output.txt"******/   
int writeoutput(char *name,int energy,int irow,int sequencesize);

printf("enter in the sequence size");
scanf("%d",&sequencesize);

       /***** read in sequence array: file expected is "input.dat"*/
testnum = fgetints("input.dat", sequencesize);

        /****** read in instructions for generating the conformation.
         ******* this file is called "inputinst.dat"*******/
numperms = getinstructions("inputinst.dat"); 

/* end of initial data input */
/* make the initial conformation:
 * this puts the first sequence number in the middle position of the 
 * conformation array. the first one is always fixed this way to give
 * us an origin. posy and posx keep track of the current position 
 * where we are. they are incremented appropriately depending on the
 * conformation building instructions, given by buildconform[][]
 */

for (i=0;i<numperms;i++)  /*loop through each permutation set*/
 { posy = MAXNUM;
   posx = MAXNUM; /*must reset initial positions for each conformation!*/
     for (ii=0;ii<MAXNUMTWO;ii++)         /*must reset array too!*/       
        for (jj=0;jj<MAXNUMTWO;jj++)
           conform[ii][jj] = 0;    
   conform[posy][posx] = sequence[0]; /*initial position input*/
   FLAG=0;                  /* flag=0=valid conform. flag=1, badconform.*/

   for (j=0; j<permlength; j++) /*loop through the instruction set */
   {   
/********** 0=place left, 1=place up, 2=place right 3=down *********/

     switch (buildconform[i][j])
     {
      case LEFT:
	  if (conform[posy][posx-1] == 0)
	   {
	     conform[posy][posx-1] = sequence[j+1];
             posx = posx-1;    /* update position reference */
             break;
           }
	  else
           {FLAG=1;	
	   break;}
	
      case UP:
          if (conform[posy-1][posx] == 0)	
	   {
	   conform[posy-1][posx] = sequence[j+1];
           posy = posy-1;
           break;
	   }
	 else 
          {FLAG=1;
	  break;}

     case DOWN:	
         if (conform[posy+1][posx] == 0)	
	   {
	   conform[posy+1][posx] = sequence[j+1];
           posy = posy+1;
           break;
	   }
	 else 
          {FLAG=1;
	  break;}
	

      case RIGHT:
	 if (conform[posy][posx+1]==0)
	    {
            conform[posy][posx+1] = sequence[j+1];
	    posx = posx+1;
            break;
            }
	else 
        {FLAG=1;
         break;}
     
    default: printf("error in switch \n");
         break;

     }/*end switch*/
    } /*end j-loop for*/

/* now, inside the i-loop we call the energy function to find the energy
 * of this conformation, if it was successful. we also write the sucessful
 * conformation directions to a file with the energy which we will re-calc.
 * later with a mathmatica script for visuals. */

if (FLAG == 0)
  { totalenergy = calcenergy(i, sequencesize);
     check=writeoutput("output.txt",totalenergy,i,sequencesize);     
   }
} /* end of i loop*/
  
} /*end of main*/



int fgetints(char *name, int max)
{
int count;
FILE *fp  = fopen(name, "r");

if (fp !=  NULL)
{
   for (count=0; count < max; count++)
      if (fscanf(fp, "%d", &sequence[count]) !=1)
         break;
  fclose(fp);
}
return count;
}



int getinstructions(char *name)
{
int countrow,countcolumn,numperms;
FILE *fp = fopen(name, "r");

if (fp != NULL)
{
/* note first two numbers of the array buildconform will be the number of
  permuations and the size. */

fscanf(fp, "%d", &numperms);
fscanf(fp, "%d", &permlength);

   for (countrow=0; countrow<numperms; countrow++)
     for (countcolumn = 0; countcolumn<permlength;countcolumn++)
		fscanf(fp, "%d", &buildconform[countrow][countcolumn]);
fclose(fp);
}
return numperms;
}


int writeoutput(char *name,int totalenergy, int irow,int sequencesize)
{
  int count,check;
  FILE *fp  = fopen(name, "a");

  if (fp !=  NULL)
  {   check = 1;
      fprintf(fp, "%d", totalenergy );
      fprintf(fp," ");      
for (count=0; count < sequencesize-1; count++)

/* note: here we are using that the conformation instruction size is always
   less than the sequence size by 1. */

   {
           fprintf(fp, "%d", buildconform[irow][count]);
 }            
  
fprintf(fp,"\n");
fclose(fp);
}
return check;
}

int calcenergy(int irow, int sequencesize)
{
int i,j,k,check,initrow,initcol,totalenergy;
int gettable(char *name); 

/* first: read in energy table. it is formatted like an array: 0-5 on top,
 * 0-5 on the side, such that the interatctions are where the array lines
 * cross: for example, energy[0][0] is the empty-empty lattice site
 * energy, energy [1][1] is the 1-1 interatction, which are the
 * Hydrophobic interatctions (=-1 for the Dill Model). Note that
 * [01]=[10] , etc.
 */
     
check = gettable("energies.dat");
initrow = MAXNUM;
initcol = MAXNUM;
totalenergy = 0;
j=0;
for (i=0;i<sequencesize;i++)
{
   switch (conform[initrow][initcol])
   {
    case EMPTY:
    printf ("error! initial one should not be zero. check input file.\n");
    break;

    case H:
       switch (conform[initrow][initcol+1]) /*compare w/the right one*/
        {
         case EMPTY:
          totalenergy=totalenergy + energy[0][1];
          break; 
         case H:
          totalenergy=totalenergy + energy[1][1];
          break; 
         case P:
          totalenergy=totalenergy + energy[2][1];
          break; 
         case HH:
          totalenergy=totalenergy + energy[3][1];
          break; 
         case M:
          totalenergy=totalenergy + energy[4][1];
          break; 
	 default:
          printf("something wrong in energy case switch\n");
	  break;
        } /*end inner switch for H:right*/

switch (conform[initrow][initcol-1]) /*compare w/the left one*/
        {
         case EMPTY:
          totalenergy=totalenergy + energy[0][1];
          break; 
         case H:
          totalenergy=totalenergy + energy[1][1];
          break; 
         case P:
          totalenergy=totalenergy + energy[2][1];
          break; 
         case HH:
          totalenergy=totalenergy + energy[3][1];
          break; 
         case M:
          totalenergy=totalenergy + energy[4][1];
          break; 
	 default:
          printf("something wrong in energy case switch\n");
	  break;
        } /*end inner switch for left*/
switch (conform[initrow-1][initcol]) /*compare w/the upper one*/
        {
         case EMPTY:
          totalenergy=totalenergy + energy[0][1];
          break; 
         case H:
          totalenergy=totalenergy + energy[1][1];
          break; 
         case P:
          totalenergy=totalenergy + energy[2][1];
          break; 
         case HH:
          totalenergy=totalenergy + energy[3][1];
          break; 
         case M:
          totalenergy=totalenergy + energy[4][1];
          break; 
	 default:
          printf("something wrong in energy case switch\n");
	  break;
        } /*end inner switch for H: up*/
switch (conform[initrow+1][initcol]) /*compare w/the lowerone*/
        {
         case EMPTY:
          totalenergy=totalenergy + energy[0][1];
          break; 
         case H:
          totalenergy=totalenergy + energy[1][1];
          break; 
         case P:
          totalenergy=totalenergy + energy[2][1];
          break; 
         case HH:
          totalenergy=totalenergy + energy[3][1];
          break; 
         case M:
          totalenergy=totalenergy + energy[4][1];
          break; 
	 default:
          printf("something wrong in energy case switch\n");
	  break;
        } /*end inner switch for lower H*/



case P:
       switch (conform[initrow][initcol+1]) /*compare w/the right one*/
        {
         case EMPTY:
          totalenergy=totalenergy + energy[0][2];
          break; 
         case H:
          totalenergy=totalenergy + energy[1][2];
          break; 
         case P:
          totalenergy=totalenergy + energy[2][2];
          break; 
         case HH:
          totalenergy=totalenergy + energy[3][2];
          break; 
         case M:
          totalenergy=totalenergy + energy[4][2];
          break; 
	 default:
          printf("something wrong in energy case switch\n");
	  break;
        } /*end inner switch for right: P*/
switch (conform[initrow][initcol-1]) /*compare w/the left one*/
        {
         case EMPTY:
          totalenergy=totalenergy + energy[0][2];
          break; 
         case H:
          totalenergy=totalenergy + energy[1][2];
          break; 
         case P:
          totalenergy=totalenergy + energy[2][2];
          break; 
         case HH:
          totalenergy=totalenergy + energy[3][2];
          break; 
         case M:
          totalenergy=totalenergy + energy[4][2];
          break; 
	 default:
          printf("something wrong in energy case switch\n");
	  break;
        } /*end inner switch for P: left*/
switch (conform[initrow-1][initcol]) /*compare w/the upper one*/
        {
         case EMPTY:
          totalenergy=totalenergy + energy[0][2];
          break; 
         case H:
          totalenergy=totalenergy + energy[1][2];
          break; 
         case P:
          totalenergy=totalenergy + energy[2][2];
          break; 
         case HH:
          totalenergy=totalenergy + energy[3][2];
          break; 
         case M:
          totalenergy=totalenergy + energy[4][2];
          break; 
	 default:
          printf("something wrong in energy case switch\n");
	  break;
        } /*end inner switch for P:up*/
switch (conform[initrow+1][initcol]) /*compare w/the lower one*/
        {
         case EMPTY:
          totalenergy=totalenergy + energy[0][2];
          break; 
         case H:
          totalenergy=totalenergy + energy[1][2];
          break; 
         case P:
          totalenergy=totalenergy + energy[2][2];
          break; 
         case HH:
          totalenergy=totalenergy + energy[3][2];
          break; 
         case M:
          totalenergy=totalenergy + energy[4][2];
          break; 
	 default:
          printf("something wrong in energy case switch\n");
	  break;
        } /*end inner switch for lower*/




case HH:
       switch (conform[initrow][initcol+1]) /*compare w/the right one*/
        {
         case EMPTY:
          totalenergy=totalenergy + energy[0][3];
          break; 
         case H:
          totalenergy=totalenergy + energy[1][3];
          break; 
         case P:
          totalenergy=totalenergy + energy[2][3];
          break; 
         case HH:
          totalenergy=totalenergy + energy[3][3];
          break; 
         case M:
          totalenergy=totalenergy + energy[4][3];
          break; 
	 default:
          printf("something wrong in energy case switch\n");
	  break;
        } /*end inner switch for HH;right*/
switch (conform[initrow][initcol-1]) /*compare w/the left one*/
        {
         case EMPTY:
          totalenergy=totalenergy + energy[0][3];
          break; 
         case H:
          totalenergy=totalenergy + energy[1][3];
          break; 
         case P:
          totalenergy=totalenergy + energy[2][3];
          break; 
         case HH:
          totalenergy=totalenergy + energy[3][3];
          break; 
         case M:
          totalenergy=totalenergy + energy[4][3];
          break; 
	 default:
          printf("something wrong in energy case switch\n");
	  break;
        } /*end inner switch for HH-left*/
switch (conform[initrow-1][initcol]) /*compare w/the upper one*/
        {
         case EMPTY:
          totalenergy=totalenergy + energy[0][3];
          break; 
         case H:
          totalenergy=totalenergy + energy[1][3];
          break; 
         case P:
          totalenergy=totalenergy + energy[2][3];
          break; 
         case HH:
          totalenergy=totalenergy + energy[3][3];
          break; 
         case M:
          totalenergy=totalenergy + energy[4][3];
          break; 
	 default:
          printf("something wrong in energy case switch\n");
	  break;
        } /*end inner switch for HH-upper*/
switch (conform[initrow+1][initcol]) /*compare w/the lower one*/
        {
         case EMPTY:
          totalenergy=totalenergy + energy[0][3];
          break; 
         case H:
          totalenergy=totalenergy + energy[1][3];
          break; 
         case P:
          totalenergy=totalenergy + energy[2][3];
          break; 
         case HH:
          totalenergy=totalenergy + energy[3][3];
          break; 
         case M:
          totalenergy=totalenergy + energy[4][3];
          break; 
	 default:
          printf("something wrong in energy case switch\n");
	  break;
        } /*end inner switch for HH;lower*/



case M:
       switch (conform[initrow][initcol+1]) /*compare w/the right one*/
        {
         case EMPTY:
          totalenergy=totalenergy + energy[0][4];
          break; 
         case H:
          totalenergy=totalenergy + energy[1][4];
          break; 
         case P:
          totalenergy=totalenergy + energy[2][4];
          break; 
         case HH:
          totalenergy=totalenergy + energy[3][4];
          break; 
         case M:
          totalenergy=totalenergy + energy[4][4];
          break; 
	 default:
          printf("something wrong in energy case switch\n");
	  break;
        } /*end inner switch for M;right*/
switch (conform[initrow][initcol-1]) /*compare w/the left one*/
        {
         case EMPTY:
          totalenergy=totalenergy + energy[0][4];
          break; 
         case H:
          totalenergy=totalenergy + energy[1][4];
          break; 
         case P:
          totalenergy=totalenergy + energy[2][4];
          break; 
         case HH:
          totalenergy=totalenergy + energy[3][4];
          break; 
         case M:
          totalenergy=totalenergy + energy[4][4];
          break; 
	 default:
          printf("something wrong in energy case switch\n");
	  break;
        } /*end inner switch for M;left*/
switch (conform[initrow-1][initcol]) /*compare w/the up one*/
        {
         case EMPTY:
          totalenergy=totalenergy + energy[0][4];
          break; 
         case H:
          totalenergy=totalenergy + energy[1][4];
          break; 
         case P:
          totalenergy=totalenergy + energy[2][4];
          break; 
         case HH:
          totalenergy=totalenergy + energy[3][4];
          break; 
         case M:
          totalenergy=totalenergy + energy[4][4];
          break; 
	 default:
          printf("something wrong in energy case switch\n");
	  break;
        } /*end inner switch for M;up*/

switch (conform[initrow+1][initcol+1]) /*compare w/the down one*/
        {
         case EMPTY:
          totalenergy=totalenergy + energy[0][4];
          break; 
         case H:
          totalenergy=totalenergy + energy[1][4];
          break; 
         case P:
          totalenergy=totalenergy + energy[2][4];
          break; 
         case HH:
          totalenergy=totalenergy + energy[3][4];
          break; 
         case M:
          totalenergy=totalenergy + energy[4][4];
          break; 
	 default:
          printf("something wrong in energy case switch\n");
	  break;
        } /*end inner switch for M-down*/


  }/*end outer switch for energy adding*/
/*now we must figure out how to increment the conform array to move
 * where we want for the next comparison. Recall 0=left, 1=up, 2=right, 
 * 3=down
 */
      switch (buildconform [irow][j])
      { 
        case LEFT:
          initcol = initcol - 1;
          break;
       case UP:
          initrow = initrow - 1;
          break;
       case RIGHT:
          initcol = initcol + 1;
          break;   
      case DOWN:
         initrow = initrow +1;
         break;
      default:
         printf("error in energy switch increment\n");
    }  /* end incrementing switch */
j = j ++;  /*next time, look at the next step of the conformation
               * instruction line */
} /*end for loop over the entire sequence*/

/* note in the above, all energies are counted twice. thus we can safely
 * do integer division by 2 to get rid of this, since there energy
 * number will always be divisible by 2. */

totalenergy = totalenergy/2;

/*now we must subtract the energy in the linked parts of
 * the monomer */

for (k=0;k<sequencesize;k++)
   totalenergy = totalenergy-energy[sequence[k]][sequence[k+1]];           

return totalenergy;
} /*end function*/





int gettable(char *name)
{
int countrow,countcolumn,check;
FILE *fp = fopen(name, "r");

if (fp != NULL)
{  check = 1;
   for (countrow=0; countrow<5; countrow++)
     for (countcolumn = 0; countcolumn<5;countcolumn++)
                fscanf(fp, "%d", &energy[countrow][countcolumn]);
fclose(fp); 
}
return check;
}