#include <stdio.h>
#include <math.h>
/*.proc files #include at end*/

#define LEVELS 20
#define TRUE 1
#define FALSE 0
#define DLEVELS 10
#define DRND 0.1
#define DALPHA 0.9524

struct   pi {  
  float  Jmin;
  float  denom[3][2][8];
  float  prob;
  float  value[8];
  short  Opti;
  short  Optj;
  short  prt[3];
  short  optflag;

  struct pi *child[3][2][8];         
  struct pi *left;                    
  struct pi *right;
  struct pi *parent;
};
typedef struct pi *list;
short levels=LEVELS,
      constrain=FALSE,
      ivalue,
      jvalue;

long int nnodes=0;

float alpha=DALPHA,
      rnd=DRND;

FILE  *fsens,
      *fdata,
      *fyear;

main() {
  void  make_sens(),                    get_constants(),
	sort_Jmin(),                    determine_optimum(),
	init_to_zero(),                 pick_mtrx(),
	set_p_mtrx(),                   init_Jmin(),
	open_files(),                   get_constraints(),
	print_tree(),                   set_ca_mtrx(),
	combine_nodes(),                set_cu_mtrx(),
	set_cx_mtrx(),                  unalloc(),
	merge_nodes(),                 	set_q_mtrx(),
        graph_year(),                   graph_data(),
        get_pi_val(),                   close_files(),
        write_pi(),                     round(),
        timestamp(),                    pr_mtrx();

  float assign_value(),
        compare_nodes(),
        p[3][8][8], 
        q[2][8][8], 
        ca[3][8], 
        cu[8], 
        cx[2][7];

  int   make_eight(),
        expand_tree(),
        alloc(),
        make_fortyeight(),
        rndmax(),
        pchoice, 
        qchoice,
        madelevel;
  
  list  side[LEVELS+1];
  short z;
  
  system("date > .time");
  system("date");

  init_to_zero(p, 3);
  init_to_zero(q, 2);

  pick_mtrx(&pchoice, &qchoice);
  get_constants(&rnd, &alpha, &levels);
  get_constraints(&constrain, &ivalue, &jvalue);

  set_p_mtrx(p, pchoice);
  set_q_mtrx(q, qchoice);
  set_ca_mtrx(ca);
  set_cu_mtrx(cu);
  set_cx_mtrx(cx);

  get_pi_val(side);

  printf("\nMaking level %d... ",z+1);
  if (constrain)
    make_eight(side, p, q, 0, ivalue, jvalue);
  else
    make_fortyeight(side, p, q, 0); 
  printf("No. of nodes made:%d\n",nnodes);

  printf("  Combining level %d... ",z+1);
  combine_nodes(side[1], 1);
  printf("No. of nodes after merging:%d\n\n",nnodes);

  for(z=1;z<levels;z++) {
    printf("Making level %d... ", z+1);
    madelevel=expand_tree(side, p, q, z);
    printf("No. of nodes made:%d\n",nnodes);
    if(madelevel) {
      printf("  Combining level %d... ",z+1);
      combine_nodes(side[z+1], z+1); 
      printf("No. of nodes after merging:%d\n\n",nnodes);
    }
  }
  init_Jmin(side[levels]);
  determine_optimum(side, ca, cu, cx, p, qchoice);  

  sort_Jmin(side);
  
  open_files();

  graph_year(side[0],pchoice, qchoice, alpha, rnd);

  for(z=0;z<levels;z++) {
    make_sens (side[0], z, pchoice, qchoice);
    print_tree (side[z], z);
  }
  close_files();
  system("date");
  timestamp("/mit/smmadana/graph.data");

}
/****************************************************************************/
int expand_tree(side, p, q, z)    
list side[];                  
float p[][8][8], q[][8][8]; 
short z; {
  int outofmem=FALSE, made48;
  list tempside;

  tempside=side[z];
  while(side[z]!=NULL && outofmem==FALSE) {
    if (constrain)
      outofmem=make_eight(side, p, q, z, ivalue, jvalue);
    else
      outofmem=make_fortyeight(side, p, q, z);     
    side[z]=side[z]->left; 
  }
  side[z]=tempside;
  if (outofmem==TRUE)
    unalloc(tempside);
  made48= (outofmem) ? FALSE : TRUE;
  return(made48);          
}
/***************************************************************************/
int make_fortyeight(side, p, q, z) 
list side[];
float p[][8][8], q[][8][8];
short z; { 
  short i, j, k, l, m, n, outofmem;
  float total, dummy, dummy2;
  for(i=0;i<3;i++) {
    for(j=0;j<2;j++)  
      for(l=0;l<8;l++){
        if(alloc(side[z],i,j,l,z)==FALSE)
          return(outofmem=TRUE);
	dummy2=0.0;
        for(m=0;m<8;m++) {
          dummy=0.0;
          for(n=0;n<8;n++) {
            dummy=dummy +(p[i][n][m]*(side[z]->value[n]));
          }
          dummy2=dummy2+(q[j][m][l]*dummy);
        }
        total=0.0;
        for(k=0;k<8;k++) {
        if(dummy2!=0)
            side[z]->child[i][j][l]->value[k]=
     		       assign_value(side[z], p[i], k, q[j][k][l],dummy2);
        else side[z]->child[i][j][l]->value[k]=0.0;
          total+=side[z]->child[i][j][l]->value[k];
        }
        if(total==0.0)
          free((char *)(side[z]->child[i][j][l]));
        else{
          nnodes++;
          side[z]->child[i][j][l]->parent=side[z];
          side[z]->child[i][j][l]->left=side[z+1];
          if(side[z+1]!=NULL)
            side[z+1]->right=side[z]->child[i][j][l];
          side[z+1]=side[z]->child[i][j][l];
          side[z+1]->prt[0]=i;
          side[z+1]->prt[1]=j;
          side[z+1]->prt[2]=l;             
          side[z+1]->optflag=FALSE;
          side[z]->denom[i][j][l]=dummy2;

          side[z+1]->prob=side[z]->denom[i][j][l]
                                *side[z]->prob;
        }
      }
    }
  return(outofmem=FALSE);
}
/*****************************************************************************/
void combine_nodes(side,z) 
list side; 
short z; { 
  list tempside, preside,tempside2;

  tempside=preside=side;
  while(tempside!=NULL) {
    side=tempside->left;
    while(side!=NULL) {
      tempside2=side->left;
      if((compare_nodes(side, tempside)) < 0.1) {
        merge_nodes(side, tempside);
        free((char *)(side));
        nnodes--;
      }
      side=tempside2;
    }
    tempside=tempside->left;
  }
  side=preside;
}
/****************************************************************************/
float compare_nodes(side, tempside)
list side, tempside; {
  short l;
  float diff=0.0;
  for(l=0;l<8;l++)
    diff+= fabs((side->value[l])-(tempside->value[l]));
  return(diff);
}
/*****************************************************************************/
void merge_nodes(side, tempnode)
list side, tempnode; {
  side->parent->child[side->prt[0]][side->prt[1]][side->prt[2]]=tempnode;
  tempnode->prob+=side->prob;
  if(side->left!=NULL){
    side->left->right=side->right;
    side->right->left=side->left;
  }
  else 
    side->right->left=NULL;
}
/*****************************************************************************/
float assign_value(head, p, k, q, dummy2)
list head;
short k;      
float p[][8], q,dummy2; {
  float dummy1=0.0;
  short n;

  for(n=0;n<8;n++)
    dummy1+=(p[n][k]*(head->value[n]));
  if(dummy2==0.0)
    return(0.0);
  else
    return(q*dummy1/dummy2);
} 
/*****************************************************************************/
void init_Jmin(side) 
list side; {                            
  list tempside;

  tempside=side;
  while(tempside!=NULL) {
    tempside->Jmin=0.0;
    tempside=tempside->left;
} }
/*****************************************************************************/
void determine_optimum(side, ca, cu, cx, p,qchoice) 
list side[]; 
int qchoice;                          
float ca[][8], cu[], cx[2][7], p[][8][8]; {
  list tempside;
  float cofm, cofu, cofa, future;
  short z;
  int i, j, l, m, k;

  if (constrain) {
    for(z=levels-1;z>(-1);z--) {
      tempside=side[z];
      i=ivalue;
      j=jvalue;
      while(tempside!=NULL) {
        tempside->Jmin=10000000.0;
	cofa=cofu=0.0;
	for(k=0;k<8;k++) {
	  cofa+=(tempside->value[k]*ca[i][k]);
	  for(m=0;m<8;m++) 
	    cofu+=(tempside->value[k]*p[i][k][m]*cu[m]);
        } 
	cofu*=alpha;
	cofm=alpha*cx[j][qchoice];
	future=0.0;
	for(l=0;l<8;l++) {
	  future+=(tempside->denom[i][j][l] * 
		   tempside->child[i][j][l]->Jmin);
	}
	future=(alpha*future)+cofa+cofu+cofm;
	if(future<tempside->Jmin)
	   tempside->Jmin=future;
	tempside->Opti=i;
	tempside->Optj=j;
        tempside=tempside->left;
  } } } 
  else {
    for(z=levels-1;z>(-1);z--) {
      tempside=side[z];
      while(tempside!=NULL) {
	tempside->Jmin=10000000.0;
	for(i=0;i<3;i++) {
	  cofa=cofu=0.0;
	  for(k=0;k<8;k++) {
	    cofa+=(tempside->value[k]*ca[i][k]);
	    for(m=0;m<8;m++) 
	      cofu+=(tempside->value[k]*p[i][k][m]*cu[m]);
	  } 
	  cofu*=alpha;
	  for(j=0;j<2;j++) {
	    cofm=alpha*cx[j][qchoice];
	    future=0.0;
	    for(l=0;l<8;l++) {
	      future+=(tempside->denom[i][j][l] * 
		     tempside->child[i][j][l]->Jmin);
	    }
	    future=(alpha*future)+cofa+cofu+cofm;
	    if(future<tempside->Jmin) {
	      tempside->Jmin=future;
	      tempside->Opti=i;
	      tempside->Optj=j;
        } } }
      tempside=tempside->left;
  } } }
  printf("Optimum nodes determined... ");
}    
/*****************************************************************************/
void sort_Jmin(side)
list side[]; {
  short z,l; 
  list tempnode;
  for(z=0;z<levels;z++){
    tempnode=side[z];
    while(side[z]!=NULL) {
      if(side[z]->optflag==TRUE) {
	for(l=0;l<8;l++) {
          if(side[z]->child[side[z]->Opti][side[z]->Optj][l]!=NULL)
	     side[z]->child[side[z]->Opti][side[z]->Optj][l]->optflag=TRUE;
      } }
      side[z]=side[z]->left;
    }
    side[z]=tempnode;
  }
  printf("Optimum nodes sorted.\n");
}
/*****************************************************************************/
void init_to_zero(a,i)
short i;
float a[][8][8];{
  short loop1, loop2, loop3;

  for(loop1=0;loop1<i;loop1++)  
    for(loop2=0;loop2<8;loop2++)
	for(loop3=0;loop3<8;loop3++)
	  a[loop1][loop2][loop3]=0.0;
}
/*****************************************************************************/
void pick_mtrx(pchoice, qchoice)
int *pchoice,*qchoice; {
  printf("\n\nEnter choices of matrixes...\n");
  printf("Choice of p matrix:");
  scanf("%d", pchoice);
  printf("Choice of q matrix:");
  scanf("%d", qchoice);
  printf("\n");
}
/*****************************************************************************/
void get_pi_val(side)
list  side[]; {
  int k;

  printf("\nEnter intial pi values...\n\n");
  side[0]=(list) malloc(sizeof(struct pi));
  for(k=0; k<8; k++){
    printf("  Pi value[%d]=",k);
    scanf("%f,",&side[0]->value[k]);
  }
  side[0]->prob=1.0;
  side[0]->optflag=TRUE;
}

/*****************************************************************************/
int alloc(side,i,j,l,z)
list side;
int i,j,l,z; {
  int finished=TRUE;

  side->child[i][j][l]=(list) malloc(sizeof(struct pi));
  if(side->child[i][j][l]==NULL) {
    finished=FALSE;
    printf("out of memory...adjusting level value\n");
    printf("Maximum no. of nodes reached:%d\n",nnodes);
    levels=z;
    printf("level adjusted to %d\n",levels);
  } 
  return(finished);
}
/*****************************************************************************/
void unalloc (side)
list side; {
  short i,j,l;
  list tempside;

  tempside=side;
  while (side!=NULL){
    for(i=0;i<3;i++)
      for(j=0;j<2;j++)
        for(l=0;l<8;l++){
          if(side->child[i][j][l]!=NULL)
            free((char *)(side->child[i][j][l]));
    }
    side=side->left;
  }
  side=tempside;
}
/*****************************************************************************/
void get_constants(rnd,alpha,levels)
short *levels;
float *rnd, *alpha; {
  printf("\nTo assign the default value to the constants, enter a zero.\n\n");
  printf("Enter value for levels (default is 10, maximum is 20):");
  scanf("%d",levels);
  if (*levels==0)
    *levels=DLEVELS;
  printf("\nEnter value for the interval (default is 0.1, minimum is 0.005):");
  scanf("%f",rnd);
  if (*rnd==0.0)
    *rnd=DRND;
  printf("\nEnter value for alpha (default is .9524):");
  scanf("%f",alpha);
  if (*alpha==0.0)
    *alpha=DALPHA;
}
/*****************************************************************************/
void get_constraints(constrain, ivalue, jvalue)
short *constrain, *ivalue, *jvalue; {
  printf("\nTo constrain i and j, enter a '1',\n");
  printf("otherwise enter a '0': ");
  scanf("%d", constrain);
  if (*constrain){
    printf("\nWhat is the i value? ");
    scanf("%d",ivalue);
    printf("\nWhat is the j value? ");
    scanf("%d", jvalue);
  }
  printf("\n");
}


/*****************************************************************************/
int make_eight(side, p, q, z, i, j) 
list side[];
float p[][8][8], q[][8][8];
short i, j, z; { 
  short k, l, m, n, outofmem;
  float total, dummy, dummy2;
      for(l=0;l<8;l++){
        if(alloc(side[z],i,j,l,z)==FALSE)
          return(outofmem=TRUE);
	dummy2=0.0;
        for(m=0;m<8;m++) {
          dummy=0.0;
          for(n=0;n<8;n++) {
            dummy=dummy +(p[i][n][m]*(side[z]->value[n]));
          }
          dummy2=dummy2+(q[j][m][l]*dummy);
        }
        total=0.0;
        for(k=0;k<8;k++) {
        if(dummy2!=0)
            side[z]->child[i][j][l]->value[k]=
     		       assign_value(side[z], p[i], k, q[j][k][l],dummy2);
        else side[z]->child[i][j][l]->value[k]=0.0;
          total+=side[z]->child[i][j][l]->value[k];
        }
        if(total==0.0)
          free((char *)(side[z]->child[i][j][l]));
        else{
          side[z]->child[i][j][l]->parent=side[z];
          side[z]->child[i][j][l]->left=side[z+1];
          if(side[z+1]!=NULL)
            side[z+1]->right=side[z]->child[i][j][l];
          side[z+1]=side[z]->child[i][j][l];
          side[z+1]->prt[0]=i;
          side[z+1]->prt[1]=j;
          side[z+1]->prt[2]=l;             
          side[z+1]->optflag=FALSE;
          side[z]->denom[i][j][l]=dummy2;

          side[z+1]->prob=side[z]->denom[i][j][l]
                                *side[z]->prob;
        }
      }
  return(outofmem=FALSE);
}
/*****************************************************************************/
void pr_mtrx (mtrx)
float mtrx [8][8]; {
  int i,j;
 
  for(i=7;i>(-1);i--) {
    for(j=7;j>(-1);j--)
      printf("%8.2f",mtrx[i][j]);
    printf("\n");
  }
}
/*****************************************************************************/
#include "proc/timestamp.proc"
#include "proc/set_mtrx.proc"
#include "proc/outrnd.proc"
#include "proc/round.proc"