//// This program generates the CPLEX input file
//// to solve the 2dHP protein folding problem
//// for a sequence of n 'h', with n an integer 
//// specified as parameter of the program.

//// Written by Alessandro Dal Palu'. 
//// Web: www2.unipr.it/~dalpalu


#include <stdio.h>
#include <string.h>
#include <math.h>

int size;

int get_point(int x, int y, int l)
{
	return x+y*size+l*size*size;
}

int main(char* argc, char** argv)
{
	
	
	FILE* f;
	char prim[100];
	
	for (int n=4;n<=26;n++)
	{
	
	for (int z=0;z<n;z++)
	 prim[z]='h';
	prim[n]=0;
	
	char fn[100];
	
	fn[0]='0'+n/10;
	fn[1]='0'+n%10;
	fn[2]=0;	

	strcat(fn,"h.txt");
	f=fopen(fn,"w+");
	
	printf("Length %d\n",n);
	
	int delta=(int)ceil(sqrt(n));
	size=2*delta+1;
	
   fprintf(f,"Maximize\n");
	for (int i=0;i<size-1;i++)	
	for (int j=0;j<size-1;j++)	
	{
	 fprintf(f,"+x%d ",get_point(i,j,n));
	 fprintf(f,"+x%d ",get_point(i,j,n+1));
	}



   fprintf(f,"\nSubject to\n");

	fprintf(f," x%d = 1\n",get_point(delta,delta,0));
	fprintf(f," x%d = 1\n",get_point(delta,delta+1,1));	
		
	/// flat layer
	for (int l=0;l<n;l++)	
	{
	 for (int i=0;i<size;i++)	
	 for (int j=0;j<size;j++)	
	   fprintf(f,"+ x%d ",get_point(i,j,l));			
	 fprintf(f," = 1\n");
	}	
	
	// no self loops
	for (int i=0;i<size;i++)	
	for (int j=0;j<size;j++)	
	{
	  for (int l=0;l<n;l++)	
	    fprintf(f,"+ x%d ",get_point(i,j,l));			
	  fprintf(f," <= 1\n");
	}	
	
	// next 
	for (int i=1;i<size-1;i++)	
	for (int j=1;j<size-1;j++)	
	{
	  for (int l=0;l<n-1;l++)	
	   {
		   fprintf(f,"x%d + x%d + x%d + x%d - x%d >= 0\n",
		         get_point(i+1,j,l+1),
		         get_point(i-1,j,l+1),
		         get_point(i,j+1,l+1),
		         get_point(i,j-1,l+1),
		         get_point(i,j,l));					   
		   fprintf(f,"x%d + x%d + x%d + x%d - x%d >= 0\n",
		         get_point(i+1,j,l),
		         get_point(i-1,j,l),
		         get_point(i,j+1,l),
		         get_point(i,j-1,l),
		         get_point(i,j,l+1));					   
		}		
	}	
	

	// energy -> map h on OR plane
	for (int i=0;i<size-1;i++)	
	for (int j=0;j<size-1;j++)	
	{
	  	
// hor
	  for (int l=0;l<n;l++)	
	  if (prim[l]=='h')
	    fprintf(f,"+ x%d",get_point(i,j,l));
	  for (int l=0;l<n;l++)	
	  if (prim[l]=='h')
	    fprintf(f,"+ x%d",get_point(i+1,j,l));
	    
	  fprintf(f,"- x%d <= 1\n",get_point(i,j,n));  ///// new var on next layer (i i+1 on n, j j+1 on n+1)
	  
	  fprintf(f," x%d ",get_point(i,j,n));  
	  for (int l=0;l<n;l++)	
	  if (prim[l]=='h')
	    fprintf(f,"- x%d",get_point(i,j,l));
	  fprintf(f, " <=0\n");
	  
	  fprintf(f," x%d ",get_point(i,j,n));  
	  for (int l=0;l<n;l++)	
	  if (prim[l]=='h')
	    fprintf(f,"- x%d",get_point(i+1,j,l));
	  fprintf(f, " <=0\n");
	  
//// vert
	  for (int l=0;l<n;l++)	
	  if (prim[l]=='h')
	    fprintf(f,"+ x%d",get_point(i,j,l));
	  for (int l=0;l<n;l++)	
	  if (prim[l]=='h')
	    fprintf(f,"+ x%d",get_point(i,j+1,l));
	    
	  fprintf(f,"- x%d <= 1\n",get_point(i,j,n+1));  
	  
	  fprintf(f," x%d ",get_point(i,j,n+1));  
	  for (int l=0;l<n;l++)	
	  if (prim[l]=='h')
	    fprintf(f,"- x%d",get_point(i,j,l));
	  fprintf(f, " <=0\n");
	  
	  fprintf(f," x%d ",get_point(i,j,n+1));  
	  for (int l=0;l<n;l++)	
	  if (prim[l]=='h')
	    fprintf(f,"- x%d",get_point(i,j+1,l));
	  fprintf(f, " <=0\n");

   }

   fprintf(f,"Binary\n");
   
	for (int l=0;l<n;l++)	
	 for (int i=0;i<size;i++)	
	 for (int j=0;j<size;j++)	
	  fprintf(f, "x%d\n",get_point(i,j,l));

		
   fprintf(f,"\nEnd\n");

	fclose(f);	
	
	}
}