statistc.cpp

一ＯＣＲ的相关资料。．希望对研究ＯＣＲ的朋友有所帮助．

          if (count > new_mode) {
            new_mode = count;
            new_centre = entry + rangemin;
          }
        }
      }
    }
                                 //need new and room
    if (new_mode > 0 && cluster_count < max_clusters) {
      cluster_count++;
      new_cluster = TRUE;
      if (!clusters[cluster_count].set_range (rangemin, rangemax))
        return 0;
      centres[cluster_count] = (float) new_centre;
      clusters[cluster_count].add (new_centre, new_mode);
      clusters[0].add (new_centre, new_mode);
      for (entry = new_centre - 1; centres[cluster_count] - entry < lower
        && entry >= rangemin
      && pile_count (entry) <= pile_count (entry + 1); entry--) {
        count = pile_count (entry) - clusters[0].pile_count (entry);
        if (count > 0) {
          clusters[cluster_count].add (entry, count);
          clusters[0].add (entry, count);
        }
      }
      for (entry = new_centre + 1; entry - centres[cluster_count] < lower
        && entry < rangemax
      && pile_count (entry) <= pile_count (entry - 1); entry++) {
        count = pile_count (entry) - clusters[0].pile_count (entry);
        if (count > 0) {
          clusters[cluster_count].add (entry, count);
          clusters[0].add (entry, count);
        }
      }
      centres[cluster_count] =
        (float) clusters[cluster_count].ile ((float) 0.5);
    }
  }
  while (new_cluster && cluster_count < max_clusters);
  free_mem(centres); 
  return cluster_count;
}


/**********************************************************************
 * STATS::local_min
 *
 * Return TRUE if this point is a local min.
 **********************************************************************/

BOOL8 STATS::local_min(         //test minness
                       INT32 x  //of x
                      ) {
  INT32 index;                   //table index

  if (buckets == NULL) {
    /*		err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
            ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
            "Empty stats");*/
    return FALSE;
  }
  if (x < rangemin)
    x = rangemin;
  if (x >= rangemax)
    x = rangemax - 1;
  x -= rangemin;
  if (buckets[x] == 0)
    return TRUE;
  for (index = x - 1; index >= 0 && buckets[index] == buckets[x]; index--);
  if (index >= 0 && buckets[index] < buckets[x])
    return FALSE;
  for (index = x + 1; index < rangemax - rangemin
    && buckets[index] == buckets[x]; index++);
  if (index < rangemax - rangemin && buckets[index] < buckets[x])
    return FALSE;
  else
    return TRUE;
}


/**********************************************************************
 * STATS::print
 *
 * Print a summary of the stats and optionally a dump of the table.
 **********************************************************************/

void STATS::print(            //print stats table
                  FILE *,     //Now uses tprintf instead
                  BOOL8 dump  //dump full table
                 ) {
  INT32 index;                   //table index

  if (buckets == NULL) {
    /*     err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
       ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
       "Empty stats"); */
    return;
  }
  if (dump) {
    for (index = 0; index < rangemax - rangemin; index++) {
      tprintf ("%4d:%-3d ", rangemin + index, buckets[index]);
      if (index % 8 == 7)
        tprintf ("\n");
    }
    tprintf ("\n");
  }

  tprintf ("Total count=%d\n", total_count);
  tprintf ("Min=%d\n", (INT32) (ile ((float) 0.0)));
  tprintf ("Lower quartile=%.2f\n", ile ((float) 0.25));
  tprintf ("Median=%.2f\n", ile ((float) 0.5));
  tprintf ("Upper quartile=%.2f\n", ile ((float) 0.75));
  tprintf ("Max=%d\n", (INT32) (ile ((float) 0.99999)));
  tprintf ("Mean= %.2f\n", mean ());
  tprintf ("SD= %.2f\n", sd ());
}


/**********************************************************************
 * STATS::min_bucket
 *
 * Find REAL minimum bucket - ile(0.0) isnt necessarily correct
 **********************************************************************/

INT32 STATS::min_bucket() {  //Find min
  INT32 min;

  if (buckets == NULL) {
    /*		err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
            ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
            "Empty stats");*/
    return rangemin;
  }

  for (min = 0; (min < rangemax - rangemin) && (buckets[min] == 0); min++);
  return rangemin + min;
}


/**********************************************************************
 * STATS::max_bucket
 *
 * Find REAL maximum bucket - ile(1.0) isnt necessarily correct
 **********************************************************************/

INT32 STATS::max_bucket() {  //Find max
  INT32 max;

  if (buckets == NULL) {
    /*		err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
            ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
            "Empty stats");*/
    return rangemin;
  }

  for (max = rangemax - rangemin - 1;
    (max > 0) && (buckets[max] == 0); max--);
  return rangemin + max;
}


/**********************************************************************
 * STATS::short_print
 *
 * Print a summary of the stats and optionally a dump of the table.
 * ( BUT ONLY THE PART OF THE TABLE BETWEEN MIN AND MAX)
 **********************************************************************/

void STATS::short_print(            //print stats table
                        FILE *,     //Now uses tprintf instead
                        BOOL8 dump  //dump full table
                       ) {
  INT32 index;                   //table index
  INT32 min = min_bucket ();
  INT32 max = max_bucket ();

  if (buckets == NULL) {
    /*     err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
       ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
       "Empty stats"); */
    return;
  }
  if (dump) {
    for (index = min; index <= max; index++) {
      tprintf ("%4d:%-3d ", rangemin + index, buckets[index]);
      if ((index - min) % 8 == 7)
        tprintf ("\n");
    }
    tprintf ("\n");
  }

  tprintf ("Total count=%d\n", total_count);
  tprintf ("Min=%d Really=%d\n", (INT32) (ile ((float) 0.0)), min);
  tprintf ("Max=%d Really=%d\n", (INT32) (ile ((float) 1.1)), max);
  tprintf ("Range=%d\n", max + 1 - min);
  tprintf ("Lower quartile=%.2f\n", ile ((float) 0.25));
  tprintf ("Median=%.2f\n", ile ((float) 0.5));
  tprintf ("Upper quartile=%.2f\n", ile ((float) 0.75));
  tprintf ("Mean= %.2f\n", mean ());
  tprintf ("SD= %.2f\n", sd ());
}


/**********************************************************************
 * STATS::plot
 *
 * Draw a histogram of the stats table.
 **********************************************************************/

#ifndef GRAPHICS_DISABLED
void STATS::plot(                //plot stats table
                 WINDOW window,  //to draw in
                 float xorigin,  //bottom left
                 float yorigin,
                 float xscale,   //one x unit
                 float yscale,   //one y unit
                 COLOUR colour   //colour to draw in
                ) {
  INT32 index;                   //table index

  if (buckets == NULL) {
    /*		err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
            ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
            "Empty stats");*/
    return;
  }
  interior_style (window, INT_HOLLOW, 1);
  perimeter_color_index(window, colour); 

  for (index = 0; index < rangemax - rangemin; index++) {
    rectangle (window, xorigin + xscale * index, yorigin,
      xorigin + xscale * (index + 1),
      yorigin + yscale * buckets[index]);
  }
}
#endif


/**********************************************************************
 * STATS::plotline
 *
 * Draw a histogram of the stats table. (Line only
 **********************************************************************/

#ifndef GRAPHICS_DISABLED
void STATS::plotline(                //plot stats table
                     WINDOW window,  //to draw in
                     float xorigin,  //bottom left
                     float yorigin,
                     float xscale,   //one x unit
                     float yscale,   //one y unit
                     COLOUR colour   //colour to draw in
                    ) {
  INT32 index;                   //table index

  if (buckets == NULL) {
    /*     err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
       ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
       "Empty stats"); */
    return;
  }
  line_color_index(window, colour); 
  line_type(window, SOLID); 

  move2d (window, xorigin, yorigin + yscale * buckets[0]);
  for (index = 0; index < rangemax - rangemin; index++) {
    draw2d (window, xorigin + xscale * index,
      yorigin + yscale * buckets[index]);
  }
}
#endif


/**********************************************************************
 * choose_nth_item
 *
 * Returns the index of what would b the nth item in the array
 * if the members were sorted, without actually sorting.
 **********************************************************************/

DLLSYM INT32 choose_nth_item(               //fast median
                             INT32 index,   //index to choose
                             float *array,  //array of items
                             INT32 count    //no of items
                            ) {
  static UINT16 seeds[3] = { SEED1, SEED2, SEED3 };
  //for nrand
  INT32 next_sample;             //next one to do
  INT32 next_lesser;             //space for new
  INT32 prev_greater;            //last one saved
  INT32 equal_count;             //no of equal ones
  float pivot;                   //proposed median
  float sample;                  //current sample

  if (count <= 1)
    return 0;
  if (count == 2) {
    if (array[0] < array[1]) {
      return index >= 1 ? 1 : 0;
    }
    else {
      return index >= 1 ? 0 : 1;
    }
  }
  else {
    if (index < 0)
      index = 0;                 //ensure lergal
    else if (index >= count)
      index = count - 1;
    #ifdef __UNIX__
    equal_count = (INT32) (nrand48 (seeds) % count);
    #else
    equal_count = (INT32) (rand () % count);
    #endif
    pivot = array[equal_count];
                                 //fill gap
    array[equal_count] = array[0];
    next_lesser = 0;
    prev_greater = count;
    equal_count = 1;
    for (next_sample = 1; next_sample < prev_greater;) {
      sample = array[next_sample];
      if (sample < pivot) {
                                 //shuffle
        array[next_lesser++] = sample;
        next_sample++;
      }
      else if (sample > pivot) {
        prev_greater--;
                                 //juggle
        array[next_sample] = array[prev_greater];
        array[prev_greater] = sample;
      }
      else {
        equal_count++;
        next_sample++;
      }
    }
    for (next_sample = next_lesser; next_sample < prev_greater;)
      array[next_sample++] = pivot;
    if (index < next_lesser)
      return choose_nth_item (index, array, next_lesser);
    else if (index < prev_greater)
      return next_lesser;        //in equal bracket
    else
      return choose_nth_item (index - prev_greater,
        array + prev_greater,
        count - prev_greater) + prev_greater;
  }
}


/**********************************************************************
 * choose_nth_item
 *
 * Returns the index of what would b the nth item in the array
 * if the members were sorted, without actually sorting.
 **********************************************************************/

DLLSYM INT32
choose_nth_item (                //fast median
INT32 index,                     //index to choose
void *array,                     //array of items
INT32 count,                     //no of items
size_t size,                     //element size
                                 //comparator
int (*compar) (const void *, const void *)
) {
  static UINT16 seeds[3] = { SEED1, SEED2, SEED3 };
  //for nrand
  int result;                    //of compar
  INT32 next_sample;             //next one to do
  INT32 next_lesser;             //space for new
  INT32 prev_greater;            //last one saved
  INT32 equal_count;             //no of equal ones
  INT32 pivot;                   //proposed median

  if (count <= 1)
    return 0;
  if (count == 2) {
    if (compar (array, (char *) array + size) < 0) {
      return index >= 1 ? 1 : 0;
    }
    else {
      return index >= 1 ? 0 : 1;
    }
  }
  if (index < 0)
    index = 0;                   //ensure lergal
  else if (index >= count)
    index = count - 1;
  #ifdef __UNIX__
  pivot = (INT32) (nrand48 (seeds) % count);
  #else
  pivot = (INT32) (rand () % count);
  #endif
  swap_entries (array, size, pivot, 0);
  next_lesser = 0;
  prev_greater = count;
  equal_count = 1;
  for (next_sample = 1; next_sample < prev_greater;) {
    result =
      compar ((char *) array + size * next_sample,
      (char *) array + size * next_lesser);
    if (result < 0) {
      swap_entries (array, size, next_lesser++, next_sample++);
      //shuffle
    }
    else if (result > 0) {
      prev_greater--;
      swap_entries(array, size, prev_greater, next_sample); 
    }
    else {
      equal_count++;
      next_sample++;
    }
  }
  if (index < next_lesser)
    return choose_nth_item (index, array, next_lesser, size, compar);
  else if (index < prev_greater)
    return next_lesser;          //in equal bracket
  else
    return choose_nth_item (index - prev_greater,
      (char *) array + size * prev_greater,
      count - prev_greater, size,
      compar) + prev_greater;
}


/**********************************************************************
 * swap_entries
 *
 * Swap 2 entries of abitrary size in-place in a table.
 **********************************************************************/

void swap_entries(               //swap in place
                  void *array,   //array of entries
                  size_t size,   //size of entry
                  INT32 index1,  //entries to swap
                  INT32 index2) {
  char tmp;
  char *ptr1;                    //to entries
  char *ptr2;
  size_t count;                  //of bytes

  ptr1 = (char *) array + index1 * size;
  ptr2 = (char *) array + index2 * size;
  for (count = 0; count < size; count++) {
    tmp = *ptr1;
    *ptr1++ = *ptr2;
    *ptr2++ = tmp;               //tedious!
  }
}