/*
  Returns a list of local extrema of the vector, organized thus:
  
   [ [:min, idmin1], [:max, idmax1], ...]
  
  The values are pushed in the order in which they are found. It works
  thus: it scans the vector and looks around the current point in a
  given window. If the current point is the maximum or the minimum, it
  is considered as a local maximum/minimum. Control over which extrema
  are included is given to the used through threshold mechanisms.

  The _options_ hash controls how the peaks are detected:
  * _window_: the number of elements on which we look on
    both sides (default 5, ie the local maximum is over 11 points)
  * _threshold_: the minimum amplitude the extrema must have to
    be considered (default 0)
  * _dthreshold_: how much over/under the average an extremum must be
    (default 0) 
  * _or_: whether the _threshold_ and _dthreshold_ tests are both
    necessary or if only one is (default false: both tests are
    necessary)

    *Note:* beware of NANs ! They *will* screw up peak detection, as
  they are neither bigger nor smaller than anything...  
*/
static VALUE dvector_extrema(int argc, VALUE *argv, VALUE self)
{
  long window = 5;
  double threshold = 0;
  double dthreshold = 0;
  int inclusive = 1;
  
  if(argc == 1) {
    VALUE t;
    t = rb_hash_aref(argv[0], rb_str_new2("window"));
    if(RTEST(t)) {
      window = FIX2LONG(t);
    }
    t = rb_hash_aref(argv[0], rb_str_new2("threshold"));
    if(RTEST(t)) {
      threshold = rb_num2dbl(t);
    }
    t = rb_hash_aref(argv[0], rb_str_new2("dthreshold"));
    if(RTEST(t)) {
      dthreshold = rb_num2dbl(t);
    }
    
    t = rb_hash_aref(argv[0], rb_str_new2("or"));
    inclusive = ! RTEST(t);
  } else if(argc > 1)
    rb_raise(rb_eArgError, "Dvector.extrema only takes 0 or 1 argument");

  /* Handling of the vector */
  long len, i,j;
  double * data = Dvector_Data_for_Read(self, &len);
  VALUE s_min = ID2SYM(rb_intern("min"));
  VALUE s_max = ID2SYM(rb_intern("max"));

  

  VALUE ret = rb_ary_new();
                       
  for(i = 0; i < len; i++) {

    /* This is stupid and will need decent optimization when I have
       time */
    long first = i > window ? i - window : 0;
    double cur_min = data[first];
    long cur_min_idx = first;
    double cur_max = data[first];
    long cur_max_idx = first;
    double average = 0;
    long nb = 0;
    
    for(j = first; (j < i+window) && (j < len); j++,nb++) {
      average += data[j];
      if(data[j] <= cur_min) {
        cur_min = data[j];
        cur_min_idx = j;
      }
      if(data[j] >= cur_max) {
        cur_max = data[j];
        cur_max_idx = j;
      }
    }
    average /= nb;

    if(cur_min_idx == i) {
      /* This is a potential minimum */
      if((inclusive && 
          (fabs(cur_min) >= threshold) && 
          (fabs(cur_min - average) >= dthreshold))
         || (!inclusive && 
             ((fabs(cur_min) >= threshold) ||
              (fabs(cur_min - average) >= dthreshold))
             )) {
        VALUE min = rb_ary_new();
        rb_ary_push(min, s_min);
        rb_ary_push(min, LONG2FIX(i));
        rb_ary_push(ret, min);
      }
    }
    else if(cur_max_idx == i) {
      /* A potential maximum */
      if((inclusive && 
          (fabs(cur_max) >= threshold) && 
          (fabs(cur_max - average) >= dthreshold))
         || (!inclusive && 
             ((fabs(cur_max) >= threshold) ||
              (fabs(cur_max - average) >= dthreshold))
             )) {
        VALUE max = rb_ary_new();
        rb_ary_push(max, s_max);
        rb_ary_push(max, LONG2FIX(i));
        rb_ary_push(ret, max);
      }
    }
  }
  return ret;
}