1 /**********************************************************************
2 * File: statistc.c (Formerly stats.c)
3 * Description: Simple statistical package for integer values.
4 * Author: Ray Smith
5 * Created: Mon Feb 04 16:56:05 GMT 1991
6 *
7 * (C) Copyright 1991, Hewlett-Packard Ltd.
8 ** Licensed under the Apache License, Version 2.0 (the "License");
9 ** you may not use this file except in compliance with the License.
10 ** You may obtain a copy of the License at
11 ** http://www.apache.org/licenses/LICENSE-2.0
12 ** Unless required by applicable law or agreed to in writing, software
13 ** distributed under the License is distributed on an "AS IS" BASIS,
14 ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15 ** See the License for the specific language governing permissions and
16 ** limitations under the License.
17 *
18 **********************************************************************/
19
20 #include "mfcpch.h" //precompiled headers
21 #include <string.h>
22 #include <math.h>
23 #include <stdlib.h>
24 #include "memry.h"
25 //#include "ipeerr.h"
26 #include "tprintf.h"
27 #include "statistc.h"
28
29 #define SEED1 0x1234 //default seeds
30 #define SEED2 0x5678
31 #define SEED3 0x9abc
32
33 /**********************************************************************
34 * STATS::STATS
35 *
36 * Construct a new stats element by allocating and zeroing the memory.
37 **********************************************************************/
38
STATS(inT32 min,inT32 max)39 STATS::STATS( //constructor
40 inT32 min, //min of range
41 inT32 max //max of range
42 ) {
43
44 if (max <= min) {
45 /* err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
46 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
47 "Illegal range for stats, Min=%d, Max=%d",min,max);*/
48 min = 0;
49 max = 1;
50 }
51 rangemin = min; //setup
52 rangemax = max;
53 buckets = (inT32 *) alloc_mem ((max - min) * sizeof (inT32));
54 if (buckets != NULL)
55 this->clear (); //zero it
56 /* else
57 err.log(RESULT_NO_MEMORY,E_LOC,ERR_PRIMITIVES,
58 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
59 "No memory for stats, Min=%d, Max=%d",min,max); */
60 }
61
62
STATS()63 STATS::STATS() { //constructor
64 rangemax = 0; //empty
65 rangemin = 0;
66 buckets = NULL;
67 }
68
69
70 /**********************************************************************
71 * STATS::set_range
72 *
73 * Alter the range on an existing stats element.
74 **********************************************************************/
75
set_range(inT32 min,inT32 max)76 bool STATS::set_range( //constructor
77 inT32 min, //min of range
78 inT32 max //max of range
79 ) {
80
81 if (max <= min) {
82 return false;
83 }
84 rangemin = min; //setup
85 rangemax = max;
86 if (buckets != NULL)
87 free_mem(buckets); //no longer want it
88 buckets = (inT32 *) alloc_mem ((max - min) * sizeof (inT32));
89 /* if (buckets==NULL)
90 return err.log(RESULT_NO_MEMORY,E_LOC,ERR_PRIMITIVES,
91 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
92 "No memory for stats, Min=%d, Max=%d",min,max);*/
93
94 this->clear (); //zero it
95 return true;
96 }
97
98
99 /**********************************************************************
100 * STATS::clear
101 *
102 * Clear out the STATS class by zeroing all the buckets.
103 **********************************************************************/
104
clear()105 void STATS::clear() { //clear out buckets
106 total_count = 0;
107 if (buckets != NULL)
108 memset (buckets, 0, (rangemax - rangemin) * sizeof (inT32));
109 //zero it
110 }
111
112
113 /**********************************************************************
114 * STATS::~STATS
115 *
116 * Destructor for a stats class.
117 **********************************************************************/
118
~STATS()119 STATS::~STATS ( //destructor
120 ) {
121 if (buckets != NULL) {
122 free_mem(buckets);
123 buckets = NULL;
124 }
125 }
126
127
128 /**********************************************************************
129 * STATS::add
130 *
131 * Add a set of samples to (or delete from) a pile.
132 **********************************************************************/
133
add(inT32 value,inT32 count)134 void STATS::add( //add sample
135 inT32 value, //bucket
136 inT32 count //no to add
137 ) {
138 if (buckets == NULL) {
139 /* err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
140 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
141 "Empty stats");*/
142 return;
143 }
144 if (value <= rangemin)
145 buckets[0] += count; //silently clip to range
146 else if (value >= rangemax)
147 buckets[rangemax - rangemin - 1] += count;
148 else
149 //add count to cell
150 buckets[value - rangemin] += count;
151 total_count += count; //keep count of total
152 }
153
154
155 /**********************************************************************
156 * STATS::mode
157 *
158 * Find the mode of a stats class.
159 **********************************************************************/
160
mode()161 inT32 STATS::mode() { //get mode of samples
162 inT32 index; //current index
163 inT32 max; //max cell count
164 inT32 maxindex; //index of max
165
166 if (buckets == NULL) {
167 /* err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
168 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
169 "Empty stats");*/
170 return rangemin;
171 }
172 for (max = 0, maxindex = 0, index = rangemax - rangemin - 1; index >= 0;
173 index--) {
174 if (buckets[index] > max) {
175 max = buckets[index]; //find biggest
176 maxindex = index;
177 }
178 }
179 return maxindex + rangemin; //index of biggest
180 }
181
182
183 /**********************************************************************
184 * STATS::mean
185 *
186 * Find the mean of a stats class.
187 **********************************************************************/
188
mean()189 float STATS::mean() { //get mean of samples
190 inT32 index; //current index
191 inT32 sum; //sum of cells
192
193 if (buckets == NULL) {
194 /* err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
195 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
196 "Empty stats");*/
197 return (float) rangemin;
198 }
199 for (sum = 0, index = rangemax - rangemin - 1; index >= 0; index--) {
200 //sum all buckets
201 sum += index * buckets[index];
202 }
203 if (total_count > 0)
204 //mean value
205 return (float) sum / total_count + rangemin;
206 else
207 return (float) rangemin; //no mean
208 }
209
210
211 /**********************************************************************
212 * STATS::sd
213 *
214 * Find the standard deviation of a stats class.
215 **********************************************************************/
216
sd()217 float STATS::sd() { //standard deviation
218 inT32 index; //current index
219 inT32 sum; //sum of cells
220 inT32 sqsum; //sum of squares
221 float variance;
222
223 if (buckets == NULL) {
224 /* err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
225 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
226 "Empty stats"); */
227 return (float) 0.0;
228 }
229 for (sum = 0, sqsum = 0, index = rangemax - rangemin - 1; index >= 0;
230 index--) {
231 //sum all buckets
232 sum += index * buckets[index];
233 //and squares
234 sqsum += index * index * buckets[index];
235 }
236 if (total_count > 0) {
237 variance = sum / ((float) total_count);
238 variance = sqsum / ((float) total_count) - variance * variance;
239 return (float) sqrt (variance);
240 }
241 else
242 return (float) 0.0;
243 }
244
245
246 /**********************************************************************
247 * STATS::ile
248 *
249 * Find an arbitrary %ile of a stats class.
250 **********************************************************************/
251
ile(float frac)252 float STATS::ile( //percentile
253 float frac //fraction to find
254 ) {
255 inT32 index; //current index
256 inT32 sum; //sum of cells
257 float target; //target value
258
259 if (buckets == NULL) {
260 /* err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
261 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
262 "Empty stats"); */
263 return (float) rangemin;
264 }
265 target = frac * total_count;
266 if (target <= 0)
267 target = (float) 1;
268 if (target > total_count)
269 target = (float) total_count;
270 for (sum = 0, index = 0; index < rangemax - rangemin
271 && sum < target; sum += buckets[index], index++);
272 if (index > 0)
273 return rangemin + index - (sum - target) / buckets[index - 1];
274 //better than just ints
275 else
276 return (float) rangemin;
277 }
278
279
280 /**********************************************************************
281 * STATS::median
282 *
283 * Finds a more usefule estimate of median than ile(0.5).
284 *
285 * Overcomes a problem with ile() - if the samples are, for example,
286 * 6,6,13,14 ile(0.5) return 7.0 - when a more useful value would be midway
287 * between 6 and 13 = 9.5
288 **********************************************************************/
289
median()290 float STATS::median() { //get median
291 float median;
292 inT32 min_pile;
293 inT32 median_pile;
294 inT32 max_pile;
295
296 if (buckets == NULL) {
297 /* err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
298 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
299 "Empty stats");*/
300 return (float) rangemin;
301 }
302 median = (float) ile ((float) 0.5);
303 median_pile = (inT32) floor (median);
304 if ((total_count > 1) && (pile_count (median_pile) == 0)) {
305 /* Find preceeding non zero pile */
306 for (min_pile = median_pile; pile_count (min_pile) == 0; min_pile--);
307 /* Find following non zero pile */
308 for (max_pile = median_pile; pile_count (max_pile) == 0; max_pile++);
309 median = (float) ((min_pile + max_pile) / 2.0);
310 }
311 return median;
312 }
313
314
315 /**********************************************************************
316 * STATS::smooth
317 *
318 * Apply a triangular smoothing filter to the stats.
319 * This makes the modes a bit more useful.
320 * The factor gives the height of the triangle, i.e. the weight of the
321 * centre.
322 **********************************************************************/
323
smooth(inT32 factor)324 void STATS::smooth( //smooth samples
325 inT32 factor //size of triangle
326 ) {
327 inT32 entry; //bucket index
328 inT32 offset; //from entry
329 inT32 entrycount; //no of entries
330 inT32 bucket; //new smoothed pile
331 //output stats
332 STATS result(rangemin, rangemax);
333
334 if (buckets == NULL) {
335 /* err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
336 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
337 "Empty stats"); */
338 return;
339 }
340 if (factor < 2)
341 return; //is a no-op
342 entrycount = rangemax - rangemin;
343 for (entry = 0; entry < entrycount; entry++) {
344 //centre weight
345 bucket = buckets[entry] * factor;
346 for (offset = 1; offset < factor; offset++) {
347 if (entry - offset >= 0)
348 bucket += buckets[entry - offset] * (factor - offset);
349 if (entry + offset < entrycount)
350 bucket += buckets[entry + offset] * (factor - offset);
351 }
352 result.add (entry + rangemin, bucket);
353 }
354 total_count = result.total_count;
355 memcpy (buckets, result.buckets, entrycount * sizeof (inT32));
356 }
357
358
359 /**********************************************************************
360 * STATS::cluster
361 *
362 * Cluster the samples into max_cluster clusters.
363 * Each call runs one iteration. The array of clusters must be
364 * max_clusters+1 in size as cluster 0 is used to indicate which samples
365 * have been used.
366 * The return value is the current number of clusters.
367 **********************************************************************/
368
cluster(float lower,float upper,float multiple,inT32 max_clusters,STATS * clusters)369 inT32 STATS::cluster( //cluster samples
370 float lower, //thresholds
371 float upper,
372 float multiple, //distance threshold
373 inT32 max_clusters, //max no to make
374 STATS *clusters //array of clusters
375 ) {
376 BOOL8 new_cluster; //added one
377 float *centres; //cluster centres
378 inT32 entry; //bucket index
379 inT32 cluster; //cluster index
380 inT32 best_cluster; //one to assign to
381 inT32 new_centre = 0; //residual mode
382 inT32 new_mode; //pile count of new_centre
383 inT32 count; //pile to place
384 float dist; //from cluster
385 float min_dist; //from best_cluster
386 inT32 cluster_count; //no of clusters
387
388 if (max_clusters < 1)
389 return 0;
390 if (buckets == NULL) {
391 /* err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
392 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
393 "Empty stats");*/
394 return 0;
395 }
396 centres = (float *) alloc_mem ((max_clusters + 1) * sizeof (float));
397 if (centres == NULL) {
398 /* err.log(RESULT_NO_MEMORY,E_LOC,ERR_PRIMITIVES,
399 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
400 "No memory for centres"); */
401 return 0;
402 }
403 for (cluster_count = 1; cluster_count <= max_clusters
404 && clusters[cluster_count].buckets != NULL
405 && clusters[cluster_count].total_count > 0; cluster_count++) {
406 centres[cluster_count] =
407 (float) clusters[cluster_count].ile ((float) 0.5);
408 new_centre = clusters[cluster_count].mode ();
409 for (entry = new_centre - 1; centres[cluster_count] - entry < lower
410 && entry >= rangemin
411 && pile_count (entry) <= pile_count (entry + 1); entry--) {
412 count = pile_count (entry) - clusters[0].pile_count (entry);
413 if (count > 0) {
414 clusters[cluster_count].add (entry, count);
415 clusters[0].add (entry, count);
416 }
417 }
418 for (entry = new_centre + 1; entry - centres[cluster_count] < lower
419 && entry < rangemax
420 && pile_count (entry) <= pile_count (entry - 1); entry++) {
421 count = pile_count (entry) - clusters[0].pile_count (entry);
422 if (count > 0) {
423 clusters[cluster_count].add (entry, count);
424 clusters[0].add (entry, count);
425 }
426 }
427 }
428 cluster_count--;
429
430 if (cluster_count == 0) {
431 clusters[0].set_range (rangemin, rangemax);
432 }
433 do {
434 new_cluster = FALSE;
435 new_mode = 0;
436 for (entry = 0; entry < rangemax - rangemin; entry++) {
437 count = buckets[entry] - clusters[0].buckets[entry];
438 //remaining pile
439 if (count > 0) { //any to handle
440 min_dist = (float) MAX_INT32;
441 best_cluster = 0;
442 for (cluster = 1; cluster <= cluster_count; cluster++) {
443 dist = entry + rangemin - centres[cluster];
444 //find distance
445 if (dist < 0)
446 dist = -dist;
447 if (dist < min_dist) {
448 min_dist = dist; //find least
449 best_cluster = cluster;
450 }
451 }
452 if (min_dist > upper //far enough for new
453 && (best_cluster == 0
454 || entry + rangemin > centres[best_cluster] * multiple
455 || entry + rangemin < centres[best_cluster] / multiple)) {
456 if (count > new_mode) {
457 new_mode = count;
458 new_centre = entry + rangemin;
459 }
460 }
461 }
462 }
463 //need new and room
464 if (new_mode > 0 && cluster_count < max_clusters) {
465 cluster_count++;
466 new_cluster = TRUE;
467 if (!clusters[cluster_count].set_range (rangemin, rangemax))
468 return 0;
469 centres[cluster_count] = (float) new_centre;
470 clusters[cluster_count].add (new_centre, new_mode);
471 clusters[0].add (new_centre, new_mode);
472 for (entry = new_centre - 1; centres[cluster_count] - entry < lower
473 && entry >= rangemin
474 && pile_count (entry) <= pile_count (entry + 1); entry--) {
475 count = pile_count (entry) - clusters[0].pile_count (entry);
476 if (count > 0) {
477 clusters[cluster_count].add (entry, count);
478 clusters[0].add (entry, count);
479 }
480 }
481 for (entry = new_centre + 1; entry - centres[cluster_count] < lower
482 && entry < rangemax
483 && pile_count (entry) <= pile_count (entry - 1); entry++) {
484 count = pile_count (entry) - clusters[0].pile_count (entry);
485 if (count > 0) {
486 clusters[cluster_count].add (entry, count);
487 clusters[0].add (entry, count);
488 }
489 }
490 centres[cluster_count] =
491 (float) clusters[cluster_count].ile ((float) 0.5);
492 }
493 }
494 while (new_cluster && cluster_count < max_clusters);
495 free_mem(centres);
496 return cluster_count;
497 }
498
499
500 /**********************************************************************
501 * STATS::local_min
502 *
503 * Return TRUE if this point is a local min.
504 **********************************************************************/
505
local_min(inT32 x)506 BOOL8 STATS::local_min( //test minness
507 inT32 x //of x
508 ) {
509 inT32 index; //table index
510
511 if (buckets == NULL) {
512 /* err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
513 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
514 "Empty stats");*/
515 return FALSE;
516 }
517 if (x < rangemin)
518 x = rangemin;
519 if (x >= rangemax)
520 x = rangemax - 1;
521 x -= rangemin;
522 if (buckets[x] == 0)
523 return TRUE;
524 for (index = x - 1; index >= 0 && buckets[index] == buckets[x]; index--);
525 if (index >= 0 && buckets[index] < buckets[x])
526 return FALSE;
527 for (index = x + 1; index < rangemax - rangemin
528 && buckets[index] == buckets[x]; index++);
529 if (index < rangemax - rangemin && buckets[index] < buckets[x])
530 return FALSE;
531 else
532 return TRUE;
533 }
534
535
536 /**********************************************************************
537 * STATS::print
538 *
539 * Print a summary of the stats and optionally a dump of the table.
540 **********************************************************************/
541
print(FILE *,BOOL8 dump)542 void STATS::print( //print stats table
543 FILE *, //Now uses tprintf instead
544 BOOL8 dump //dump full table
545 ) {
546 inT32 index; //table index
547
548 if (buckets == NULL) {
549 /* err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
550 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
551 "Empty stats"); */
552 return;
553 }
554 if (dump) {
555 for (index = 0; index < rangemax - rangemin; index++) {
556 tprintf ("%4d:%-3d ", rangemin + index, buckets[index]);
557 if (index % 8 == 7)
558 tprintf ("\n");
559 }
560 tprintf ("\n");
561 }
562
563 tprintf ("Total count=%d\n", total_count);
564 tprintf ("Min=%d\n", (inT32) (ile ((float) 0.0)));
565 tprintf ("Lower quartile=%.2f\n", ile ((float) 0.25));
566 tprintf ("Median=%.2f\n", ile ((float) 0.5));
567 tprintf ("Upper quartile=%.2f\n", ile ((float) 0.75));
568 tprintf ("Max=%d\n", (inT32) (ile ((float) 0.99999)));
569 tprintf ("Mean= %.2f\n", mean ());
570 tprintf ("SD= %.2f\n", sd ());
571 }
572
573
574 /**********************************************************************
575 * STATS::min_bucket
576 *
577 * Find REAL minimum bucket - ile(0.0) isnt necessarily correct
578 **********************************************************************/
579
min_bucket()580 inT32 STATS::min_bucket() { //Find min
581 inT32 min;
582
583 if (buckets == NULL) {
584 /* err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
585 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
586 "Empty stats");*/
587 return rangemin;
588 }
589
590 for (min = 0; (min < rangemax - rangemin) && (buckets[min] == 0); min++);
591 return rangemin + min;
592 }
593
594
595 /**********************************************************************
596 * STATS::max_bucket
597 *
598 * Find REAL maximum bucket - ile(1.0) isnt necessarily correct
599 **********************************************************************/
600
max_bucket()601 inT32 STATS::max_bucket() { //Find max
602 inT32 max;
603
604 if (buckets == NULL) {
605 /* err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
606 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
607 "Empty stats");*/
608 return rangemin;
609 }
610
611 for (max = rangemax - rangemin - 1;
612 (max > 0) && (buckets[max] == 0); max--);
613 return rangemin + max;
614 }
615
616
617 /**********************************************************************
618 * STATS::short_print
619 *
620 * Print a summary of the stats and optionally a dump of the table.
621 * ( BUT ONLY THE PART OF THE TABLE BETWEEN MIN AND MAX)
622 **********************************************************************/
623
short_print(FILE *,BOOL8 dump)624 void STATS::short_print( //print stats table
625 FILE *, //Now uses tprintf instead
626 BOOL8 dump //dump full table
627 ) {
628 inT32 index; //table index
629 inT32 min = min_bucket ();
630 inT32 max = max_bucket ();
631
632 if (buckets == NULL) {
633 /* err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
634 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
635 "Empty stats"); */
636 return;
637 }
638 if (dump) {
639 for (index = min; index <= max; index++) {
640 tprintf ("%4d:%-3d ", rangemin + index, buckets[index]);
641 if ((index - min) % 8 == 7)
642 tprintf ("\n");
643 }
644 tprintf ("\n");
645 }
646
647 tprintf ("Total count=%d\n", total_count);
648 tprintf ("Min=%d Really=%d\n", (inT32) (ile ((float) 0.0)), min);
649 tprintf ("Max=%d Really=%d\n", (inT32) (ile ((float) 1.1)), max);
650 tprintf ("Range=%d\n", max + 1 - min);
651 tprintf ("Lower quartile=%.2f\n", ile ((float) 0.25));
652 tprintf ("Median=%.2f\n", ile ((float) 0.5));
653 tprintf ("Upper quartile=%.2f\n", ile ((float) 0.75));
654 tprintf ("Mean= %.2f\n", mean ());
655 tprintf ("SD= %.2f\n", sd ());
656 }
657
658
659 /**********************************************************************
660 * STATS::plot
661 *
662 * Draw a histogram of the stats table.
663 **********************************************************************/
664
plot(ScrollView * window,float xorigin,float yorigin,float xscale,float yscale,ScrollView::Color colour)665 void STATS::plot( //plot stats table
666 ScrollView* window, //to draw in
667 float xorigin, //bottom left
668 float yorigin,
669 float xscale, //one x unit
670 float yscale, //one y unit
671 ScrollView::Color colour //colour to draw in
672 ) {
673 #ifndef GRAPHICS_DISABLED
674 inT32 index; //table index
675
676 if (buckets == NULL) {
677 /* err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
678 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
679 "Empty stats");*/
680 return;
681 }
682 window->Pen(colour);
683
684 for (index = 0; index < rangemax - rangemin; index++) {
685 window->Rectangle( xorigin + xscale * index, yorigin,
686 xorigin + xscale * (index + 1),
687 yorigin + yscale * buckets[index]);
688 }
689 #endif
690 }
691
692
693 /**********************************************************************
694 * STATS::plotline
695 *
696 * Draw a histogram of the stats table. (Line only
697 **********************************************************************/
698
plotline(ScrollView * window,float xorigin,float yorigin,float xscale,float yscale,ScrollView::Color colour)699 void STATS::plotline( //plot stats table
700 ScrollView* window, //to draw in
701 float xorigin, //bottom left
702 float yorigin,
703 float xscale, //one x unit
704 float yscale, //one y unit
705 ScrollView::Color colour //colour to draw in
706 ) {
707 #ifndef GRAPHICS_DISABLED
708 inT32 index; //table index
709
710 if (buckets == NULL) {
711 /* err.log(RESULT_LOGICAL_ERROR,E_LOC,ERR_PRIMITIVES,
712 ERR_SCROLLING,ERR_CONTINUE,ERR_ERROR,
713 "Empty stats"); */
714 return;
715 }
716 window->Pen(colour);
717
718 window->SetCursor(xorigin, yorigin + yscale * buckets[0]);
719 for (index = 0; index < rangemax - rangemin; index++) {
720 window->DrawTo(xorigin + xscale * index, yorigin + yscale * buckets[index]);
721 }
722 #endif
723 }
724
725
726 /**********************************************************************
727 * choose_nth_item
728 *
729 * Returns the index of what would b the nth item in the array
730 * if the members were sorted, without actually sorting.
731 **********************************************************************/
732
choose_nth_item(inT32 index,float * array,inT32 count)733 DLLSYM inT32 choose_nth_item( //fast median
734 inT32 index, //index to choose
735 float *array, //array of items
736 inT32 count //no of items
737 ) {
738 static uinT16 seeds[3] = { SEED1, SEED2, SEED3 };
739 //for nrand
740 inT32 next_sample; //next one to do
741 inT32 next_lesser; //space for new
742 inT32 prev_greater; //last one saved
743 inT32 equal_count; //no of equal ones
744 float pivot; //proposed median
745 float sample; //current sample
746
747 if (count <= 1)
748 return 0;
749 if (count == 2) {
750 if (array[0] < array[1]) {
751 return index >= 1 ? 1 : 0;
752 }
753 else {
754 return index >= 1 ? 0 : 1;
755 }
756 }
757 else {
758 if (index < 0)
759 index = 0; //ensure lergal
760 else if (index >= count)
761 index = count - 1;
762 #ifdef __UNIX__
763 equal_count = (inT32) (nrand48 (seeds) % count);
764 #else
765 equal_count = (inT32) (rand () % count);
766 #endif
767 pivot = array[equal_count];
768 //fill gap
769 array[equal_count] = array[0];
770 next_lesser = 0;
771 prev_greater = count;
772 equal_count = 1;
773 for (next_sample = 1; next_sample < prev_greater;) {
774 sample = array[next_sample];
775 if (sample < pivot) {
776 //shuffle
777 array[next_lesser++] = sample;
778 next_sample++;
779 }
780 else if (sample > pivot) {
781 prev_greater--;
782 //juggle
783 array[next_sample] = array[prev_greater];
784 array[prev_greater] = sample;
785 }
786 else {
787 equal_count++;
788 next_sample++;
789 }
790 }
791 for (next_sample = next_lesser; next_sample < prev_greater;)
792 array[next_sample++] = pivot;
793 if (index < next_lesser)
794 return choose_nth_item (index, array, next_lesser);
795 else if (index < prev_greater)
796 return next_lesser; //in equal bracket
797 else
798 return choose_nth_item (index - prev_greater,
799 array + prev_greater,
800 count - prev_greater) + prev_greater;
801 }
802 }
803
804
805 /**********************************************************************
806 * choose_nth_item
807 *
808 * Returns the index of what would b the nth item in the array
809 * if the members were sorted, without actually sorting.
810 **********************************************************************/
811
812 DLLSYM inT32
choose_nth_item(inT32 index,void * array,inT32 count,size_t size,int (* compar)(const void *,const void *))813 choose_nth_item ( //fast median
814 inT32 index, //index to choose
815 void *array, //array of items
816 inT32 count, //no of items
817 size_t size, //element size
818 //comparator
819 int (*compar) (const void *, const void *)
820 ) {
821 static uinT16 seeds[3] = { SEED1, SEED2, SEED3 };
822 //for nrand
823 int result; //of compar
824 inT32 next_sample; //next one to do
825 inT32 next_lesser; //space for new
826 inT32 prev_greater; //last one saved
827 inT32 equal_count; //no of equal ones
828 inT32 pivot; //proposed median
829
830 if (count <= 1)
831 return 0;
832 if (count == 2) {
833 if (compar (array, (char *) array + size) < 0) {
834 return index >= 1 ? 1 : 0;
835 }
836 else {
837 return index >= 1 ? 0 : 1;
838 }
839 }
840 if (index < 0)
841 index = 0; //ensure lergal
842 else if (index >= count)
843 index = count - 1;
844 #ifdef __UNIX__
845 pivot = (inT32) (nrand48 (seeds) % count);
846 #else
847 pivot = (inT32) (rand () % count);
848 #endif
849 swap_entries (array, size, pivot, 0);
850 next_lesser = 0;
851 prev_greater = count;
852 equal_count = 1;
853 for (next_sample = 1; next_sample < prev_greater;) {
854 result =
855 compar ((char *) array + size * next_sample,
856 (char *) array + size * next_lesser);
857 if (result < 0) {
858 swap_entries (array, size, next_lesser++, next_sample++);
859 //shuffle
860 }
861 else if (result > 0) {
862 prev_greater--;
863 swap_entries(array, size, prev_greater, next_sample);
864 }
865 else {
866 equal_count++;
867 next_sample++;
868 }
869 }
870 if (index < next_lesser)
871 return choose_nth_item (index, array, next_lesser, size, compar);
872 else if (index < prev_greater)
873 return next_lesser; //in equal bracket
874 else
875 return choose_nth_item (index - prev_greater,
876 (char *) array + size * prev_greater,
877 count - prev_greater, size,
878 compar) + prev_greater;
879 }
880
881
882 /**********************************************************************
883 * swap_entries
884 *
885 * Swap 2 entries of abitrary size in-place in a table.
886 **********************************************************************/
887
swap_entries(void * array,size_t size,inT32 index1,inT32 index2)888 void swap_entries( //swap in place
889 void *array, //array of entries
890 size_t size, //size of entry
891 inT32 index1, //entries to swap
892 inT32 index2) {
893 char tmp;
894 char *ptr1; //to entries
895 char *ptr2;
896 size_t count; //of bytes
897
898 ptr1 = (char *) array + index1 * size;
899 ptr2 = (char *) array + index2 * size;
900 for (count = 0; count < size; count++) {
901 tmp = *ptr1;
902 *ptr1++ = *ptr2;
903 *ptr2++ = tmp; //tedious!
904 }
905 }
906