// This version tries to trace the path and find the next maximum
// So far does not work
// Second maximum correct. Path not correct
#include <stdio.h> 
#include <stdlib.h>
#include <time.h>
//#define rows 15
//#define cols 15
#define rows 4
#define cols 8
struct cand{int v; int i; int j; int k; int dir;};
struct cand max, record[10000];
struct cand wpath[10000];
struct cand wpath1[10000];
struct cand npath[10000];
int aw[rows+1][rows+1][cols+1]; // max W gain
int an[rows+1][rows+1][cols+1]; // max N gain
int bw[rows+1][rows+1][cols+1]; int bn[rows+1][rows+1][cols+1]; 
int pref[rows+1][cols+1];
int fw(int i, int j, int k); int fn(int i, int j, int k);
int a[rows+1][cols+1]; // original matrix
int region[100][100]; int ri, rj, rk; int ii,jj,kk;
int mii, iii,jjj,kkk, i1;
int max1;
int count=0, maxcount=0;;

int main(void) { FILE *fp,*fp1;
        int cc, ccc;
        int i,m,n,l,trows, tcols,k=0,t,j,x, sum;
        int wn, maxwn=-1000000, im=0, jm=0, km=0;
//        fp=fopen("250x250.txt", "r");
        fp=fopen("4x8.txt", "r");
          if (fp==NULL) { printf("Cant open the file"); 
          exit(EXIT_FAILURE);
        }
        fscanf(fp,"%d",&trows); fscanf(fp,"%d",&tcols);
        for(m=1;m<=rows;m++){
          for(n=1;n<=cols;n++){ fscanf(fp, "%d", &i); a[m][n]=i;
          };
        };
        fclose(fp);
        fp1 = fopen("output.txt", "ab");
        if (fp1==NULL){ printf ("cant open the file"); exit(EXIT_FAILURE);
        }
  printf("here"); getchar();
/********** New addition prefix sum *********************/
                                          //          j
{int i, j, k;                             //   |     | |    |
  for(i=1; i<=rows; i++)                  //   |    k| |    |
    for(j=1; j<=cols; j++){               //  i|     | |    |
      pref[0][j]=0;                       //   |     | |    |
      for(k=1; k<=i; k++)
        pref[k][j]=pref[k-1][j]+a[k][j];
    }
}
/********************************************************/
    //W Calculation
    for(k=0;k<cols;k++){ // 1 --> 0
        for(t=0; t<=rows-1;t++){
            for(i=1; i<=rows-t; i++){               
              j=i+t; fw(i,j,k);
            };
        };
    };
    //N Calculation
    for(k=cols;k>=1;k--){
        for(t=0; t<=rows-1;t++){
            for(i=1; i<=rows-t; i++){               
                j=i+t; fn(i,j,k);
            };
        };
    };
    //WN Calculation
    //fprintf(fp1,"\n WN Calculation \n");
    for(k=1;k<=cols;k++){
        for(t=0; t<=rows-1;t++){
            for(i=1; i<=rows-t; i++){               
                j=i+t;
                sum=pref[j][k]-pref[i-1][k];
                wn=aw[i][j][k]+an[i][j][k]-sum;
                count++;
                record[count].v=wn;
                record[count].i=i;
                record[count].j=j;
                record[count].k=k;
                record[count].dir=0;
                if (wn>=maxwn){ maxwn = wn; im=i; jm=j; km=k; maxcount=count;
                }
            };
        };
    };     
   cc=count;
//    for(iii=1; iii<=rows+1; iii++)
 //   for(kkk=1; kkk<=cols; kkk++)aw[iii][iii-1][kkk]=-9999;
//    for(kkk=1; kkk<=cols; kkk++)aw[iii][iii-1][kkk]=0;
//    for(iii=1; iii<=rows; iii++)
//    for(jj=1; jj<=cols; jj++){ aw[iii][jj][0]=-9999; }
   printf("MaxWN[%d][%d][%d] =  %d \n", im,jm,km,maxwn);
   printf("aw =  %d \n", aw[im][jm][km]);
   printf("an =  %d \n", maxwn-aw[im][jm][km]);
   original_a(jm, km);
   kk=km; ii=im; jj=jm;
  contour(ii,jj,kk);
  getchar();
/*** Path finding ***/
   im=ii; jm=jj; km=kk;
   count=0;
   while((km>0)&&(im<=jm)){ 
      count++;
      wpath[count].v=aw[im][jm][km];
      wpath[count].i=im;
      wpath[count].j=jm;
      wpath[count].k=km;
      wpath[count].dir=bw[im][jm][km];
      if(bw[im][jm][km]==1)km--; else
      if(bw[im][jm][km]==2)im++; else
      if(bw[im][jm][km]==3)jm--;
   }
  printf("count %d ", count); getchar();
/*
  for(i1=1;i1<=count;i1++)printf("aw, i, j, k %3d %d %d %d\n", 
     aw[wpath[i1].i][wpath[i1].j][wpath[i1].k],
     wpath[i1].i, wpath[i1].j, wpath[i1].k);
 getchar();
*/
/*** second maximum candidate***/
   max.v=-999;
   for(i=1; i<=count; i++)if(record[i].v>max.v){
     max=record[i]; mii=i;
   }
   im=max.i; jm=max.j; km=max.k;
   record[maxcount].v=-999;
   max.v=-999;
   for(i=1; i<=cc; i++)if(record[i].v>max.v){
     max=record[i]; mii=i;
   }
  printf("second by global scanning = %d, anchor column %d %d %d", 
         max.v, max.i, max.j, max.k); getchar();

   im=max.i; jm=max.j; km=max.k;
//   dump(count);
/*** second maximum ***/
   max1=next(1,ii,jj,kk);;
  printf("max1= %d", max1); getchar();
  printf("second by trace back = %d, anchor column %d %d %d", 
    max1+an[ii][jj][kk]-(pref[jj][kk]-pref[ii-1][kk]), ii, jj, kk); 
  getchar();
   ccc=find(1);
   ii=wpath[ccc].i; jj=wpath[ccc].j; kk=wpath[ccc].k;
   if(wpath1[ccc].dir==1)kk--;
   if(wpath1[ccc].dir==2)ii++;
   if(wpath1[ccc].dir==3)jj--;
   contour(ii,jj,kk);
//   dump(count);
   fclose(fp1);
}  // end main
int next(int c, int i, int j, int k){
  int f1, f2, f3, sum, max1, i1, j1, k1, dir;
  int x;
//  if(c<=count){
  if(c<count){
    sum=pref[j][k]-pref[i-1][k];          
    f1=aw[i][j][k-1]+sum; 
    f2=aw[i+1][j][k]+a[i][k]; 
    f3=aw[i][j-1][k]+a[j][k];
    x=maximum(f1,f2,f3);
    if(bw[i][j][k]==1)k1=k-1; else k1=k;
    if(bw[i][j][k]==2)i1=i+1; else i1=i;
    if(bw[i][j][k]==3)j1=j-1; else j1=j;
//  printf("x=   %5d i, j, k %d %d %d | %d %d %d", x,i,j,k,f1,f2,f3); getchar();
    max1=next(c+1,i1,j1,k1); 
    if(bw[i][j][k]==1)f1=max1+sum;    // Second maximum not yet
    if(bw[i][j][k]==2)f2=max1+a[i][k]; 
    if(bw[i][j][k]==3)f3=max1+a[j][k];
    x=maximum(f1,f2,f3);
//  printf("max= %5d i, j, k %d %d %d | %d %d %d", x,i,j,k,f1,f2,f3); getchar();
    if(k==1)x=-9999;
    if(x==f1)wpath1[c].dir=1;
    if(x==f2)wpath1[c].dir=2;
    if(x==f3)wpath1[c].dir=3;
    return x;
  } else return -9999;
}
int find(int c){
while(wpath[c].dir==wpath1[c].dir){c++; printf("here"); getchar();}
return c;
}
/*** Use the same wpath data structure repeatedly ***/
int maximum(int x, int y, int z){
  if((x>=y)&&(x>=z)){ return x;}
  if((y>=x)&&(y>=z)){ return y;}
  if((z>=x)&&(z>=y)){ return z;}
}
dump(int cc){int i;
for(i=1;i<=cc;i++)printf("v %6d i %2d j %2d k %2d dir %2d\n", wpath[i].v, 
                              wpath[i].i,
                              wpath[i].j, wpath[i].k, wpath[i].dir);
}
original_a(int j, int k){ // Output original matrix
int ii,jj;
for(ii=1;ii<=rows; ii++){
  printf("%2d  ", ii);
  for(jj=1;jj<=cols;jj++)printf("%3d ", a[ii][jj]); printf("\n");
}
}
contour(int im, int jm, int km){
   {int i, j;
      for(i=1;i<=rows; i++){
        for(j=1;j<=cols;j++)region[i][j]=0;
      }
   }
   while((km>0)&&(aw[im][jm][km]>0)){ // Record position (im,km) and (jm,km) in region
      region[im][km]=1; region[jm][km]=1;
      if(bw[im][jm][km]==1)km--; else
      if(bw[im][jm][km]==2)im++; else
      if(bw[im][jm][km]==3)jm--;
   }
   im=ii; jm=jj; km=kk;
   while(km<=cols){
      region[im][km]=1; region[jm][km]=1;
      if(bn[im][jm][km]==1)km++; else
      if(bn[im][jm][km]==2)im++; else
      if(bn[im][jm][km]==3)jm--;
   }
   {int i, j;
      for(i=1;i<=rows; i++){
        printf("%2d  ", i);
        for(j=1;j<=cols;j++)printf("%3d ", region[i][j]); printf("\n");
      }
   }
}
int fw(int i, int j, int k){
int f1,f2,f3,sum=0,l,pp;
    if (k<=0) aw[i][j][k]=0;  
    if(k==1 && i==j) aw[i][j][k]=a[i][k];
    sum=pref[j][k]-pref[i-1][k];  // New line
    f1=aw[i][j][k-1]+sum;
    f2=aw[i+1][j][k]+a[i][k];
    f3=aw[i][j-1][k]+a[j][k];
    //return the maximum of f1,f2,f3
    if(i<=j){  // new
    if     (f1>=f2 && f1>=f3){ aw[i][j][k]=f1; bw[i][j][k]=1; }
    else{if(f2>=f1 && f2>=f3){ aw[i][j][k]=f2; bw[i][j][k]=2; }
    else                     { aw[i][j][k]=f3; bw[i][j][k]=3; }
    }
    }  // new
}
int fn(int i, int j, int k){
    int f1,f2,f3,sum=0,l,pp;
    //IF the colimn number is 0 then return 0
    if (k>=cols)an[i][j][k]=0;
    if (j<i) an[i][j][k]=0;
    //initialize fn(1,1,n), fn(2,2,n), fn(3,3,n) etc.
    if(k==cols-1 && i==j) an[i][j][k]=a[i][k];
    sum=pref[j][k]-pref[i-1][k];  // New line
    // calculate f1, f2, f3
    f1=an[i][j][k+1]+sum;
    f2=an[i+1][j][k] + a[i][k];;
    f3=an[i][j-1][k]+a[j][k];
    //return the maximum of f1,f2,f3
    if (f1>=f2 && f1>=f3)      {an[i][j][k]=f1; bn[i][j][k]=1;}
    else{ if(f2>=f1 && f2>=f3) {an[i][j][k]=f2; bn[i][j][k]=2;} 
        else                   {an[i][j][k]=f3; bn[i][j][k]=3;} 
    }
}