[gnuastro-commits] master 619b5be: Flag initialized every time in MakeCatalog's upperlimit

[Mohammad Akhlaghi]

branch: master
commit 619b5be17b8c8df8f7603ce109776ee459655f29
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Flag initialized every time in MakeCatalog's upperlimit
    
    To speed things up and not take too much memory and CPU resources, we ask
    the sigma-clipping of the values measured during the upper-limit
    calculations to be done in place. But we weren't actually initializing the
    sorting flag (necessary for the current implementation of sigma-clipping).
    Therefore, the sorting flag would be set for the random values of the first
    object in a thread, and the flag would persit for all the next objects
    (which were written in the same data container). This would create
    un-reliable values.
    
    With this commit, we initialize the sorting flag to zero before starting
    each object's upper-limit magnitude measurement. In the process some
    cleaning up was also done in other parts of the code for easy readability.
    
    This fixes bug #54188.
---
 NEWS                       |  1 +
 bin/mkcatalog/mkcatalog.c  | 21 ++++++++++++++++-----
 bin/mkcatalog/upperlimit.c |  6 ++++--
 doc/gnuastro.texi          |  6 ++++--
 lib/gnuastro/statistics.h  | 11 +----------
 lib/statistics.c           |  1 +
 6 files changed, 27 insertions(+), 19 deletions(-)

diff --git a/NEWS b/NEWS
index ef6d63d..28fbc74 100644
--- a/NEWS
+++ b/NEWS
@@ -21,6 +21,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
   bug #54057: Building failure due to not finding gsl_interp_steffen.
   bug #54063: Match tests in make check fail randomly.
   bug #54186: MakeCatalog's --checkupperlimit not keeping output's name.
+  bug #54188: MakeCatalog's Upperlimit not being sigma-clipped properly.
 
 
 
diff --git a/bin/mkcatalog/mkcatalog.c b/bin/mkcatalog/mkcatalog.c
index 724c430..fd41bee 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -125,18 +125,29 @@ mkcatalog_single_object(void *in_prm)
 
   /* If we have upper-limit mode, then allocate the container to keep the
      values to calculate the standard deviation. */
-  pp.up_vals = ( p->upperlimit
-                 ? gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &p->upnum,
+  if(p->upperlimit)
+    {
+      /* Allocate the space to keep the upper-limit values. */
+      pp.up_vals = gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &p->upnum,
                                   NULL, 0, p->cp.minmapsize, NULL, NULL,
-                                  NULL)
-                 : NULL );
+                                  NULL);
+
+      /* Set the blank checked flag to 1. By definition, this dataset won't
+         have any blank values. Also `flag' is initialized to `0'. So we
+         just have to set the checked flag (`GAL_DATA_FLAG_BLANK_CH') to
+         one to inform later steps that there are no blank values. */
+      pp.up_vals->flag |= GAL_DATA_FLAG_BLANK_CH;
+    }
+  else
+    pp.up_vals=NULL;
+
 
   /* Fill the desired columns for all the objects given to this thread. */
   for(i=0; tprm->indexs[i]!=GAL_BLANK_SIZE_T; ++i)
     {
       /* For easy reading. Note that the object IDs start from one while
          the array positions start from 0. */
-      pp.ci=NULL;
+      pp.ci     = NULL;
       pp.object = tprm->indexs[i] + 1;
       pp.tile   = &p->tiles[ tprm->indexs[i] ];
 
diff --git a/bin/mkcatalog/upperlimit.c b/bin/mkcatalog/upperlimit.c
index 406a66e..3645c2f 100644
--- a/bin/mkcatalog/upperlimit.c
+++ b/bin/mkcatalog/upperlimit.c
@@ -464,7 +464,8 @@ upperlimit_measure(struct mkcatalog_passparams *pp, int32_t 
clumplab,
             case UI_KEY_UPPERLIMITSIGMA:
             case UI_KEY_UPPERLIMITONESIGMA:
 
-              /* We only need to do this once. */
+              /* We only need to do this once, but the columns can be
+                 requested in any order. */
               if(sigclip==NULL)
                 {
                   /* Calculate the sigma-clipped standard deviation. Since
@@ -474,7 +475,7 @@ upperlimit_measure(struct mkcatalog_passparams *pp, int32_t 
clumplab,
                   init_size=pp->up_vals->size;
                   sigclip=gal_statistics_sigma_clip(pp->up_vals,
                                                     p->upsigmaclip[0],
-                                                    p->upsigmaclip[1], 1, 1);
+                                                    p->upsigmaclip[1],1,1);
                   pp->up_vals->size=pp->up_vals->dsize[0]=init_size;
                   scarr=sigclip->array;
 
@@ -573,6 +574,7 @@ upperlimit_one_tile(struct mkcatalog_passparams *pp, 
gal_data_t *tile,
   /* Initializations. */
   tarray=tile->array;
   gsl_rng_set(pp->rng, seed);
+  pp->up_vals->flag &= ~GAL_DATA_FLAG_SORT_CH;
 
 
   /* Set the range of random values for this tile. */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 45c734f..a28d000 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -26540,9 +26540,11 @@ is larger than @code{1} (one), it must be an integer 
and will be
 interpreted as the number clips to do. If it is less than @code{1} (one),
 it is interpreted as the tolerance level to stop the iteration.
 
-The output dataset has the following elements with a
address@hidden type:
+The returned dataset (let's call it @code{out}) contains a four-element
+array with type @code{GAL_TYPE_FLOAT32}. The final number of clips is
+stored in the @code{out->status}.
 @example
+float *array=out->array;
 array[0]: Number of points used.
 array[1]: Median.
 array[2]: Mean.
diff --git a/lib/gnuastro/statistics.h b/lib/gnuastro/statistics.h
index 92b7868..ad96d01 100644
--- a/lib/gnuastro/statistics.h
+++ b/lib/gnuastro/statistics.h
@@ -167,7 +167,7 @@ gal_statistics_cfp(gal_data_t *data, gal_data_t *bins, int 
normalize);
 
 
 /****************************************************************
- *****************        Sigma clip         ********************
+ *****************         Outliers          ********************
  ****************************************************************/
 gal_data_t *
 gal_statistics_sigma_clip(gal_data_t *input, float multip, float param,
@@ -177,15 +177,6 @@ gal_statistics_sigma_clip(gal_data_t *input, float multip, 
float param,
 
 
 
-/****************************************************************/
-/*************         Identify outliers         ****************/
-/****************************************************************/
-void
-gal_statistics_remove_outliers_flat_cdf(float *sorted, size_t *outsize);
-
-
-
-
 __END_C_DECLS    /* From C++ preparations */
 
 #endif           /* __GAL_STATISTICS_H__ */
diff --git a/lib/statistics.c b/lib/statistics.c
index a29303e..7d6e2f4 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -2068,6 +2068,7 @@ gal_statistics_sigma_clip(gal_data_t *input, float 
multip, float param,
 
 
 
+
 #if 0
 /* Using the cumulative distribution function this function will
    remove outliers from a dataset. */