/* This was an attempt at a depth-first search, but runs out of memory.  Should be 
   fixed by making it non-recursive.  */
#include <stdio.h>
#include <math.h>
#include <sys/time.h>

#define TRUE 1
#define FALSE 0
#define EMPTY (-1)

#define LEVELS 20
#define DLEVELS 10
#define DRND 0.1
#define DALPHA 0.9524

struct   pi {
  float  Jmin,
         denom[3][2][8],
         value[8],
         cube;

  short  child[3][2][8];
};  
         
typedef struct pi *list;             
             
short levels=LEVELS,
      mllc=0, pisize=0;

float alpha=DALPHA,
      rnd=DRND;

main() {
     
  list node;

  void  get_constants(),     
        determine_optimum(),
	init_to_zero(),                 pick_mtrx(),
	set_p_mtrx(),                   
        set_cu_mtrx(),
	set_ca_mtrx(),
	set_cx_mtrx(),	                set_q_mtrx(),                   
        get_pi_val(),                   init_tree(),
        make_fortyeight(),
        init_node(),                    init_block();

  float assign_value(),
        round(),
        p[3][8][8],
        q[2][8][8],
        ca[3][8], 
        cu[8],                       
        cx[2][7];                    
             
  int   rndmax(),        pchoice, qchoice;
  
  pisize=sizeof(struct pi); 

  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);       
            
  node=(list)malloc(pisize);             
  init_tree(node);                  
                      
  make_fortyeight(node, p, q, 0, ca, cu, cx, qchoice);

  printf("Jmin: %f\n",node->Jmin);
}            
/********************************************************************************/
void make_fortyeight(dad, p, q, z, ca, cu, cx, qchoice)
int qchoice;
list dad;
float p[][8][8], q[][8][8], ca[][8], cu[], cx[2][7];
short z; {   
  short i, j, k=0, l, next=0, merged;
  int kid=0;
  list node[48];

  for(k=0;k<48;k++) {
    node[k]=(list)malloc(pisize);
    mllc++;
    if(node[k]==NULL)
      printf("error in allocation.\n");
  }

  for(i=0;i<3;i++) 
    for(j=0;j<2;j++) { 
      for(l=0;l<8;l++) {
        node[kid]->cube=assign_value(dad, node[kid], i, j, l, p, q);
        if(node[kid]->cube==0.0)
          dad->child[i][j][l]=EMPTY;
        else {
          dad->child[i][j][l]=kid;
          init_node(node[kid]);
          merged=FALSE;
          if(z<(levels-1)) {
            k=0;
            while((k<kid)&&(merged==FALSE)) {
              if(node[k]->cube==node[kid]->cube) {
                node[kid]->Jmin=node[k]->Jmin;
                merged=TRUE;
              }
              k++;
            }
            if(!merged)
              make_fortyeight(node[kid], p, q, z+1, ca, cu, cx, qchoice);
          }
          determine_optimum(node, dad, ca, cu, cx, p, qchoice);
          kid++;
        }
      }
    }
  for(k=0;k<48;k++) {
    free((char *) node[kid]);
    mllc--;
  }
  printf("Freed to:%d  ", mllc);
  printf("Jmin for lev %d: %05.2f\n",z, dad->Jmin);
}
/*****************************************************************************/
void determine_optimum(node, dad, ca, cu, cx, p, qchoice)        
int qchoice; 
float ca[][8], cu[], cx[2][7], p[][8][8]; 
list dad, node[48]; {                  
  float cofm, cofu, cofa, future=0.0, cost;
  int i, j, l, m, k;                 

  dad->Jmin=10000.0; 
  for(i=0;i<3;i++) {             
    cofa=cofu=0.0;               
    for(k=0;k<8;k++) {           
      cofa+=(dad->value[k]*ca[i][k]);               
      for(m=0;m<8;m++)           
        cofu+=(dad->value[k]*p[i][k][m]*cu[m]);     
    }    
    cofu*=alpha;                 
    for(j=0;j<2;j++) {        
      cofm=alpha*cx[j][qchoice];
      future=0.0;                
      cost=cofa+cofu+cofm;       
      for(l=0;l<8;l++) {         
        if(dad->child[i][j][l]!=EMPTY)              
	  future+=(dad->denom[i][j][l]*node[dad->child[i][j][l]]->Jmin);
      }  
      future=(alpha*future)+cost;
      if(future<dad->Jmin) 
        dad->Jmin=future;   
    }    
  } 
}
/*****************************************************************************/
void init_node(side)                 
list side; { 
  int i,j,l;
  
  for(i=0;i<3;i++)
    for(j=0;j<2;j++)
      for(l=0;l<8;l++) {
        side->child[i][j][l]=EMPTY;
        side->denom[i][j][l]=0.0;
      }
   side->Jmin=0.0;
}
/*****************************************************************************/
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(node) 
list node; {                  
  int k;

  printf("\nEnter intial pi values...\n\n");
  for(k=0; k<8; k++){
    printf("  Pi value[%d]=",k);
    scanf("%f,",&node->value[k]);
  }
  printf("\n");
}
/*****************************************************************************/
void init_tree(node)
list node; {
             
  init_node(node);                
             
  get_pi_val(node);                      
             
  node->cube=round(node->value,rnd);                   
}            
/*****************************************************************************/
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 %.4f):",DALPHA);
  scanf("%f",alpha);
  if (*alpha==0.0)
    *alpha=DALPHA;
  printf("  Enter value for levels (default is %d, maximum is %d):",
         DLEVELS,LEVELS);
  scanf("%d",levels);
  if (*levels==0)
    *levels=DLEVELS;
  printf("  Enter value for the interval (default is %.3f, minimum is 0.005):",
         DRND);
  scanf("%f",rnd);
  if (*rnd==0.0)
    *rnd=DRND;
}
/*****************************************************************************/
float assign_value(ndad, nkid,i,j,l,p,q)
list ndad, nkid;
int i, j, l; 
float p[][8][8],q[][8][8]; {
  int k, n, m, biggest, rndmax(), round[8];
  float sum=0.0, tmp, cube=0.0, delta[8],
        dummy, dummy1, dummy2, pi[8];

  dummy2=0.0;
  for(m=0;m<8;m++) {
    dummy=0.0;
    for(n=0;n<8;n++)
      dummy+=(p[i][n][m]*(ndad->value[n]));
    dummy2+=(q[j][m][l]*dummy);
  }

  if(dummy2==0.0) {
    ndad->denom[i][j][l]=0.0;
    return (0.0);
  }
  else {
    ndad->denom[i][j][l]=dummy2;
    for(k=0;k<8;k++) {
      dummy1=0.0;
      for(n=0;n<8;n++)
        dummy1+=p[i][n][k]*ndad->value[n];
      tmp=(q[j][k][l]*dummy1/dummy2);    
      pi[k]=((float)(rint(tmp/rnd))*rnd);
      delta[k]=fabs(pi[k]-tmp);          
      if (pi[k]>tmp)                     
        round[k]=1;
      else         
        round[k]=0;
      sum+=pi[k];  
    }              
    if (sum>(1.0+rnd-0.001)) {           
      while (sum>1.0) {                  
        biggest=rndmax(delta,round,1);   
        pi[biggest]-=rnd;                
        delta[biggest]=0.0;              
        sum-=rnd;  
      }            
    }              
    if (sum<(1.0-rnd+0.001)) {           
      while (sum<1.0) {                  
        biggest=rndmax(delta,round,0);
        pi[biggest]+=rnd;
        delta[biggest]=0.0;
        sum+=rnd;
      }
    }
    j=1;
    for(i=0;i<8;i++) {                   
      cube+=(j*j*j)*pi[i]*pi[i]*pi[i];   
      nkid->value[i]=pi[i];
      j++;
    }
    return(cube);
  }
}
/******************************************************************************/
float round(pi,rnd)
float pi[8], rnd; {
  int i,biggest, rndmax();                  
  register int j=1;
  static int round[8];
  float sum=0.0, tmp,cube=0.0;
  static float delta[8];
  
  for(i=0;i<8;i++) {
    tmp=(float)(rint(pi[i]/rnd)*rnd);
    delta[i]=fabs(pi[i]-tmp);
    if (pi[i]<tmp)
      round[i]=1;
    else 
      round[i]=0;
    pi[i]=tmp;
    sum+=pi[i];

  }
  if (sum>(1.0+rnd-0.001)) {
    while (sum>1.0) {
      biggest=rndmax(delta,round,1);
      pi[biggest]-=rnd;
      delta[biggest]=0.0;
      sum-=rnd;
    }
  }
  if (sum<(1.0-rnd+0.001)) {
    while (sum<1.0) {
      biggest=rndmax(delta,round,0);
      pi[biggest]+=rnd;
      delta[biggest]=0.0;
      sum+=rnd;
    }
  }
  for(i=0;i<8;i++) {
    cube+=(j*j*j)*pi[i]*pi[i]*pi[i];
    j++;
  }
  return(cube);
}
/*****************************************************************************/
int rndmax(delta, round, val)
float delta[];
int round[],val; {

  int i, choice=0;
  float diff=0.0;

  for(i=0;i<8;i++) {
    if(round[i]==val) {
      if(delta[i]>diff) {
        diff=delta[i];
        choice=i;
      }
    }
  }
  return(choice);
}
#include "set_mtrx.proc"