#include <stdio.h>
#include <math.h>
#include <sys/time.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;
  float  cube;

  struct pi *child[3][2][8];         
  struct pi *left;                    
  struct pi *right;
  struct pi *parent;
};
typedef struct pi *list;



short levels=LEVELS,
      ivalue,
      jvalue,
      constrain=FALSE;

int   nnod[LEVELS];

long int comp=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(), 
	print_tree(),                   set_ca_mtrx(),
        set_cu_mtrx(),
	set_cx_mtrx(),                  unalloc(),
	set_q_mtrx(),
        graph_year(),                   graph_data(),
        get_pi_val(),                   close_files(),
        write_pi(),                   
        expand_tree(),                  data_head(),
        clock();

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

  int   alloc(),
        make_fortyeight(),
        rndmax(),
        insert_node(),
        pchoice, 
        qchoice;

  
  list  side[LEVELS+1], straighten(), head;
  short z, time[3];


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

  pick_mtrx(&pchoice, &qchoice);
  get_constants(&rnd, &alpha, &levels);

  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);
 
  clock(time);
  setbuf(stdout,NULL);
  head=(list) malloc(sizeof(struct pi));

  printf("\nLevel  1...\n");
  make_fortyeight(side, p, q, 0); 
  side[1]=straighten(head,side[1],0);
  printf("\n");

  for(z=1;z<levels;z++) {
    expand_tree(side, p, q, z);
    side[z+1]=straighten(head,side[z+1],0);
    printf("\n");
    printf("(nodes:%04d, ",nnod[z]);
    printf("ratio:%.2f).\n",(((float)nnod[z])/(float)nnod[z-1]));
  }
  init_Jmin(side[levels]);
  determine_optimum(side, ca, cu, cx, p, qchoice);  

  sort_Jmin(side);
  
  open_files();

  data_head(alpha, rnd, pchoice, qchoice, time, side);
  graph_year(pchoice, qchoice);
  for(z=0;z<levels;z++) {
    make_sens (side[0], z, pchoice, qchoice);
    print_tree (side[z], z);
  }
  close_files();
  printf("\nRun Time: %d:%02d:%02d",time[0],time[1],time[2]);
  printf("                                         Jmin:%f.\n\n",side[0]->Jmin);
  printf("comparisons:%d\n",comp);
}
/****************************************************************************/
void expand_tree(side, p, q, z)    
list side[];                  
float p[][8][8], q[][8][8]; 
short z; {
  int outofmem=FALSE,dots=0;
  list tempside;
  float space,count=1.0;

  space=40.0/(float)nnod[z-1];
  printf("Level %2d", z+1);
  nnod[z]=0;
  tempside=side[z];
  while(side[z]!=NULL && outofmem==FALSE) {
    outofmem=make_fortyeight(side, p, q, z);     
    side[z]=side[z]->left;
    count-=space;
    while(count<=0) {
      count+=1.0;
      printf(".");
      dots++;
    }
  }
  side[z]=tempside;
  if (outofmem==TRUE)
    unalloc(tempside);
  if(dots<40) 
    printf(".");
}
/***************************************************************************/
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, combined;
  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(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]));
          side[z]->denom[i][j][l]=0.0;
        }
        else{
          side[z]->child[i][j][l]->cube=
                                  round(side[z]->child[i][j][l]->value ,rnd);
          side[z]->denom[i][j][l]=dummy2;
          combined=FALSE;
          if(side[z+1]!=NULL)
            combined=insert_node(side,i,j,l,z);
          else {
            side[z+1]=side[z]->child[i][j][l];
            side[z+1]->left=side[z+1]->right=NULL;
          }
          if(!combined) {
            nnod[z]++;
            side[z]->child[i][j][l]->parent=side[z];
            side[z]->child[i][j][l]->prt[0]=i;
            side[z]->child[i][j][l]->prt[1]=j;
            side[z]->child[i][j][l]->prt[2]=l;             
            side[z]->child[i][j][l]->optflag=FALSE;
            side[z]->child[i][j][l]->prob=side[z]->denom[i][j][l]*side[z]->prob;
            side[z]->child[i][j][l]->left=NULL;
            side[z]->child[i][j][l]->right=NULL;
          }
        }
      }
  return(FALSE);
}
/*****************************************************************************/
int insert_node(side,i,j,l,z) 
list side[];
short i,j,l,z; {
  list node,kid;

  kid=side[z]->child[i][j][l];
  node=side[z+1];
  while(node!=NULL) {
      comp++;
    if(node->cube==kid->cube)
      if((compare_nodes(node,kid))<(rnd-0.001)) {
        free((char *)(side[z]->child[i][j][l]));
        side[z]->child[i][j][l]=node;
        return(TRUE);
    }
      comp++;
    if(kid->cube < node->cube) {
      if(node->left==NULL) {
        node->left=kid;
        return(FALSE);
      }
      node=node->left;
    }
    else {
      if(node->right==NULL) {
         node->right=kid;
         return(FALSE);
       }
      node=node->right;
    }
  }
  return(FALSE);
}

/*****************************************************************************/
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]));
    comp++;
    if (diff>(rnd-0.001)) return(diff);
    }
  return(diff);
}
/*****************************************************************************/
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;

    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;
    } }
}    
/*****************************************************************************/
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;
  }
}
/*****************************************************************************/
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("\nEnter choices of matrixes...\n\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;
  side[0]->left=side[0]->right=NULL;
}

/*****************************************************************************/
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");
    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 alpha (default is .9524):");
  scanf("%f",alpha);
  if (*alpha==0.0)
    *alpha=DALPHA;
  printf("\nEnter 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;
}
/*****************************************************************************/
list straighten(head,tree,depth)
list head, tree; 
int depth; {
  list templ,tempr;
  static int max=0;

  depth++;
  templ=tree->left;
  tempr=tree->right;

  if(templ!=NULL)
    head=straighten(head,templ,depth);
  if(tempr!=NULL)
    head=straighten(head,tempr,depth);
  if (head!=NULL)
    head->right=tree;
  tree->left=head;
  if(depth>max) {
    printf("%d ",depth);
    max=depth;
  }
  depth--;
  return(tree);
}


#include "proc/clock.proc"
#include "proc/set_mtrx.proc"
#include "proc/outfiles.proc"
#include "proc/round.proc"