/*
   File:        hype.c
   Author:      Andrew W. Moore
   Created:     Sat Sep 12 14:41:19 EDT 1992
   Description: Hyperrectangles and operations

   Copyright (C) 1992, Andrew W. Moore
*/

#include <stdio.h>
#include <math.h>
#include "ambs.h"      /* Very basic operations */
#include "amma.h"      /* Fast, non-fragmenting, memory management */
#include "amar.h"      /* Obvious operations on 1-d arrays */
#include "maxdim.h"    /* The MAX_DIM declaration */
#include "gpro.h"      /* Graphics projections kd->2d space */
#include "hype.h"      /* Hyper-rectangles from ../kdtr */


void set_limit_infinite(hlim,size)
hlimit *hlim;
int size;
{
  set_floats_constant(hlim->values,size,-77.77);
  set_bools_constant(hlim->finite,size,FALSE);
}

void set_hype_infinite(hy)
hype *hy;
/* hype size has already been set up, and its internal arrays allocated.
   Makes an infinite-in-all-directions hype
*/
{
  set_limit_infinite(&hy->lo,hy->size);
  set_limit_infinite(&hy->hi,hy->size);
}

void init_limit(hlim,size)
hlimit *hlim;
int size;
{
  hlim -> values = AM_MALLOC_ARRAY(float,size);
  hlim -> finite = AM_MALLOC_ARRAY(bool,size);
}

hype *create_hype(size)
int size;
/* Makes an infinite-in-all-directions hype */
/*
   mallocs the memory for the internal arrays
*/
{
  hype *res = AM_MALLOC(hype);
  res -> size = size;
  init_limit(&res->lo,size);
  init_limit(&res->hi,size);
  set_hype_infinite(res);
  return(res);
}

void free_limit(hlim,size)
hlimit *hlim;
int size;
{
  am_free((char *) hlim -> values,size * sizeof(float));
  am_free((char *) hlim -> finite,size * sizeof(bool));
}

void free_hype(hy)
hype *hy;
/*
   Frees the internal arrays
*/
{
  free_limit(&hy->lo,hy->size);
  free_limit(&hy->hi,hy->size);
  am_free((char *) hy,sizeof(hype));
}

void check_hype(hy,i,opname)
hype *hy;
int i;
char *opname;
{
  if ( hy == NULL )
    my_error("check_hype NULL");
  else if ( i < 0 || i >= hy -> size )
  {
    printf("%s(hy,%d) when hy->size = %d\n",opname,i,hy->size);
    my_error("hype index error");
  }
}

hlimit *hlimit_of(hy,limit)
hype *hy;
int limit;
{
  if ( hy == NULL || (limit != LO && limit != HI) )
    my_error("hype NULL or limit not LO or HI in hlimit_of()");
  else
    return((limit == LO) ? &hy->lo : &hy->hi);
}

void hype_update(hy,comp,limit,value)
hype *hy;
int comp;
int limit;
float value;
{
  hlimit *lim = hlimit_of(hy,limit);
  check_hype(hy,comp,"hype_update");
  if ( (limit == LO && hy->hi.finite[comp] && hy->hi.values[comp] < value) ||
       (limit == HI && hy->lo.finite[comp] && hy->lo.values[comp] > value)
     )
  {
    fprintf(stderr,"*** Can't turn hyper-rect inside out\n");
    fprintf(stderr,"comp = %d\n",comp);
    fprintf(stderr,"limit = %d\n",limit);
    fprintf(stderr,"value = %g\n",value);
    fprintf_hype(stderr,"hy",hy,"\n");
    my_error("hype_update()");
  }

  lim->values[comp] = value;
  lim->finite[comp] = TRUE;
}

void hype_update_infinite(hy,comp,limit)
hype *hy;
int comp;
int limit;
{
  hlimit *lim = hlimit_of(hy,limit);
  check_hype(hy,comp,"hype_update_infinite");
  lim->values[comp] = -88.88;
  lim->finite[comp] = FALSE;
}

bool hype_has_limit(hy,comp,limit)
hype *hy;
int comp;
int limit;
{
  hlimit *lim = hlimit_of(hy,limit);
  check_hype(hy,comp,"hype_update_infinite");
  return(lim->finite[comp]);
}

float hype_ref(hy,comp,limit)
hype *hy;
int comp;
int limit;
{
  hlimit *lim = hlimit_of(hy,limit);
  check_hype(hy,comp,"hype_update_infinite");
  if ( !lim->finite[comp] )
  {
    fprintf_hype(stderr,"hype = ",hy,"\n");
    printf("comp = %d, limit = %d, illegal\n",comp,limit);
    my_error("hype_ref");
  }
  else
    return(lim->values[comp]);
}

void fprintf_hype_comp(s,hy,i,limit)
FILE *s;
hype *hy;
int i;
int limit;
{
  if ( hype_has_limit(hy,i,limit) )
    fprintf(s,"%g",hype_ref(hy,i,limit));
  else
    fprintf(s,"%sinfinity",(limit == LO) ? "-" : "+");
}

void fprintf_hype(s,m1,hy,m2)
FILE *s;
char *m1;
hype *hy;
char *m2;
{
  int i;
  fprintf(s,"%s (",m1);
  for ( i = 0 ; i < hy -> size ; i++ )
  {
    fprintf(s,"[");
    fprintf_hype_comp(s,hy,i,LO);
    fprintf(s,",");
    fprintf_hype_comp(s,hy,i,HI);
    fprintf(s,"]%s",(i==hy->size-1) ? "" : " ");
  }
  fprintf(s,") %s",m2);
}

bool is_inside_hype(hy,farr)
hype *hy;
float *farr;
{
  bool result = TRUE;
  int i;
  for ( i = 0 ; result && i < hy->size ; i++ )
    result = (!hy->lo.finite[i] || hy->lo.values[i] <= farr[i]) &&
             (!hy->hi.finite[i] || hy->hi.values[i] >= farr[i]);
  return(result);
}

bool hype_high_limit_is_above(hy,comp,val)
hype *hy;
int comp;
float val;
/*
    returns TRUE iff the HI value of the hyperrectangle in component "comp"
    is strictly above "val".
*/
{
  return( !hy->hi.finite[comp] || hy->hi.values[comp] > val );
}

bool hype_low_limit_is_below(hy,comp,val)
hype *hy;
int comp;
float val;
/*
    returns TRUE iff the LO value of the hyperrectangle in component "comp"
    is strictly below "val".
*/
{
  return( !hy->lo.finite[comp] || hy->lo.values[comp] < val );
}

bool comp_inside_hype(hy,comp,val)
hype *hy;
int comp;
float val;
/*
   returns TRUE iff val is above or equal to the LO value of hy and
                    val is below or equal to the HI value of hy.

   Note
     comp_inside_hype(hy,comp,val) =
        (hy->lo[comp] == val || hype_low_limit_is_below(hy,comp,val)) &&
        (hy->hi[comp] == val || hype_high_limit_is_above(hy,comp,val)) &&
   is always true.
*/
{
  return( ( !hy->hi.finite[comp] || hy->hi.values[comp] >= val ) &&
          ( !hy->lo.finite[comp] || hy->lo.values[comp] <= val )
        );
}

float clipped_below_value(hy,comp,val)
hype *hy;
int comp;
float val;
/*
   Is the "comp"th component of "hy" below "val"?
      If it is, we return "val"
      Otherwise we return hy[comp]
*/
{
  if ( !hy->lo.finite[comp] )
    return(val);
  else
    return(real_max(hy->lo.values[comp],val));
}

float clipped_above_value(hy,comp,val)
hype *hy;
int comp;
float val;
/*
   Is the "comp"th component of "hy" above "val"?
      If it is, we return "val"
      Otherwise we return hy[comp]
*/
{
  if ( !hy->hi.finite[comp] )
    return(val);
  else
    return(real_min(hy->hi.values[comp],val));
}

void gp_draw_maybe_filled_hype(gp,hy,is_filled)
gproj *gp;
hype *hy;
bool is_filled;
{
  int i;
  bool draw_me = TRUE;
  for ( i = 0 ; draw_me && i < hy->size ; i++ )
    draw_me = ( i == gp->x.comp ||
                i == gp->y.comp ||
                comp_inside_hype(hy,i,gp->base[i])
              );

  if ( draw_me )
  {
    float xlo = clipped_below_value(hy,gp->x.comp,gp->x.lo);
    float ylo = clipped_below_value(hy,gp->y.comp,gp->y.lo);
    float xhi = clipped_above_value(hy,gp->x.comp,gp->x.hi);
    float yhi = clipped_above_value(hy,gp->y.comp,gp->y.hi);
    if ( is_filled )
      gp_draw_box(gp,xlo,ylo,xhi,yhi);
    else
      gp_draw_rectangle(gp,xlo,ylo,xhi,yhi);
  }
}

void gp_draw_hype(gp,hy)
gproj *gp;
hype *hy;
{
  gp_draw_maybe_filled_hype(gp,hy,TRUE);
}

void gp_draw_hype_border(gp,hy)
gproj *gp;
hype *hy;
{
  gp_draw_maybe_filled_hype(gp,hy,FALSE);
}

void m_and_s_from_hype(hy,middle,scales)
hype *hy;
float *middle,*scales;
/*
  hype must be completely finite.
  middle and scales must have at least hy->size array-cells allocated.
*/
{
  int i;
  for ( i = 0 ; i < hy -> size ; i++ )
  {
    float lo = hype_ref(hy,i,LO);
    float hi = hype_ref(hy,i,HI);
    middle[i] = (lo + hi)/2.0;
    scales[i] = (hi - lo)/2.0;
  }
}

void shrink_hype_comp(hy,factor,middle,comp,limit)
hype *hy;
float factor;
float *middle;
int comp,limit;
{
  if ( hype_has_limit(hy,comp,limit) )
  {
    hype_update(hy,comp,limit,factor * hype_ref(hy,comp,limit) +
                              (1.0 - factor) * middle[comp]
               );
  }
}

void shrink_hype(hy,factor,middle)
hype *hy;
float factor;
float *middle;
/*
   Alters hy to shrink the distance from middle to each hype
   extreme by a factor of factor. (IE 0.5 halves it).

   An infinite component is unaffected.
*/
{
  int i;
  for ( i = 0 ; i < hy -> size ; i++ )
  {
    shrink_hype_comp(hy,factor,middle,i,LO);
    shrink_hype_comp(hy,factor,middle,i,HI);
  }
}

void clip_hype_with_m_and_s(hy,middle,scales)
hype *hy;
float *middle,*scales;
{
  int i;
  for ( i = 0 ; i < hy -> size ; i++ )
  {
    if (!hype_has_limit(hy,i,LO) || hype_ref(hy,i,LO) < middle[i] - scales[i])
      hype_update(hy,i,LO,middle[i]-scales[i]);
    
    if (!hype_has_limit(hy,i,HI) || hype_ref(hy,i,HI) > middle[i] + scales[i])
      hype_update(hy,i,HI,middle[i]+scales[i]);
  }
}

float hype_volume(hy,scales)
hype *hy;
float *scales;
/*
    Returns normalized volume of hype: normalized by scales.
*/
{
  int i;
  float result = 1.0;
  
  for ( i = 0 ; i < hy -> size ; i++ )
    result *= (hype_ref(hy,i,HI) - hype_ref(hy,i,LO)) / scales[i];

  return(result);
}

void copy_hlimit(h1,h2,size)
hlimit *h1,*h2;
int size;
{
  copy_floats(h1->values,h2->values,size);
  copy_bools(h1->finite,h2->finite,size);
}

void copy_hype(h1,h2)
hype *h1,*h2;
{
  if ( h1 -> size != h2 -> size )
    my_error("copy_hype().. the bloody sizes are different");
  else
  {
    copy_hlimit(&h1->lo,&h2->lo,h1->size);
    copy_hlimit(&h1->hi,&h2->hi,h1->size);
  }
}

int longest_scaled_comp(hy,scales)
hype *hy;
float *scales;
/*
   Find longest dimension relative to scales (i.e. if
     width[i] / width[j] = scales[i] / scales[j] its a draw).
*/
{
  int result = -1;
  float longest_sw = 0.0;
  int i;

  for ( i = 0 ; i < hy->size ; i++ )
  {
    float sw = ( hype_ref(hy,i,HI) - hype_ref(hy,i,LO) ) / scales[i];
    if ( result < 0 || sw > longest_sw )
    {
      result = i;
      longest_sw = sw;
    }
  }

  return(result);
}