1 /* -*- mode: C; mode: fold; -*- */
2
3 /*
4 Copyright (c) 2003-2007 Massachusetts Institute of Technology
5 Copyright (c) 2013-2017,2018 John E. Davis <jed@jedsoft.org>
6
7 This software was developed by the MIT Center for Space Research
8 under contract SV1-61010 from the Smithsonian Institution.
9
10 Permission to use, copy, modify, distribute, and sell this software
11 and its documentation for any purpose is hereby granted without fee,
12 provided that the above copyright notice appear in all copies and
13 that both that copyright notice and this permission notice appear in
14 the supporting documentation, and that the name of the Massachusetts
15 Institute of Technology not be used in advertising or publicity
16 pertaining to distribution of the software without specific, written
17 prior permission. The Massachusetts Institute of Technology makes
18 no representations about the suitability of this software for any
19 purpose. It is provided "as is" without express or implied warranty.
20
21 THE MASSACHUSETTS INSTITUTE OF TECHNOLOGY DISCLAIMS ALL WARRANTIES
22 WITH REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF
23 MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE MASSACHUSETTS
24 INSTITUTE OF TECHNOLOGY BE LIABLE FOR ANY SPECIAL, INDIRECT OR
25 CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
26 OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
27 NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
28 WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
29 */
30
31 /* Author: John E. Davis (davis@space.mit.edu) */
32
33 #include "config.h"
34
35 #include <stdio.h>
36 #include <math.h>
37 #include <slang.h>
38 #include <string.h>
39
40 #ifdef __cplusplus
41 extern "C"
42 {
43 #endif
44 SLANG_MODULE(histogram);
45 #ifdef __cplusplus
46 }
47 #endif
48
49 #define MODULE_VERSION_STRING "0.4.0"
50 #define MODULE_VERSION_NUMBER 400
51
52 typedef unsigned int HistData_Type;
53 #define HISTDATA_TYPE SLANG_UINT_TYPE
54
55 static const char *Module_Version_String = MODULE_VERSION_STRING;
56
57 #ifndef HAVE_ISNAN
58 # define isnan(x) (x != x)
59 #endif
60
check_grid(double * x,SLuindex_Type n)61 static int check_grid (double *x, SLuindex_Type n)
62 {
63 double xlo;
64 SLuindex_Type i;
65
66 if (n == 0)
67 return 0;
68
69 xlo = x[0];
70
71 if (isnan(xlo))
72 goto return_error;
73
74 for (i = 0; i < n; i++)
75 {
76 if (isnan(x[i]))
77 goto return_error;
78
79 if (x[i] < xlo)
80 goto return_error;
81
82 xlo = x[i];
83 }
84 return 0;
85
86 return_error:
87 SLang_verror (SL_INVALID_PARM, "Invalid grid: Expecting one in increasing order");
88 return -1;
89 }
90
double_binary_search(double x,double * xp,SLuindex_Type n)91 static SLuindex_Type double_binary_search (double x, double *xp, SLuindex_Type n)
92 {
93 SLuindex_Type n0, n1, n2;
94 /* SLuindex_Type j; */
95 double xp0, xp1;
96
97 if (n < 2)
98 return 0;
99
100 if (x < (xp0 = xp[0]))
101 return 0;
102 if (x >= (xp1 = xp[n-1]))
103 return n-1;
104 #if 0
105 j = (SLuindex_Type) (((x - xp0)/(xp1 - xp0))*(n-1));
106 if (j < n)
107 {
108 if (x >= xp[j])
109 {
110 if ((j == n-1) || (xp[j+1] < x))
111 return j;
112 n0 = j;
113 n1 = n;
114 }
115 else
116 {
117 n0 = 0;
118 n1 = j;
119 }
120 }
121 else
122 {
123 n0 = 0;
124 n1 = n;
125 }
126 #else
127 n0 = 0;
128 n1 = n;
129 #endif
130
131 while (n1 > n0 + 1)
132 {
133 /* Note: since n1 > n0+1, then n1 >= n0+2, and n0+n1 >= 2*n0+2
134 * ==> n2 = (n0+n1)/2 >= n0+1.
135 * ==> n2 > n0
136 * Also, n1 >= n0+2 ==> n1-2 >= n0 ==> 2*n1 - 2 >= (n0+n1)
137 * ==> n2 = (n0+n1)/2 <= n1-1 ==> n2+1 <= n1
138 * ==> n2 < n1
139 * Hence, n0 < n1 < n2. This means that the interval determined
140 * by the following code will become smaller.
141 */
142 n2 = (n0 + n1) / 2;
143 if (xp[n2] > x)
144 n1 = n2;
145 else
146 n0 = n2;
147 }
148 return n0;
149 }
150
151 #define BIN_EDGES_TYPE double
152 #define BINARY_SEARCH double_binary_search
153 #define PTS_TYPE unsigned char
154 #define HISTOGRAM_1D uc_histogram_1d
155 #define HISTOGRAM_2D uc_histogram_2d
156 #define CHECK_NANS 0
157 #include "histogram-module.inc"
158
159 #define BIN_EDGES_TYPE double
160 #define BINARY_SEARCH double_binary_search
161 #define PTS_TYPE int
162 #define HISTOGRAM_1D i_histogram_1d
163 #define HISTOGRAM_2D i_histogram_2d
164 #define CHECK_NANS 0
165 #include "histogram-module.inc"
166
167 #define BIN_EDGES_TYPE double
168 #define BINARY_SEARCH double_binary_search
169 #define PTS_TYPE float
170 #define HISTOGRAM_1D f_histogram_1d
171 #define HISTOGRAM_2D f_histogram_2d
172 #define CHECK_NANS 1
173 #include "histogram-module.inc"
174
175 #define BIN_EDGES_TYPE double
176 #define BINARY_SEARCH double_binary_search
177 #define PTS_TYPE double
178 #define HISTOGRAM_1D d_histogram_1d
179 #define HISTOGRAM_2D d_histogram_2d
180 #define HISTOGRAM_REBIN d_hist1d_rebin
181 #define CHECK_NANS 1
182 #include "histogram-module.inc"
183
uc_fast_hist_1d(unsigned char * pts,SLuindex_Type npts,double * bin_edges,SLuindex_Type nbins,HistData_Type * histogram)184 static int uc_fast_hist_1d (unsigned char *pts, SLuindex_Type npts,
185 double *bin_edges, SLuindex_Type nbins,
186 HistData_Type *histogram)
187 {
188 HistData_Type h[256];
189 SLuindex_Type i;
190 SLuindex_Type nbins_m1;
191
192 if (nbins == 0)
193 return 0;
194
195 if (-1 == check_grid (bin_edges, nbins))
196 return -1;
197
198 for (i = 0; i < 256; i++)
199 h[i] = 0.0;
200
201 for (i = 0; i < npts; i++)
202 h[pts[i]]++;
203
204 nbins_m1 = nbins - 1;
205 for (i = 0; i < nbins_m1; i++)
206 {
207 double x0, x1;
208 SLuindex_Type ix0;
209
210 x1 = bin_edges[i+1];
211 if (x1 <= 0.0)
212 continue;
213
214 x0 = bin_edges[i];
215 if (x0 < 0.0)
216 x0 = 0.0;
217 ix0 = (SLuindex_Type) ceil (x0);
218
219 while (i < nbins_m1)
220 {
221 SLuindex_Type ix1;
222 SLuindex_Type j, jmax;
223
224 x1 = bin_edges[i+1];
225 ix1 = (SLuindex_Type) ceil (x1);
226 jmax = ix1;
227 if (jmax > 256) jmax = 256;
228
229 for (j = ix0; j < jmax; j++)
230 histogram[i] += h[j];
231
232 ix0 = ix1;
233 if (ix0 > 255)
234 break;
235 i++;
236 }
237 }
238
239 /* Handle overflow bin */
240
241 if (bin_edges[nbins_m1] < 0)
242 i = 0;
243 else
244 i = (SLuindex_Type) ceil (bin_edges[nbins_m1]);
245
246 while (i < 256)
247 {
248 histogram[nbins_m1] += h[i];
249 i++;
250 }
251 return 0;
252 }
253
convert_reverse_indices(SLindex_Type * r,SLuindex_Type num_r,SLang_Array_Type * h)254 static SLang_Array_Type *convert_reverse_indices (SLindex_Type *r, SLuindex_Type num_r, SLang_Array_Type *h)
255 {
256 SLang_Array_Type *new_r;
257 SLang_Array_Type **new_r_data;
258 SLuindex_Type i, num_h;
259 SLindex_Type *lens;
260
261 if (NULL == (new_r = SLang_create_array (SLANG_ARRAY_TYPE, 0, NULL, h->dims, h->num_dims)))
262 return NULL;
263
264 num_h = h->num_elements;
265
266 if (NULL == (lens = (SLindex_Type *)SLmalloc (num_h * sizeof (SLindex_Type))))
267 {
268 SLang_free_array (new_r);
269 return NULL;
270 }
271 memset ((char *)lens, 0, num_h*sizeof(SLindex_Type));
272
273 for (i = 0; i < num_r; i++)
274 {
275 SLindex_Type r_i = r[i];
276
277 if (r_i >= 0)
278 lens[r_i]++;
279 }
280
281 new_r_data = (SLang_Array_Type **) new_r->data;
282 for (i = 0; i < num_h; i++)
283 {
284 if (NULL == (new_r_data[i] = SLang_create_array (SLANG_ARRAY_INDEX_TYPE, 0, NULL, &lens[i], 1)))
285 goto return_error;
286
287 lens[i] = 0;
288 }
289
290 for (i = 0; i < num_r; i++)
291 {
292 SLang_Array_Type *at;
293 SLindex_Type r_i = r[i];
294
295 if (r_i < 0)
296 continue;
297
298 at = new_r_data[r_i];
299
300 ((SLindex_Type *)at->data)[lens[r_i]] = i;
301 lens[r_i]++;
302 }
303
304 SLfree ((char *)lens);
305 return new_r;
306
307 return_error:
308 SLfree ((char *) lens);
309 SLang_free_array (new_r);
310 return NULL;
311 }
312
map_to_best_type(int type,int * mtypep)313 static int map_to_best_type (int type, int *mtypep)
314 {
315 switch (type)
316 {
317 case SLANG_UCHAR_TYPE:
318 *mtypep = SLANG_UCHAR_TYPE;
319 break;
320
321 case SLANG_CHAR_TYPE:
322 case SLANG_SHORT_TYPE:
323 case SLANG_INT_TYPE:
324 *mtypep = SLANG_INT_TYPE;
325 break;
326
327 case SLANG_FLOAT_TYPE:
328 *mtypep = SLANG_FLOAT_TYPE;
329 break;
330
331 default:
332 *mtypep = SLANG_DOUBLE_TYPE;
333 break;
334 }
335 return 0;
336 }
337
pop_1d_array_of_type(SLang_Array_Type ** atp,int type)338 static int pop_1d_array_of_type (SLang_Array_Type **atp, int type)
339 {
340 SLang_Array_Type *at;
341
342 if (-1 == SLang_pop_array_of_type (&at, type))
343 return -1;
344
345 if (at->num_dims != 1)
346 {
347 SLang_verror (SL_INVALID_PARM, "Expecting a 1-d array");
348 SLang_free_array (at);
349 return -1;
350 }
351 *atp = at;
352 return 0;
353 }
354
pop_1d_double_arrays(SLang_Array_Type ** ap,SLang_Array_Type ** bp)355 static int pop_1d_double_arrays (SLang_Array_Type **ap, SLang_Array_Type **bp)
356 {
357 SLang_Array_Type *a, *b;
358
359 *ap = *bp = NULL;
360
361 if (-1 == pop_1d_array_of_type (&b, SLANG_DOUBLE_TYPE))
362 return -1;
363
364 if (-1 == pop_1d_array_of_type (&a, SLANG_DOUBLE_TYPE))
365 {
366 SLang_free_array (b);
367 return -1;
368 }
369
370 if (a->num_elements != b->num_elements)
371 {
372 SLang_verror (SL_INVALID_PARM, "Arrays do not match in size");
373 SLang_free_array (b);
374 SLang_free_array (a);
375 return -1;
376 }
377
378 *ap = a;
379 *bp = b;
380 return 0;
381 }
382
pop_hist1d_pts_array(SLang_Array_Type ** at)383 static int pop_hist1d_pts_array (SLang_Array_Type **at)
384 {
385 int type;
386
387 if (-1 == map_to_best_type (SLang_peek_at_stack1 (), &type))
388 return -1;
389 #if 1
390 if (-1 == SLang_pop_array_of_type (at, type))
391 return -1;
392 #else
393 if (-1 == pop_1d_array_of_type (at, type))
394 return -1;
395 #endif
396 return 0;
397 }
398
pop_hist2d_pts_array(SLang_Array_Type ** atxp,SLang_Array_Type ** atyp)399 static int pop_hist2d_pts_array (SLang_Array_Type **atxp,
400 SLang_Array_Type **atyp)
401 {
402 int ytype, xtype;
403 SLang_Array_Type *atx, *aty;
404
405 ytype = SLang_peek_at_stack1 ();
406 if (-1 == SLroll_stack (2))
407 return -1;
408 xtype = SLang_peek_at_stack1 ();
409
410 if (-1 == map_to_best_type (ytype, &ytype))
411 return -1;
412 if (-1 == map_to_best_type (xtype, &xtype))
413 return -1;
414
415 if (xtype != ytype)
416 ytype = xtype = SLANG_DOUBLE_TYPE;
417
418 if (-1 == pop_1d_array_of_type (&atx, xtype))
419 return -1;
420
421 if (-1 == pop_1d_array_of_type (&aty, ytype))
422 {
423 SLang_free_array (atx);
424 return -1;
425 }
426
427 if (atx->num_elements != aty->num_elements)
428 {
429 SLang_verror (SL_INVALID_PARM, "hist2d: x and y points arrays must match in size");
430 SLang_free_array (aty);
431 SLang_free_array (atx);
432 return -1;
433 }
434
435 *atxp = atx;
436 *atyp = aty;
437 return 0;
438 }
439
alloc_reverse_indices(SLuindex_Type num)440 static SLindex_Type *alloc_reverse_indices (SLuindex_Type num)
441 {
442 SLuindex_Type i;
443 SLindex_Type *r;
444
445 if (NULL == (r = (SLindex_Type *) SLmalloc ((num + 1) * sizeof(SLindex_Type))))
446 return NULL;
447
448 for (i = 0; i < num; i++)
449 r[i] = -1;
450
451 return r;
452 }
453
hist1d(void)454 static void hist1d (void)
455 {
456 SLang_Array_Type *edges_at, *hist_at, *pts_at, *indices_at;
457 SLang_Ref_Type *ref;
458 SLindex_Type *rev_indices;
459 int status, type, has_hist = 0;
460
461 ref = NULL;
462 switch (SLang_Num_Function_Args)
463 {
464 case 4:
465 has_hist = 1;
466 if (SLang_peek_at_stack () == SLANG_NULL_TYPE)
467 (void) SLdo_pop ();
468 else if (-1 == SLang_pop_ref (&ref))
469 return;
470 break;
471
472 case 3:
473 type = SLang_peek_at_stack ();
474 if (type == SLANG_REF_TYPE)
475 {
476 if (-1 == SLang_pop_ref (&ref))
477 return;
478 }
479 else if (type == SLANG_NULL_TYPE)
480 (void) SLdo_pop ();
481 else
482 has_hist = 1;
483 break;
484
485 case 2:
486 break;
487
488 default:
489 SLang_verror (SL_USAGE_ERROR, "h = hist1d ([h,] points, bins [,&reverse-indices])");
490 return;
491 }
492
493 if (-1 == SLang_pop_array_of_type (&edges_at, SLANG_DOUBLE_TYPE))
494 {
495 SLang_free_ref (ref);
496 return;
497 }
498
499 pts_at = NULL;
500 hist_at = NULL;
501 indices_at = NULL;
502 rev_indices = NULL;
503
504 if (-1 == pop_hist1d_pts_array (&pts_at))
505 goto free_and_return;
506
507 if (has_hist && (SLang_peek_at_stack () == SLANG_NULL_TYPE))
508 {
509 (void) SLdo_pop ();
510 has_hist = 0;
511 }
512
513 if (has_hist)
514 {
515 if (-1 == SLang_pop_array_of_type (&hist_at, HISTDATA_TYPE))
516 goto free_and_return;
517 if (hist_at->flags & SLARR_DATA_VALUE_IS_READ_ONLY)
518 {
519 SLang_verror (SL_InvalidParm_Error, "Input histogram array is read-only");
520 goto free_and_return;
521 }
522 if ((hist_at->num_dims != 1)
523 || (hist_at->dims[0] != edges_at->dims[0]))
524 {
525 SLang_verror (SL_InvalidParm_Error, "Input histogram array is not 1d or does not match grid");
526 goto free_and_return;
527 }
528 }
529 else
530 {
531 if (NULL == (hist_at = SLang_create_array (HISTDATA_TYPE, 0, NULL, edges_at->dims, 1)))
532 goto free_and_return;
533 /* hist_at->data is already initialized to 0 */
534 }
535
536 if ((ref != NULL)
537 && (NULL == (rev_indices = alloc_reverse_indices (pts_at->num_elements))))
538 goto free_and_return;
539
540 switch (pts_at->data_type)
541 {
542 case SLANG_UCHAR_TYPE:
543 if (rev_indices == NULL)
544 status = uc_fast_hist_1d ((unsigned char *)pts_at->data, pts_at->num_elements,
545 (double *)edges_at->data, edges_at->num_elements,
546 (HistData_Type *) hist_at->data);
547 else status = uc_histogram_1d ((unsigned char *)pts_at->data, pts_at->num_elements,
548 (double *)edges_at->data, edges_at->num_elements,
549 (HistData_Type *) hist_at->data, rev_indices);
550 break;
551
552 case SLANG_INT_TYPE:
553 status = i_histogram_1d ((int *)pts_at->data, pts_at->num_elements,
554 (double *)edges_at->data, edges_at->num_elements,
555 (HistData_Type *) hist_at->data, rev_indices);
556 break;
557
558 case SLANG_FLOAT_TYPE:
559 status = f_histogram_1d ((float *)pts_at->data, pts_at->num_elements,
560 (double *)edges_at->data, edges_at->num_elements,
561 (HistData_Type *) hist_at->data, rev_indices);
562 break;
563
564 case SLANG_DOUBLE_TYPE:
565 status = d_histogram_1d ((double *)pts_at->data, pts_at->num_elements,
566 (double *)edges_at->data, edges_at->num_elements,
567 (HistData_Type *) hist_at->data, rev_indices);
568 break;
569
570 default:
571 SLang_verror (SL_INTERNAL_ERROR, "Error in hist1d: array not supported");
572 status = -1;
573 }
574
575 if (status == -1)
576 goto free_and_return;
577
578 if (ref != NULL)
579 {
580 if (NULL == (indices_at = convert_reverse_indices (rev_indices, pts_at->num_elements, hist_at)))
581 goto free_and_return;
582
583 if (-1 == SLang_assign_to_ref (ref, SLANG_ARRAY_TYPE, (VOID_STAR) &indices_at))
584 goto free_and_return;
585 }
586
587 (void) SLang_push_array (hist_at, 0);
588
589 /* NULLs ok below */
590 free_and_return:
591 if (rev_indices != NULL) SLfree ((char *) rev_indices);
592 SLang_free_ref (ref);
593 SLang_free_array (indices_at);
594 SLang_free_array (edges_at);
595 SLang_free_array (pts_at);
596 SLang_free_array (hist_at);
597 }
598
hist2d(void)599 static void hist2d (void)
600 {
601 SLang_Array_Type *xedges_at, *yedges_at, *xpts_at, *ypts_at;
602 SLang_Array_Type *hist_at, *indices_at;
603 SLang_Ref_Type *ref;
604 SLindex_Type *rev_indices;
605 int has_hist = 0;
606 int type;
607 SLindex_Type dims[2];
608
609 ref = NULL;
610 switch (SLang_Num_Function_Args)
611 {
612 case 6:
613 has_hist = 1;
614 if (SLang_peek_at_stack () == SLANG_NULL_TYPE)
615 (void) SLdo_pop ();
616 else if (-1 == SLang_pop_ref (&ref))
617 return;
618 break;
619
620 case 5:
621 type = SLang_peek_at_stack ();
622 if (type == SLANG_REF_TYPE)
623 {
624 if (-1 == SLang_pop_ref (&ref))
625 return;
626 }
627 else if (type == SLANG_NULL_TYPE)
628 (void) SLdo_pop ();
629 else
630 has_hist = 1;
631 break;
632
633 case 4:
634 break;
635
636 default:
637 SLang_verror (SL_USAGE_ERROR, "h = hist2d ([h,] xpnts, ypnts, xbins, ybins [,&reverse-indices])");
638 return;
639 }
640
641 if (-1 == SLang_pop_array_of_type (&yedges_at, SLANG_DOUBLE_TYPE))
642 {
643 SLang_free_ref (ref);
644 return;
645 }
646
647 if (-1 == SLang_pop_array_of_type (&xedges_at, SLANG_DOUBLE_TYPE))
648 {
649 SLang_free_array (yedges_at);
650 SLang_free_ref (ref);
651 return;
652 }
653
654 xpts_at = ypts_at = NULL;
655 hist_at = NULL;
656 indices_at = NULL;
657 rev_indices = NULL;
658
659 if (-1 == pop_hist2d_pts_array (&xpts_at, &ypts_at))
660 goto free_and_return;
661
662 if (has_hist && (SLang_peek_at_stack () == SLANG_NULL_TYPE))
663 {
664 (void) SLdo_pop ();
665 has_hist = 0;
666 }
667
668 dims[0] = xedges_at->num_elements;
669 dims[1] = yedges_at->num_elements;
670
671 if (has_hist)
672 {
673 if (-1 == SLang_pop_array_of_type (&hist_at, HISTDATA_TYPE))
674 goto free_and_return;
675 if (hist_at->flags & SLARR_DATA_VALUE_IS_READ_ONLY)
676 {
677 SLang_verror (SL_InvalidParm_Error, "Input histogram array is read-only");
678 goto free_and_return;
679 }
680 if ((hist_at->num_dims != 2)
681 || (hist_at->dims[0] != dims[0])
682 || (hist_at->dims[1] != dims[1]))
683 {
684 SLang_verror (SL_InvalidParm_Error, "Input histogram array is not 2d or does not match grids");
685 goto free_and_return;
686 }
687 }
688 else if (NULL == (hist_at = SLang_create_array (HISTDATA_TYPE, 0, NULL, dims, 2)))
689 goto free_and_return;
690 /* hist_at->data is already initialized to 0 */
691
692 if ((ref != NULL)
693 && (NULL == (rev_indices = alloc_reverse_indices (xpts_at->num_elements))))
694 goto free_and_return;
695
696 switch (xpts_at->data_type)
697 {
698 case SLANG_UCHAR_TYPE:
699 if (-1 == uc_histogram_2d ((unsigned char *)xpts_at->data, (unsigned char *)ypts_at->data, ypts_at->num_elements,
700 (double *)xedges_at->data, xedges_at->num_elements,
701 (double *)yedges_at->data, yedges_at->num_elements,
702 (HistData_Type *) hist_at->data, rev_indices))
703 goto free_and_return;
704 break;
705
706 case SLANG_INT_TYPE:
707 if (-1 == i_histogram_2d ((int *)xpts_at->data, (int *)ypts_at->data, ypts_at->num_elements,
708 (double *)xedges_at->data, xedges_at->num_elements,
709 (double *)yedges_at->data, yedges_at->num_elements,
710 (HistData_Type *) hist_at->data, rev_indices))
711 goto free_and_return;
712 break;
713
714 case SLANG_FLOAT_TYPE:
715 if (-1 == f_histogram_2d ((float *)xpts_at->data, (float *)ypts_at->data, ypts_at->num_elements,
716 (double *)xedges_at->data, xedges_at->num_elements,
717 (double *)yedges_at->data, yedges_at->num_elements,
718 (HistData_Type *) hist_at->data, rev_indices))
719 goto free_and_return;
720 break;
721
722 case SLANG_DOUBLE_TYPE:
723 if (-1 == d_histogram_2d ((double *)xpts_at->data, (double *)ypts_at->data, ypts_at->num_elements,
724 (double *)xedges_at->data, xedges_at->num_elements,
725 (double *)yedges_at->data, yedges_at->num_elements,
726 (HistData_Type *) hist_at->data, rev_indices))
727 goto free_and_return;
728 break;
729
730 default:
731 SLang_verror (SL_INTERNAL_ERROR, "Error in hist2d: array not supported");
732 goto free_and_return;
733 }
734
735 if (ref != NULL)
736 {
737 if (NULL == (indices_at = convert_reverse_indices (rev_indices, xpts_at->num_elements, hist_at)))
738 goto free_and_return;
739
740 if (-1 == SLang_assign_to_ref (ref, SLANG_ARRAY_TYPE, (VOID_STAR) &indices_at))
741 goto free_and_return;
742 }
743
744 (void) SLang_push_array (hist_at, 0);
745
746 /* NULLs ok below */
747 free_and_return:
748 if (rev_indices != NULL) SLfree ((char *) rev_indices);
749 SLang_free_ref (ref);
750 SLang_free_array (indices_at);
751 SLang_free_array (yedges_at);
752 SLang_free_array (xedges_at);
753 SLang_free_array (xpts_at);
754 SLang_free_array (ypts_at);
755 SLang_free_array (hist_at);
756 }
757
binary_search_intrin(void)758 static void binary_search_intrin (void)
759 {
760 SLang_Array_Type *at;
761 SLang_Array_Type *xt;
762 SLang_Array_Type *at_ind;
763 double x, *xp, *yp;
764 SLuindex_Type i, num_indices, num;
765 SLindex_Type ind, *indices;
766
767 if (SLang_Num_Function_Args != 2)
768 {
769 SLang_verror (SL_USAGE_ERROR, "i = hist_bsearch (x, a); %% a[i]<=x<a[i+1]");
770 return;
771 }
772
773 if (-1 == SLang_pop_array_of_type (&at, SLANG_DOUBLE_TYPE))
774 return;
775
776 if (SLang_peek_at_stack () == SLANG_ARRAY_TYPE)
777 {
778 if (-1 == SLang_pop_array_of_type (&xt, SLANG_DOUBLE_TYPE))
779 {
780 SLang_free_array (at);
781 return;
782 }
783
784 if (NULL == (at_ind = SLang_create_array (SLANG_INT_TYPE, 0, NULL, xt->dims, xt->num_dims)))
785 {
786 SLang_free_array (at);
787 SLang_free_array (xt);
788 return;
789 }
790 xp = (double *) xt->data;
791 num_indices = xt->num_elements;
792 indices = (SLindex_Type *)at_ind->data;
793 }
794 #if SLANG_VERSION < 20000
795 else if (0 == SLang_pop_double (&x, NULL, NULL))
796 #else
797 else if (0 == SLang_pop_double (&x))
798 #endif
799 {
800 xt = NULL;
801 at_ind = NULL;
802 xp = &x;
803 num_indices = 1;
804 indices = &ind;
805 }
806 else
807 {
808 SLang_free_array (at);
809 return;
810 }
811
812 num = at->num_elements;
813 yp = (double *)at->data;
814
815 if (-1 == check_grid (yp, num))
816 {
817 SLang_free_array (at);
818 SLang_free_array (xt);
819 return;
820 }
821
822 for (i = 0; i < num_indices; i++)
823 indices[i] = double_binary_search (xp[i], yp, num);
824
825 SLang_free_array (at);
826 SLang_free_array (xt);
827
828 if (at_ind != NULL)
829 (void) SLang_push_array (at_ind, 1);
830 else
831 (void) SLang_push_array_index (indices[0]);
832 }
833
834 /* Usage: h_new = rebin (new_grid, old_grid, h_old); */
hist1d_rebin(void)835 static void hist1d_rebin (void)
836 {
837 SLang_Array_Type *grid_old, *h_old, *grid_new, *h_new;
838 SLindex_Type old_num_bins, new_num_bins;
839
840 if (SLang_Num_Function_Args != 3)
841 {
842 SLang_verror (SL_USAGE_ERROR, "h_new = hist1d_rebin (new_grid, old_grid, h_old)");
843 return;
844 }
845
846 if (-1 == pop_1d_double_arrays (&grid_old, &h_old))
847 return;
848
849 if (-1 == pop_1d_array_of_type (&grid_new, SLANG_DOUBLE_TYPE))
850 {
851 SLang_free_array (h_old);
852 SLang_free_array (grid_old);
853 return;
854 }
855
856 old_num_bins = grid_old->num_elements;
857 new_num_bins = grid_new->num_elements;
858
859 if (NULL == (h_new = SLang_create_array (SLANG_DOUBLE_TYPE, 0, NULL, &new_num_bins, 1)))
860 {
861 SLang_free_array (grid_new);
862 SLang_free_array (grid_old);
863 SLang_free_array (h_old);
864 return;
865 }
866
867 if (0 == d_hist1d_rebin ((double *)grid_new->data, (SLuindex_Type) new_num_bins,
868 (double *)grid_old->data, (double *)h_old->data, (SLuindex_Type) old_num_bins,
869 (double *)h_new->data))
870 (void) SLang_push_array (h_new, 0);
871
872 SLang_free_array (h_new);
873 SLang_free_array (grid_new);
874 SLang_free_array (grid_old);
875 SLang_free_array (h_old);
876 }
877
878 #define A SLANG_ARRAY_TYPE
879 #define V SLANG_VOID_TYPE
880
881 static SLang_Intrin_Fun_Type Module_Intrinsics [] =
882 {
883 MAKE_INTRINSIC_0("hist1d", hist1d, V),
884 MAKE_INTRINSIC_0("hist2d", hist2d, V),
885 MAKE_INTRINSIC_0("hist_bsearch", binary_search_intrin, V),
886 MAKE_INTRINSIC_0("hist1d_rebin", hist1d_rebin, V),
887 SLANG_END_INTRIN_FUN_TABLE
888 };
889
890 static SLang_Intrin_Var_Type Module_Variables [] =
891 {
892 MAKE_VARIABLE("_histogram_module_version_string", &Module_Version_String, SLANG_STRING_TYPE, 1),
893 SLANG_END_INTRIN_VAR_TABLE
894 };
895
896 static SLang_IConstant_Type Module_IConstants [] =
897 {
898 MAKE_ICONSTANT("_histogram_module_version", MODULE_VERSION_NUMBER),
899 SLANG_END_ICONST_TABLE
900 };
901
init_histogram_module_ns(char * ns_name)902 int init_histogram_module_ns (char *ns_name)
903 {
904 SLang_NameSpace_Type *ns = SLns_create_namespace (ns_name);
905 if (ns == NULL)
906 return -1;
907
908 if (
909 (-1 == SLns_add_intrin_var_table (ns, Module_Variables, NULL))
910 || (-1 == SLns_add_intrin_fun_table (ns, Module_Intrinsics, NULL))
911 || (-1 == SLns_add_iconstant_table (ns, Module_IConstants, NULL))
912 )
913 return -1;
914
915 return 0;
916 }
917
918 /* This function is optional */
deinit_histogram_module(void)919 void deinit_histogram_module (void)
920 {
921 }
922