[gnuastro-commits] master 014057b2: Arithmetic: correcting collapse-number to not use median

[Mohammad Akhlaghi]

branch: master
commit 014057b23cb93b3ccdfcf3a04331d93d41e2f941
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Arithmetic: correcting collapse-number to not use median
    
    Until now, when the Arithmetic program was given the 'collapse-number'
    operator, it would mistakenly call the internal macro that is used for
    'collapse-median'! This originated in Commit 3aa904eca8c (from September
    2022 where the sigma-clipping collapse operators were added). Furthermore,
    while checking the sort-based collapse operators ('collapse-median' and all
    the clip-based collapses), I noticed they are also giving segmentation
    faults on an input which had a full row of NaN values.
    
    With this commit, the first problem is fixed by fixing the typo (and giving
    the correct macro is for 'collapse-number'). The second problem was caused
    by the new feature of the blank-checking operator (that would free the
    space when all elements were NaN). To avoid the problem, we now added a
    check immediately after the sort-based operators and if the array is freed,
    we'll allocate a new one for use in later iterations.
    
    The original problem (with 'collapse-number') was reported by Sepideh
    Eskandarlou.
    
    This fixes bug #65743.
---
 NEWS                        |  3 +++
 bin/arithmetic/arithmetic.c |  2 +-
 lib/data.c                  |  4 ++--
 lib/dimension.c             | 25 ++++++++++++++++++-------
 lib/statistics.c            |  6 +++---
 5 files changed, 27 insertions(+), 13 deletions(-)

diff --git a/NEWS b/NEWS
index 87566844..132347ef 100644
--- a/NEWS
+++ b/NEWS
@@ -181,6 +181,9 @@ See the end of the file for license conditions.
     even when the output is not FITS (plain-text or on the standard
     output). Reported by Sepideh Eskandarlou.
 
+  - bug #65743: Arithmetic's collapse-number mistakenly calculating
+    collapse-median instead. Reported by Sepideh Eskandarlou.
+
 
 
 
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 973e1564..64b8f04f 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -1513,7 +1513,7 @@ arithmetic_set_operator(char *string, size_t 
*num_operands, int *inlib)
       else if (!strcmp(string, "collapse-mean"))
         { op=ARITHMETIC_OP_COLLAPSE_MEAN;         *num_operands=0; }
       else if (!strcmp(string, "collapse-number"))
-        { op=ARITHMETIC_OP_COLLAPSE_MEDIAN;       *num_operands=0; }
+        { op=ARITHMETIC_OP_COLLAPSE_NUMBER;       *num_operands=0; }
       else if (!strcmp(string, "collapse-median"))
         { op=ARITHMETIC_OP_COLLAPSE_MEDIAN;       *num_operands=0; }
       else if (!strcmp(string, "collapse-madclip-mad"))
diff --git a/lib/data.c b/lib/data.c
index 3d7f1bae..8b62193b 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -81,8 +81,8 @@ gal_data_alloc(void *array, uint8_t type, size_t ndim, size_t 
*dsize,
           __func__, sizeof *out);
 
   /* Initialize the allocated array. */
-  gal_data_initialize(out, array, type, ndim, dsize, wcs, clear, minmapsize,
-                      quietmmap, name, unit, comment);
+  gal_data_initialize(out, array, type, ndim, dsize, wcs, clear,
+                      minmapsize, quietmmap, name, unit, comment);
 
   /* Return the final structure. */
   return out;
diff --git a/lib/dimension.c b/lib/dimension.c
index a5a13153..46898146 100644
--- a/lib/dimension.c
+++ b/lib/dimension.c
@@ -918,7 +918,7 @@ dimension_csb_copy(gal_data_t *in, size_t from, gal_data_t 
*work,
 static gal_data_t *
 dimension_collapse_sortbased_operation(struct dimension_sortbased_p *p,
                                        gal_data_t *work, uint8_t clipflags,
-                                       uint8_t isfill)
+                                       uint8_t isfill, size_t wdsize)
 {
   gal_data_t *out=NULL;
 
@@ -928,6 +928,7 @@ dimension_collapse_sortbased_operation(struct 
dimension_sortbased_p *p,
     case DIMENSION_COLLAPSE_MEDIAN:
       out=gal_statistics_median(work, 1);
       break;
+
     case DIMENSION_COLLAPSE_MADCLIP_MAD:
     case DIMENSION_COLLAPSE_MADCLIP_STD:
     case DIMENSION_COLLAPSE_MADCLIP_MEAN:
@@ -941,8 +942,8 @@ dimension_collapse_sortbased_operation(struct 
dimension_sortbased_p *p,
       /* When we are dealing with a clip-fill operator, the order of
          the elements matter, so we don't want it to be done inplace.*/
       out=gal_statistics_clip_mad(work, p->sclipmultip,
-                                   p->sclipparam, clipflags,
-                                   isfill?0:1, 1);
+                                  p->sclipparam, clipflags,
+                                  isfill?0:1, 1);
       break;
 
     case DIMENSION_COLLAPSE_SIGCLIP_MAD:
@@ -967,6 +968,13 @@ dimension_collapse_sortbased_operation(struct 
dimension_sortbased_p *p,
             "recognized", __func__, PACKAGE_BUGREPORT, p->operator);
     }
 
+  /* All the statistical operators that are based on sorting will free the
+     input array if it doesn't have any elements! But we are re-using the
+     'work' array many times. So we need to re-allocate the array here. */
+  if(work->array==NULL)
+    work->array=gal_pointer_allocate(work->type, wdsize, 0, __func__,
+                                     "work->array");
+
   /* Return the output. */
   return out;
 }
@@ -979,7 +987,7 @@ static gal_data_t *
 dimension_collapse_sortbased_fill(struct dimension_sortbased_p *p,
                                   gal_data_t *fstat, gal_data_t *work,
                                   gal_data_t *conv, uint8_t clipflags,
-                                  int check)
+                                  size_t wdsize, int check)
 {
   size_t one=1;
   int std1_mad0=0;
@@ -1087,7 +1095,8 @@ dimension_collapse_sortbased_fill(struct 
dimension_sortbased_p *p,
   /* Apply the operation: note that 'isfill' is zero here because we don't
      care about the order of elements any more ('isfill' is only used to do
      the operation inplace or not). */
-  out=dimension_collapse_sortbased_operation(p, work, clipflags, 0);
+  out=dimension_collapse_sortbased_operation(p, work, clipflags, 0,
+                                             wdsize);
 
   /* Clean up and return. */
   gal_data_free(tmp);
@@ -1289,13 +1298,15 @@ dimension_collapse_sortbased_worker(void *in_prm)
 
       /* Do the necessary statistical operation. */
       stat=dimension_collapse_sortbased_operation(p, conv?conv:work,
-                                                  clipflags, isfill);
+                                                  clipflags, isfill,
+                                                  wdsize);
 
       /* If this is a "filling" operation, then repeat the operation with
          the fill. */
       if(isfill)
         stat=dimension_collapse_sortbased_fill(p, stat, work, conv,
-                                               clipflags, index==-1);
+                                               clipflags, wdsize,
+                                               index==-1);
 
       /* Set the index in the output 'stat' array. These can't be set in
          the main operation 'switch' because the functions are different,
diff --git a/lib/statistics.c b/lib/statistics.c
index bd3bbbfc..acefdd28 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -1741,7 +1741,6 @@ gal_statistics_no_blank_sorted(gal_data_t *input, int 
inplace)
         }
       else contig=input;
 
-
       /* Make sure there are no blanks in the array that will be
          used. After this step, we won't be dealing with 'input' any more,
          but with 'noblank'. */
@@ -2829,11 +2828,12 @@ statistics_clip(gal_data_t *input, float multip, float 
param,
   /* Measure and report the remaining elements if requested. */
   if(extrastats) statistics_clip_stats_extra(nbs, oa, extrastats);
 
-  /* Clean up and return. */
+  /* Fix the 'array' pointer, clean up and return. */
   nbs->array=nbs_array;
   gal_data_free(center_i);
   gal_data_free(spread_i);
-  if(nbs!=input) gal_data_free(nbs);
+  if(nbs==input) input->array=nbs->array;
+  else           gal_data_free(nbs);
   return out;
 }