// multires peak cache object

#include <stdlib.h>
#include <math.h>
#include "debug.h"
#include "macros.h"
#include "mrpc.h"




/* get nb of pairs at level */
int mrpc::nb_pairs(int level){
    ASSERT (level >= this->top_level);    // not represented (to fine)
    ASSERT (level <= this->bottom_level); // .. (to coarse)
    return (1 << (this->bottom_level - level));
}

/* get first cache block at level */
pair_t * mrpc::level(int level){
    ASSERT (level >= this->top_level);    // not represented (to fine)
    ASSERT (level <= this->bottom_level); // .. (to coarse)
    int index = this->total_nb_pairs - (this->total_nb_pairs >> (level - this->top_level));
    ASSERT (index >= 0);
    ASSERT (index < this->total_nb_pairs);
    return this->multi + index;
}

/* utility functions for multires stuff */

int mrpc::in_vector(float *f){
  int i = f - this->vec;
  return ((i >= 0) && (i < this->elements));
}

/* scan a part of the vector */
void mrpc::scan_bound(float *f, int elements, int stride, pair_t *pair){
    pair->min = 0.0f;
    pair->max = 0.0f;
    if (!elements) return;

    if (this->in_vector(f)){
        pair->min = f[0];
        pair->max = f[0];
        elements--;
        f += stride;
    }

    while (elements--){
        float x = this->in_vector(f) ? f[0] : 0.0f;
        pair->min = dmin(pair->min, x);
        pair->max = dmax(pair->max, x);
        f += stride;
    }
}


/* compute all lower levels from top level */
void mrpc::refine(int block){
    int l;
    for (l = this->top_level; l < this->bottom_level; l++){
        int src_pairs = this->nb_pairs(l);
        int dst_pairs = this->nb_pairs(l+1);
        pair_t *src_pair = this->level(l)   + block * src_pairs;
        pair_t *dst_pair = this->level(l+1) + block * dst_pairs;
        
        while (dst_pairs--){
            dst_pair->min = dmin(src_pair[0].min, src_pair[1].min);
            dst_pair->max = dmax(src_pair[0].max, src_pair[1].max);
            dst_pair++;
            src_pair+=2;
        }
    }
}


/* build top level 
   the only sour point here is that the last block does
   not have proper size, so we workaround it by NOT computing it until it works (FXME)
*/
void mrpc::build_top(int block){
    int point = block * this->samples_per_block;
    ASSERT(point >= 0);
    
    int i = this->nb_top_pairs;
    int chunk = 1 << this->top_level; // sample chunk size that will be used for top level pair
    float *vec = this->vec + point;
    pair_t *p = this->level(this->top_level) + block * this->nb_top_pairs;

    while (i--){
        this->scan_bound(vec, chunk, 1, p);
        vec += chunk;
        p++;
    }
}


// update cache for given range
void mrpc::update(int start, int endx){
  if (endx < start) {return;}

  int block_start = start / this->samples_per_block;
  int block_endx = ((endx-1) / this->samples_per_block) + 1;

  this->build(clip(0, block_start, this->blocks),
	      clip(0, block_endx, this->blocks));
  

}

void mrpc::build(int block_start, int block_endx){
    int i;
    /* for each block */
    for(i=block_start; i<block_endx; i++){
        
        /* generate top level */
        this->build_top(i);

        /* refine */
        this->refine(i);
    }
}


mrpc::mrpc(float *vec, int elements, int top_level, int bottom_level){

    ASSERT(vec);
    ASSERT(elements);

    /* TODO determine bottom level from array size */

    /* init state */
    this->vec = vec;
    this->elements = elements;
    this->top_level = top_level;
    this->bottom_level = bottom_level;
    this->nb_top_pairs = 1 << (bottom_level - top_level);
    this->samples_per_block = 1 << bottom_level;
    this->blocks = ((elements - 1) / this->samples_per_block) + 1;
    this->total_nb_pairs = this->blocks * this->nb_top_pairs * 2;  // 2 = 1 + 1/2 + 1/4 + ...
    this->multi = (pair_t *)calloc(this->total_nb_pairs, sizeof(pair_t));


    /* build cache */
    this->build(0, this->blocks);


}

mrpc::~mrpc(){
  free(this->multi);
}

/* note: floats are not deep enough for using incremental updates
   recomputing everything from ints would be a better approach,
   but it is slower. not that it matters that much though.
 */
void mrpc::peaks_direct(pair_t *p, int nb_pairs, 
			double fpoint, double fsamples_per_pixel){

    int samples_per_pixel = (int)fsamples_per_pixel;
    while (nb_pairs--){
	int point = (int)fpoint;
	this->scan_bound(this->vec + point, samples_per_pixel, 1, p++); // get pair
	fpoint += fsamples_per_pixel; // move point
    }
}

/* note: if fpairs_per_pixel >= 1.0f subsequent steps do not introduce aliasing
   instead of respecting this, we allow aliasing to occur for very small window size
   and full view. */
void mrpc::peaks_multires(pair_t *p, int nb_pairs, double fpoint, 
                         double fsamples_per_pixel, int level){


    /* prepare for traversing cache level */
    double subsamp =  ((double)(1 << level));     // subsampling factor for this level
    double fpairs_per_pixel = fsamples_per_pixel / subsamp; // get pair increment

    double fpair = fpoint / subsamp; // get pair point
    int level_start = this->level(level) - this->multi; // start index of this level's pairs
    int level_endx = level_start + this->blocks * this->nb_pairs(level); // end index
    
    /* fill peak array */
    while (nb_pairs--){
	int pair = ((int)fpair) + level_start;
	if ((pair >= level_start) && (pair < level_endx)){
	    p->min = this->multi[pair].min;
	    p->max = this->multi[pair].max;
	}
	else {
	    // no data == 0.0
	    p->min = 0.0f;
	    p->max = 0.0f;
	}
	p++;
	fpair += fpairs_per_pixel;
    }
}


/* fill an array of peaks */
void mrpc::peaks(int start, int endx, pair_t *p, int nb_pairs)
{

  //fprintf(stderr, "mrpc::peaks start=%d endx=%d\n", start, endx);

    ASSERT(start <= endx);
    ASSERT(nb_pairs > 0);

    double fsamples_per_pixel = dmax(1.0f, ((double)(endx - start)) / ((double)nb_pairs));

    // maybe this is not necessary. however, the entire widget respects this
    ASSERT(fsamples_per_pixel >= 1.0); 

    double flevel = log(fsamples_per_pixel) / log(2);
    double fpoint = start;

    /* determine log samplingfactor */
    int level = (int)ceil(flevel);               // round up to lessen aliasing

    /* use cache if zoomlevel coarse enough */
    if (level >= this->top_level + 1){
      pair_t p2[nb_pairs];
      this->peaks_multires(p, nb_pairs, fpoint, fsamples_per_pixel, level); // get hires array
      this->peaks_multires(p2, nb_pairs, fpoint, fsamples_per_pixel, level+1); // get lowres array
      
      // (0->lowres version, 1->hires version)
      float res = ((double)level - flevel);
      int i;
      /*  interpolate hires and lowres array (mipmap)
	  the result is roughly speaking that the corner frequency of the envelope follower is at nyquist/2
	  
	  the transition between cache and direct is noticable becase direct has a nyqyust envelope follower,
	  (1 pixel wide) to ensure that full zoom (sample = pixel) is accurate. at high zoom levels aliasing
	  is less important because the data is usually of low frequency nature.

       */
      for (i=0; i<nb_pairs; i++){
	p[i].min *= res;
	p[i].min += (1.0f - res) * p2[i].min;
	p[i].max *= res;
	p[i].max += (1.0f - res) * p2[i].max;
      }
    }
    /* dont interpolate, just use coarsest level */
    else if (level >= this->top_level){
      this->peaks_multires(p, nb_pairs, fpoint, fsamples_per_pixel, level); // get hires array
    }

    /* if zoomlevel too fine, compute directly from array */
    else {
      this->peaks_direct(p, nb_pairs, fpoint, fsamples_per_pixel);
    }




}