[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 60cd76e 090/125: Statistics mostly complete
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 60cd76e 090/125: Statistics mostly complete |
|
Date: |
Sun, 23 Apr 2017 22:36:45 -0400 (EDT) |
branch: master
commit 60cd76e3c75fb815cce6212f8232a1f8be5958af
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Statistics mostly complete
The statistics program has been mostly complete, the mode finding algorithm
hasn't been applied yet though. Generally, Statistics has been mostly
re-written with mainly new options. The manual has also been updated to
fully describe the new added features.
General to all the programs:
- The `option_keys_enum' has been moved to the `ui.h' file instead of
`ui.c'. This was done because in some cases (like in Statistics), the
keys might be useful for other source files in the program too.
- A new function in `table.h' has been written to directly write the data
into a log file without the program having to worry about some details.
- The `Image manipulation' and `Image analysis' chapter titles of the book
have been changed to `Data manipulation' and `Data analysis'.
- The statistical operations returning one value have been moved from
`lib/arithmetic.c' to `lib/statistics.c' as separate functions. They are
still callable within the Arithmetic library and program, but the main
functions are now defined in `lib/statistics.c'.
- The blank functions were also cleaned up with better use of macros so
the `case' of each type now fully fits in one line, making them easier
to inspect/debug/improve.
---
bin/convertt/ui.c | 27 --
bin/convertt/ui.h | 33 ++
bin/convolve/ui.c | 28 --
bin/convolve/ui.h | 34 ++
bin/cosmiccal/ui.c | 23 -
bin/cosmiccal/ui.h | 28 ++
bin/crop/crop.c | 18 +-
bin/crop/ui.c | 38 --
bin/crop/ui.h | 43 ++
bin/header/ui.c | 25 -
bin/header/ui.h | 30 ++
bin/mknoise/ui.c | 19 -
bin/mknoise/ui.h | 24 +
bin/mkprof/mkprof.c | 12 +-
bin/mkprof/ui.c | 67 +--
bin/mkprof/ui.h | 58 +++
bin/statistics/args.h | 88 ++--
bin/statistics/aststatistics.conf | 6 +
bin/statistics/main.h | 12 +-
bin/statistics/statistics.c | 678 +++++++++++++++++---------
bin/statistics/ui.c | 355 +++++++++-----
bin/statistics/ui.h | 46 ++
bin/table/ui.c | 31 +-
bin/table/ui.h | 22 +
bin/warp/ui.c | 29 --
bin/warp/ui.h | 34 ++
doc/gnuastro.texi | 836 ++++++++++++++++++--------------
lib/arithmetic.c | 315 ++----------
lib/blank.c | 464 ++++++++----------
lib/commonopts.h | 2 +-
lib/data.c | 6 +
lib/fits.c | 55 +--
lib/gnuastro/fits.h | 7 +-
lib/gnuastro/statistics.h | 54 +--
lib/gnuastro/table.h | 15 +-
lib/gnuastro/txt.h | 4 +-
lib/mode.c | 2 +-
lib/options.c | 111 ++---
lib/options.h | 5 -
lib/statistics.c | 984 ++++++++++++++++++++++++++++----------
lib/table.c | 82 +++-
lib/txt.c | 12 +-
tests/statistics/basicstats.sh | 4 +-
43 files changed, 2761 insertions(+), 2005 deletions(-)
diff --git a/bin/convertt/ui.c b/bin/convertt/ui.c
index 9d6ec6a..5a2298c 100644
--- a/bin/convertt/ui.c
+++ b/bin/convertt/ui.c
@@ -87,33 +87,6 @@ enum program_args_groups
-/* Available letters for short options:
-
- d e f g j k l n p r s t v y z
- A B E F G I J M O Q R T U W X Y Z */
-enum option_keys_enum
-{
- /* With short-option version. */
- ARGS_OPTION_KEY_QUALITY = 'u',
- ARGS_OPTION_KEY_WIDTHINCM = 'w',
- ARGS_OPTION_KEY_BORDERWIDTH = 'b',
- ARGS_OPTION_KEY_HEX = 'x',
- ARGS_OPTION_KEY_FLUXLOW = 'L',
- ARGS_OPTION_KEY_FLUXHIGH = 'H',
- ARGS_OPTION_KEY_HIGH = 'H',
- ARGS_OPTION_KEY_MAXBYTE = 'm',
- ARGS_OPTION_KEY_FLMINBYTE = 'a',
- ARGS_OPTION_KEY_FHMAXBYTE = 'b',
- ARGS_OPTION_KEY_CHANGE = 'c',
- ARGS_OPTION_KEY_CHANGEAFTERTRUNC = 'C',
- ARGS_OPTION_KEY_INVERT = 'i',
-
- /* Only with long version (start with a value 1000, the rest will be set
- automatically). */
-};
-
-
-
diff --git a/bin/convertt/ui.h b/bin/convertt/ui.h
index 442770c..0cf6df5 100644
--- a/bin/convertt/ui.h
+++ b/bin/convertt/ui.h
@@ -23,6 +23,39 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#ifndef UI_H
#define UI_H
+
+
+
+
+/* Available letters for short options:
+
+ d e f g j k l n p r s t v y z
+ A B E F G I J M O Q R T U W X Y Z */
+enum option_keys_enum
+{
+ /* With short-option version. */
+ ARGS_OPTION_KEY_QUALITY = 'u',
+ ARGS_OPTION_KEY_WIDTHINCM = 'w',
+ ARGS_OPTION_KEY_BORDERWIDTH = 'b',
+ ARGS_OPTION_KEY_HEX = 'x',
+ ARGS_OPTION_KEY_FLUXLOW = 'L',
+ ARGS_OPTION_KEY_FLUXHIGH = 'H',
+ ARGS_OPTION_KEY_HIGH = 'H',
+ ARGS_OPTION_KEY_MAXBYTE = 'm',
+ ARGS_OPTION_KEY_FLMINBYTE = 'a',
+ ARGS_OPTION_KEY_FHMAXBYTE = 'b',
+ ARGS_OPTION_KEY_CHANGE = 'c',
+ ARGS_OPTION_KEY_CHANGEAFTERTRUNC = 'C',
+ ARGS_OPTION_KEY_INVERT = 'i',
+
+ /* Only with long version (start with a value 1000, the rest will be set
+ automatically). */
+};
+
+
+
+
+
void
ui_read_check_inputs_setup(int argc, char *argv[], struct converttparams *p);
diff --git a/bin/convolve/ui.c b/bin/convolve/ui.c
index 12b4edd..3f1d54e 100644
--- a/bin/convolve/ui.c
+++ b/bin/convolve/ui.c
@@ -89,34 +89,6 @@ enum program_args_groups
-/* Available letters for short options:
-
- e f g i j l n r p t u v w x y z
- A B E F G H I J M O Q R T W X Y Z */
-enum option_keys_enum
-{
- /* With short-option version. */
- ARGS_OPTION_KEY_KERNEL = 'k',
- ARGS_OPTION_KEY_KHDU = 'U',
- ARGS_OPTION_KEY_MINSHARPSPEC = 'c',
- ARGS_OPTION_KEY_CHECKFREQSTEPS = 'C',
- ARGS_OPTION_KEY_MESHSIZE = 'c',
- ARGS_OPTION_KEY_NCH1 = 'a',
- ARGS_OPTION_KEY_NCH2 = 'b',
- ARGS_OPTION_KEY_LASTMESHFRAC = 'L',
- ARGS_OPTION_KEY_DOMAIN = 'd',
- ARGS_OPTION_KEY_MAKEKERNEL = 'm',
-
- /* Only with long version (start with a value 1000, the rest will be set
- automatically). */
- ARGS_OPTION_KEY_NOKERNELFLIP = 1000,
- ARGS_OPTION_KEY_NOKERNELNORM,
- ARGS_OPTION_KEY_CHECKMESH,
- ARGS_OPTION_KEY_FULLCONVOLUTION,
-};
-
-
-
diff --git a/bin/convolve/ui.h b/bin/convolve/ui.h
index 937df9c..f5f5616 100644
--- a/bin/convolve/ui.h
+++ b/bin/convolve/ui.h
@@ -23,6 +23,40 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#ifndef UI_H
#define UI_H
+
+
+
+
+/* Available letters for short options:
+
+ e f g i j l n r p t u v w x y z
+ A B E F G H I J M O Q R T W X Y Z */
+enum option_keys_enum
+{
+ /* With short-option version. */
+ ARGS_OPTION_KEY_KERNEL = 'k',
+ ARGS_OPTION_KEY_KHDU = 'U',
+ ARGS_OPTION_KEY_MINSHARPSPEC = 'c',
+ ARGS_OPTION_KEY_CHECKFREQSTEPS = 'C',
+ ARGS_OPTION_KEY_MESHSIZE = 'c',
+ ARGS_OPTION_KEY_NCH1 = 'a',
+ ARGS_OPTION_KEY_NCH2 = 'b',
+ ARGS_OPTION_KEY_LASTMESHFRAC = 'L',
+ ARGS_OPTION_KEY_DOMAIN = 'd',
+ ARGS_OPTION_KEY_MAKEKERNEL = 'm',
+
+ /* Only with long version (start with a value 1000, the rest will be set
+ automatically). */
+ ARGS_OPTION_KEY_NOKERNELFLIP = 1000,
+ ARGS_OPTION_KEY_NOKERNELNORM,
+ ARGS_OPTION_KEY_CHECKMESH,
+ ARGS_OPTION_KEY_FULLCONVOLUTION,
+};
+
+
+
+
+
void
ui_read_check_inputs_setup(int argc, char *argv[], struct convolveparams *p);
diff --git a/bin/cosmiccal/ui.c b/bin/cosmiccal/ui.c
index fe4aa75..92b17ab 100644
--- a/bin/cosmiccal/ui.c
+++ b/bin/cosmiccal/ui.c
@@ -74,29 +74,6 @@ doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" will do
cosmological "
-/* Available letters for short options:
-
- a b d e f g j k l m n p r s t u v w x y z
- A B C E F G H J L M O Q R T U W X Y Z */
-enum option_keys_enum
-{
- /* With short-option version. */
- ARGS_OPTION_KEY_REDSHIFT = 'z',
- ARGS_OPTION_KEY_H0 = 'H',
- ARGS_OPTION_KEY_OLAMBDA = 'l',
- ARGS_OPTION_KEY_OMATTER = 'm',
- ARGS_OPTION_KEY_ORADIATION = 'r',
- ARGS_OPTION_KEY_ONLYVOLUME = 'v',
- ARGS_OPTION_KEY_ONLYABSMAGCONV = 'a',
-
-
- /* Only with long version (start with a value 1000, the rest will be set
- automatically). */
-};
-
-
-
-
diff --git a/bin/cosmiccal/ui.h b/bin/cosmiccal/ui.h
index 3667e11..42ff358 100644
--- a/bin/cosmiccal/ui.h
+++ b/bin/cosmiccal/ui.h
@@ -23,6 +23,34 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#ifndef UI_H
#define UI_H
+
+
+
+
+/* Available letters for short options:
+
+ a b d e f g j k l m n p r s t u v w x y z
+ A B C E F G H J L M O Q R T U W X Y Z */
+enum option_keys_enum
+{
+ /* With short-option version. */
+ ARGS_OPTION_KEY_REDSHIFT = 'z',
+ ARGS_OPTION_KEY_H0 = 'H',
+ ARGS_OPTION_KEY_OLAMBDA = 'l',
+ ARGS_OPTION_KEY_OMATTER = 'm',
+ ARGS_OPTION_KEY_ORADIATION = 'r',
+ ARGS_OPTION_KEY_ONLYVOLUME = 'v',
+ ARGS_OPTION_KEY_ONLYABSMAGCONV = 'a',
+
+
+ /* Only with long version (start with a value 1000, the rest will be set
+ automatically). */
+};
+
+
+
+
+
void
ui_read_check_inputs_setup(int argc, char *argv[],
struct cosmiccalparams *p);
diff --git a/bin/crop/crop.c b/bin/crop/crop.c
index 7226c7d..0c04a83 100644
--- a/bin/crop/crop.c
+++ b/bin/crop/crop.c
@@ -376,7 +376,7 @@ void
crop(struct cropparams *p)
{
int err=0;
- char *comments;
+ char *tmp;
pthread_t t; /* We don't use the thread id, so all are saved here. */
pthread_attr_t attr;
pthread_barrier_t b;
@@ -384,6 +384,7 @@ crop(struct cropparams *p)
size_t i, *indexs, thrdcols;
size_t nt=p->cp.numthreads, nb;
void *(*modefunction)(void *)=NULL;
+ struct gal_linkedlist_stll *comments=NULL;
/* Set the function to run: */
@@ -446,13 +447,14 @@ crop(struct cropparams *p)
if(p->cp.log)
{
if(p->checkcenter)
- asprintf(&comments, "# Width of central check box: %zu\n#",
- p->checkcenter);
- else
- comments=NULL;
- gal_options_print_log(p->log, PROGRAM_STRING, &p->rawtime, comments,
- LOGFILENAME, &p->cp);
- if(comments) free(comments);
+ {
+ asprintf(&tmp, "Width of central check box: %zu",
+ p->checkcenter);
+ gal_linkedlist_add_to_stll(&comments, tmp, 0);
+ }
+ gal_table_write_log(p->log, PROGRAM_STRING, &p->rawtime, comments,
+ LOGFILENAME, p->cp.dontdelete, p->cp.quiet);
+ gal_linkedlist_free_stll(comments, 1);
}
/* Print the final verbose info, save log, and clean up: */
diff --git a/bin/crop/ui.c b/bin/crop/ui.c
index da5038f..021445c 100644
--- a/bin/crop/ui.c
+++ b/bin/crop/ui.c
@@ -94,44 +94,6 @@ enum program_args_groups
-/* Available letters for short options:
-
- e k m t u v
- A B E F G H I J L O Q R T U W X Y Z */
-enum option_keys_enum
-{
- /* With short-option version. */
- ARGS_OPTION_KEY_CATALOG = 'C',
- ARGS_OPTION_KEY_NOBLANK = 'b',
- ARGS_OPTION_KEY_CHECKCENTER = 'c',
- ARGS_OPTION_KEY_SUFFIX = 'p',
- ARGS_OPTION_KEY_NAMECOL = 'n',
- ARGS_OPTION_KEY_RACOL = 'f',
- ARGS_OPTION_KEY_DECCOL = 'g',
- ARGS_OPTION_KEY_RA = 'r',
- ARGS_OPTION_KEY_DEC = 'd',
- ARGS_OPTION_KEY_XCOL = 'i',
- ARGS_OPTION_KEY_YCOL = 'j',
- ARGS_OPTION_KEY_XC = 'x',
- ARGS_OPTION_KEY_YC = 'y',
- ARGS_OPTION_KEY_IWIDTH = 'a',
- ARGS_OPTION_KEY_WWIDTH = 'w',
- ARGS_OPTION_KEY_SECTION = 's',
- ARGS_OPTION_KEY_POLYGON = 'l',
- ARGS_OPTION_KEY_ZEROISNOTBLANK = 'z',
- ARGS_OPTION_KEY_MODE = 'M',
-
- /* Only with long version (start with a value 1000, the rest will be set
- automatically). */
- ARGS_OPTION_KEY_CATHDU = 1000,
- ARGS_OPTION_KEY_HSTARTWCS,
- ARGS_OPTION_KEY_HENDWCS,
- ARGS_OPTION_KEY_OUTPOLYGON,
-};
-
-
-
-
diff --git a/bin/crop/ui.h b/bin/crop/ui.h
index 96cf34a..e4c4226 100644
--- a/bin/crop/ui.h
+++ b/bin/crop/ui.h
@@ -23,6 +23,49 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#ifndef UI_H
#define UI_H
+
+
+
+
+/* Available letters for short options:
+
+ e k m t u v
+ A B E F G H I J L O Q R T U W X Y Z */
+enum option_keys_enum
+{
+ /* With short-option version. */
+ ARGS_OPTION_KEY_CATALOG = 'C',
+ ARGS_OPTION_KEY_NOBLANK = 'b',
+ ARGS_OPTION_KEY_CHECKCENTER = 'c',
+ ARGS_OPTION_KEY_SUFFIX = 'p',
+ ARGS_OPTION_KEY_NAMECOL = 'n',
+ ARGS_OPTION_KEY_RACOL = 'f',
+ ARGS_OPTION_KEY_DECCOL = 'g',
+ ARGS_OPTION_KEY_RA = 'r',
+ ARGS_OPTION_KEY_DEC = 'd',
+ ARGS_OPTION_KEY_XCOL = 'i',
+ ARGS_OPTION_KEY_YCOL = 'j',
+ ARGS_OPTION_KEY_XC = 'x',
+ ARGS_OPTION_KEY_YC = 'y',
+ ARGS_OPTION_KEY_IWIDTH = 'a',
+ ARGS_OPTION_KEY_WWIDTH = 'w',
+ ARGS_OPTION_KEY_SECTION = 's',
+ ARGS_OPTION_KEY_POLYGON = 'l',
+ ARGS_OPTION_KEY_ZEROISNOTBLANK = 'z',
+ ARGS_OPTION_KEY_MODE = 'M',
+
+ /* Only with long version (start with a value 1000, the rest will be set
+ automatically). */
+ ARGS_OPTION_KEY_CATHDU = 1000,
+ ARGS_OPTION_KEY_HSTARTWCS,
+ ARGS_OPTION_KEY_HENDWCS,
+ ARGS_OPTION_KEY_OUTPOLYGON,
+};
+
+
+
+
+
void
ui_read_check_inputs_setup(int argc, char *argv[], struct cropparams *p);
diff --git a/bin/header/ui.c b/bin/header/ui.c
index eaf5dce..7c1f4c0 100644
--- a/bin/header/ui.c
+++ b/bin/header/ui.c
@@ -73,31 +73,6 @@ doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" print the
header "
-/* Available letters for short options:
-
- a b d e f g j k l m n p r s t u v w x y z
- A B C E F G H J L M O Q R T U W X Y Z */
-enum option_keys_enum
-{
- /* With short-option version. */
- ARGS_OPTION_KEY_ASIS = 'a',
- ARGS_OPTION_KEY_DELETE = 'd',
- ARGS_OPTION_KEY_RENAME = 'r',
- ARGS_OPTION_KEY_UPDATE = 'u',
- ARGS_OPTION_KEY_WRITE = 'w',
- ARGS_OPTION_KEY_COMMENT = 'c',
- ARGS_OPTION_KEY_HISTORY = 'h',
- ARGS_OPTION_KEY_DATE = 't',
- ARGS_OPTION_KEY_QUITONERROR = 'Q',
-
-
- /* Only with long version (start with a value 1000, the rest will be set
- automatically). */
-};
-
-
-
-
diff --git a/bin/header/ui.h b/bin/header/ui.h
index 2fef561..6936149 100644
--- a/bin/header/ui.h
+++ b/bin/header/ui.h
@@ -23,6 +23,36 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#ifndef UI_H
#define UI_H
+
+
+
+
+/* Available letters for short options:
+
+ a b d e f g j k l m n p r s t u v w x y z
+ A B C E F G H J L M O Q R T U W X Y Z */
+enum option_keys_enum
+{
+ /* With short-option version. */
+ ARGS_OPTION_KEY_ASIS = 'a',
+ ARGS_OPTION_KEY_DELETE = 'd',
+ ARGS_OPTION_KEY_RENAME = 'r',
+ ARGS_OPTION_KEY_UPDATE = 'u',
+ ARGS_OPTION_KEY_WRITE = 'w',
+ ARGS_OPTION_KEY_COMMENT = 'c',
+ ARGS_OPTION_KEY_HISTORY = 'h',
+ ARGS_OPTION_KEY_DATE = 't',
+ ARGS_OPTION_KEY_QUITONERROR = 'Q',
+
+
+ /* Only with long version (start with a value 1000, the rest will be set
+ automatically). */
+};
+
+
+
+
+
void
ui_read_check_inputs_setup(int argc, char *argv[],
struct headerparams *p);
diff --git a/bin/mknoise/ui.c b/bin/mknoise/ui.c
index 6221182..077622a 100644
--- a/bin/mknoise/ui.c
+++ b/bin/mknoise/ui.c
@@ -74,25 +74,6 @@ doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" will add
noise to all the "
-/* Available letters for short options:
-
- a c f g i j k l m n p r t u v w x y
- A B C E F G H I J L M O Q R T U W X Y Z */
-enum option_keys_enum
-{
- /* With short-option version. */
- ARGS_OPTION_KEY_STDADD = 's',
- ARGS_OPTION_KEY_BACKGROUND = 'b',
- ARGS_OPTION_KEY_ZEROPOINT = 'z',
- ARGS_OPTION_KEY_ENVSEED = 'e',
-
- /* Only with long version (start with a value 1000, the rest will be set
- automatically). */
-};
-
-
-
-
diff --git a/bin/mknoise/ui.h b/bin/mknoise/ui.h
index 1e8731c..070237d 100644
--- a/bin/mknoise/ui.h
+++ b/bin/mknoise/ui.h
@@ -23,6 +23,30 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#ifndef UI_H
#define UI_H
+
+
+
+
+/* Available letters for short options:
+
+ a c f g i j k l m n p r t u v w x y
+ A B C E F G H I J L M O Q R T U W X Y Z */
+enum option_keys_enum
+{
+ /* With short-option version. */
+ ARGS_OPTION_KEY_STDADD = 's',
+ ARGS_OPTION_KEY_BACKGROUND = 'b',
+ ARGS_OPTION_KEY_ZEROPOINT = 'z',
+ ARGS_OPTION_KEY_ENVSEED = 'e',
+
+ /* Only with long version (start with a value 1000, the rest will be set
+ automatically). */
+};
+
+
+
+
+
void
ui_read_check_inputs_setup(int argc, char *argv[], struct mknoiseparams *p);
diff --git a/bin/mkprof/mkprof.c b/bin/mkprof/mkprof.c
index ecdaf44..6b5d77b 100644
--- a/bin/mkprof/mkprof.c
+++ b/bin/mkprof/mkprof.c
@@ -568,14 +568,15 @@ void
mkprof(struct mkprofparams *p)
{
int err;
+ char *tmp;
pthread_t t; /* Thread id not used, all are saved here. */
- char *comments;
pthread_attr_t attr;
pthread_barrier_t b;
struct mkonthread *mkp;
size_t i, *indexs, thrdcols;
size_t nt=p->cp.numthreads, nb;
long onaxes[2], os=p->oversample;
+ struct gal_linkedlist_stll *comments=NULL;
/* Allocate the arrays to keep the thread and parameters for each
@@ -642,10 +643,11 @@ mkprof(struct mkprofparams *p)
/* Write the log file. */
if(p->cp.log)
{
- asprintf(&comments, "# Zeropoint: %.4f", p->zeropoint);
- gal_options_print_log(p->log, PROGRAM_STRING, &p->rawtime, comments,
- LOGFILENAME, &p->cp);
- free(comments);
+ asprintf(&tmp, "Zeropoint: %g", p->zeropoint);
+ gal_linkedlist_add_to_stll(&comments, tmp, 0);
+ gal_table_write_log(p->log, PROGRAM_STRING, &p->rawtime, comments,
+ LOGFILENAME, p->cp.dontdelete, p->cp.quiet);
+ gal_linkedlist_free_stll(comments, 1);
}
/* If numthreads>1, then wait for all the jobs to finish and destroy
diff --git a/bin/mkprof/ui.c b/bin/mkprof/ui.c
index f5fb997..f763128 100644
--- a/bin/mkprof/ui.c
+++ b/bin/mkprof/ui.c
@@ -91,60 +91,6 @@ enum program_args_groups
-/* Keys for each option.
-
- Available letters (-V which is used by GNU is also removed):
-
- a d f g j l u v
- A E G H I J L M O Q U W Z */
-enum option_keys_enum
-{
- /* With short-option version. */
- ARGS_OPTION_KEY_BACKGROUND = 'k',
- ARGS_OPTION_KEY_BACKHDU = 'B',
- ARGS_OPTION_KEY_NAXIS1 = 'x',
- ARGS_OPTION_KEY_NAXIS2 = 'y',
- ARGS_OPTION_KEY_CLEARCANVAS = 'C',
- ARGS_OPTION_KEY_OVERSAMPLE = 's',
- ARGS_OPTION_KEY_INDIVIDUAL = 'i',
- ARGS_OPTION_KEY_NOMERGED = 'm',
- ARGS_OPTION_KEY_NUMRANDOM = 'r',
- ARGS_OPTION_KEY_TOLERANCE = 't',
- ARGS_OPTION_KEY_TUNITINP = 'p',
- ARGS_OPTION_KEY_XSHIFT = 'X',
- ARGS_OPTION_KEY_YSHIFT = 'Y',
- ARGS_OPTION_KEY_PREPFORCONV = 'c',
- ARGS_OPTION_KEY_ZEROPOINT = 'z',
- ARGS_OPTION_KEY_CIRCUMWIDTH = 'w',
- ARGS_OPTION_KEY_REPLACE = 'R',
- ARGS_OPTION_KEY_ENVSEED = 'e',
- ARGS_OPTION_KEY_MFORFLATPIX = 'F',
-
- /* Only with long version. */
- ARGS_OPTION_KEY_PSFINIMG = 1000,
- ARGS_OPTION_KEY_MAGATPEAK,
- ARGS_OPTION_KEY_XCOL,
- ARGS_OPTION_KEY_YCOL,
- ARGS_OPTION_KEY_RACOL,
- ARGS_OPTION_KEY_DECCOL,
- ARGS_OPTION_KEY_FCOL,
- ARGS_OPTION_KEY_RCOL,
- ARGS_OPTION_KEY_NCOL,
- ARGS_OPTION_KEY_PCOL,
- ARGS_OPTION_KEY_QCOL,
- ARGS_OPTION_KEY_MCOL,
- ARGS_OPTION_KEY_TCOL,
- ARGS_OPTION_KEY_CRPIX1,
- ARGS_OPTION_KEY_CRPIX2,
- ARGS_OPTION_KEY_CRVAL1,
- ARGS_OPTION_KEY_CRVAL2,
- ARGS_OPTION_KEY_RESOLUTION,
-};
-
-
-
-
-
@@ -846,7 +792,7 @@ ui_finalize_coordinates(struct mkprofparams *p)
static void
ui_make_log(struct mkprofparams *p)
{
- char *comment;
+ char *name, *comment;
/* Return if no long file is to be created. */
if(p->cp.log==0) return;
@@ -873,15 +819,12 @@ ui_make_log(struct mkprofparams *p)
"Magnitude of profile's overlap with merged image.");
/* Row number in input catalog. */
- if(gal_fits_name_is_fits(p->catname))
- asprintf(&comment, "Row number of profile in %s (hdu: %s).",
- p->catname, p->cp.hdu);
- else
- asprintf(&comment, "Row number of profile in %s.", p->catname);
+ name=gal_fits_name_save_as_string(p->catname, p->cp.hdu);
+ asprintf(&comment, "Row number of profile in %s.", name);
gal_data_add_to_ll(&p->log, NULL, GAL_DATA_TYPE_UINT64, 1, &p->num, NULL,
- 1, p->cp.minmapsize, "INPUT_ROW_NO", "count",
- "Row number of profile in ");
+ 1, p->cp.minmapsize, "INPUT_ROW_NO", "count", comment);
free(comment);
+ free(name);
}
diff --git a/bin/mkprof/ui.h b/bin/mkprof/ui.h
index c986541..b166599 100644
--- a/bin/mkprof/ui.h
+++ b/bin/mkprof/ui.h
@@ -23,6 +23,64 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#ifndef UI_H
#define UI_H
+
+
+
+
+/* Keys for each option.
+
+ Available letters (-V which is used by GNU is also removed):
+
+ a d f g j l u v
+ A E G H I J L M O Q U W Z */
+enum option_keys_enum
+{
+ /* With short-option version. */
+ ARGS_OPTION_KEY_BACKGROUND = 'k',
+ ARGS_OPTION_KEY_BACKHDU = 'B',
+ ARGS_OPTION_KEY_NAXIS1 = 'x',
+ ARGS_OPTION_KEY_NAXIS2 = 'y',
+ ARGS_OPTION_KEY_CLEARCANVAS = 'C',
+ ARGS_OPTION_KEY_OVERSAMPLE = 's',
+ ARGS_OPTION_KEY_INDIVIDUAL = 'i',
+ ARGS_OPTION_KEY_NOMERGED = 'm',
+ ARGS_OPTION_KEY_NUMRANDOM = 'r',
+ ARGS_OPTION_KEY_TOLERANCE = 't',
+ ARGS_OPTION_KEY_TUNITINP = 'p',
+ ARGS_OPTION_KEY_XSHIFT = 'X',
+ ARGS_OPTION_KEY_YSHIFT = 'Y',
+ ARGS_OPTION_KEY_PREPFORCONV = 'c',
+ ARGS_OPTION_KEY_ZEROPOINT = 'z',
+ ARGS_OPTION_KEY_CIRCUMWIDTH = 'w',
+ ARGS_OPTION_KEY_REPLACE = 'R',
+ ARGS_OPTION_KEY_ENVSEED = 'e',
+ ARGS_OPTION_KEY_MFORFLATPIX = 'F',
+
+ /* Only with long version. */
+ ARGS_OPTION_KEY_PSFINIMG = 1000,
+ ARGS_OPTION_KEY_MAGATPEAK,
+ ARGS_OPTION_KEY_XCOL,
+ ARGS_OPTION_KEY_YCOL,
+ ARGS_OPTION_KEY_RACOL,
+ ARGS_OPTION_KEY_DECCOL,
+ ARGS_OPTION_KEY_FCOL,
+ ARGS_OPTION_KEY_RCOL,
+ ARGS_OPTION_KEY_NCOL,
+ ARGS_OPTION_KEY_PCOL,
+ ARGS_OPTION_KEY_QCOL,
+ ARGS_OPTION_KEY_MCOL,
+ ARGS_OPTION_KEY_TCOL,
+ ARGS_OPTION_KEY_CRPIX1,
+ ARGS_OPTION_KEY_CRPIX2,
+ ARGS_OPTION_KEY_CRVAL1,
+ ARGS_OPTION_KEY_CRVAL2,
+ ARGS_OPTION_KEY_RESOLUTION,
+};
+
+
+
+
+
void
ui_read_check_inputs_setup(int argc, char *argv[], struct mkprofparams *p);
diff --git a/bin/statistics/args.h b/bin/statistics/args.h
index 8b0e848..7ae18c3 100644
--- a/bin/statistics/args.h
+++ b/bin/statistics/args.h
@@ -71,17 +71,18 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
- "quantrange",
- ARGS_OPTION_KEY_QUANTRANGE,
- "FLT",
+ "qrange",
+ ARGS_OPTION_KEY_QRANGE,
+ "FLT[,FLT]",
0,
- "Quantile (Q) range: from Q to 1-Q.",
+ "Quantile range: one (from Q to 1-Q) or two.",
GAL_OPTIONS_GROUP_INPUT,
- &p->quantilerange,
- GAL_DATA_TYPE_FLOAT32,
+ NULL,
+ GAL_DATA_TYPE_STRING,
GAL_OPTIONS_RANGE_ANY,
GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
+ GAL_OPTIONS_NOT_SET,
+ ui_parse_numbers
},
@@ -90,7 +91,7 @@ struct argp_option program_options[] =
{
0, 0, 0, 0,
- "Values to print in one row (on command-line)",
+ "Single value measurements (to print in one row)",
ARGS_GROUP_IN_ONE_ROW
},
{
@@ -259,7 +260,7 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_CUMULATIVE,
0,
0,
- "Save the cumulative frequency plot.",
+ "Save the cumulative frequency plot in output.",
ARGS_GROUP_PARTICULAR_STAT,
&p->cumulative,
GAL_OPTIONS_NO_ARG_TYPE,
@@ -272,26 +273,14 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_SIGMACLIP,
"FLT,FLT",
0,
- "Multiple of sigma and tolerance or number.",
+ "Multiple of sigma and tolerance/number.",
ARGS_GROUP_PARTICULAR_STAT,
- &p->sigclipstr,
+ NULL,
GAL_DATA_TYPE_STRING,
GAL_OPTIONS_RANGE_ANY,
GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
- {
- "mirror",
- ARGS_OPTION_KEY_MIRROR,
- "FLT",
- 0,
- "Quantile of mirror distribution to save.",
- ARGS_GROUP_PARTICULAR_STAT,
- &p->mirrorquant,
- GAL_DATA_TYPE_FLOAT32,
- GAL_OPTIONS_RANGE_GT_0_LT_1,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
+ GAL_OPTIONS_NOT_SET,
+ ui_parse_numbers
},
@@ -317,14 +306,27 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
- "lowerbin",
- ARGS_OPTION_KEY_LOWERBIN,
+ "numasciibins",
+ ARGS_OPTION_KEY_NUMASCIIBINS,
+ "INT",
0,
+ "Number of bins in ASCII histogram or CFP plots.",
+ ARGS_GROUP_HIST_CFP,
+ &p->numasciibins,
+ GAL_DATA_TYPE_SIZE_T,
+ GAL_OPTIONS_RANGE_GT_0,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "asciiheight",
+ ARGS_OPTION_KEY_ASCIIHEIGHT,
+ "INT",
0,
- "Save interval lower limit, not center",
+ "Height of ASCII histogram or CFP plots.",
ARGS_GROUP_HIST_CFP,
- &p->lowerbin,
- GAL_OPTIONS_NO_ARG_TYPE,
+ &p->asciiheight,
+ GAL_DATA_TYPE_SIZE_T,
GAL_OPTIONS_RANGE_GT_0,
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET
@@ -343,19 +345,6 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
- "onebinstart",
- ARGS_OPTION_KEY_ONEBINSTART,
- "FLT",
- 0,
- "Shift bins so one bin starts on this value.",
- ARGS_GROUP_HIST_CFP,
- &p->onebinstart,
- GAL_DATA_TYPE_FLOAT32,
- GAL_OPTIONS_RANGE_ANY,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
- {
"maxbinone",
ARGS_OPTION_KEY_MAXBINONE,
0,
@@ -368,6 +357,19 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET
},
+ {
+ "onebinstart",
+ ARGS_OPTION_KEY_ONEBINSTART,
+ "FLT",
+ 0,
+ "Shift bins so one bin starts on this value.",
+ ARGS_GROUP_HIST_CFP,
+ &p->onebinstart,
+ GAL_DATA_TYPE_FLOAT32,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
{0}
diff --git a/bin/statistics/aststatistics.conf
b/bin/statistics/aststatistics.conf
index dd22c54..f6a6c46 100644
--- a/bin/statistics/aststatistics.conf
+++ b/bin/statistics/aststatistics.conf
@@ -20,6 +20,12 @@
# Input image:
hdu 0
+# Histogram and CFP settings
+ numasciibins 70
+ asciiheight 10
+ numbins 100
+
+
# Common options
minmapsize 1000000000
tableformat fits-binary
diff --git a/bin/statistics/main.h b/bin/statistics/main.h
index 3e26f67..af87ccc 100644
--- a/bin/statistics/main.h
+++ b/bin/statistics/main.h
@@ -49,24 +49,26 @@ struct statisticsparams
char *column; /* Column name or number if input is table. */
float greaterequal; /* Only use values >= this value. */
float lessthan; /* Only use values < this value. */
- float quantilerange; /* Quantile (Q) range: from Q to 1-Q. */
+ float quantmin; /* Quantile min or range: from Q to 1-Q. */
+ float quantmax; /* Quantile maximum. */
uint8_t asciihist; /* Print an ASCII histogram. */
uint8_t asciicfp; /* Print an ASCII cumulative frequency plot.*/
uint8_t histogram; /* Save histogram in output. */
uint8_t cumulative; /* Save cumulative distibution in output. */
- char *sigclipstr; /* Multiple of sigma, and tolerance or num. */
- float mirrorquant; /* Quantile of mirror distribution to save. */
+ float sigclipmultip; /* Multiple of sigma in sigma clipping. */
+ float sigclipparam; /* Tolerance to stop or number of clips. */
size_t numbins; /* Number of bins in histogram or CFP. */
- uint8_t lowerbin; /* Save interval lower limit, not center. */
+ size_t numasciibins; /* Number of bins in ASCII plots. */
+ size_t asciiheight; /* Height of ASCII histogram or CFP plots. */
uint8_t normalize; /* set the sum of all bins to 1. */
float onebinstart; /* Shift bins to start at this value. */
uint8_t maxbinone; /* Set the maximum bin to 1. */
/* Internal */
uint8_t needssort; /* If sorting is needed. */
- uint8_t basicinfo; /* Print the basic information. */
gal_data_t *input; /* Input data structure. */
+ gal_data_t *sorted; /* Sorted input data structure. */
int isfits; /* Input is a FITS file. */
int hdu_type; /* Type of HDU (image or table). */
time_t rawtime; /* Starting time of the program. */
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index c08b0a5..cc5abc2 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -30,67 +30,32 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <string.h>
#include <stdlib.h>
+#include <gnuastro/fits.h>
#include <gnuastro/arithmetic.h>
#include <gnuastro/statistics.h>
#include <timing.h>
+#include <checkset.h>
#include "main.h"
+
+#include "ui.h"
#include "statistics.h"
-/* This function will report the simple immediate statistics of the
- data. For the average and standard deviation, the unsorted data is
- used so we don't suddenly encounter rounding errors. */
-void
-reportsimplestats(struct statisticsparams *p)
+
+
+
+/*******************************************************************/
+/************** Print in one row ***************/
+/*******************************************************************/
+static void
+ui_print_one_number(gal_data_t *out)
{
- printf("\n... in reportsimplestats ...\n");
- exit(1);
-
-#if 0
-
- double sum;
- size_t modeindex;
- float modequant, symvalue;
- float ave, std, med, modesym;
-
- sum=gal_statistics_float_sum(p->img, p->size);
- gal_statistics_f_ave_std(p->img, p->size, &ave, &std, NULL);
- med=p->sorted[gal_statistics_index_from_quantile(p->size, 0.5f)];
-
- /* Very simple and basic: */
- printf(SNAMEVAL FNAMEVAL FNAMEVAL, "Number of points", p->size,
- "Minimum", p->sorted[0], "Maximum", p->sorted[p->size-1]);
- printf(FNAMEVAL FNAMEVAL FNAMEVAL FNAMEVAL, "Sum", sum, "Mean", ave,
- "Standard deviation", std, "Median", med);
-
- /* The mode: */
- gal_statistics_mode_index_in_sorted(p->sorted, p->size, p->mirrordist,
- &modeindex, &modesym);
- modequant=(float)(modeindex)/(float)(p->size);
-
- /* Report the values: */
- printf(" -- %-45s%.4f %g\n", "Mode (quantile, value)",
- modequant, p->sorted[modeindex]);
- symvalue=gal_statistics_mode_value_from_sym(p->sorted, p->size,
- modeindex, modesym);
- printf(" -- %-45s%.4f %g\n", "Mode symmetricity and its cutoff"
- " value", modesym, symvalue);
- if(modesym<GAL_STATISTICS_MODE_SYM_GOOD)
- printf(" ## MODE SYMMETRICITY IS TOO LOW ##\n");
-
- /* Save the mode histogram and cumulative frequency plot. Note
- that if the histograms are to be built, then
- mhistname!=NULL. */
- if(p->mhistname)
- gal_statistics_mode_mirror_plots(p->sorted, p->size, modeindex,
- p->histmin, p->histmax,
- p->histnumbins, p->mhistname,
- p->mcfpname, (p->histrangeformirror
- ? 0.0f
- : p->mirrorplotdist) );
-#endif
+ char *toprint=gal_data_write_to_string(out->array, out->type, 0);
+ printf("%s ", toprint);
+ gal_data_free(out);
+ free(toprint);
}
@@ -98,37 +63,103 @@ reportsimplestats(struct statisticsparams *p)
static void
-print_ascii_plot(gal_data_t *plot, gal_data_t *bins, int h1_c0)
+ui_print_one_row(struct statisticsparams *p)
+{
+ struct gal_linkedlist_ill *tmp;
+
+ /* Print the numbers. */
+ for(tmp=p->toprint; tmp!=NULL; tmp=tmp->next)
+ switch(tmp->v)
+ {
+ case ARGS_OPTION_KEY_NUMBER:
+ ui_print_one_number( gal_statistics_number(p->input) ); break;
+ case ARGS_OPTION_KEY_MINIMUM:
+ ui_print_one_number( gal_statistics_minimum(p->input) ); break;
+ case ARGS_OPTION_KEY_MAXIMUM:
+ ui_print_one_number( gal_statistics_maximum(p->input) ); break;
+ case ARGS_OPTION_KEY_SUM:
+ ui_print_one_number( gal_statistics_sum(p->input) ); break;
+ case ARGS_OPTION_KEY_MEAN:
+ ui_print_one_number( gal_statistics_mean(p->input) ); break;
+ case ARGS_OPTION_KEY_STD:
+ ui_print_one_number( gal_statistics_std(p->input) ); break;
+ case ARGS_OPTION_KEY_MEDIAN:
+ ui_print_one_number( gal_statistics_median(p->sorted, 0) ); break;
+ case ARGS_OPTION_KEY_MODE:
+ error(EXIT_FAILURE, 0, "mode isn't implemented yet!");
+ break;
+ default:
+ error(EXIT_FAILURE, 0, "A bug! Operation code %d not recognized in "
+ "`ui_print_row'. Please contact us at %s so we can address "
+ "the problem", tmp->v, PACKAGE_BUGREPORT);
+ }
+
+ /* Print a new line. */
+ printf("\n");
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*******************************************************************/
+/************** ASCII plots ***************/
+/*******************************************************************/
+static void
+print_ascii_plot(struct statisticsparams *p, gal_data_t *plot,
+ gal_data_t *bins, int h1_c0, int printinfo)
{
int i, j;
- size_t height=10;
- float *p, *b, *f, *ff, halfbinwidth;
+ size_t *s, *sf, max=0;
+ double *b, v, halfbinwidth, correction;
- /* Print the range so the user knows. */
- b=bins->array;
- halfbinwidth = (b[1]-b[0])/2;
- printf("\n%s:\nX: (linear: %g -- %g)\nY: (linear)\n\n",
- ( h1_c0 ? "Histogram" : "Cumulative frequency plot"),
- b[0]-halfbinwidth, b[ bins->size - 1 ] + halfbinwidth);
+ /* Find the maximum of the plot. */
+ sf=(s=plot->array)+plot->size; do max = *s>max ? *s : max; while(++s<sf);
- /* The maximum values are set to one. So, multiply all the pixels by the
- desired height in character-height-units. */
- ff=(f=plot->array)+plot->size; do *f++ *= height; while(f<ff);
+ /* Print the range so the user knows. */
+ if(printinfo)
+ {
+ b=bins->array;
+ halfbinwidth = (b[1]-b[0])/2;
+ printf("\nASCII %s:\n", ( h1_c0 ? "Histogram" :
+ "Cumulative frequency plot") );
+ if(h1_c0) printf("Number: %zu\n", p->input->size);
+ printf("Y: (linear: 0 to %zu)\n", max);
+ printf("X: (linear: %g -- %g, in %zu bins)\n", b[0]-halfbinwidth,
+ b[ bins->size - 1 ] + halfbinwidth, bins->size);
+ }
/* Print the ASCII plot: */
- p=plot->array;
- for(i=height;i>=0;--i)
+ s=plot->array;
+ correction = (double)(p->asciiheight) / (double)max;
+ for(i=p->asciiheight;i>=0;--i)
{
- printf(" |");
- for(j=0;j<bins->size;++j)
+ printf(" |");
+ for(j=0;j<plot->size;++j)
{
- if(p[j]>=((float)i-0.5f) && p[j]>0.0f) printf("*");
+ v = (double)s[j] * correction;
+ if( v >= ((double)i-0.5f) && v > 0.0f ) printf("*");
else printf(" ");
}
printf("\n");
}
- printf(" |");
- for(j=0;j<bins->size;++j) printf("-");
+ printf(" |");
+ for(j=0;j<plot->size;++j) printf("-");
printf("\n\n");
}
@@ -138,21 +169,20 @@ print_ascii_plot(gal_data_t *plot, gal_data_t *bins, int
h1_c0)
static void
ascii_plots(struct statisticsparams *p)
{
- size_t numbins=70;
gal_data_t *bins, *hist, *cfp;
/* Make the bins and the respective plot. */
- bins=gal_statistics_regular_bins(p->input, NULL, numbins, NAN);
- hist=gal_statistics_histogram(p->input, bins, 0, 1);
+ bins=gal_statistics_regular_bins(p->input, NULL, p->numasciibins, NAN);
+ hist=gal_statistics_histogram(p->input, bins, 0, 0);
if(p->asciicfp)
{
bins->next=hist;
- cfp=gal_statistics_cfp(p->input, bins, 1);
+ cfp=gal_statistics_cfp(p->input, bins, 0);
}
/* Print the plots. */
- if(p->asciihist) print_ascii_plot(hist, bins, 1);
- if(p->asciicfp) print_ascii_plot(cfp, bins, 0);
+ if(p->asciihist) print_ascii_plot(p, hist, bins, 1, 1);
+ if(p->asciicfp) print_ascii_plot(p, cfp, bins, 0, 1);
/* Clean up.*/
gal_data_free(bins);
@@ -163,193 +193,379 @@ ascii_plots(struct statisticsparams *p)
-void
-printhistcfp(struct statisticsparams *p, float *bins, size_t numbins,
- char *filename, char *outputtype)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*******************************************************************/
+/******* Histogram and cumulative frequency tables ***********/
+/*******************************************************************/
+static void
+save_hist_and_or_cfp(struct statisticsparams *p)
{
- printf("\n... in printhistcfp ...\n");
- exit(1);
-
-#if 0
- float d;
- size_t i;
- FILE *out;
- time_t rawtime;
- int int0float1=1;
-
- /* Open the file: */
- errno=0;
- out=fopen(filename, "w");
- if(out==NULL)
- error(EXIT_FAILURE, errno, "couldn't open file %s", filename);
-
- /* Get the time to print on the report. */
- time(&rawtime);
- fprintf(out, "# %s \n# %s, created on %s", SPACK_STRING, outputtype,
- ctime(&rawtime));
- fprintf(out, "# Input (hdu): %s (%s)\n", p->up.inputname, p->cp.hdu);
- if(p->up.masknameset)
- fprintf(out, "# Mask (hdu): %s (%s)\n", p->up.maskname, p->up.mhdu);
-
- if(p->lowerbin)
- fprintf(out, "# Column 1: Flux of lower value of each bin\n");
+ gal_data_t *bins, *hist, *cfp=NULL;
+ struct gal_linkedlist_stll *comments=NULL;
+ char *tmp, *contents, *output, *suf, *fix, *suffix;
+
+
+ /* Set the bins and make the histogram, this is necessary for both the
+ histogram and CFP (recall that the CFP is built from the
+ histogram). */
+ bins=gal_statistics_regular_bins(p->input, NULL, p->numbins,
+ p->onebinstart);
+ hist=gal_statistics_histogram(p->input, bins, p->normalize, p->maxbinone);
+
+
+ /* Set the histogram as the next pointer of bins. This is again necessary
+ in both cases: when only a histogram is requested, it is used for the
+ plotting. When only a CFP is desired, it is used as input into
+ `gal_statistics_cfp'. */
+ bins->next=hist;
+
+
+ /* Make the cumulative frequency plot if the user wanted it. Make the
+ CFP, note that for the CFP, `maxbinone' and `normalize' are the same:
+ the last bin (largest value) must be one. So if any of them are given,
+ then set the last argument to 1.*/
+ if(p->cumulative)
+ cfp=gal_statistics_cfp(p->input, bins, p->normalize || p->maxbinone);
+
+
+ /* FITS tables don't accept `uint64_t', so to be consistent, we'll conver
+ the histogram and CFP to `uint32_t'.*/
+ if(hist->type==GAL_DATA_TYPE_UINT64)
+ hist=gal_data_copy_to_new_type_free(hist, GAL_DATA_TYPE_UINT32);
+ if(cfp && cfp->type==GAL_DATA_TYPE_UINT64)
+ cfp=gal_data_copy_to_new_type_free(cfp, GAL_DATA_TYPE_UINT32);
+
+
+ /* Finalize the next pointers. */
+ bins->next=hist;
+ hist->next=cfp;
+
+
+ /* Prepare the contents. */
+ if(p->histogram && p->cumulative)
+ { suf="_hist_cfp"; contents="Histogram and cumulative frequency plot"; }
+ else if(p->histogram)
+ { suf="_hist"; contents="Histogram"; }
else
- fprintf(out, "# Column 1: Flux in the middle of each bin\n");
+ { suf="_cfp"; contents="Cumulative frequency plot"; }
- if(strcmp(outputtype, CFPSTRING)==0)
- {
- fprintf(out, "# Column 2: Average of the sorted index of all points "
- "in this bin");
- if(p->normcfp)
- fprintf(out, " (normalized).\n");
- else if (p->maxcfpeqmaxhist)
- fprintf(out, " (Scaled to the histogram).\n");
- else
- {
- fprintf(out, ".\n");
- int0float1=0;
- }
- }
- else if (strcmp(outputtype, HISTSTRING)==0)
+
+ /* Set the basename and output suffix. */
+ fix = ( p->cp.output
+ ? gal_fits_name_is_fits(p->cp.output) ? "fits" : "txt"
+ : "txt" );
+ asprintf(&suffix, "%s.%s", suf, fix);
+
+
+ /* If another file is to be created or output name isn't set, we'll need
+ to use Automatic output. */
+ output = ( p->cp.output
+ ? p->cp.output
+ : gal_checkset_automatic_output(&p->cp, p->inputname, suffix) );
+
+
+ /* Write the comments, NOTE: we are writing the first two in reverse of
+ the order we want them. They will later be freed as part of the list's
+ freeing.*/
+ tmp=gal_fits_name_save_as_string(p->inputname, p->cp.hdu);
+ gal_linkedlist_add_to_stll(&comments, tmp, 0);
+
+ asprintf(&tmp, "%s created from:", contents);
+ gal_linkedlist_add_to_stll(&comments, tmp, 0);
+
+ if(strcmp(fix, "fits")) /* The intro info will be in FITS files anyway.*/
+ gal_table_comments_add_intro(&comments, PROGRAM_STRING, &p->rawtime);
+
+
+ /* Write the table. */
+ gal_table_write(bins, comments, p->cp.tableformat, output,
+ p->cp.dontdelete);
+
+
+ /* Let the user know, if we aren't in quiet mode. */
+ if(!p->cp.quiet)
+ printf("%s created.\n", output);
+
+
+ /* Clean up. */
+ free(suffix);
+ if(output!=p->cp.output) free(output);
+ gal_linkedlist_free_stll(comments, 1);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*******************************************************************/
+/************** Basic information ***************/
+/*******************************************************************/
+/* To keep things in `print_basics' clean, we'll define the input data
+ here, then only print the values there. */
+void
+print_input_info(struct statisticsparams *p)
+{
+ char *str, *name, *col=NULL;
+
+ /* Print the program name and version. */
+ printf("%s\n", PROGRAM_STRING);
+
+ /* Print the input information, if the input was a table, we also need to
+ give the column information. When the column has a name, it will be
+ printed, when it doesn't, we'll use the same string the user gave. */
+ printf("-------\n");
+ name=gal_fits_name_save_as_string(p->inputname, p->cp.hdu);
+ printf("Input: %s\n", name);
+
+ /* If a table was given, print the column. */
+ if(p->column) printf("Column: %s\n", p->column);
+
+ /* Range. */
+ str=NULL;
+ if( !isnan(p->greaterequal) && !isnan(p->lessthan) )
+ asprintf(&str, "from (inclusive) %g, upto (exclusive) %g",
+ p->greaterequal, p->lessthan);
+ else if( !isnan(p->greaterequal) )
+ asprintf(&str, "from (inclusive) %g", p->greaterequal);
+ else if( !isnan(p->lessthan) )
+ asprintf(&str, "upto (exclusive) %g", p->lessthan);
+ if(str)
{
- if(p->normhist)
- fprintf(out, "# Column 2: Fraction of points in this bin. \n");
- else if(p->maxhistone)
- fprintf(out, "# Column 2: Histogram if the maximum bin is "
- "set to 1.\n");
- else
- {
- fprintf(out, "# Column 2: Number of points in this bin. \n");
- int0float1=0;
- }
+ printf("Range: %s.\n", str);
+ free(str);
}
- /* Put the data in the file: */
- if(p->lowerbin==0) d=(bins[2]-bins[0])/2;
- else d=0.0f;
- if(int0float1)
- for(i=0;i<numbins;++i)
- fprintf(out, "%-20.6f"PRINTFLT, bins[i*2]+d, bins[i*2+1]);
- else
- for(i=0;i<numbins;++i)
- fprintf(out, "%-20.6f"PRINTINT, bins[i*2]+d, bins[i*2+1]);
- fclose(out);
-#endif
+ /* Units. */
+ if(p->input->unit) printf("Unit: %s\n", p->input->unit);
+
+ /* Clean up. */
+ if(col) free(col);
+ free(name);
+ printf("-------\n");
}
+/* This function will report the simple immediate statistics of the
+ data. For the average and standard deviation, the unsorted data is
+ used so we don't suddenly encounter rounding errors. */
void
-statistics(struct statisticsparams *p)
+print_basics(struct statisticsparams *p)
{
+ char *str;
+ double mean, std;
+ int namewidth=40;
+ gal_data_t *tmp, *bins, *hist;
+
+ /* Define the input dataset. */
+ print_input_info(p);
+
+ /* Print the number: */
+ printf(" %-*s %zu\n", namewidth, "Number of elements:", p->input->size);
+
+ /* Minimum: */
+ tmp=gal_statistics_minimum(p->input);
+ str=gal_data_write_to_string(tmp->array, tmp->type, 0);
+ printf(" %-*s %s\n", namewidth, "Minimum:", str);
+ gal_data_free(tmp);
+ free(str);
+
+ /* Maximum: */
+ tmp=gal_statistics_maximum(p->input);
+ str=gal_data_write_to_string(tmp->array, tmp->type, 0);
+ printf(" %-*s %s\n", namewidth, "Maximum:", str);
+ gal_data_free(tmp);
+ free(str);
+
+ /* Find the mean and standard deviation, but don't print them, see
+ explanations under median. */
+ tmp=gal_statistics_mean_std(p->input);
+ mean = ((double *)(tmp->array))[0];
+ std = ((double *)(tmp->array))[1];
+ gal_data_free(tmp);
+
+ /* Find and print the median: we want the median to be found in place to
+ save time/memory. But having a sorted array can decrease the floating
+ point accuracy of the standard deviation. So we'll do the median
+ calculation in the end.*/
+ tmp=gal_statistics_median(p->input, 1);
+ str=gal_data_write_to_string(tmp->array, tmp->type, 0);
+ printf(" %-*s %s\n", namewidth, "Median:", str);
+ gal_data_free(tmp);
+ free(str);
+
+ /* Print the mean and standard deviation. */
+ printf(" %-*s %g\n", namewidth, "Mean:", mean);
+ printf(" %-*s %g\n", namewidth, "Standard deviation:", std);
+
+ /* Ascii histogram. Note that we don't want to force the user to have the
+ plotting parameters. */
+ printf("-------\nHistogram:\n");
+ p->asciiheight = p->asciiheight ? p->asciiheight : 10;
+ p->numasciibins = p->numasciibins ? p->numasciibins : 70;
+ bins=gal_statistics_regular_bins(p->input, NULL, p->numasciibins, NAN);
+ hist=gal_statistics_histogram(p->input, bins, 0, 0);
+ print_ascii_plot(p, hist, bins, 1, 0);
+ gal_data_free(bins);
+ gal_data_free(hist);
+}
+
+
+
+
+
+
+
+
+
+
- /* If none of the processes here are requested, return. */
- if( p->asciihist==0
- && p->asciicfp==0
- && p->histogram==0
- && p->cumulative==0
- && p->sigclipstr==NULL
- && isnan(p->mirrorquant) )
- return;
- /* For the main body of the program, we will assume the data is in
- float32. */
- p->input=gal_data_copy_to_new_type_free(p->input, GAL_DATA_TYPE_FLOAT32);
- /* Print the ASCII histogram if requested. */
- if(p->asciihist || p->asciicfp) ascii_plots(p);
-#if 0
- int r;
- size_t i;
- float quant=-1.0f; /* The quantile was already */
- float ave, std, med;
- float maxhist=-FLT_MAX, *bins=NULL; /* taken into affect in ui.c. */
- /* Report the simple statistics: */
- if(p->cp.verb)
+
+
+
+/*******************************************************************/
+/************** Sigma clipping ***************/
+/*******************************************************************/
+void
+print_sigma_clip(struct statisticsparams *p)
+{
+ float *a;
+ char *mode;
+ int namewidth=40;
+ gal_data_t *sigclip;
+
+ /* Set the mode for printing: */
+ if( p->sigclipparam>=1.0f )
+ asprintf(&mode, "for %g clips", p->sigclipparam);
+ else
+ asprintf(&mode, "until relative change in STD is less than %g",
+ p->sigclipparam);
+
+ /* Report the status */
+ if(!p->cp.quiet)
{
- reportsimplestats(p);
- if(p->asciihist) printasciihist(p);
+ print_input_info(p);
+ printf("%g-sigma clipping steps %s:\n\n", p->sigclipmultip, mode);
}
- /* Make the histogram: */
- if(p->histname)
+ /* Do the Sigma clipping: */
+ sigclip=gal_statistics_sigma_clip(p->sorted, p->sigclipmultip,
+ p->sigclipparam, p->cp.quiet);
+ a=sigclip->array;
+
+ /* Finish the introduction. */
+ if(!p->cp.quiet)
+ printf("-------\nSummary:\n");
+ else
+ printf("%g-sigma clipped %s:\n", p->sigclipmultip, mode);
+
+ /* Print the final results: */
+ printf(" %-*s %zu\n", namewidth, "Number of input elements:",
+ p->input->size);
+ if( p->sigclipparam < 1.0f )
+ printf(" %-*s %d\n", namewidth, "Number of clips:", sigclip->status);
+ printf(" %-*s %g\n", namewidth, "Final number of elements:", a[0]);
+ printf(" %-*s %g\n", namewidth, "Median:", a[1]);
+ printf(" %-*s %g\n", namewidth, "Mean:", a[2]);
+ printf(" %-*s %g\n", namewidth, "Standard deviation:", a[3]);
+
+ /* Clean up. */
+ free(mode);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*******************************************************************/
+/************** Main function ***************/
+/*******************************************************************/
+void
+statistics(struct statisticsparams *p)
+{
+ int print_basic_info=1;
+
+ /* Print the one-row numbers if the user asked for them. */
+ if(p->toprint)
{
- /* Make the actual data histogram and save it. */
- gal_statistics_set_bins(p->sorted, p->size, p->histnumbins, p->histmin,
- p->histmax, p->onebinvalue, quant, &bins);
- gal_statistics_histogram(p->sorted, p->size, bins, p->histnumbins,
- p->normhist, p->maxhistone);
- printhistcfp(p, bins, p->histnumbins, p->histname, HISTSTRING);
-
- /* Get the hisogram maximum value if it is needed for the
- cumulative frequency plot. */
- if(p->maxcfpeqmaxhist)
- for(i=0;i<p->histnumbins;++i)
- if(bins[i*2+1]>maxhist)
- maxhist=bins[i*2+1];
+ print_basic_info=0;
+ ui_print_one_row(p);
}
- /* Make the cumulative distribution function: */
- if(p->cfpname)
+ /* Print the ASCII plots if requested. */
+ if(p->asciihist || p->asciicfp)
{
- if(p->cfpsimhist)
- {
- p->cfpnum=p->histnumbins;
- for(i=0;i<p->cfpnum;++i)
- bins[i*2+1]=0.0f;
- }
- else
- {
- if(p->histname) free(bins);
- gal_statistics_set_bins(p->sorted, p->size, p->cfpnum, p->cfpmin,
- p->cfpmax, p->onebinvalue, quant, &bins);
- }
- gal_statistics_cumulative_fp(p->sorted, p->size, bins,
- p->cfpnum, p->normcfp);
-
- if(p->maxcfpeqmaxhist)
- for(i=0;i<p->cfpnum;++i)
- bins[i*2+1]*=(maxhist/p->size);
+ ascii_plots(p);
+ print_basic_info=0;
+ }
- printhistcfp(p, bins, p->cfpnum, p->cfpname, CFPSTRING);
+ /* Save the histogram and CFP as tables if requested. */
+ if(p->histogram || p->cumulative)
+ {
+ print_basic_info=0;
+ save_hist_and_or_cfp(p);
}
- /* Make the mirror distribution if asked for: */
- if(isnan(p->mirror)==0)
- gal_statistics_mode_mirror_plots(p->sorted, p->size,
- gal_statistics_index_from_quantile(p->size,
- p->mirror),
- p->histmin, p->histmax, p->histnumbins,
- p->mirrorhist, p->mirrorcfp,
- (p->histrangeformirror ? 0.0f
- : p->mirrorplotdist));
-
- /* Print out the Sigma clippings: */
- if(p->sigclip && p->cp.verb)
+ /* Print the sigma-clipped results. */
+ if( !isnan(p->sigclipmultip ) )
{
- printf(" - Sigma clipping results (Median, Mean, STD, Number):\n");
- printf(" - %.2f times sigma by convergence (tolerance: %.4f):\n",
- p->sigclipmultip, p->sigcliptolerance);
- r=gal_statistics_sigma_clip_converge(p->sorted, 1, p->size,
- p->sigclipmultip,
- p->sigcliptolerance, &ave,
- &med, &std, 1);
- if(r==0)
- printf(" #### Could not converge\n");
- printf(" - %.2f sigma-clipping %zu times:\n",
- p->sigclipmultip, p->sigclipnum);
- gal_statistics_sigma_clip_certain_num(p->sorted, 1, p->size,
- p->sigclipmultip, p->sigclipnum,
- &ave, &med, &std, 1);
+ print_basic_info=0;
+ print_sigma_clip(p);
}
- /* Free the allocated arrays: */
- if(p->histname || p->cfpname) free(bins);
-#endif
+ /* If nothing was requested print the simple statistics. */
+ if(print_basic_info)
+ print_basics(p);
+
}
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index 9fe149b..1a9b7a7 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -65,9 +65,14 @@ static char
args_doc[] = "ASTRdata";
const char
-doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" will print the basic "
- "statistics of the input image pixel flux distribution. All blank pixels "
- "or pixels specified by a mask image will be ignored.\n"
+doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" will do statistical "
+ "analysis on the input dataset (table column or image). All blank "
+ "pixels or pixels outside of the given range are ignored. You can "
+ "either directly ask for certain statistics in one line/row as shown "
+ "below with the same order as requested, or get tables of different "
+ "statistical measures like the histogram, cumulative frequency style "
+ "and etc. If no particular statistic is requested, some basic "
+ "information about the dataset is printed on the command-line.\n"
GAL_STRINGS_MORE_HELP_INFO
/* After the list of options: */
"\v"
@@ -88,48 +93,6 @@ enum program_args_groups
-/* Available letters for short options:
-
- a b c d e f i j k p r u v w x y z
- B E F G J L R W X Y Z
-*/
-enum option_keys_enum
-{
- /* With short-option version. */
- ARGS_OPTION_KEY_COLUMN = 'c',
- ARGS_OPTION_KEY_GREATEREQUAL = 'g',
- ARGS_OPTION_KEY_LESSTHAN = 'l',
- ARGS_OPTION_KEY_QUANTRANGE = 'Q',
- ARGS_OPTION_KEY_MEAN = 'm',
- ARGS_OPTION_KEY_STD = 't',
- ARGS_OPTION_KEY_MEDIAN = 'M',
- ARGS_OPTION_KEY_MODE = 'O',
- ARGS_OPTION_KEY_ASCIIHIST = 'A',
- ARGS_OPTION_KEY_HISTOGRAM = 'H',
- ARGS_OPTION_KEY_CUMULATIVE = 'C',
- ARGS_OPTION_KEY_SIGMACLIP = 's',
- ARGS_OPTION_KEY_NORMALIZE = 'n',
-
- /* Only with long version (start with a value 1000, the rest will be set
- automatically). */
- ARGS_OPTION_KEY_NUMBER = 1000,
- ARGS_OPTION_KEY_MINIMUM,
- ARGS_OPTION_KEY_MAXIMUM,
- ARGS_OPTION_KEY_SUM,
- ARGS_OPTION_KEY_ASCIICFP,
- ARGS_OPTION_KEY_MIRROR,
- ARGS_OPTION_KEY_NUMBINS,
- ARGS_OPTION_KEY_LOWERBIN,
- ARGS_OPTION_KEY_ONEBINSTART,
- ARGS_OPTION_KEY_MAXBINONE,
-};
-
-
-
-
-
-
-
@@ -165,9 +128,11 @@ ui_initialize_options(struct statisticsparams *p,
/* Program-specific initializers */
p->lessthan = NAN;
p->onebinstart = NAN;
- p->mirrorquant = NAN;
p->greaterequal = NAN;
- p->quantilerange = NAN;
+ p->quantmin = NAN;
+ p->quantmax = NAN;
+ p->sigclipparam = NAN;
+ p->sigclipmultip = NAN;
/* Set the mandatory common options. */
for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
@@ -242,6 +207,14 @@ ui_add_to_print_in_row(struct argp_option *option, char
*arg,
{
struct statisticsparams *p=(struct statisticsparams *)params;
+ if(lineno==-1)
+ error(EXIT_FAILURE, 0, "currently the options to be printed in one row "
+ "(like `--number', `--mean', and etc) do not support printing "
+ "with the `--printparams' (`-P'), or writing into configuration "
+ "files due to lack of time when implementing these features. "
+ "Please get in touch with us at `%s', so we can implement it if "
+ "it is possible now, thank you", PACKAGE_BUGREPORT);
+
/* If this option is given in a configuration file, then `arg' will not
be NULL and we don't want to do anything if it is `0'. */
if( arg && arg[1]!='\0' && *arg!='0' && *arg!='1' )
@@ -263,6 +236,135 @@ ui_add_to_print_in_row(struct argp_option *option, char
*arg,
+static void *
+ui_parse_numbers(struct argp_option *option, char *arg,
+ char *filename, size_t lineno, void *params)
+{
+ char *str;
+ gal_data_t *in;
+ struct statisticsparams *p=(struct statisticsparams *)params;
+
+ /* For the `--printparams' (`-P') option:*/
+ if(lineno==-1)
+ {
+ switch(option->key)
+ {
+ case ARGS_OPTION_KEY_SIGMACLIP:
+ asprintf(&str, "%g,%g", p->sigclipmultip, p->sigclipparam);
+ break;
+ case ARGS_OPTION_KEY_QRANGE:
+ if( isnan(p->quantmax) ) asprintf(&str, "%g", p->quantmin);
+ else asprintf(&str, "%g,%g", p->quantmin, p->quantmax);
+ break;
+ default:
+ error(EXIT_FAILURE, 0, "a bug! option `%s' not recognized "
+ "in `ui_parse_numbers' (when called for printing). Please "
+ "contact us at %s to fix the problem", option->name,
+ PACKAGE_BUGREPORT);
+ }
+ return str;
+ }
+
+ /* Parse the inputs. */
+ in=gal_options_parse_list_of_numbers(arg, filename, lineno);
+
+ /* Checks depending on the option. */
+ switch(option->key)
+ {
+ case ARGS_OPTION_KEY_SIGMACLIP:
+ /* Check if there was only two numbers. */
+ if(in->size!=2)
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "the `--%s' "
+ "option takes two values (separated by a comma) for "
+ "defining the sigma-clip. However, %zu numbers were "
+ "read in the string `%s' (value to this option).\n\n"
+ "The first number is the multiple of sigma, and the "
+ "second is either the tolerance (if its is less than "
+ "1.0), or a specific number of times to clip (if it "
+ "is equal or larger than 1.0).", option->name, in->size,
+ arg);
+
+ /* Read the values in. */
+ p->sigclipparam = ((double *)(in->array))[1];
+ p->sigclipmultip = ((double *)(in->array))[0];
+
+ /* Some sanity checks. */
+ if( p->sigclipmultip<=0 )
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "the first value to "
+ "the `--%s' option (multiple of sigma), must be "
+ "greater than zero. From the string `%s' (value to "
+ "this option), you have given a value of %g for the "
+ "first value", option->name, arg, p->sigclipmultip);
+
+ /* If the second value must also be positive. */
+ if( p->sigclipparam<=0 )
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "the second value "
+ "to the `--%s' option (tolerance to stop clipping or "
+ "number of clips), must be greater than zero. From "
+ "the string `%s' (value to this option), you have "
+ "given a value of %g for the second value",
+ option->name, arg, p->sigclipparam);
+
+ /* if the second value is larger or equal to 1.0, it must be an
+ integer. */
+ if( p->sigclipparam >= 1.0f && ceil(p->sigclipparam) != p->sigclipparam)
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "when the second "
+ "value to the `--%s' option is >=1, it is interpretted "
+ "as an absolute number of clips. So it must be an "
+ "integer. However, your second value is a floating "
+ "point number: %g (parsed from `%s')", option->name,
+ p->sigclipparam, arg);
+ break;
+
+ case ARGS_OPTION_KEY_QRANGE:
+ /* Check if there was only two numbers. */
+ if(in->size!=1 && in->size!=2)
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "the `--%s' "
+ "option takes one or two values values (separated by "
+ "a comma) to define the range of used values with "
+ "quantiles. However, %zu numbers were read in the "
+ "string `%s' (value to this option).\n\n"
+ "If there is only one number as input, it will be "
+ "interpretted as the lower quantile (Q) range. The "
+ "higher range will be set to the quantile (1-Q). "
+ "When two numbers are given, they will be used as the "
+ "lower and higher quantile range respectively",
+ option->name, in->size, arg);
+
+ /* Read the values in. */
+ p->quantmin = ((double *)(in->array))[0];
+ if(in->size==2) p->quantmax = ((double *)(in->array))[1];
+
+ /* Make sure the values are between 0 and 1. */
+ if( (p->quantmin<0 || p->quantmin>1)
+ || ( !isnan(p->quantmax) && (p->quantmax<0 || p->quantmax>1) ) )
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "values to the "
+ "`--quantrange' option must be between 0 and 1 "
+ "(inclusive). Your input was: `%s'", arg);
+ break;
+
+ /* When only one value is given, make sure it is less than 0.5. */
+ if( !isnan(p->quantmax) && p->quantmin>0.5 )
+ error(EXIT_FAILURE, 0, "%g>=0.5! When only one value is given to the "
+ "`--%s' option, the range is defined as Q and 1-Q. Thus, the "
+ "value must be less than 0.5", p->quantmin, option->name);
+
+ default:
+ error(EXIT_FAILURE, 0, "a bug! option `%s' not recognized "
+ "in `ui_parse_numbers' (when called for printing). Please "
+ "contact us at %s to fix the problem", option->name,
+ PACKAGE_BUGREPORT);
+ }
+
+
+ /* No return value necessary for this function. */
+ return NULL;
+}
+
+
+
+
+
@@ -303,23 +405,30 @@ ui_read_check_only_options(struct statisticsparams *p)
/* Less-than and greater-equal cannot be called together with
quantrange. */
if( ( !isnan(p->lessthan) || !isnan(p->greaterequal) )
- && !isnan(p->quantilerange) )
+ && !isnan(p->quantmin) )
error(EXIT_FAILURE, 0, "`--lessthan' and/or `--greaterequal' cannot "
"be called together with `--quantrange'");
- /* The quantile range only makes sense with value less than 0.5. */
- if( !isnan(p->quantilerange) && p->quantilerange>=0.5)
- error(EXIT_FAILURE, 0, "%g>=0.5! The quantile range must be less than "
- "0.5. Recall the the quantile (Q) range is defined to be: Q to "
- "1-Q", p->quantilerange);
-
-
/* When binned outputs are requested, make sure that `numbins' is set. */
- if(p->histogram || p->cumulative)
+ if( (p->histogram || p->cumulative) && p->numbins==0)
error(EXIT_FAILURE, 0, "`--numbins' isn't set. When the histogram or "
"cumulative frequency plots are requested, the number of bins "
"(`--numbins') is necessary");
+
+
+ /* If an ascii plot is requested, check if the ascii number of bins and
+ height are given. */
+ if( (p->asciihist || p->asciicfp)
+ && (p->numasciibins==0 || p->asciiheight==0) )
+ error(EXIT_FAILURE, 0, "when an ascii plot is requested, "
+ "`--numasciibins' and `--asciiheight' are mandatory, but atleast "
+ "one of these has not been given");
+
+
+ /* Reverse the list of statistics to print in one row, so it has the same
+ order the user wanted. */
+ gal_linkedlist_reverse_ill(&p->toprint);
}
@@ -401,62 +510,6 @@ ui_check_options_and_arguments(struct statisticsparams *p)
/*************** Preparations *******************/
/**************************************************************/
static void
-ui_print_one_number(gal_data_t *out)
-{
- char *toprint=gal_data_write_to_string(out->array, out->type, 0);
- printf("%s ", toprint);
- gal_data_free(out);
- free(toprint);
-}
-
-
-
-
-
-void
-ui_print_one_row(struct statisticsparams *p)
-{
- struct gal_linkedlist_ill *tmp;
-
- /* Reverse the list to have the same order the user wanted. */
- gal_linkedlist_reverse_ill(&p->toprint);
-
- /* Print the numbers. */
- for(tmp=p->toprint; tmp!=NULL; tmp=tmp->next)
- switch(tmp->v)
- {
- case ARGS_OPTION_KEY_NUMBER:
- ui_print_one_number( gal_statistics_number(p->input) ); break;
- case ARGS_OPTION_KEY_MINIMUM:
- ui_print_one_number( gal_statistics_minimum(p->input) ); break;
- case ARGS_OPTION_KEY_MAXIMUM:
- ui_print_one_number( gal_statistics_maximum(p->input) ); break;
- case ARGS_OPTION_KEY_SUM:
- ui_print_one_number( gal_statistics_sum(p->input) ); break;
- case ARGS_OPTION_KEY_MEAN:
- ui_print_one_number( gal_statistics_mean(p->input) ); break;
- case ARGS_OPTION_KEY_STD:
- ui_print_one_number( gal_statistics_std(p->input) ); break;
- case ARGS_OPTION_KEY_MEDIAN:
- ui_print_one_number( gal_statistics_median(p->input) ); break;
- case ARGS_OPTION_KEY_MODE:
- error(EXIT_FAILURE, 0, "mode isn't implemented yet!");
- break;
- default:
- error(EXIT_FAILURE, 0, "A bug! Operation code %d not recognized in "
- "`ui_print_row'. Please contact us at %s so we can address "
- "the problem", tmp->v, PACKAGE_BUGREPORT);
- }
-
- /* Print a new line. */
- printf("\n");
-}
-
-
-
-
-
-static void
ui_out_of_range_to_blank(struct statisticsparams *p)
{
size_t one=1;
@@ -466,6 +519,24 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
| GAL_ARITHMETIC_INPLACE
| GAL_ARITHMETIC_NUMOK );
+ /* If the user has given a quantile range, then set the `greaterequal'
+ and `lessthan' values. */
+ if( !isnan(p->quantmin) )
+ {
+ /* If only one value was given, set the maximum quantile range. */
+ if( isnan(p->quantmax) ) p->quantmax = 1 - p->quantmin;
+
+ /* Set the greater-equal value. */
+ tmp=gal_statistsics_quantile(p->input, p->quantmin, 1);
+ tmp=gal_data_copy_to_new_type_free(tmp, GAL_DATA_TYPE_FLOAT32);
+ p->greaterequal=*((float *)(tmp->array));
+
+ /* Set the lower-than value. */
+ tmp=gal_statistsics_quantile(p->input, p->quantmax, 1);
+ tmp=gal_data_copy_to_new_type_free(tmp, GAL_DATA_TYPE_FLOAT32);
+ p->lessthan=*((float *)(tmp->array));
+ }
+
/* Set the condition. Note that the `greaterequal' name is for the data
we want. So we will set the condition based on those that are
less-than */
@@ -506,7 +577,9 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
0, -1, NULL, NULL, NULL);
*((float *)(tmp->array)) = NAN;
- /* Set all the pixels that satisfy the condition to blank. */
+ /* Set all the pixels that satisfy the condition to blank. Note that a
+ blank value will be used in the proper type of the input in the
+ `where' operator.*/
gal_arithmetic(GAL_ARITHMETIC_OP_WHERE, flagsor, p->input, cond, tmp);
}
@@ -514,6 +587,46 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
+/* Check if a sorted array is necessary and if so, then make a sorted
+ array. */
+static void
+ui_make_sorted_if_necessary(struct statisticsparams *p)
+{
+ int is_necessary=0;
+ struct gal_linkedlist_ill *tmp;
+
+ /* Check in the one-row outputs. */
+ for(tmp=p->toprint; tmp!=NULL; tmp=tmp->next)
+ switch(tmp->v)
+ {
+ case ARGS_OPTION_KEY_MEDIAN:
+ case ARGS_OPTION_KEY_MODE:
+ is_necessary=1;
+ break;
+ }
+
+ /* Check in the rest of the outputs. */
+ if( is_necessary==0 && !isnan(p->sigclipmultip) )
+ is_necessary=1;
+
+
+ /* Do the sorting, we will keep the sorted array in a separate space,
+ since the unsorted nature of the original dataset will help decrease
+ floating point errors. If the input is already sorted, we'll just
+ point it to the input.*/
+ if( gal_statistics_is_sorted(p->input) )
+ p->sorted=p->input;
+ else
+ {
+ p->sorted=gal_data_copy(p->input);
+ gal_statistics_sort_increasing(p->sorted);
+ }
+}
+
+
+
+
+
void
ui_preparations(struct statisticsparams *p)
{
@@ -551,11 +664,14 @@ ui_preparations(struct statisticsparams *p)
/* Only keep the numbers we want. */
gal_blank_remove(p->input);
- /* Print the one-row numbers if the user asked for them. We want this to
- be done before converting the array to a float, since these operations
- can be done on any type and if the user has just asked for these
- operations, we don't want to waste their time for nothing. */
- if(p->toprint) ui_print_one_row(p);
+ /* Make sure there is any data remaining: */
+ if(p->input->size==0)
+ error(EXIT_FAILURE, 0, "%s: no data, maybe the `--greaterequal' or "
+ "`--lessthan' options need to be adjusted",
+ gal_fits_name_save_as_string(p->inputname, p->cp.hdu) );
+
+ /* Make the sorted array if necessary. */
+ ui_make_sorted_if_necessary(p);
}
@@ -660,4 +776,7 @@ ui_free_report(struct statisticsparams *p)
/* Free the allocated arrays: */
free(p->cp.hdu);
free(p->cp.output);
+ if(p->sorted!=p->input)
+ gal_data_free(p->sorted);
+ gal_data_free(p->input);
}
diff --git a/bin/statistics/ui.h b/bin/statistics/ui.h
index 88bd088..5162519 100644
--- a/bin/statistics/ui.h
+++ b/bin/statistics/ui.h
@@ -23,6 +23,52 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#ifndef UI_H
#define UI_H
+
+
+
+
+/* Available letters for short options:
+
+ a b c d e f i j k p r u v w x y z
+ B E F G J L R W X Y Z
+*/
+enum option_keys_enum
+{
+ /* With short-option version. */
+ ARGS_OPTION_KEY_COLUMN = 'c',
+ ARGS_OPTION_KEY_GREATEREQUAL = 'g',
+ ARGS_OPTION_KEY_LESSTHAN = 'l',
+ ARGS_OPTION_KEY_QRANGE = 'Q',
+ ARGS_OPTION_KEY_MEAN = 'm',
+ ARGS_OPTION_KEY_STD = 't',
+ ARGS_OPTION_KEY_MEDIAN = 'M',
+ ARGS_OPTION_KEY_MODE = 'O',
+ ARGS_OPTION_KEY_ASCIIHIST = 'A',
+ ARGS_OPTION_KEY_HISTOGRAM = 'H',
+ ARGS_OPTION_KEY_CUMULATIVE = 'C',
+ ARGS_OPTION_KEY_SIGMACLIP = 's',
+ ARGS_OPTION_KEY_NORMALIZE = 'n',
+
+ /* Only with long version (start with a value 1000, the rest will be set
+ automatically). */
+ ARGS_OPTION_KEY_NUMBER = 1000,
+ ARGS_OPTION_KEY_MINIMUM,
+ ARGS_OPTION_KEY_MAXIMUM,
+ ARGS_OPTION_KEY_SUM,
+ ARGS_OPTION_KEY_ASCIICFP,
+ ARGS_OPTION_KEY_NUMBINS,
+ ARGS_OPTION_KEY_NUMASCIIBINS,
+ ARGS_OPTION_KEY_ASCIIHEIGHT,
+ ARGS_OPTION_KEY_LOWERBIN,
+ ARGS_OPTION_KEY_ONEBINSTART,
+ ARGS_OPTION_KEY_MAXBINONE,
+};
+
+
+
+
+
+/* Functions */
void
ui_read_check_inputs_setup(int argc, char *argv[],
struct statisticsparams *p);
diff --git a/bin/table/ui.c b/bin/table/ui.c
index 9278a63..6bd2ce0 100644
--- a/bin/table/ui.c
+++ b/bin/table/ui.c
@@ -78,23 +78,6 @@ doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" can be used
to view the "
-/* Available letters for short options:
-
- a b d e f g j k l m n p r s t u v w x y z
- A B C E F G H J L M O Q R T U W X Y Z */
-enum option_keys_enum
-{
- /* With short-option version. */
- ARGS_OPTION_KEY_COLUMN = 'c',
- ARGS_OPTION_KEY_INFORMATION = 'i',
-
- /* Only with long version (start with a value 1000, the rest will be set
- automatically). */
-};
-
-
-
-
@@ -270,7 +253,7 @@ ui_check_options_and_arguments(struct tableparams *p)
void
ui_preparations(struct tableparams *p)
{
- char *numstr;
+ char *tmp;
int tableformat;
gal_data_t *allcols;
size_t i, numcols, numrows;
@@ -294,11 +277,9 @@ ui_preparations(struct tableparams *p)
{
/* Print the file information. */
printf("--------\n");
- printf("%s", p->filename);
- if(gal_fits_name_is_fits(p->filename))
- printf(" (hdu: %s)\n", cp->hdu);
- else
- printf("\n");
+ tmp=gal_fits_name_save_as_string(p->filename, p->cp.hdu);
+ printf("%s", tmp);
+ free(tmp);
/* Print each column's information. */
gal_table_print_info(allcols, numcols, numrows);
@@ -323,8 +304,8 @@ ui_preparations(struct tableparams *p)
else
for(i=1;i<=numcols;++i)
{
- asprintf(&numstr, "%zu", i);
- gal_linkedlist_add_to_stll(&p->columns, numstr, 0);
+ asprintf(&tmp, "%zu", i);
+ gal_linkedlist_add_to_stll(&p->columns, tmp, 0);
}
}
diff --git a/bin/table/ui.h b/bin/table/ui.h
index a0f5420..57e4a85 100644
--- a/bin/table/ui.h
+++ b/bin/table/ui.h
@@ -23,6 +23,28 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#ifndef UI_H
#define UI_H
+
+
+
+
+/* Available letters for short options:
+
+ a b d e f g j k l m n p r s t u v w x y z
+ A B C E F G H J L M O Q R T U W X Y Z */
+enum option_keys_enum
+{
+ /* With short-option version. */
+ ARGS_OPTION_KEY_COLUMN = 'c',
+ ARGS_OPTION_KEY_INFORMATION = 'i',
+
+ /* Only with long version (start with a value 1000, the rest will be set
+ automatically). */
+};
+
+
+
+
+
void
ui_read_check_inputs_setup(int argc, char *argv[], struct tableparams *p);
diff --git a/bin/warp/ui.c b/bin/warp/ui.c
index 771ff5c..4d46e3a 100644
--- a/bin/warp/ui.c
+++ b/bin/warp/ui.c
@@ -86,35 +86,6 @@ enum program_args_groups
-/* Available letters for short options:
-
- b g i j l n u v w x y
- A B E F G H I J L M O Q R T U W X Y Z */
-enum option_keys_enum
-{
- /* With short-option version. */
- ARGS_OPTION_KEY_KEEPWCS = 'k',
- ARGS_OPTION_KEY_COVEREDFRAC = 'C',
- ARGS_OPTION_KEY_TYPE = 'T',
- ARGS_OPTION_KEY_ALIGN = 'a',
- ARGS_OPTION_KEY_ROTATE = 'r',
- ARGS_OPTION_KEY_SCALE = 's',
- ARGS_OPTION_KEY_FLIP = 'f',
- ARGS_OPTION_KEY_SHEAR = 'e',
- ARGS_OPTION_KEY_TRANSLATE = 't',
- ARGS_OPTION_KEY_PROJECT = 'p',
- ARGS_OPTION_KEY_MATRIX = 'm',
- ARGS_OPTION_KEY_CENTERONCORNER = 'c',
-
- /* Only with long version (start with a value 1000, the rest will be set
- automatically). */
- ARGS_OPTION_KEY_HSTARTWCS = 1000,
- ARGS_OPTION_KEY_HENDWCS,
-};
-
-
-
-
diff --git a/bin/warp/ui.h b/bin/warp/ui.h
index 39e37c3..e49a73c 100644
--- a/bin/warp/ui.h
+++ b/bin/warp/ui.h
@@ -23,6 +23,40 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#ifndef UI_H
#define UI_H
+
+
+
+
+/* Available letters for short options:
+
+ b g i j l n u v w x y
+ A B E F G H I J L M O Q R T U W X Y Z */
+enum option_keys_enum
+{
+ /* With short-option version. */
+ ARGS_OPTION_KEY_KEEPWCS = 'k',
+ ARGS_OPTION_KEY_COVEREDFRAC = 'C',
+ ARGS_OPTION_KEY_TYPE = 'T',
+ ARGS_OPTION_KEY_ALIGN = 'a',
+ ARGS_OPTION_KEY_ROTATE = 'r',
+ ARGS_OPTION_KEY_SCALE = 's',
+ ARGS_OPTION_KEY_FLIP = 'f',
+ ARGS_OPTION_KEY_SHEAR = 'e',
+ ARGS_OPTION_KEY_TRANSLATE = 't',
+ ARGS_OPTION_KEY_PROJECT = 'p',
+ ARGS_OPTION_KEY_MATRIX = 'm',
+ ARGS_OPTION_KEY_CENTERONCORNER = 'c',
+
+ /* Only with long version (start with a value 1000, the rest will be set
+ automatically). */
+ ARGS_OPTION_KEY_HSTARTWCS = 1000,
+ ARGS_OPTION_KEY_HENDWCS,
+};
+
+
+
+
+
/* Functions */
void
ui_read_check_inputs_setup(int argc, char *argv[], struct warpparams *p);
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index f495497..f360814 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -192,8 +192,8 @@ sub-component to a title is present.
* Installation:: Requirements and installation.
* Common program behavior:: Common behavior in all programs.
* Extensions and Tables:: Tools to operate on extensions and tables.
-* Image manipulation:: Tools for basic image manipulation.
-* Image analysis:: Analyze images.
+* Data manipulation:: Tools for basic image manipulation.
+* Data analysis:: Analyze images.
* Modeling and fittings:: Make and fit models.
* High-level calculations:: Physical calculations.
* Libraries:: Use Gnuastro in your own code.
@@ -351,7 +351,7 @@ Table
* Invoking asttable:: Options and arguments to Table.
-Image manipulation
+Data manipulation
* Crop:: Crop region(s) from a dataset.
* Arithmetic:: Arithmetic on input data.
@@ -428,7 +428,7 @@ Tiling an image
* Checking grid values:: Check the values on the grid.
* Mesh grid options:: Common options related to the mesh grid.
-Image analysis
+Data analysis
* Statistics:: Calculate dataset statistics.
* NoiseChisel:: Detect objects in an image.
@@ -6073,7 +6073,7 @@ END
address@hidden Extensions and Tables, Image manipulation, Common program
behavior, Top
address@hidden Extensions and Tables, Data manipulation, Common program
behavior, Top
@chapter Extensions and Tables
@cindex File operations
@@ -7002,8 +7002,8 @@ it is possible to output one column more than once.
address@hidden Image manipulation, Image analysis, Extensions and Tables, Top
address@hidden Image manipulation
address@hidden Data manipulation, Data analysis, Extensions and Tables, Top
address@hidden Data manipulation
Images are one of the major formats of data that is used in astronomy. The
functions in this chapter explain the GNU Astronomy Utilities which are
@@ -7019,7 +7019,7 @@ transformation to it.
* SubtractSky:: Find the sky and subtract it from image.
@end menu
address@hidden Crop, Arithmetic, Image manipulation, Image manipulation
address@hidden Crop, Arithmetic, Data manipulation, Data manipulation
@section Crop
@cindex Section of an image
@@ -7717,7 +7717,7 @@ created (unless @option{--individual} is called).
address@hidden Arithmetic, Convolve, Crop, Image manipulation
address@hidden Arithmetic, Convolve, Crop, Data manipulation
@section Arithmetic
It is commonly necessary to do operations on some or all of the elements of
@@ -7881,29 +7881,32 @@ Minimum (non-blank) value in the top operand on the
stack, so
image onto the stack. Therefore this operator is mainly intended for data
(for example images), if the top operand is a number, this operator just
returns it without any change. So note that when this operator acts on a
-single image, the output will no longer be an image, but a number.
+single image, the output will no longer be an image, but a number. The
+output of this operand is in the same type as the input.
@item maxvalue
-Maximum (non-blank) value of first operand, similar to @command{minvalue}.
+Maximum (non-blank) value of first operand in the same type, similar to
address@hidden
@item numvalue
-Number of non-blank elements in first operand, similar to
address@hidden
+Number of non-blank elements in first operand in the @code{uint64} type,
+similar to @command{minvalue}.
@item sumvalue
-Sum of non-blank elements in first operand, similar to @command{minvalue}.
+Sum of non-blank elements in first operand in the @code{float32} type,
+similar to @command{minvalue}.
@item meanvalue
-Mean value of non-blank elements in first operand, similar to
address@hidden
+Mean value of non-blank elements in first operand in the @code{float32}
+type, similar to @command{minvalue}.
@item stdvalue
-Standard deviation of non-blank elements in first operand, similar to
address@hidden
+Standard deviation of non-blank elements in first operand in the
address@hidden type, similar to @command{minvalue}.
@item medianvalue
-Median of non-blank elements in first operand, similar to
address@hidden
+Median of non-blank elements in first operand with the same type, similar
+to @command{minvalue}.
@cindex NaN
@item min
@@ -8356,7 +8359,7 @@ $ echo "" | awk '@{print (10.32-3.84)address@hidden'
address@hidden Convolve, Warp, Arithmetic, Image manipulation
address@hidden Convolve, Warp, Arithmetic, Data manipulation
@section Convolve
@cindex Convolution
@@ -9752,7 +9755,7 @@ explanations under the @option{--makekernel} option for
more information.
address@hidden Warp, SubtractSky, Convolve, Image manipulation
address@hidden Warp, SubtractSky, Convolve, Data manipulation
@section Warp
Image warping is the process of mapping the pixels of one image onto a new
pixel grid. This process is sometimes known as transformation, however
@@ -10338,7 +10341,7 @@ of the pixels in the input and output images will be
the same).
address@hidden SubtractSky, , Warp, Image manipulation
address@hidden SubtractSky, , Warp, Data manipulation
@section SubtractSky
@cindex Atmosphere
@@ -11057,13 +11060,15 @@ how the standard deviation on each mesh is finally
found too.
address@hidden Image analysis, Modeling and fittings, Image manipulation, Top
address@hidden Image analysis
address@hidden Data analysis, Modeling and fittings, Data manipulation, Top
address@hidden Data analysis
-Astronomical images contain very valuable information, the tools in
-this section can help in extracting and quantifying that
-information. For example calculating image statistics, or finding the
-sky value or detecting objects within an image.
+Astronomical datasets (images or tables) contain very valuable information,
+the tools in this section can help in analysing, extracting, and
+quantifying that information. For example getting general or specific
+statistics of the dataset (with @ref{Statistics}), detecting signal within
+a noisy dataset (with @ref{NoiseChisel}), or creating a catalog from an
+input dataset (with @ref{MakeCatalog}).
@menu
* Statistics:: Calculate dataset statistics.
@@ -11071,23 +11076,27 @@ sky value or detecting objects within an image.
* MakeCatalog:: Catalog from input and labeled images.
@end menu
address@hidden Statistics, NoiseChisel, Image analysis, Image analysis
address@hidden Statistics, NoiseChisel, Data analysis, Data analysis
@section Statistics
-The distribution of pixel values in an image can give us valuable
-information about the image, for example if it is a positively skewed
-distribution, we can see that there is significant data in the
-image. If the distribution is roughly symmetric, we can tell that
-there is no significant data in the image.
+The distribution of values in a dataset can provide valuable information
+about it. For example, in an image, if it is a positively skewed
+distribution, we can see that there is significant data in the image. If
+the distribution is roughly symmetric, we can tell that there is no
+significant data in the image. In a table, when we need to select a sample
+of objects, it is important to first get a general view of the whole
+sample.
-On the other hand, in some measurements that we do on the image, we
-might need to know the certain statistical parameters of the
-image. For example, if we have run a detection algorithm on an image,
+On the other hand, you might need to know certain statistical parameters of
+the dataset. For example, if we have run a detection algorithm on an image,
and we want to see how accurate it was, one method is to calculate the
-average of the undetected pixels and see how reasonable it is (if
-detection is done correctly, the average of undetected pixels should
-be approximately equal to the background value, see @ref{Sky
-value}). Statistics is built for precisely such situations.
+average of the undetected pixels and see how reasonable it is (if detection
+is done correctly, the average of undetected pixels should be approximately
+equal to the background value, see @ref{Sky value}). In a table, you might
+have calcuated the magnitudes of a certain class of objects and want to get
+some general characteristics of the distribution immediately on the
+command-line (very fast!), to possibly change some parameters. The
+Statistics program is designed for such situations.
@menu
* Histogram and Cumulative Frequency Plot:: Basic definitions.
@@ -11100,137 +11109,169 @@ value}). Statistics is built for precisely such
situations.
@node Histogram and Cumulative Frequency Plot, Sigma clipping, Statistics,
Statistics
@subsection Histogram and Cumulative Frequency Plot
-Histograms and the cumulative frequency plots are both used to study the
-distribution of data. The histogram is mainly easier to understand for the
-untrained eye, while the cumulative frequency plot is more accurate, but
-needs a good level of experience for interpretation.
-
@cindex Histogram
+Histograms and the cumulative frequency plots are both used to visually
+study the distribution of a dataset. A histogram shows the number of data
+points which lie within pre-defined intervals (bins). So on the horizontal
+axis we have the bin centers and on the vertical, the number of points that
+are in that bin. You can use it to get a general view of the distribution:
+which values have been repeated the most? how close/far are the most
+significant bins? Are there more values in the larger part of the range of
+the dataset, or in the lower part? Similarly, many very important
+properties about the dataset can be deduced from a visual inspection of the
+histogram. In the Statistics program, the histogram can be either output to
+a table to plot with your favorite plotting address@hidden recommend
address@hidden://pgfplots.sourceforge.net/,PGFPlots} which generates your plots
+directly within @TeX{} (the same tool that generates your document).}, or
+it can be shown with ASCII characters on the command-line, which is very
+crude, but good enough for a fast and on-the-go analysis, see the example
+in @ref{Invoking aststatistics}.
+
@cindex Intervals, histogram
@cindex Bin width, histogram
@cindex Normalizing histogram
address@hidden Probability distribution function
-A histogram shows the number of data points which lie within
-pre-defined intervals (bins). It is used to get a general view of the
-distribution and its shape. The width of the bins is one of the most
-important parameters for a histogram. In the limiting case that the
-bin-widths tend to zero (and assuming there is data for each bin),
-then the normalized histogram would show the probability distribution
-function of the distribution. Normalizing a histogram means to divide
-the number of data points in each bin by the total number of data.
address@hidden Probability density function
+The width of the bins is only necessary parameter for a histogram. In the
+limiting case that the bin-widths tend to zero (while assuming the number
+of points in the dataset tend to infinity), then the histogram will tend to
+the @url{https://en.wikipedia.org/wiki/Probability_density_function,
+probability density function} of the distribution. When the absolute number
+of points in each bin is not relevant to the study (only the shape of the
+histogram is important), you can @emph{normalize} a histogram so like the
+probability density function, the sum of all its bins will be one.
@cindex Cumulative Frequency Plot
-In the cumulative frequency plot of a distribution, the x axis is the
-sorted data values and the y axis is the index of each data in the
-sorted distribution. Unlike a histogram, a cumulative frequency plot
-does not involve intervals or bins. This makes it less prone to any
-sort of bias or error that a given bin-width would have on the
-analysis. When a larger number of the data points have roughly the
-same value, then the cumulative frequency plot will become steep in
-that vicinity. This occurs because on the x axis (data values), there
-is little change while on the y axis the indexes constantly
-increase. Normalizing a cumulative frequency plot means to divide each
-index (y axis) by the total number of data points.
+In the cumulative frequency plot of a distribution, the horizontal axis is
+the sorted data values and the y axis is the index of each data in the
+sorted distribution. Unlike a histogram, a cumulative frequency plot does
+not involve intervals or bins. This makes it less prone to any sort of bias
+or error that a given bin-width would have on the analysis. When a larger
+number of the data points have roughly the same value, then the cumulative
+frequency plot will become steep in that vicinity. This occurs because on
+the horizontal axis, there is little change while on the vertical axis, the
+indexes constantly increase. Normalizing a cumulative frequency plot means
+to divide each index (y axis) by the total number of data points (or the
+last value).
Unlike the histogram which has a limited number of bins, ideally the
cumulative frequency plot should have one point for every data
-point. Even in small images (for example a @mymath{200\times200}) this
-will result in an unreasonably larger number of points to plot
-(40000)! So when the cumulative frequency plot of an image is stored
-in a text file, it is best to only store its value on a certain number
-of points (intervals) rather than the whole data. The number of points
-to use for the final plot can be specified with the @option{--cfpnum}
-option.
-
-Note that the interval's lower value is considered to be part of each
-interval, but its larger value is not. Formally, an interval between a
-and b is represented by [a, b). This is true for all the intervals
-except the last one. The last interval is closed or [a, b].
+element. Even in small datasets (for example a @mymath{200\times200} image)
+this will result in an unreasonably large number of points to plot (40000)!
+As a result, for practical reasons, it is common to only store its value on
+a certain number of points (intervals) in the input range rather than the
+whole dataset, so you should determine the number of bins you want when
+asking for a cumulative frequency plot. In Gnuastro (and thus the
+Statistics program), the number reported for each bin is the total number
+of datapoints until the larger interval value for that bin. You can see an
+example histogram and cumulative frequency plot of a single dataset under
+the @option{--asciihist} and @option{--asciicfp} options of @ref{Invoking
+aststatistics}.
+
+So as a summary, both the histogram and cumulative frequency plot in
+Statistics will work with bins. Within each bin/interval, the lower value
+is considered to be within then bin (it is inclusive), but its larger value
+is not (it is exclusive). Formally, an interval/bin between a and b is
+represented by [a, b). When the over-all range of the dataset is specified
+(with the @option{--greaterequal}, @option{--lessthan}, or
address@hidden options), the acceptable values of the dataset are also
+defined with a similar inclusive-exclusive manner. But when the range is
+determined from the actual dataset (none of these options is called), the
+last element in the dataset is included in the last bin's count.
@node Sigma clipping, Mirror distribution, Histogram and Cumulative Frequency
Plot, Statistics
@subsection Sigma clipping
Let's assume that you have pure noise (centered on zero) with a clear
-Gaussian distribution, see @ref{Photon counting noise}. Now let's assume
-you add very bright objects (signal) on the image which have a very sharp
-boundary. By a sharp boundary, we mean that there is a clear cutoff at the
-place the objects finish. In other words, at their boundaries, the objects
-do not fade away into the noise. In such a case, when you plot the
-histogram (see @ref{Histogram and Cumulative Frequency Plot}) of the
-distribution, the pixels relating to those objects will be clearly separate
-from pixels that belong to parts of the image that did not have data. In
-the cumulative frequency plot, you would observe a long flat region were
-for a certain range of data (x axis), there is no increase in the index (y
-axis).
address@hidden://en.wikipedia.org/wiki/Normal_distribution,Gaussian
+distribution}, or see @ref{Photon counting noise}. Now let's assume you add
+very bright objects (signal) on the image which have a very sharp
+boundary. By a sharp boundary, we mean that there is a clear cutoff (from
+the noise) at the pixels the objects finish. In other words, at their
+boundaries, the objects do not fade away into the noise. In such a case,
+when you plot the histogram (see @ref{Histogram and Cumulative Frequency
+Plot}) of the distribution, the pixels relating to those objects will be
+clearly separate from pixels that belong to parts of the image that did not
+have any signal (were just noise). In the cumulative frequency plot, after
+a steady rise (due to the noise), you would observe a long flat region were
+for a certain range of data (horizontal axis), there is no increase in the
+index (vertical axis).
@cindex Blurring
@cindex Cosmic rays
@cindex Aperture blurring
@cindex Atmosphere blurring
-In cases like the above, @mymath{\sigma}-clipping is a very useful
-tool to remove the effect (bias) of such forms of signal from the
-distribution. In astronomical applications, cosmic rays (when they
-collide at a near normal incidence angle) are a very good example of
-such signals. The tracks they leave behind in the image are perfectly
-immune to the blurring caused by the atmosphere and the aperture. They
-are also very energetic and so their borders are usually clearly
-separated from the surrounding noise. So @mymath{\sigma}-clipping is
-very useful in removing their effect on the data. See Figure 15 in
-Akhlaghi and Ichikawa (2015).
-
address@hidden is the very simple iteration below. In each
-iteration, the range of input data might decrease and so when the data
-in the image have the conditions above, they will be removed. The
-criteria to stop the iteration will be discussed below.
+Outliers like the example above can significantly bias the measurement of
+noise statistics. @mymath{\sigma}-clipping is defined as a way to avoid the
+effect of such outliers. In astronomical applications, cosmic rays (when
+they collide at a near normal incidence angle) are a very good example of
+such outliers. The tracks they leave behind in the image are perfectly
+immune to the blurring caused by the atmosphere and the aperture. They are
+also very energetic and so their borders are usually clearly separated from
+the surrounding noise. So @mymath{\sigma}-clipping is very useful in
+removing their effect on the data. See Figure 15 in Akhlaghi and Ichikawa,
address@hidden://arxiv.org/abs/1505.01664,2015}.
+
address@hidden is defined as the very simple iteration below. In
+each iteration, the range of input data might decrease and so when the
+outliers have the conditions above, the outliers will be removed through
+this iteration. The exit criteria will be discussed below.
@enumerate
@item
-Calculate the mean, standard deviation (@mymath{\sigma}) and median
-(@mymath{m}) of a distribution.
+Calculate the standard deviation (@mymath{\sigma}) and median (@mymath{m})
+of a distribution.
@item
Remove all points that are smaller or larger than
@mymath{m\pm\alpha\sigma}.
address@hidden
+Go back to step 1, unless the selected exit criteria is reached.
@end enumerate
@noindent
-The reason the median is used as a reference and not the mean is that
-the mean is too significantly affected by the presence of signal,
-while the median is less affected, see @ref{Finding the sky value}. As
-you can tell from this algorithm, besides the condition above (that
-the signal have clear high signal to noise boundaries)
address@hidden is only useful when the signal does not cover
-more than half of the full data set. If they do, then the median will
-lie over the signal and sigma clipping might remove the pixels with no
-signal.
-
-In the literature researchers use either one of two criteria to stop
-this iteration above:
+The reason the median is used as a reference and not the mean is that the
+mean is too significantly affected by the presence of outliers, while the
+median is less affected, see @ref{Finding the sky value}. As you can tell
+from this algorithm, besides the condition above (that the signal have
+clear high signal to noise boundaries) @mymath{\sigma}-clipping is only
+useful when the signal does not cover more than half of the full data
+set. If they do, then the median will lie over the outliers and
address@hidden might remove the pixels with no signal.
+
+There are commonly two exit criteria to stop the @mymath{\sigma}-clipping
+iteration:
@itemize
@item
-When a certain number of iterations has taken place.
+When a certain number of iterations has taken place (second value to the
address@hidden option is larger than 1).
@item
When the new measured standard deviation is within a certain tolerance
-level of the old one. The tolerance level is defined by:
+level of the old one (second value to the @option{--sigclip} option is less
+than 1). The tolerance level is defined by:
@dispmath{\sigma_{old}-\sigma_{new} \over \sigma_{new}}
-The standard deviation is heavily influenced by the presence of
-outliers. Therefore the fact that it stops changing between two
-iterations is a sign that we have successfully removed outliers. Note
-that in each clipping, the dispersion in the distribution is either
-less or equal. So @mymath{\sigma_{old}\geq\sigma_{new}}.
+The standard deviation is used because it is heavily influenced by the
+presence of outliers. Therefore the fact that it stops changing between two
+iterations is a sign that we have successfully removed outliers. Note that
+in each clipping, the dispersion in the distribution is either less or
+equal. So @mymath{\sigma_{old}\geq\sigma_{new}}.
@end itemize
@cartouche
@noindent
-Other objects in astronomical data analysis like galaxies and stars
-are blurred by the atmosphere and the telescope aperture. Galaxies in
-particular do not appear to have a clear high signal to noise cutoff
-at all. Therefore @mymath{\sigma}-clipping will not be useful in
-removing their effect on the data. In astronomy, it is only useful for
-removing the effect of Cosmic rays.
+When working on astronomical images, objects like galaxies and stars are
+blurred by the atmosphere and the telescope aperture, therefore their
+signal sinks into the noise very gradually. Galaxies in particular do not
+appear to have a clear high signal to noise cutoff at all. Therefore
address@hidden will not be useful in removing their effect on the
+data.
+
+To guage if @mymath{\sigma}-clipping will be useful for your dataset, look
+at the histogram (see @ref{Histogram and Cumulative Frequency Plot}). The
+ASCII histogram that is printed on the command-line with
address@hidden is good enough in most cases.
@end cartouche
@node Mirror distribution, Invoking aststatistics, Sigma clipping, Statistics
@@ -11331,8 +11372,8 @@ is called, the given ranges will also shift
accordingly).
@node Invoking aststatistics, , Mirror distribution, Statistics
@subsection Invoking Statistics
-Statistics will print the major statistical measures of an image's pixel
-value distribution. The executable name is @file{aststatistics} with the
+Statistics will print statistical measures of an input dataset (table
+column or image). The executable name is @file{aststatistics} with the
following general template
@example
@@ -11343,278 +11384,347 @@ $ aststatistics [OPTION ...] InputImage.fits
One line examples:
@example
-$ aststatistics input.fits
-$ aststatistics animage.fits --ignoremin --nohist
+$ aststatistics image.fits
+$ aststatistics tab.fits -c/MAG/ --histogram
+$ aststatistics catalog.fits -h1 --column=MAG_F160W
+$ aststatistics tab.txt -cMAG_F160W --median --std
+$ aststatistics cat.fits -cMAG -g26 -l27 --sigmaclip=3,0.2
@end example
@noindent
-If Statistics is to do any data processing, an input image should be
-provided with the recognized extensions as input data, see
address@hidden See @ref{Common options} for the list of options that are
-shared by all programs. All the main statistical operations have their
-specific set of options. If a string is given to the @option{--output}
-option, it is used as the base name for the generated files. Without this
-option, the input image name is used as the name-base.
-
-Some of the options are necessary and if they are not included in the
-configuration file, Statistics will not run, see @ref{Configuration
-files}. However, for some others this is not so: @option{--histmin},
address@hidden, @option{--histquant}, @option{--cfpmin},
address@hidden, @option{--cfpquant}. These are options to do with the
-range of values in the histogram and cumulative frequency plots. If no
-value is given for these options when Statistics is about to start
-processing the data, then the full data range will be used. Such that the
-minimum image value will be set for the minimums and the maximum image
-value will be used for the maximum.
+An input image or table is necessary when processing is to be done. If a
+name is given to the @option{--output} option, it is used as the base name
+for the generated files (when applicable). Without @option{--output}, the
+input name will be used to generate an output name, see @ref{Automatic
+output}. The options described below are particular to Statistics, but for
+general operations, it shares a large collection of options with the other
+Gnuastro programs, see @ref{Common options}. Options can also be given in
+configuration files, for more, please see @ref{Configuration files}.
+
+The input dataset may have blank values (see @ref{Blank pixels}), in this
+case, all blank pixels are ignored during the calculation. Initially, the
+full dataset will be read, but it is possible to select a specific range of
+data elements to use in the analysis of each run. You can either directly
+specify a minimum and maximum value for the range of data elements to use
+(with @option{--greaterequal} or @option{--lessthan}), or specify the range
+using quantiles (with @option{--qrange}). If a range is specified all
+pixels outside of it are ignored before any processing.
@cindex ASCII plot
-By default, in verbose mode @footnote{If the @option{-q} option is not
-called then all programs will operate in verbose mode, see @ref{Common
-options}.}, along with a short summary of the basic data statistics, a
-simple ASCII histogram will also be printed. This can be useful for a very
-quick and general view of the distribution. An example verbose output of
-Statistics on one of the @command{$ make check} outputs can be seen below:
+When no operation is requested, Statistics will print some general basic
+properties of the input dataset on the command-line like the example below
+(ran on one of the output images of @command{make check}). It includes the
+basic information about the input dataset: filename, FITS extension (if a
+FITS file), the column name (if it was a table) or number if it the column
+doesn't have a name, along with the range and units of the data (if they
+were present in the input). Following that, the basic statistics of the
+dataset like number of points, minimum, maximum, and etc are printed. The
+outputs are finalized with an ASCII histogram to give you a general feeling
+of how the data are distributed in the given range.
@example
-$ aststatistics ./tests/convolve_spatial_warped_noised.fits \
- --histquant=0.05
-Statistics started on AAA BBB CC DD:EE:FF GGGG
- - Input read: ./tests/convolve_spatial_warped_noised.fits (hdu: 0)
- -- Number of points 10000
- -- Minimum -38.2066
- -- Maximum 1268.72
- -- Sum 154927
- -- Mean 15.4927
- -- Standard deviation 60.5407
- -- Median 4.82691
- -- Mode (quantile, value) 0.4335, 2.90301
- -- Mode symmetricity and its cutoff value 0.5909, 17.2507
- -- ASCII histogram in the range: -13.912957 -- 65.058487:
- | * *
- | ********
- | ********* **
- | **************
- | ******************
- | ********************
- | ***********************
- |***************************
- |******************************
- |*****************************************
- |************************************************************
- |------------------------------------------------------------
-
- - Sigma clipping results (Median, Mean, STD, Number):
- - 4.00 times sigma by convergence (tolerance: 0.2000):
- 1: 4.826907, 15.492665, 60.540707, 10000
- 2: 4.665353, 10.117773, 27.793823, 9881
- 3: 4.433090, 7.458610, 18.510950, 9715
- 4: 4.216213, 6.176818, 15.231371, 9575
- 5: 4.088199, 5.679162, 14.183748, 9502
- - 4.00 sigma-clipping 5 times:
- 1: 4.826907, 15.492665, 60.540707, 10000
- 2: 4.665353, 10.117773, 27.793823, 9881
- 3: 4.433090, 7.458610, 18.510950, 9715
- 4: 4.216213, 6.176818, 15.231371, 9575
- 5: 4.088199, 5.679162, 14.183748, 9502
-Statistics finished in: 0.006964 (seconds)
+Statistics (GNU Astronomy Utilities) X.X
+-------
+Input: convolve_spatial_scaled_noised.fits (hdu: 0)
+Range: from (inclusive) 9500, upto (exclusive) 11000.
+Unit: Brightness
+-------
+ Number of elements: 9074
+ Minimum: 9622.35
+ Maximum: 10999.8
+ Median: 10093.7
+ Mean: 10144
+ Standard deviation: 221.81
+-------
+Histogram:
+ | **
+ | ******
+ | *******
+ | *********
+ | *************
+ | **************
+ | ******************
+ | ********************
+ | *************************** *
+ | ***************************************** ***
+ |* **************************************************************
+ |-----------------------------------------------------------------
@end example
+The following set of options are for specifying the input/outputs of
+Statistics. Recall that besides these options that are only meaningful in
+Statistics, there are many other options that are common to all Gnuastro
+programs including Statistics, see @ref{Common options} for those.
+
@table @option
address@hidden -r
address@hidden --ignoremin
-Ignore all data elements that have a value equal to the minimum value in
-the image. In practice this is like setting them as @ref{Blank pixels},
-their values will not be used.
address@hidden -c STR/INT
address@hidden --column=STR/INT
+The column selector when the input file is a table, plain text, FITS ASCII,
+and FITS binary tables are all acceptable. See @ref{Selecting table
+columns} for a full description of how to use this option. For more on how
+tables are read in Gnuastro, please see @ref{Tables in Gnuastro}.
address@hidden -l
address@hidden --lowerbin
-Set the first column of the histogram and cumulative frequency plots
-to the lower interval boundary. By default (without calling this
-option), the central interval value is used.
address@hidden -g FLT
address@hidden --greaterequal=FLT
+Limit the range of inputs into those with values greater and equal to what
+is given to this option. None of the values below this value will be used
+in any of the processing steps below.
address@hidden --onebinvalue=FLT
-Make sure that one bin starts with the value to this option. In practice,
-this will shift the bins used to find the histogram and cumulative
-frequency plot such that one bin's lower interval becomes this value. For
-example when the histogram range includes negative values, but the data
-doesn't. If zero is somewhere between one bin, then the viewers of the
-plot(s) will think negative data is also present. By setting
address@hidden, you can make sure that the viewers of the
-histogram will not be confused.
address@hidden -l FLT
address@hidden --lessthan=FLT
+Limit the range of inputs into those with values less-than what is given to
+this option. None of the values greater or equal to this value will be used
+in any of the processing steps below.
+
address@hidden -Q FLT[,FLT]
address@hidden --qrange=FLT[,FLT]
+Specify the range of usable inputs using the quantile. This option can take
+one or two quantiles to specify the range. When only one number is input
+(let's call it @mymath{Q}), the range will be those values in the quantile
+range @mymath{Q} to @mymath{1-Q}. So when only one value is given, it must
+be less than 0.5. When two values are given, the first is used as the lower
+quantile range and the second is used as the larger quantile range.
address@hidden NaN
-Note that by default, the first row of the histogram and cumulative
-frequency plot show the central values of each bin. So in the example
-above you will not see the 0.000. To see it, add the
address@hidden option to show the lower value of each bin. If you
-don't care about the bin positions within the specified range you can
-set the value to this option to a Not-a-Number (NaN) value on the
-command-line (@option{--onebinvalue=nan}) or in the configuration
-files with a @code{nan} following the option name. If the value is not
-within the specified bin range, it will be ignored.
-
address@hidden --noasciihist
-Do not show an ASCII plot on the command-line.
-
address@hidden --mirrorquant=FLT
-Quantile to put the mirror. A value between 0 and 1. See @ref{Mirror
-distribution} for a complete explanation. Outputs two files with suffixes
address@hidden and @file{_mirrorcfp.txt}.
-
address@hidden --checkmode
-The mode of the data is found by comparing the input data distribution with
-its mirror distribution. If this option is called, the mirror
-distribution's histogram and cumulative frequency plots will be saved in to
-plain text files ending with @file{_modehist.txt} and
address@hidden See the explanation for @ref{Mirror distribution} for
-more details about these two files and how you can easily plot the
-outputs. This option only works when when Statistics is in verbose
-mode. Since otherwise the mode is not calculated.
-
-To draw the plots you can use the script in @ref{Mirror
-distribution}. Just change the appended suffixes in the two calls to
address@hidden in the Python script.
-
address@hidden --histrangeformirror
-Use @option{--histmin} and @option{--histmax} for the range of the
-mirror distributions (which are produced with the
address@hidden and @option{--checkmode} options).
address@hidden Quantile
+The quantile of a given element in a dataset is defined by the fraction of
+its index to the total number of values in the sorted input array. So the
+smallest and largest values in the dataset have a quantile of 0.0 and
+1.0. The quantile is a very useful non-parametric (making no assumptions
+about the input) relative measure to specify a range. It can best be
+understood in terms of the cumulative frequency plot, see @ref{Histogram
+and Cumulative Frequency Plot}. The quantile of each horizontal axis value
+in the cumulative frequency plot is the vertical axis value associate with
+it.
@end table
address@hidden
address@hidden:} The stored histogram is stored in a text file
-ending with @file{_hist.txt}.
address@hidden AWK
address@hidden GNU AWK
+When you only need a single valued statistical measurement, you can use the
+group of options below. They will print only the requested value as one
+field in a line/row, like the @option{--mean}, @option{--median}
+options. These options can be called any number of times and in any
+order. The outputs of all such options will be printed on one line
+following each other (with a space character between them). This feature
+makes these options very useful in scripts, or to redirect into programs
+like GNU AWK for further immediate higher-level processing, or when you
+don't want the clutter of all the information when no option is
+specified. These are some of the most basic measures, Gnuastro is still
+under heavy development, if you want another statistical parameter, please
+contact us and we will do out best to add it to this list.
@table @option
address@hidden --nohist
-Do not calculate or save the histogram.
address@hidden -n
address@hidden --number
+Print the number of all used elements.
address@hidden --normhist
-Make a normalized histogram, see @ref{Histogram and Cumulative Frequency
-Plot}.
address@hidden --minimum
+Print the minimum value of all used elements.
address@hidden --maxhistone
-Divide all histogram bins by the number in the bin with the most data
-points. This is very useful if you want to plot the histogram along
-with a normalized cumulative frequency plot in one plot. Note that if
-the histogram numbers are important in showing along with the
-cumulative frequency plot, you can use @option{--maxcfpeqmaxhist}, see
-below.
address@hidden --maximum
+Print the maximum value of all used elements.
address@hidden -n INT
address@hidden --histnumbins=INT
-The number of bins in the histogram. Note that in practice, this is also
-equivalent to the number of rows in the output text file.
-
address@hidden -i FLT
address@hidden --histmin=FLT
-The minimum value to use in the histogram. If @option{--histquant} is
-given, then any value given @option{--histmin} or @option{--histmax} is
-ignored.
address@hidden --sum
+Print the sum of all used elements.
address@hidden -x FLT
address@hidden --histmax=FLT
-The maximum value to use in the histogram. Similar to @option{--histmin}.
address@hidden -m
address@hidden --mean
+Print the mean (average) of all used elements.
address@hidden -Q FLT
address@hidden --histquant=FLT
-Set the range of the histogram based on the image quantile. So
address@hidden is given, all the data from the 0.05 quantile to
-0.95 quantile will be used in the histogram. This is useful when there is a
-small number of outliers in the image. Note that if this option is given,
-any (possible) value given to @option{--histmin} or @option{--histmax} are
-ignored.
address@hidden -t
address@hidden --std
+Print the standard deviation of all used elements.
address@hidden table
address@hidden -M
address@hidden --median
+Print the median of all used elements.
address@hidden
address@hidden Frequency Plot:} The cumulative frequency plot will
-be stored in a text file ending with @file{_cfp.txt}. To be more
-realistic, the average of the indexs in each interval is used as the
-second column, see @ref{Histogram and Cumulative Frequency Plot}.
address@hidden -O
address@hidden --mode
+Print the mode of all used elements.
address@hidden @option
address@hidden table
address@hidden --nocfp
-Do not calculate or store the cumulative frequency plot.
+The list of options below are for those statistical operations that output
+more than one value. So while they can be called together in one run, their
+outputs will be distinct (each one's output will usually take in more than
+one line).
address@hidden --normcfp
-Normalize the cumulative frequency plot, see @ref{Histogram and Cumulative
-Frequency Plot}.
address@hidden @option
address@hidden --maxcfpeqmaxhist
-Set the maximum cumulative frequency plot value to the maximum value
-in the histogram (if it is to be created). This is a useful in
-plotting the histogram and cumulative frequency plots together when
-the histogram numbers are important.
address@hidden -A
address@hidden --asciihist
+Print an ASCII histogram of the usable values within the input dataset
+along with some basic information like the example below (from the UVUDF
address@hidden@url{https://asd.gsfc.nasa.gov/UVUDF/uvudf_rafelski_2015.fits.gz}}).
The
+width and height of the histogram (in units of character widths and heights
+on your command-line terminal) can be set with the @option{--numasciibins}
+(for the width) and @option{--asciiheight} options.
+
+For a full description of the histogram, please see @ref{Histogram and
+Cumulative Frequency Plot}. An ASCII plot is certainly very crude and
+cannot be used in any publication, but it is very useful for getting a
+general feeling of the input dataset very fast and easily on the
+command-line without having to take your hands off the keyboard (which is a
+major distraction!). If you want to try it out, you can write it all in one
+line and ignore the @key{\} and extra spaces.
address@hidden --cfpsimhist
-Set the range of the cumulative frequency plot and the number of
-points to store to the same range as the histogram. If the two are to
-be plotted together, this is very useful, since the first axis
-(column) of the two will become identical.
address@hidden
+$ aststatistics uvudf_rafelski_2015.fits.gz --hdu=1 \
+ --column=MAG_F160W --lessthan=40 \
+ --asciihist --numasciibins=55
+ASCII Histogram:
+Number: 8593
+Y: (linear: 0 to 660)
+X: (linear: 17.7735 -- 31.4679, in 55 bins)
+ | ****
+ | *****
+ | ******
+ | ********
+ | *********
+ | ***********
+ | **************
+ | *****************
+ | ***********************
+ | ********************************
+ |*** ***************************************************
+ |-------------------------------------------------------
address@hidden example
address@hidden -p INT
address@hidden --cfpnum=INT
-The number of points to store the cumulative frequency plot. They will be
-evenly distributed between the range of pixel values.
address@hidden --asciicfp
+Print the cumulative frequency plot of the usable elements in the input
+dataset. Please see descriptions under @option{--asciihist} for more, the
+example below is from the same input table as that example. To better
+understand the cumulative frequency plot, please see @ref{Histogram and
+Cumulative Frequency Plot}.
address@hidden -a FLT
address@hidden --cfpmin=FLT
-The minimum value to use for the cumulative pixel value range. If
address@hidden is given, then any value given @option{--cfpmin} or
address@hidden is ignored.
address@hidden
+$ aststatistics uvudf_rafelski_2015.fits.gz --hdu=1 \
+ --column=MAG_F160W --lessthan=40 \
+ --asciicfp --numasciibins=55
+ASCII Cumulative frequency plot:
+Y: (linear: 0 to 8593)
+X: (linear: 17.7735 -- 31.4679, in 55 bins)
+ | *******
+ | **********
+ | ***********
+ | *************
+ | **************
+ | ***************
+ | *****************
+ | *******************
+ | ***********************
+ | ******************************
+ |*******************************************************
+ |-------------------------------------------------------
address@hidden example
address@hidden -n FLT
address@hidden --cfpmax=FLT
-The maximum value to use for the cumulative pixel value range. Similar to
address@hidden
address@hidden -H
address@hidden --histogram
+Save the histogram of the usable values in the input dataset into a
+table. The first column is the value at the center of the bin and the
+second is the number of points in that bin. If the @option{--cumulative}
+option is also called with this option in a run, then the table will have
+three columns (the third is the cumulative frequency plot). Through the
address@hidden and @option{--lowerbin} you can modify the first column
+values and with @option{--normalize} and @option{--maxbinone} you can
+modify the second columns. See below for the description of each.
+
+By default (when no @option{--output} is specified) a plain text table will
+be created, see @ref{Gnuastro text table format}. If a FITS name is
+specified, you can use the common option @option{--tableformat} to have it
+as a FITS ASCII or FITS binary format, see @ref{Common options}. This table
+can then be fed into your favorite plotting tool and get a much more clean
+and nice histogram than what the raw command-line can offer you (with the
address@hidden option).
address@hidden -U FLT
address@hidden --cfpquant=FLT
-Similar to @option{--histquant} but for the cumulative frequency plot.
address@hidden -C
address@hidden --cumulative
+Save the cumulative frequency plot of the usable values in the input
+dataset into a table, similar to @option{--histogram}.
+
address@hidden -s FLT,FLT
address@hidden --sigmaclip=FLT,FLT
+Do @mymath{\sigma}-clipping on the usable pixels of the input dataset. See
address@hidden clipping} for a full description on @mymath{\sigma}-clipping and
+also to better understand this option.
+
+This option takes two values which are separated by a comma (@key{,}). Each
+value can either be written as a single number or as a fraction of two
+numbers (for example @code{3,1/10}). The first value to this option is the
+multiple of @mymath{\sigma} that will be clipped (@mymath{\alpha} in that
+section). The second value is the exit criteria. If it is less than 1, then
+it is interpretted as tolerance and if it is larger than one it is a
+specific number. Hence, in the latter case the value must be an integer.
@end table
address@hidden
address@hidden@mymath{\sigma}-clipping:} The result of each iteration of
address@hidden will be printed in the terminal for both
-types of sigma clipping: A certain number of times and convergence of
-the standard deviation.
+The final set of options in Statistics describe the histogram and
+cumulative frequency plot settings (for the @option{--histogram},
address@hidden, @option{--asciihist}, and @option{--asciicfp}
+options).
@table @option
address@hidden --nosigclip
-If this option is called, no @mymath{\sigma}-clipping will take place.
address@hidden --numbins
+The number of bins (rows) to use in the histogram and the cumulative
+frequency plot tables (outputs of @option{--histogram} and
address@hidden).
address@hidden -u FLT
address@hidden --sigclipmultip=FLT
-The multiple of the standard deviation above which to clip. This value is
-demonstrated by @mymath{\alpha} in @ref{Sigma clipping}.
address@hidden --numasciibins
+The number of bins (characters) to use in the ASCII plots when printing the
+histogram and the cumulative frequency plot (outputs of
address@hidden and @option{--asciicfp}).
address@hidden -t FLT
address@hidden --sigcliptolerance=FLT
-If the fractional difference of the standard deviation becomes less than
-this value, then @mymath{\sigma}-clipping will halt, see @ref{Sigma
-clipping}.
address@hidden --asciiheight
+The number of lines to use when printing the ASCII histogram and cumulative
+frequency plot on the command-line (outputs of @option{--asciihist} and
address@hidden).
+
address@hidden -n
address@hidden --normalize
+Normalize the histogram or cumulative frequency plot tables (outputs of
address@hidden and @option{--cumulative}). For a histogram, the sum
+of all bins will become one and for a cumulative frequency plot the last
+bin value will be one.
+
address@hidden --maxbinone
+Divide all the histogram values by the maximum bin value so it becomes one
+and the rest are similarly scaled. In some situations (for example if you
+want to plot the histogram and cumulative frequency plot in one plot) this
+can be very useful.
+
address@hidden --onebinvalue=FLT
+Make sure that one bin starts with the value to this option. In practice,
+this will shift the bins used to find the histogram and cumulative
+frequency plot such that one bin's lower interval becomes this value. For
+example when the histogram range includes negative and positive values and
+zero has a special significance in your analysis, then zero will be
+somewhere in one bin and will mix counts of positive and negative. By
+setting @option{--onebinvalue=0}, you can make sure that the viewers of the
+histogram will not be confused without doing the math of setting a range
+and number of bins.
+
address@hidden NaN
+Note that by default, the first row of the histogram and cumulative
+frequency plot show the central values of each bin. So in the example above
+you will not see the 0.000 in the first column, you will see two symmetric
+values. If the value is not within the usable input range, this option will
+be ignored.
address@hidden -g INT
address@hidden --sigclipnum=INT
-The number of iterations for the case where the @mymath{\sigma}-clipping
-iteration stops after a certain number of runs.
@end table
address@hidden NoiseChisel, MakeCatalog, Statistics, Image analysis
+
address@hidden NoiseChisel, MakeCatalog, Statistics, Data analysis
@section NoiseChisel
Once raw data have gone through the initial reduction process (through the
-programs in @ref{Image manipulation}). We are ready to derive scientific
+programs in @ref{Data manipulation}). We are ready to derive scientific
results out of them. Unfortunately in most cases, the scientifically
interesting targets are deeply drowned in a sea of noise. NoiseChisel is
Gnuastro's tool to detect signal in noise. In fact, NoiseChisel was the
@@ -12230,7 +12340,7 @@ pixel.
address@hidden MakeCatalog, , NoiseChisel, Image analysis
address@hidden MakeCatalog, , NoiseChisel, Data analysis
@section MakeCatalog
Detecting and segmenting signal over an image results in labeled
@@ -13302,7 +13412,7 @@ the actual columns.
address@hidden Modeling and fittings, High-level calculations, Image analysis,
Top
address@hidden Modeling and fittings, High-level calculations, Data analysis,
Top
@chapter Modeling and fitting
@cindex Fitting
@@ -14802,14 +14912,14 @@ were of integer types.
@node High-level calculations, Libraries, Modeling and fittings, Top
@chapter High-level calculations
-After the reduction of raw data (for example with the programs in
address@hidden manipulation}) you will have reduced images/data ready for
-processing/analyzing (for example with the programs in @ref{Image
-analysis}). But the processed/analyzed data (or catalogs) are still
-not enough to derive any scientific result. Even higher-level analysis
-is still needed to convert the observed magnitudes, sizes or volumes
-into physical quantities that we associate with each catalog entry or
-detected object which is the purpose of the tools in this section.
+After the reduction of raw data (for example with the programs in @ref{Data
+manipulation}) you will have reduced images/data ready for
+processing/analyzing (for example with the programs in @ref{Data
+analysis}). But the processed/analyzed data (or catalogs) are still not
+enough to derive any scientific result. Even higher-level analysis is still
+needed to convert the observed magnitudes, sizes or volumes into physical
+quantities that we associate with each catalog entry or detected object
+which is the purpose of the tools in this section.
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 4490a3f..68d301e 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -260,189 +260,6 @@ arithmetic_check_float_input(gal_data_t *in, int
operator, char *numstr)
/*************** Unary functions **************/
/***********************************************************************/
-/* Note that for floating point types b!=b (by definition of NaN). And in
- such cases, even if there are blank values, the smaller and larger
- condition checked will fail, therefore the final result will be what we
- want (to ignore the blank values). */
-#define UNIFUNC_MINVALUE(IT) { \
- IT *oa=o->array; \
- int blankeq = (b==b && gal_blank_present(in)) ? 1 : 0; \
- if(blankeq) \
- do if(*ia!=b) *oa = *ia < *oa ? *ia : *oa; while(++ia<iaf); \
- else \
- do *oa = *ia < *oa ? *ia : *oa; while(++ia<iaf); \
- }
-
-
-#define UNIFUNC_MAXVALUE(IT) { \
- IT *oa=o->array; \
- int blankeq = (b==b && gal_blank_present(in)) ? 1 : 0; \
- if(blankeq) \
- do if(*ia!=b) *oa = *ia > *oa ? *ia : *oa; while(++ia<iaf); \
- else \
- do *oa = *ia > *oa ? *ia : *oa; while(++ia<iaf); \
- }
-
-
-#define UNIFUNC_NUMVALUE { \
- uint64_t *oa=o->array, num=0; \
- if( gal_blank_present(in) ) \
- { \
- if(b==b) /* Is integer type. */ \
- do if(*ia!=b) ++num; while(++ia<iaf); \
- else /* Is float type with NaN blank. */ \
- do if(*ia==*ia) ++num; while(++ia<iaf); \
- } \
- else num=in->size; \
- *oa=num; \
- }
-
-
-#define UNIFUNC_SUMVALUE { \
- double sum=0.0f; \
- float *oa=o->array; \
- if( gal_blank_present(in) ) \
- { \
- if(b==b) /* Is integer type. */ \
- do if(*ia!=b) sum += *ia; while(++ia<iaf); \
- else /* Is float type with NaN blank. */ \
- do if(*ia==*ia) sum += *ia; while(++ia<iaf); \
- } \
- else \
- do sum += *ia; while(++ia<iaf); \
- *oa=sum; \
- }
-
-
-#define UNIFUNC_MEANVALUE { \
- int64_t num=0; \
- double sum=0.0f; \
- float *oa=o->array; \
- if( gal_blank_present(in) ) \
- { \
- if(b==b) /* Is integer type. */ \
- do if(*ia!=b) { ++num; sum += *ia; } while(++ia<iaf); \
- else /* Is float type with NaN blank. */ \
- do if(*ia==*ia) { ++num; sum += *ia; } while(++ia<iaf); \
- } \
- else \
- { \
- do sum += *ia; while(++ia<iaf); \
- num=in->size; \
- } \
- *oa=sum/num; \
- }
-
-
-#define UNIFUNC_STDVALUE { \
- int64_t n=0; \
- float *oa=o->array; \
- double s=0.0f, s2=0.0f; \
- if( gal_blank_present(in) ) \
- { \
- if(b==b) /* Is integer type. */ \
- do if(*ia!=b) \
- { ++n; s += *ia; s2 += *ia * *ia; } while(++ia<iaf); \
- else /* Is float type with NaN blank. */ \
- do if(*ia==*ia) \
- { ++n; s += *ia; s2 += *ia * *ia; } while(++ia<iaf); \
- } \
- else \
- { \
- do { s += *ia; s2 += *ia * *ia; } while(++ia<iaf); \
- n=in->size; \
- } \
- *oa=sqrt( (s2-s*s/n)/n ); \
- }
-
-
-#define UNIFUNC_MEDIANVALUE(IT) { \
- size_t n=0; \
- gal_data_t *noblank, *sorted; \
- IT *u=NULL, *s, *a, *oa=o->array; \
- \
- /* After this, there is no more blanks: only `noblank' is used.*/ \
- if( gal_blank_present(in) ) \
- { \
- if( flags & GAL_ARITHMETIC_FREE ) \
- { \
- gal_blank_remove(in); \
- noblank=in; \
- } \
- else \
- { \
- u=a=gal_data_malloc_array(in->type, in->size); \
- if(b==b) do if (*ia!=b) { *a++=*ia; ++n;} while(++ia<iaf); \
- else do if (*ia==*ia) { *a++=*ia; ++n;} while(++ia<iaf); \
- noblank=gal_data_alloc(u, in->type, 1, &n, NULL, 0, \
- in->minmapsize, NULL, NULL, NULL); \
- } \
- } \
- else noblank=in; \
- \
- /* After this, the array is sorted, only `sorted' is used. */ \
- if( gal_statistics_is_sorted(noblank) ) sorted=noblank; \
- else \
- { \
- if( flags & GAL_ARITHMETIC_FREE ) sorted=noblank; \
- else \
- { \
- if(u) sorted=noblank; /* New space is already allocated.*/ \
- else sorted=gal_data_copy(noblank); \
- } \
- gal_statistics_sort_increasing(sorted); \
- } \
- \
- /* Find the median. */ \
- n=sorted->size; \
- s=sorted->array; \
- *oa = (sorted->size)%2 ? s[n/2] : (s[n/2] + s[n/2-1])/2 ; \
- \
- /* Clean up. If `sorted' doesn't equal `in', then it was */ \
- /* allocated in this block and must be freed. */ \
- if( sorted!=in ) gal_data_free(sorted); \
- }
-
-
-
-
-
-#define UNIFUNC_RUN_FUNCTION_ON_ARRAY(IT){ \
- IT b, *ia=in->array, *iaf=ia + in->size; \
- gal_blank_write(&b, in->type); \
- switch(operator) \
- { \
- case GAL_ARITHMETIC_OP_MINVAL: \
- UNIFUNC_MINVALUE(IT); \
- break; \
- case GAL_ARITHMETIC_OP_MAXVAL: \
- UNIFUNC_MAXVALUE(IT); \
- break; \
- case GAL_ARITHMETIC_OP_NUMVAL: \
- UNIFUNC_NUMVALUE; \
- break; \
- case GAL_ARITHMETIC_OP_SUMVAL: \
- UNIFUNC_SUMVALUE; \
- break; \
- case GAL_ARITHMETIC_OP_MEANVAL: \
- UNIFUNC_MEANVALUE; \
- break; \
- case GAL_ARITHMETIC_OP_STDVAL: \
- UNIFUNC_STDVALUE; \
- break; \
- case GAL_ARITHMETIC_OP_MEDIANVAL: \
- UNIFUNC_MEDIANVALUE(IT); \
- break; \
- default: \
- error(EXIT_FAILURE, 0, "the operator code %d is not " \
- "recognized in UNIFUNC_RUN_FUNCTION_ON_ARRAY", operator); \
- } \
- }
-
-
-
-
-
#define UNIFUNC_RUN_FUNCTION_ON_ELEMENT(IT, OP){ \
IT *ia=in->array, *oa=o->array, *iaf=ia + in->size; \
do *oa++ = OP(*ia++); while(ia<iaf); \
@@ -494,89 +311,18 @@ arithmetic_check_float_input(gal_data_t *in, int
operator, char *numstr)
-#define UNIARY_FUNCTION_ON_ARRAY \
- switch(in->type) \
- { \
- case GAL_DATA_TYPE_UINT8: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY( uint8_t ); \
- break; \
- case GAL_DATA_TYPE_INT8: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY( int8_t ); \
- break; \
- case GAL_DATA_TYPE_UINT16: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY( uint16_t ); \
- break; \
- case GAL_DATA_TYPE_INT16: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY( int16_t ); \
- break; \
- case GAL_DATA_TYPE_UINT32: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY( uint32_t ); \
- break; \
- case GAL_DATA_TYPE_INT32: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY( int32_t ); \
- break; \
- case GAL_DATA_TYPE_UINT64: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY( uint64_t ); \
- break; \
- case GAL_DATA_TYPE_INT64: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY( int64_t ); \
- break; \
- case GAL_DATA_TYPE_FLOAT32: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY( float ); \
- break; \
- case GAL_DATA_TYPE_FLOAT64: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY( double ); \
- break; \
- default: \
- error(EXIT_FAILURE, 0, "type code %d not recognized in " \
- "`UNIARY_FUNCTION_ON_ARRAY'", in->type); \
- }
-
-
-
-
-
static gal_data_t *
arithmetic_unary_function(int operator, unsigned char flags, gal_data_t *in)
{
gal_data_t *o;
- size_t dsize=1;
/* If we want inplace output, set the output pointer to the input
pointer, for every pixel, the operation will be independent. */
- switch(operator)
- {
-
- /* Operators with only one value as output. */
- case GAL_ARITHMETIC_OP_MINVAL:
- case GAL_ARITHMETIC_OP_MAXVAL:
- case GAL_ARITHMETIC_OP_MEDIANVAL:
- o = gal_data_alloc(NULL, in->type, 1, &dsize, NULL, 0, -1,
- NULL, NULL, NULL);
- if(operator==GAL_ARITHMETIC_OP_MINVAL) /* For the min and max */
- gal_data_type_max(o->type, o->array); /* operators we need to */
- else if (operator==GAL_ARITHMETIC_OP_MAXVAL) /* start with the max */
- gal_data_type_min(o->type, o->array); /* and min values */
- break; /* respectively. */
- case GAL_ARITHMETIC_OP_NUMVAL:
- o = gal_data_alloc(NULL, GAL_DATA_TYPE_UINT64, 1, &dsize,
- NULL, 0, -1, NULL, NULL, NULL);
- break;
- case GAL_ARITHMETIC_OP_SUMVAL:
- case GAL_ARITHMETIC_OP_MEANVAL:
- case GAL_ARITHMETIC_OP_STDVAL:
- o = gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, 1, &dsize,
- NULL, 0, -1, NULL, NULL, NULL);
- break;
-
- /* The other operators */
- default:
- if(flags & GAL_ARITHMETIC_INPLACE)
- o = in;
- else
- o = gal_data_alloc(NULL, in->type, in->ndim, in->dsize, in->wcs,
- 0, in->minmapsize, NULL, NULL, NULL);
- }
+ if(flags & GAL_ARITHMETIC_INPLACE)
+ o = in;
+ else
+ o = gal_data_alloc(NULL, in->type, in->ndim, in->dsize, in->wcs,
+ 0, in->minmapsize, NULL, NULL, NULL);
/* Start setting the operator and operands. */
switch(operator)
@@ -593,16 +339,6 @@ arithmetic_unary_function(int operator, unsigned char
flags, gal_data_t *in)
UNIARY_FUNCTION_ON_ELEMENT( log10 );
break;
- case GAL_ARITHMETIC_OP_MINVAL:
- case GAL_ARITHMETIC_OP_MAXVAL:
- case GAL_ARITHMETIC_OP_NUMVAL:
- case GAL_ARITHMETIC_OP_SUMVAL:
- case GAL_ARITHMETIC_OP_MEANVAL:
- case GAL_ARITHMETIC_OP_STDVAL:
- case GAL_ARITHMETIC_OP_MEDIANVAL:
- UNIARY_FUNCTION_ON_ARRAY;
- break;
-
default:
error(EXIT_FAILURE, 0, "operator code %d not recognized in "
"arithmetic_unary_function", operator);
@@ -1565,6 +1301,7 @@ gal_arithmetic_operator_string(int operator)
case GAL_ARITHMETIC_OP_MEANVAL: return "meanvalue";
case GAL_ARITHMETIC_OP_STDVAL: return "stdvalue";
case GAL_ARITHMETIC_OP_MEDIANVAL: return "medianvalue";
+
case GAL_ARITHMETIC_OP_MIN: return "min";
case GAL_ARITHMETIC_OP_MAX: return "max";
case GAL_ARITHMETIC_OP_NUM: return "num";
@@ -1654,6 +1391,38 @@ gal_arithmetic_convert_to_compiled_type(gal_data_t *in,
unsigned char flags)
+/* Call functions in the `gnuastro/statistics' library. */
+static gal_data_t *
+arithmetic_from_statistics(int operator, unsigned char flags,
+ gal_data_t *input)
+{
+ gal_data_t *out;
+ int ip=(flags & GAL_ARITHMETIC_INPLACE) || (flags & GAL_ARITHMETIC_FREE);
+
+ switch(operator)
+ {
+ case GAL_ARITHMETIC_OP_MINVAL: out=gal_statistics_minimum(input); break;
+ case GAL_ARITHMETIC_OP_MAXVAL: out=gal_statistics_maximum(input); break;
+ case GAL_ARITHMETIC_OP_NUMVAL: out=gal_statistics_number(input); break;
+ case GAL_ARITHMETIC_OP_SUMVAL: out=gal_statistics_sum(input); break;
+ case GAL_ARITHMETIC_OP_MEANVAL: out=gal_statistics_mean(input); break;
+ case GAL_ARITHMETIC_OP_STDVAL: out=gal_statistics_std(input); break;
+ case GAL_ARITHMETIC_OP_MEDIANVAL:
+ out=gal_statistics_median(input, ip); break;
+ default:
+ error(EXIT_FAILURE, 0, "operator code %d not recognized in "
+ "`arithmetic_from_statistics'", operator);
+ }
+
+ /* If the input is to be freed, then do so and return the output. */
+ if( flags & GAL_ARITHMETIC_FREE ) gal_data_free(input);
+ return out;
+}
+
+
+
+
+
gal_data_t *
gal_arithmetic(int operator, unsigned char flags, ...)
{
@@ -1709,6 +1478,11 @@ gal_arithmetic(int operator, unsigned char flags, ...)
case GAL_ARITHMETIC_OP_SQRT:
case GAL_ARITHMETIC_OP_LOG:
case GAL_ARITHMETIC_OP_LOG10:
+ d1 = va_arg(va, gal_data_t *);
+ out=arithmetic_unary_function(operator, flags, d1);
+ break;
+
+ /* Statistical operators that return one value. */
case GAL_ARITHMETIC_OP_MINVAL:
case GAL_ARITHMETIC_OP_MAXVAL:
case GAL_ARITHMETIC_OP_NUMVAL:
@@ -1717,9 +1491,10 @@ gal_arithmetic(int operator, unsigned char flags, ...)
case GAL_ARITHMETIC_OP_STDVAL:
case GAL_ARITHMETIC_OP_MEDIANVAL:
d1 = va_arg(va, gal_data_t *);
- out=arithmetic_unary_function(operator, flags, d1);
+ out=arithmetic_from_statistics(operator, flags, d1);
break;
+ /* Absolute operator */
case GAL_ARITHMETIC_OP_ABS:
d1 = va_arg(va, gal_data_t *);
out=arithmetic_abs(flags, d1);
diff --git a/lib/blank.c b/lib/blank.c
index e4a64e6..2faed9b 100644
--- a/lib/blank.c
+++ b/lib/blank.c
@@ -36,58 +36,37 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
-/* Write the blank value of the type into an already allocate space.
-
- Note that for strings, pointer should actually be `char **'. */
+/* Write the blank value of the type into an already allocate space. Note
+ that for STRINGS, pointer should actually be `char **'. */
void
-gal_blank_write(void *pointer, uint8_t type)
+gal_blank_write(void *ptr, uint8_t type)
{
switch(type)
{
+ /* Numeric types */
+ case GAL_DATA_TYPE_UINT8: *(uint8_t *)ptr = GAL_BLANK_UINT8; break;
+ case GAL_DATA_TYPE_INT8: *(int8_t *)ptr = GAL_BLANK_INT8; break;
+ case GAL_DATA_TYPE_UINT16: *(uint16_t *)ptr = GAL_BLANK_UINT16; break;
+ case GAL_DATA_TYPE_INT16: *(int16_t *)ptr = GAL_BLANK_INT16; break;
+ case GAL_DATA_TYPE_UINT32: *(uint32_t *)ptr = GAL_BLANK_UINT32; break;
+ case GAL_DATA_TYPE_INT32: *(int32_t *)ptr = GAL_BLANK_INT32; break;
+ case GAL_DATA_TYPE_UINT64: *(uint64_t *)ptr = GAL_BLANK_UINT64; break;
+ case GAL_DATA_TYPE_INT64: *(int64_t *)ptr = GAL_BLANK_INT64; break;
+ case GAL_DATA_TYPE_FLOAT32: *(float *)ptr = GAL_BLANK_FLOAT32; break;
+ case GAL_DATA_TYPE_FLOAT64: *(double *)ptr = GAL_BLANK_FLOAT64; break;
+
+ /* String type. */
case GAL_DATA_TYPE_STRING:
- gal_checkset_allocate_copy(GAL_BLANK_STRING, pointer);
- break;
-
- case GAL_DATA_TYPE_UINT8:
- *(uint8_t *)pointer = GAL_BLANK_UINT8;
- break;
-
- case GAL_DATA_TYPE_INT8:
- *(int8_t *)pointer = GAL_BLANK_INT8;
- break;
-
- case GAL_DATA_TYPE_UINT16:
- *(uint16_t *)pointer = GAL_BLANK_UINT16;
- break;
-
- case GAL_DATA_TYPE_INT16:
- *(int16_t *)pointer = GAL_BLANK_INT16;
- break;
-
- case GAL_DATA_TYPE_UINT32:
- *(uint32_t *)pointer = GAL_BLANK_UINT32;
- break;
-
- case GAL_DATA_TYPE_INT32:
- *(int32_t *)pointer = GAL_BLANK_INT32;
- break;
-
- case GAL_DATA_TYPE_UINT64:
- *(uint64_t *)pointer = GAL_BLANK_UINT64;
- break;
-
- case GAL_DATA_TYPE_INT64:
- *(int64_t *)pointer = GAL_BLANK_INT64;
+ gal_checkset_allocate_copy(GAL_BLANK_STRING, ptr);
break;
- case GAL_DATA_TYPE_FLOAT32:
- *(float *)pointer = GAL_BLANK_FLOAT32;
- break;
-
- case GAL_DATA_TYPE_FLOAT64:
- *(double *)pointer = GAL_BLANK_FLOAT64;
- break;
+ /* Complex types */
+ case GAL_DATA_TYPE_COMPLEX32:
+ case GAL_DATA_TYPE_COMPLEX64:
+ error(EXIT_FAILURE, 0, "complex types are not yet supported in "
+ "`gal_blank_write'");
+ /* Unrecognized type. */
default:
error(EXIT_FAILURE, 0, "type code %d not recognized in "
"`gal_blank_write'", type);
@@ -119,179 +98,56 @@ gal_blank_alloc_write(uint8_t type)
-/* Print the blank value as a string. */
-char *
-gal_blank_as_string(uint8_t type, int width)
-{
- char *blank;
- switch(type)
- {
- case GAL_DATA_TYPE_BIT:
- error(EXIT_FAILURE, 0, "bit types are not implemented in "
- "`gal_data_blank_as_string' yet.");
- break;
-
- case GAL_DATA_TYPE_STRING:
- if(width)
- asprintf(&blank, "%*s", width, GAL_BLANK_STRING);
- else
- asprintf(&blank, "%s", GAL_BLANK_STRING);
- break;
-
- case GAL_DATA_TYPE_UINT8:
- if(width)
- asprintf(&blank, "%*u", width, (uint8_t)GAL_BLANK_UINT8);
- else
- asprintf(&blank, "%u", (uint8_t)GAL_BLANK_UINT8);
- break;
-
- case GAL_DATA_TYPE_INT8:
- if(width)
- asprintf(&blank, "%*d", width, (int8_t)GAL_BLANK_INT8);
- else
- asprintf(&blank, "%d", (int8_t)GAL_BLANK_INT8);
- break;
-
- case GAL_DATA_TYPE_UINT16:
- if(width)
- asprintf(&blank, "%*u", width, (uint16_t)GAL_BLANK_UINT16);
- else
- asprintf(&blank, "%u", (uint16_t)GAL_BLANK_UINT16);
- break;
-
- case GAL_DATA_TYPE_INT16:
- if(width)
- asprintf(&blank, "%*d", width, (int16_t)GAL_BLANK_INT16);
- else
- asprintf(&blank, "%d", (int16_t)GAL_BLANK_INT16);
- break;
-
- case GAL_DATA_TYPE_UINT32:
- if(width)
- asprintf(&blank, "%*u", width, (uint32_t)GAL_BLANK_UINT32);
- else
- asprintf(&blank, "%u", (uint32_t)GAL_BLANK_UINT32);
- break;
-
- case GAL_DATA_TYPE_INT32:
- if(width)
- asprintf(&blank, "%*d", width, (int32_t)GAL_BLANK_INT32);
- else
- asprintf(&blank, "%d", (int32_t)GAL_BLANK_INT32);
- break;
-
- case GAL_DATA_TYPE_UINT64:
- if(width)
- asprintf(&blank, "%*lu", width, (uint64_t)GAL_BLANK_UINT64);
- else
- asprintf(&blank, "%lu", (uint64_t)GAL_BLANK_UINT64);
- break;
-
- case GAL_DATA_TYPE_INT64:
- if(width)
- asprintf(&blank, "%*ld", width, (int64_t)GAL_BLANK_INT64);
- else
- asprintf(&blank, "%ld", (int64_t)GAL_BLANK_INT64);
- break;
-
- case GAL_DATA_TYPE_FLOAT32:
- if(width)
- asprintf(&blank, "%*f", width, (float)GAL_BLANK_FLOAT32);
- else
- asprintf(&blank, "%f", (float)GAL_BLANK_FLOAT32);
- break;
-
- case GAL_DATA_TYPE_FLOAT64:
- if(width)
- asprintf(&blank, "%*f", width, (double)GAL_BLANK_FLOAT64);
- else
- asprintf(&blank, "%f", (double)GAL_BLANK_FLOAT64);
- break;
-
- default:
- error(EXIT_FAILURE, 0, "type code %d not recognized in "
- "`gal_blank_as_string'", type);
- }
- return blank;
-}
-
-
-
-
-
-
-/* For integers a simple equality is enough. */
-#define HAS_BLANK_INT(CTYPE, BLANK) { \
- CTYPE *a=data->array, *af=a+data->size; \
- do if(*a++ == BLANK) return 1; while(a<af); \
- }
-
-/* Note that a NaN value is not equal to another NaN value, so we can't use
- the easy check for cases were the blank value is NaN. Also note that
- `isnan' is actually a macro, so it works for both float and double
- types.*/
-#define HAS_BLANK_FLT(CTYPE, BLANK, MULTIP) { \
- CTYPE *a=data->array, *af=a+(MULTIP*data->size); \
- if(isnan(BLANK)) do if(isnan(*a++)) return 1; while(a<af); \
- else do if(*a++==BLANK) return 1; while(a<af); \
- }
-
/* Return 1 if the dataset has a blank value and zero if it doesn't. */
+#define HAS_BLANK(IT) { \
+ IT b, *a=data->array, *af=a+data->size; \
+ gal_blank_write(&b, data->type); \
+ if(b==b) do if(*a==b) return 1; while(++a<af); \
+ else do if(*a!=*a) return 1; while(++a<af); \
+ }
int
gal_blank_present(gal_data_t *data)
{
char **str=data->array, **strf=str+data->size;
+ /* If there is nothing in the array (its size is zero), then return 0 (no
+ blank is present. */
+ if(data->size==0) return 0;
+
/* Go over the pixels and check: */
switch(data->type)
{
- case GAL_DATA_TYPE_BIT:
- error(EXIT_FAILURE, 0, "Currently Gnuastro doesn't support bit "
- "datatype, please get in touch with us to implement it.");
-
- case GAL_DATA_TYPE_UINT8:
- HAS_BLANK_INT(uint8_t, GAL_BLANK_UINT8); break;
-
- case GAL_DATA_TYPE_INT8:
- HAS_BLANK_INT(int8_t, GAL_BLANK_INT8); break;
-
- case GAL_DATA_TYPE_UINT16:
- HAS_BLANK_INT(uint16_t, GAL_BLANK_UINT16); break;
-
- case GAL_DATA_TYPE_INT16:
- HAS_BLANK_INT(int16_t, GAL_BLANK_INT16); break;
-
- case GAL_DATA_TYPE_UINT32:
- HAS_BLANK_INT(uint32_t, GAL_BLANK_UINT32); break;
-
- case GAL_DATA_TYPE_INT32:
- HAS_BLANK_INT(int32_t, GAL_BLANK_INT32); break;
-
- case GAL_DATA_TYPE_UINT64:
- HAS_BLANK_INT(uint64_t, GAL_BLANK_UINT64); break;
-
- case GAL_DATA_TYPE_INT64:
- HAS_BLANK_INT(int64_t, GAL_BLANK_INT64); break;
-
- case GAL_DATA_TYPE_FLOAT32:
- HAS_BLANK_FLT(float, GAL_BLANK_FLOAT32, 1); break;
-
- case GAL_DATA_TYPE_FLOAT64:
- HAS_BLANK_FLT(double, GAL_BLANK_FLOAT64, 1); break;
+ /* Numeric types */
+ case GAL_DATA_TYPE_UINT8: HAS_BLANK( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: HAS_BLANK( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: HAS_BLANK( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: HAS_BLANK( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: HAS_BLANK( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: HAS_BLANK( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: HAS_BLANK( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: HAS_BLANK( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: HAS_BLANK( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: HAS_BLANK( double ); break;
+
+ /* String. */
+ case GAL_DATA_TYPE_STRING:
+ do if(!strcmp(*str++,GAL_BLANK_STRING)) return 1; while(str<strf);
+ break;
+ /* Complex types */
case GAL_DATA_TYPE_COMPLEX32:
- HAS_BLANK_FLT(float, GAL_BLANK_FLOAT32, 2); break;
-
case GAL_DATA_TYPE_COMPLEX64:
- HAS_BLANK_FLT(double, GAL_BLANK_FLOAT64, 2); break;
+ error(EXIT_FAILURE, 0, "complex types are not yet supported in "
+ "`gal_blank_write'");
- case GAL_DATA_TYPE_STRING:
- do if(!strcmp(*str++,GAL_BLANK_STRING)) return 1; while(str<strf);
- break;
+ /* Bit */
+ case GAL_DATA_TYPE_BIT:
+ error(EXIT_FAILURE, 0, "bit type datasets are not yet supported in "
+ "`gal_blank_present'");
default:
error(EXIT_FAILURE, 0, "a bug! type value (%d) not recognized "
- "in `gal_data_has_blank'", data->type);
+ "in `gal_blank_present'", data->type);
}
/* If there was a blank value, then the function would have returned with
@@ -304,32 +160,19 @@ gal_blank_present(gal_data_t *data)
-/* For integers a simple equality is enough. */
-#define FLAG_BLANK_INT(CTYPE, BLANK) { \
- CTYPE *a=data->array; do *o = (*a==BLANK); while(++o<of); \
- }
-/* Note that a NaN value is not equal to another NaN value, so we can't use
- the easy check for cases were the blank value is NaN. Also note that
- `isnan' is actually a macro, so it works for both float and double
- types.*/
-#define FLAG_BLANK_FLT(CTYPE, BLANK) { \
- CTYPE *a=data->array; \
- if(isnan(BLANK)) do *o = isnan(*a++); while(++o<of); \
- else do *o = (*a++==BLANK); while(++o<of); \
- }
-
-#define FLAG_BLANK_COMPLEX(CTYPE, BLANK) { \
- CTYPE *a=data->array; \
- if(isnan(BLANK)) \
- do { *o=(isnan(*a) || isnan(*(a+1))); a+=2; } while(++o<of); \
- else \
- do { *o=(*a==BLANK || *(a+1)==BLANK); a+=2; } while(++o<of); \
- }
-/* Output a data-set of the the same size as the input, but with an uint8_t
+/* Create a dataset of the the same size as the input, but with an uint8_t
type that has a value of 1 for data that are blank and 0 for those that
aren't. */
+#define FLAG_BLANK(IT) { \
+ IT b, *a=data->array; \
+ gal_blank_write(&b, data->type); \
+ if(b==b) /* Blank value can be checked with the equal comparison */ \
+ do { *o = *a==b; ++a; } while(++o<of); \
+ else /* Blank value will fail with the equal comparison */ \
+ do { *o = *a!=*a; ++a; } while(++o<of); \
+ }
gal_data_t *
gal_blank_flag(gal_data_t *data)
{
@@ -348,50 +191,31 @@ gal_blank_flag(gal_data_t *data)
/* Go over the pixels and set the output values. */
switch(data->type)
{
- case GAL_DATA_TYPE_BIT:
- error(EXIT_FAILURE, 0, "Currently Gnuastro doesn't support bit "
- "datatype, please get in touch with us to implement it.");
-
- case GAL_DATA_TYPE_UINT8:
- FLAG_BLANK_INT(uint8_t, GAL_BLANK_UINT8); break;
-
- case GAL_DATA_TYPE_INT8:
- FLAG_BLANK_INT(int8_t, GAL_BLANK_INT8); break;
-
- case GAL_DATA_TYPE_UINT16:
- FLAG_BLANK_INT(uint16_t, GAL_BLANK_UINT16); break;
-
- case GAL_DATA_TYPE_INT16:
- FLAG_BLANK_INT(int16_t, GAL_BLANK_INT16); break;
-
- case GAL_DATA_TYPE_UINT32:
- FLAG_BLANK_INT(uint32_t, GAL_BLANK_UINT32); break;
-
- case GAL_DATA_TYPE_INT32:
- FLAG_BLANK_INT(int32_t, GAL_BLANK_INT32); break;
-
- case GAL_DATA_TYPE_UINT64:
- FLAG_BLANK_INT(uint64_t, GAL_BLANK_UINT64); break;
-
- case GAL_DATA_TYPE_INT64:
- FLAG_BLANK_INT(int64_t, GAL_BLANK_INT64); break;
-
- case GAL_DATA_TYPE_FLOAT32:
- FLAG_BLANK_FLT(float, GAL_BLANK_FLOAT32); break;
-
- case GAL_DATA_TYPE_FLOAT64:
- FLAG_BLANK_FLT(double, GAL_BLANK_FLOAT64); break;
-
- case GAL_DATA_TYPE_COMPLEX32:
- FLAG_BLANK_COMPLEX(float, GAL_BLANK_FLOAT32); break;
-
- case GAL_DATA_TYPE_COMPLEX64:
- FLAG_BLANK_COMPLEX(double, GAL_BLANK_FLOAT64); break;
-
+ /* Numeric types */
+ case GAL_DATA_TYPE_UINT8: FLAG_BLANK( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: FLAG_BLANK( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: FLAG_BLANK( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: FLAG_BLANK( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: FLAG_BLANK( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: FLAG_BLANK( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: FLAG_BLANK( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: FLAG_BLANK( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: FLAG_BLANK( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: FLAG_BLANK( double ); break;
+
+ /* String. */
case GAL_DATA_TYPE_STRING:
do *o++ = !strcmp(*str,GAL_BLANK_STRING); while(++str<strf);
break;
+ /* Currently unsupported types. */
+ case GAL_DATA_TYPE_BIT:
+ case GAL_DATA_TYPE_COMPLEX32:
+ case GAL_DATA_TYPE_COMPLEX64:
+ error(EXIT_FAILURE, 0, "%s type not yet supported in `gal_blank_flag'",
+ gal_data_type_as_string(data->type, 1));
+
+ /* Bad input. */
default:
error(EXIT_FAILURE, 0, "type value (%d) not recognized "
"in `gal_blank_flag'", data->type);
@@ -405,24 +229,22 @@ gal_blank_flag(gal_data_t *data)
+/* Remove blank elements from a dataset, convert it to a 1D dataset and
+ adjust the size properly. In practice this function doesn't `realloc'
+ the input array, all it does is to shift the blank eleemnts to the end
+ and adjust the size elements of the `gal_data_t'. */
#define BLANK_REMOVE(IT) { \
IT b, *a=data->array, *af=a+data->size, *o=data->array; \
if( gal_blank_present(data) ) \
{ \
gal_blank_write(&b, data->type); \
- if(b==b) /* Can be checked with equal: isn't NaN */ \
+ if(b==b) /* Blank value can be be checked with equal. */ \
do if(*a!=b) { ++num; *o++=*a; } while(++a<af); \
- else /* Blank value is NaN, so equals will fail on blank elements*/ \
+ else /* Blank value will fail on equal comparison. */ \
do if(*a==*a) { ++num; *o++=*a; } while(++a<af); \
} \
else num=data->size; \
}
-
-
-/* Remove blank elements from a dataset, convert it to a 1D dataset and
- adjust the size properly. In practice this function doesn't `realloc'
- the input array, all it does is to shift the blank eleemnts to the end
- and adjust the size elements of the `gal_data_t'. */
void
gal_blank_remove(gal_data_t *data)
{
@@ -450,3 +272,103 @@ gal_blank_remove(gal_data_t *data)
data->ndim=1;
data->dsize[0]=data->size=num;
}
+
+
+
+
+
+/* Print the blank value as a string. */
+char *
+gal_blank_as_string(uint8_t type, int width)
+{
+ char *blank;
+ switch(type)
+ {
+ case GAL_DATA_TYPE_BIT:
+ error(EXIT_FAILURE, 0, "bit types are not implemented in "
+ "`gal_data_blank_as_string' yet.");
+ break;
+
+ case GAL_DATA_TYPE_STRING:
+ if(width)
+ asprintf(&blank, "%*s", width, GAL_BLANK_STRING);
+ else
+ asprintf(&blank, "%s", GAL_BLANK_STRING);
+ break;
+
+ case GAL_DATA_TYPE_UINT8:
+ if(width)
+ asprintf(&blank, "%*u", width, (uint8_t)GAL_BLANK_UINT8);
+ else
+ asprintf(&blank, "%u", (uint8_t)GAL_BLANK_UINT8);
+ break;
+
+ case GAL_DATA_TYPE_INT8:
+ if(width)
+ asprintf(&blank, "%*d", width, (int8_t)GAL_BLANK_INT8);
+ else
+ asprintf(&blank, "%d", (int8_t)GAL_BLANK_INT8);
+ break;
+
+ case GAL_DATA_TYPE_UINT16:
+ if(width)
+ asprintf(&blank, "%*u", width, (uint16_t)GAL_BLANK_UINT16);
+ else
+ asprintf(&blank, "%u", (uint16_t)GAL_BLANK_UINT16);
+ break;
+
+ case GAL_DATA_TYPE_INT16:
+ if(width)
+ asprintf(&blank, "%*d", width, (int16_t)GAL_BLANK_INT16);
+ else
+ asprintf(&blank, "%d", (int16_t)GAL_BLANK_INT16);
+ break;
+
+ case GAL_DATA_TYPE_UINT32:
+ if(width)
+ asprintf(&blank, "%*u", width, (uint32_t)GAL_BLANK_UINT32);
+ else
+ asprintf(&blank, "%u", (uint32_t)GAL_BLANK_UINT32);
+ break;
+
+ case GAL_DATA_TYPE_INT32:
+ if(width)
+ asprintf(&blank, "%*d", width, (int32_t)GAL_BLANK_INT32);
+ else
+ asprintf(&blank, "%d", (int32_t)GAL_BLANK_INT32);
+ break;
+
+ case GAL_DATA_TYPE_UINT64:
+ if(width)
+ asprintf(&blank, "%*lu", width, (uint64_t)GAL_BLANK_UINT64);
+ else
+ asprintf(&blank, "%lu", (uint64_t)GAL_BLANK_UINT64);
+ break;
+
+ case GAL_DATA_TYPE_INT64:
+ if(width)
+ asprintf(&blank, "%*ld", width, (int64_t)GAL_BLANK_INT64);
+ else
+ asprintf(&blank, "%ld", (int64_t)GAL_BLANK_INT64);
+ break;
+
+ case GAL_DATA_TYPE_FLOAT32:
+ if(width)
+ asprintf(&blank, "%*f", width, (float)GAL_BLANK_FLOAT32);
+ else
+ asprintf(&blank, "%f", (float)GAL_BLANK_FLOAT32);
+ break;
+
+ case GAL_DATA_TYPE_FLOAT64:
+ if(width)
+ asprintf(&blank, "%*f", width, (double)GAL_BLANK_FLOAT64);
+ else
+ asprintf(&blank, "%f", (double)GAL_BLANK_FLOAT64);
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`gal_blank_as_string'", type);
+ }
+ return blank;
+}
diff --git a/lib/commonopts.h b/lib/commonopts.h
index 0c67391..521748e 100644
--- a/lib/commonopts.h
+++ b/lib/commonopts.h
@@ -123,7 +123,7 @@ struct argp_option gal_commonopts_options[] =
GAL_OPTIONS_KEY_TABLEFORMAT,
"STR",
0,
- "Output table format: `fits-ascii', `fits-binary'.",
+ "Table format: `fits-ascii', `fits-binary'.",
GAL_OPTIONS_GROUP_OUTPUT,
&cp->tableformat,
GAL_DATA_TYPE_STRING,
diff --git a/lib/data.c b/lib/data.c
index 7fd3c9f..322965e 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -1379,6 +1379,12 @@ gal_data_copy_to_new_type(gal_data_t *in, uint8_t
newtype)
out=gal_data_alloc(NULL, newtype, in->ndim, in->dsize, in->wcs,
0, in->minmapsize, in->name, in->unit, in->comment);
+ /* Set the rest of the values. */
+ out->next=in->next;
+ out->status=in->status;
+ out->disp_width=in->disp_width;
+ out->disp_precision=in->disp_precision;
+
/* For debugging.
printf("in: %d (%s)\nout: %d (%s)\n\n", in->type,
gal_data_type_as_string(in->type, 1), out->type,
diff --git a/lib/fits.c b/lib/fits.c
index 4423254..30b569e 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -134,42 +134,17 @@ gal_fits_suffix_is_fits(char *suffix)
-/* We have the name of the input file. But in most cases, the files
- that should be used (for example a mask image) are other extensions
- in the same file. So the user only has to give the HDU. The job of
- this function is to determine which is the case and set othername
- to the appropriate value. */
-void
-gal_fits_file_or_ext_name(char *inputname, char *inhdu, int othernameset,
- char **othername, char *ohdu, int ohduset,
- char *type)
+/* If the name is a FITS name, then put a `(hdu: ...)' after it and return
+ the string. If it isn't a FITS file, just print the name. Note that the
+ space is allocated. */
+char *
+gal_fits_name_save_as_string(char *filename, char *hdu)
{
- if(othernameset)
- {
- /* In some cases, for example a mask image, both the name and
- HDU are optional. So just to be safe, we will check this all
- the time. */
- if(ohduset==0)
- error(EXIT_FAILURE, 0, "a %s image was specified (%s). However, "
- "no HDU is given for it. Please add a HDU. If you regularly "
- "use the same HDU as %s, you may consider adding it to "
- "the configuration file. For more information, please see "
- "the `Configuration files' section of the %s manual by "
- "running ` info gnuastro ' on the command-line", type,
- *othername, type, PACKAGE_NAME);
- if(strcmp(inputname, *othername)==0)
- {
- if(strcmp(ohdu, inhdu)==0)
- error(EXIT_FAILURE, 0, "the specified %s name and "
- "input image name (%s) are the same while the input "
- "image hdu name and mask hdu are also identical (%s)",
- type, inputname, inhdu);
- }
- }
- else if(ohduset && strcmp(ohdu, inhdu))
- *othername=inputname;
- else
- *othername=NULL;
+ char *name;
+ if( gal_fits_name_is_fits(filename) )
+ asprintf(&name, "%s (hdu: %s)", filename, hdu);
+ else gal_checkset_allocate_copy(filename, &name);
+ return name;
}
@@ -2208,8 +2183,8 @@ fits_write_tnull_tcomm(fitsfile *fptr, gal_data_t *col,
int tabletype,
/* Write the given columns (a linked list of `gal_data_t') into a FITS
table.*/
void
-gal_fits_tab_write(gal_data_t *cols, char *comments, int tabletype,
- char *filename, int dontdelete)
+gal_fits_tab_write(gal_data_t *cols, struct gal_linkedlist_stll *comments,
+ int tabletype, char *filename, int dontdelete)
{
void *blank;
fitsfile *fptr;
@@ -2217,6 +2192,7 @@ gal_fits_tab_write(gal_data_t *cols, char *comments, int
tabletype,
size_t i, numrows=-1;
char **ttype, **tform, **tunit;
int tbltype, numcols=0, status=0;
+ struct gal_linkedlist_stll *strt;
/* Make sure all the input columns have the same number of elements */
@@ -2280,6 +2256,11 @@ gal_fits_tab_write(gal_data_t *cols, char *comments, int
tabletype,
}
+ /* Write the comments if there were any. */
+ for(strt=comments; strt!=NULL; strt=strt->next)
+ fits_write_comment(fptr, strt->v, &status);
+
+
/* Write all the headers and the version information. */
gal_fits_key_write_version(fptr, NULL, NULL);
diff --git a/lib/gnuastro/fits.h b/lib/gnuastro/fits.h
index f491876..c9b5437 100644
--- a/lib/gnuastro/fits.h
+++ b/lib/gnuastro/fits.h
@@ -115,6 +115,9 @@ gal_fits_file_or_ext_name(char *inputname, char *inhdu, int
othernameset,
char **othername, char *ohdu, int ohduset,
char *type);
+char *
+gal_fits_name_save_as_string(char *filename, char *hdu);
+
/*************************************************************
@@ -245,8 +248,8 @@ gal_fits_tab_read(char *filename, char *hdu, size_t numrows,
int minmapsize);
void
-gal_fits_tab_write(gal_data_t *cols, char *comments, int tabletype,
- char *filename, int dontdelete);
+gal_fits_tab_write(gal_data_t *cols, struct gal_linkedlist_stll *comments,
+ int tabletype, char *filename, int dontdelete);
diff --git a/lib/gnuastro/statistics.h b/lib/gnuastro/statistics.h
index 34857f4..f9fce77 100644
--- a/lib/gnuastro/statistics.h
+++ b/lib/gnuastro/statistics.h
@@ -75,28 +75,37 @@ enum bin_status
/****************************************************************
******** Simple statistics *******
- **** (wrappers for functions in `arithmetic.h') ****
****************************************************************/
+
+gal_data_t *
+gal_statistics_number(gal_data_t *input);
+
gal_data_t *
-gal_statistics_number(gal_data_t *data);
+gal_statistics_minimum(gal_data_t *input);
gal_data_t *
-gal_statistics_minimum(gal_data_t *data);
+gal_statistics_maximum(gal_data_t *input);
gal_data_t *
-gal_statistics_maximum(gal_data_t *data);
+gal_statistics_sum(gal_data_t *input);
gal_data_t *
-gal_statistics_sum(gal_data_t *data);
+gal_statistics_mean(gal_data_t *input);
gal_data_t *
-gal_statistics_mean(gal_data_t *data);
+gal_statistics_std(gal_data_t *input);
gal_data_t *
-gal_statistics_std(gal_data_t *data);
+gal_statistics_mean_std(gal_data_t *input);
gal_data_t *
-gal_statistics_median(gal_data_t *data);
+gal_statistics_median(gal_data_t *input, int inplace);
+
+gal_data_t *
+gal_statistsics_quantile(gal_data_t *input, float quantile, int inplace);
+
+size_t
+gal_statistics_quantile_index(size_t size, float quant);
@@ -115,7 +124,8 @@ gal_statistics_sort_increasing(gal_data_t *data);
void
gal_statistics_sort_decreasing(gal_data_t *data);
-
+gal_data_t *
+gal_statistics_no_blank_sorted(gal_data_t *input, int inplace);
@@ -137,32 +147,12 @@ gal_statistics_cfp(gal_data_t *data, gal_data_t *bins,
int normalize);
-
-/****************************************************************
- ***************** Quantiles ********************
- ****************************************************************/
-size_t
-gal_statistics_index_from_quantile(size_t size, float quant);
-
-
-
-
-
/****************************************************************
***************** Sigma clip ********************
****************************************************************/
-int
-gal_statistics_sigma_clip_converge(float *array, int o1_n0, size_t num_elem,
- float sigma_multiple, float accuracy,
- float *outave, float *outmed,
- float *outstd, int print);
-
-int
-gal_statistics_sigma_clip_certain_num(float *array, int o1_n0,
- size_t num_elem, float sigma_multiple,
- size_t numtimes, float *outave,
- float *outmed, float *outstd,
- int print);
+gal_data_t *
+gal_statistics_sigma_clip(gal_data_t *input, float multip, float param,
+ int quiet);
diff --git a/lib/gnuastro/table.h b/lib/gnuastro/table.h
index 4471144..60041b5 100644
--- a/lib/gnuastro/table.h
+++ b/lib/gnuastro/table.h
@@ -188,8 +188,19 @@ gal_table_read(char *filename, char *hdu, struct
gal_linkedlist_stll *cols,
/*************** Write a table ***************/
/************************************************************************/
void
-gal_table_write(gal_data_t *cols, char *comments, int tableformat,
- char *filename, int dontdelete);
+gal_table_comments_add_intro(struct gal_linkedlist_stll **comments,
+ char *program_string, time_t *rawtime);
+
+void
+gal_table_write(gal_data_t *cols, struct gal_linkedlist_stll *comments,
+ int tableformat, char *filename, int dontdelete);
+
+void
+gal_table_write_log(gal_data_t *logll, char *program_string,
+ time_t *rawtime, struct gal_linkedlist_stll *comments,
+ char *filename, int dontdelete, int quiet);
+
+
diff --git a/lib/gnuastro/txt.h b/lib/gnuastro/txt.h
index 86a54a5..447e1b5 100644
--- a/lib/gnuastro/txt.h
+++ b/lib/gnuastro/txt.h
@@ -96,8 +96,8 @@ gal_data_t *
gal_txt_image_read(char *filename, size_t minmapsize);
void
-gal_txt_write(gal_data_t *cols, char *comment, char *filename,
- int dontdelete);
+gal_txt_write(gal_data_t *cols, struct gal_linkedlist_stll *comment,
+ char *filename, int dontdelete);
diff --git a/lib/mode.c b/lib/mode.c
index 7cc5ec1..ac37f4a 100644
--- a/lib/mode.c
+++ b/lib/mode.c
@@ -354,7 +354,7 @@ modesymmetricity(float *a, size_t size, size_t mi, float
errorstdm,
mf=a[mi];
errdiff=errorstdm*sqrt(mi);
topi = 2*mi>size-1 ? size-1 : 2*mi;
- af=a[gal_statistics_index_from_quantile(2*mi+1,
+ af=a[gal_statistics_quantile_index(2*mi+1,
GAL_STATISTICS_MODE_SYMMETRICITY_LOW_QUANT)];
/* This loop is very similar to that of mirrormaxdiff(). It will
diff --git a/lib/options.c b/lib/options.c
index 3f0885c..b66681f 100644
--- a/lib/options.c
+++ b/lib/options.c
@@ -304,8 +304,15 @@ void *
gal_options_read_type(struct argp_option *option, char *arg,
char *filename, size_t lineno, void *junk)
{
+ char *str;
if(lineno==-1)
- return gal_data_type_as_string( *(uint8_t *)(option->value), 1);
+ {
+ /* Note that `gal_data_type_as_string' returns a static string. But
+ the output must be an allocated string so we can free it. */
+ gal_checkset_allocate_copy(
+ gal_data_type_as_string( *(uint8_t *)(option->value), 1), &str);
+ return str;
+ }
else
{
/* If the option is already set, just return. */
@@ -336,8 +343,15 @@ void *
gal_options_read_searchin(struct argp_option *option, char *arg,
char *filename, size_t lineno, void *junk)
{
+ char *str;
if(lineno==-1)
- return gal_table_searchin_as_string( *(uint8_t *)(option->value));
+ {
+ /* Note that `gal_data_type_as_string' returns a static string. But
+ the output must be an allocated string so we can free it. */
+ gal_checkset_allocate_copy(
+ gal_table_searchin_as_string( *(uint8_t *)(option->value)), &str);
+ return str;
+ }
else
{
/* If the option is already set, just return. */
@@ -369,8 +383,15 @@ void *
gal_options_read_tableformat(struct argp_option *option, char *arg,
char *filename, size_t lineno, void *junk)
{
+ char *str;
if(lineno==-1)
- return gal_table_format_as_string( *(uint8_t *)(option->value));
+ {
+ /* Note that `gal_data_type_as_string' returns a static string. But
+ the output must be an allocated string so we can free it. */
+ gal_checkset_allocate_copy(
+ gal_table_format_as_string( *(uint8_t *)(option->value)), &str);
+ return str;
+ }
else
{
/* If the option is already set, then you don't have to do anything. */
@@ -1288,13 +1309,14 @@ option_is_printable(struct argp_option *option)
value. */
static int
options_print_any_type(struct argp_option *option, void *ptr, int type,
- int width, FILE *fp)
+ int width, FILE *fp,
+ struct gal_options_common_params *cp)
{
char *str;
/* Write the value into a string. */
str = ( option->func
- ? option->func(option, NULL, NULL, (size_t)(-1), NULL)
+ ? option->func(option, NULL, NULL, (size_t)(-1), cp->program_struct)
: gal_data_write_to_string(ptr, type, 1) );
/* If only the width was desired, don't actually print the string, just
@@ -1304,9 +1326,8 @@ options_print_any_type(struct argp_option *option, void
*ptr, int type,
else
width=strlen(str);
- /* Free the allocated space and return. When the value was taken from a
- function, it is static, so it must not be freed. */
- if(!option->func) free(str);
+ /* Free the allocated space and return. */
+ free(str);
return width;
}
@@ -1318,7 +1339,8 @@ options_print_any_type(struct argp_option *option, void
*ptr, int type,
lengths. */
static void
options_correct_max_lengths(struct argp_option *option, int *max_nlen,
- int *max_vlen)
+ int *max_vlen,
+ struct gal_options_common_params *cp)
{
int vlen;
struct gal_linkedlist_stll *tmp;
@@ -1342,7 +1364,7 @@ options_correct_max_lengths(struct argp_option *option,
int *max_nlen,
{
/* Get the length of this node: */
vlen=options_print_any_type(option, &tmp->v, GAL_DATA_TYPE_STRING,
- 0, NULL);
+ 0, NULL, cp);
/* If its larger than the maximum length, then put it in. */
if( vlen > *max_vlen )
@@ -1352,7 +1374,7 @@ options_correct_max_lengths(struct argp_option *option,
int *max_nlen,
else
{
vlen=options_print_any_type(option, option->value, option->type,
- 0, NULL);
+ 0, NULL, cp);
if( vlen > *max_vlen )
*max_vlen=vlen;
}
@@ -1371,14 +1393,16 @@ options_correct_max_lengths(struct argp_option *option,
int *max_nlen,
and their values. */
static void
options_set_lengths(struct argp_option *poptions,
- struct argp_option *coptions, int *namelen, int *valuelen)
+ struct argp_option *coptions,
+ int *namelen, int *valuelen,
+ struct gal_options_common_params *cp)
{
int i, max_nlen=0, max_vlen=0;
/* For program specific options. */
for(i=0; !gal_options_is_last(&poptions[i]); ++i)
if(poptions[i].name && poptions[i].set)
- options_correct_max_lengths(&poptions[i], &max_nlen, &max_vlen);
+ options_correct_max_lengths(&poptions[i], &max_nlen, &max_vlen, cp);
/* For common options. Note that the options that will not be printed are
in this category, so we also need to check them. The detailed steps
@@ -1386,7 +1410,7 @@ options_set_lengths(struct argp_option *poptions,
for(i=0; !gal_options_is_last(&coptions[i]); ++i)
if( coptions[i].name && coptions[i].set
&& option_is_printable(&coptions[i]) )
- options_correct_max_lengths(&coptions[i], &max_nlen, &max_vlen);
+ options_correct_max_lengths(&coptions[i], &max_nlen, &max_vlen, cp);
/* Save the final values in the output pointers. */
*namelen=max_nlen;
@@ -1409,7 +1433,7 @@ options_print_doc(FILE *fp, const char *doc, int nvwidth)
/* The `+3' is because of the three extra spaces in this line: one before
the variable name, one after it and one after the value. */
- int i, prewidth=nvwidth+3, width=78-prewidth, cwidth;
+ int i, prewidth=nvwidth+3, width=77-prewidth, cwidth;
/* We only want the formatting when writing to stdout. */
if(len<width)
@@ -1443,7 +1467,8 @@ options_print_doc(FILE *fp, const char *doc, int nvwidth)
static void
options_print_all_in_group(struct argp_option *options, int groupint,
- int namelen, int valuelen, FILE *fp)
+ int namelen, int valuelen, FILE *fp,
+ struct gal_options_common_params *cp)
{
size_t i;
struct gal_linkedlist_stll *tmp;
@@ -1462,7 +1487,8 @@ options_print_all_in_group(struct argp_option *options,
int groupint,
{
fprintf(fp, " %-*s ", namewidth, options[i].name);
options_print_any_type(&options[i], &tmp->v,
- GAL_DATA_TYPE_STRING, valuewidth, fp);
+ GAL_DATA_TYPE_STRING, valuewidth,
+ fp, cp);
options_print_doc(fp, options[i].doc, namewidth+valuewidth);
}
@@ -1471,7 +1497,7 @@ options_print_all_in_group(struct argp_option *options,
int groupint,
{
fprintf(fp, " %-*s ", namewidth, options[i].name);
options_print_any_type(&options[i], options[i].value,
- options[i].type, valuewidth, fp);
+ options[i].type, valuewidth, fp, cp);
options_print_doc(fp, options[i].doc, namewidth+valuewidth);
}
}
@@ -1560,7 +1586,7 @@ options_print_all(struct gal_options_common_params *cp,
char *dirname,
gal_linkedlist_reverse_ill(&group);
/* Get the maximum width of names and values. */
- options_set_lengths(poptions, coptions, &namelen, &valuelen);
+ options_set_lengths(poptions, coptions, &namelen, &valuelen, cp);
/* Go over each topic and print every option that is in this group. */
while(topic)
@@ -1571,13 +1597,16 @@ options_print_all(struct gal_options_common_params *cp,
char *dirname,
/* First print the topic, */
fprintf(fp, "\n# %s\n", topicstr);
+ /*
fprintf(fp, "# ");
i=0; while(i++<strlen(topicstr)) fprintf(fp, "%c", '-');
fprintf(fp, "\n");
-
+ */
/* Then, print all the options that are in this group. */
- options_print_all_in_group(coptions, groupint, namelen, valuelen, fp);
- options_print_all_in_group(poptions, groupint, namelen, valuelen, fp);
+ options_print_all_in_group(coptions, groupint, namelen, valuelen,
+ fp, cp);
+ options_print_all_in_group(poptions, groupint, namelen, valuelen,
+ fp, cp);
}
/* Let the user know. */
@@ -1649,41 +1678,3 @@ gal_options_print_state(struct gal_options_common_params
*cp)
break;
}
}
-
-
-
-
-
-void
-gal_options_print_log(gal_data_t *logll, char *program_string,
- time_t *rawtime, char *comments, char *filename,
- struct gal_options_common_params *cp)
-{
- char *finalcomment;
- char *msg, gitdescribe[100], *gd;
-
- /* Get the Git description in the running folder. */
- gd=gal_git_describe();
- if(gd) sprintf(gitdescribe, " from %s,", gd);
- else gitdescribe[0]='\0';
-
- /* Write the top level comment. */
- asprintf(&finalcomment, "# %s\n# Created%s on %s%s",
- program_string, gitdescribe, ctime(rawtime),
- comments ? comments : "#");
-
- /* Write the log file to disk */
- gal_table_write(logll, finalcomment, GAL_TABLE_FORMAT_TXT, filename,
- cp->dontdelete);
-
- /* In verbose mode, print the information. */
- if(!cp->quiet)
- {
- asprintf(&msg, "%s created.", filename);
- gal_timing_report(NULL, msg, 1);
- free(msg);
- }
-
- /* Clean up. */
- free(finalcomment);
-}
diff --git a/lib/options.h b/lib/options.h
index ebbe7bc..21b1871 100644
--- a/lib/options.h
+++ b/lib/options.h
@@ -261,9 +261,4 @@ gal_options_read_config_set(struct
gal_options_common_params *cp);
void
gal_options_print_state(struct gal_options_common_params *cp);
-void
-gal_options_print_log(gal_data_t *logll, char *program_string,
- time_t *rawtime, char *comments, char *filename,
- struct gal_options_common_params *cp);
-
#endif
diff --git a/lib/statistics.c b/lib/statistics.c
index 1c51fe7..7604418 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <stdlib.h>
#include <gnuastro/data.h>
+#include <gnuastro/blank.h>
#include <gnuastro/qsort.h>
#include <gnuastro/arithmetic.h>
#include <gnuastro/statistics.h>
@@ -50,50 +51,431 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/****************************************************************
******** Simple statistics *******
- **** (wrappers for functions in `arithmetic.h') ****
****************************************************************/
-#define SIMP_FLAGS GAL_ARITHMETIC_NUMOK
+/* Return the number of non-blank elements in an array as a single element,
+ uint64 type data structure. */
+#define STATS_NUM(IT) { \
+ IT b, *a=input->array, *af=a+input->size; \
+ gal_blank_write(&b, input->type); \
+ if( gal_blank_present(input) ) \
+ { \
+ if(b==b) /* Blank can be found with the equal sign. */ \
+ do if(*a!=b) ++(*num); while(++a<af); \
+ else /* Blank is NaN (fails on equal). */ \
+ do if(*a==*a) ++(*num); while(++a<af); \
+ } \
+ else *num=input->size; \
+ }
+gal_data_t *
+gal_statistics_number(gal_data_t *input)
+{
+ uint64_t *num;
+ size_t dsize=1;
+ gal_data_t *out=gal_data_alloc(NULL, GAL_DATA_TYPE_UINT64, 1, &dsize,
+ NULL, 1, -1, NULL, NULL, NULL);
+ num=out->array;
+ switch(input->type)
+ {
+ case GAL_DATA_TYPE_UINT8: STATS_NUM( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: STATS_NUM( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: STATS_NUM( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: STATS_NUM( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: STATS_NUM( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: STATS_NUM( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: STATS_NUM( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: STATS_NUM( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: STATS_NUM( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: STATS_NUM( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`gal_statistics_number'", input->type);
+ }
+ return out;
+}
+
+
+
+
+/* Return the minimum (non-blank) value of a dataset in the same type as
+ the dataset. Note that a NaN (blank in floating point) will fail on any
+ comparison. So when finding the minimum or maximum, when the blank value
+ is NaN, we can safely assume there is no blank value at all. */
+#define STATS_MIN(IT) { \
+ IT b, *o=out->array, *a=input->array, *af=a+input->size; \
+ gal_blank_write(&b, input->type); \
+ if( b==b && gal_blank_present(input) ) \
+ do if(*a!=b) { *o = *a < *o ? *a : *o; ++n; } while(++a<af); \
+ else \
+ do { *o = *a < *o ? *a : *o; ++n; } while(++a<af); \
+ \
+ /* If there were no usable elements, set the output to blank. */ \
+ if(n==0) *o=b; \
+ }
gal_data_t *
-gal_statistics_number(gal_data_t *data)
+gal_statistics_minimum(gal_data_t *input)
{
- return gal_arithmetic(GAL_ARITHMETIC_OP_NUMVAL, SIMP_FLAGS, data);
+ size_t dsize=1, n=0;
+ gal_data_t *out=gal_data_alloc(NULL, input->type, 1, &dsize,
+ NULL, 1, -1, NULL, NULL, NULL);
+ gal_data_type_max(out->type, out->array);
+ switch(input->type)
+ {
+ case GAL_DATA_TYPE_UINT8: STATS_MIN( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: STATS_MIN( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: STATS_MIN( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: STATS_MIN( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: STATS_MIN( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: STATS_MIN( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: STATS_MIN( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: STATS_MIN( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: STATS_MIN( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: STATS_MIN( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`gal_statistics_minimum'", input->type);
+ }
+ return out;
}
+
+
+
+
+/* Return the maximum (non-blank) value of a dataset in the same type as
+ the dataset. See explanations of `gal_statistics_minimum'. */
+#define STATS_MAX(IT) { \
+ IT b, *o=out->array, *a=input->array, *af=a+input->size; \
+ gal_blank_write(&b, input->type); \
+ if( b==b && gal_blank_present(input) ) \
+ do if(*a!=b) { *o = *a > *o ? *a : *o; ++n; } while(++a<af); \
+ else \
+ do { *o = *a > *o ? *a : *o; ++n; } while(++a<af); \
+ \
+ /* If there were no usable elements, set the output to blank. */ \
+ if(n==0) *o=b; \
+ }
gal_data_t *
-gal_statistics_minimum(gal_data_t *data)
+gal_statistics_maximum(gal_data_t *input)
{
- return gal_arithmetic(GAL_ARITHMETIC_OP_MINVAL, SIMP_FLAGS, data);
+ size_t dsize=1, n=0;
+ gal_data_t *out=gal_data_alloc(NULL, input->type, 1, &dsize,
+ NULL, 1, -1, NULL, NULL, NULL);
+ gal_data_type_min(out->type, out->array);
+ switch(input->type)
+ {
+ case GAL_DATA_TYPE_UINT8: STATS_MAX( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: STATS_MAX( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: STATS_MAX( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: STATS_MAX( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: STATS_MAX( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: STATS_MAX( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: STATS_MAX( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: STATS_MAX( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: STATS_MAX( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: STATS_MAX( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`gal_statistics_maximum'", input->type);
+ }
+ return out;
}
+
+
+
+
+/* Return the sum of the input dataset as a single element dataset of type
+ float64. */
+#define STATS_SUM_NUM(IT) { \
+ IT b, *a=input->array, *af=a+input->size; \
+ gal_blank_write(&b, input->type); \
+ if( gal_blank_present(input) ) \
+ { \
+ if(b==b) /* Blank value can be checked with an equals. */ \
+ do if(*a!=b) { ++n; sum += *a; } while(++a<af); \
+ else /* Blank value will fail an equal comparison. */ \
+ do if(*a==*a) { ++n; sum += *a; } while(++a<af); \
+ } \
+ else \
+ { \
+ do sum += *a; while(++a<af); \
+ n = input->size; \
+ } \
+ }
gal_data_t *
-gal_statistics_maximum(gal_data_t *data)
+gal_statistics_sum(gal_data_t *input)
{
- return gal_arithmetic(GAL_ARITHMETIC_OP_MAXVAL, SIMP_FLAGS, data);
+ double sum=0.0f;
+ size_t dsize=1, n=0;
+ gal_data_t *out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &dsize,
+ NULL, 1, -1, NULL, NULL, NULL);
+ switch(input->type)
+ {
+ case GAL_DATA_TYPE_UINT8: STATS_SUM_NUM( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: STATS_SUM_NUM( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: STATS_SUM_NUM( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: STATS_SUM_NUM( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: STATS_SUM_NUM( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: STATS_SUM_NUM( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: STATS_SUM_NUM( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: STATS_SUM_NUM( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: STATS_SUM_NUM( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: STATS_SUM_NUM( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`gal_statistics_maximum'", input->type);
+ }
+ *((double *)(out->array)) = n ? sum : GAL_BLANK_FLOAT64;
+ return out;
}
+
+
+
+
+/* Return the mean of the input dataset as a float64 type single-element
+ dataset. */
gal_data_t *
-gal_statistics_sum(gal_data_t *data)
+gal_statistics_mean(gal_data_t *input)
{
- return gal_arithmetic(GAL_ARITHMETIC_OP_SUMVAL, SIMP_FLAGS, data);
+ double sum=0.0f;
+ size_t dsize=1, n=0;
+ gal_data_t *out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &dsize,
+ NULL, 1, -1, NULL, NULL, NULL);
+ switch(input->type)
+ {
+ case GAL_DATA_TYPE_UINT8: STATS_SUM_NUM( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: STATS_SUM_NUM( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: STATS_SUM_NUM( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: STATS_SUM_NUM( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: STATS_SUM_NUM( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: STATS_SUM_NUM( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: STATS_SUM_NUM( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: STATS_SUM_NUM( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: STATS_SUM_NUM( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: STATS_SUM_NUM( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`gal_statistics_maximum'", input->type);
+ }
+ *((double *)(out->array)) = n ? sum/n : GAL_BLANK_FLOAT64;
+ return out;
}
+
+
+
+
+/* Return the standard deviation of the input dataset as a single element
+ dataset of type float64. */
+#define STATS_NSS2(IT) { \
+ IT b, *a=input->array, *af=a+input->size; \
+ gal_blank_write(&b, input->type); \
+ if( gal_blank_present(input) ) \
+ { \
+ if(b==b) /* Blank value can be checked with an equals. */ \
+ do if(*a!=b) { ++n; s += *a; s2 += *a * *a; } while(++a<af); \
+ else /* Blank value will fail an equal comparison. */ \
+ do if(*a==*a){ ++n; s += *a; s2 += *a * *a; } while(++a<af); \
+ } \
+ else \
+ { \
+ do { s += *a; s2 += *a * *a; } while(++a<af); \
+ n = input->size; \
+ } \
+ }
gal_data_t *
-gal_statistics_mean(gal_data_t *data)
+gal_statistics_std(gal_data_t *input)
{
- return gal_arithmetic(GAL_ARITHMETIC_OP_MEANVAL, SIMP_FLAGS, data);
+ size_t dsize=1, n=0;
+ double s=0.0f, s2=0.0f;
+ gal_data_t *out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &dsize,
+ NULL, 1, -1, NULL, NULL, NULL);
+ switch(input->type)
+ {
+ case GAL_DATA_TYPE_UINT8: STATS_NSS2( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: STATS_NSS2( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: STATS_NSS2( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: STATS_NSS2( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: STATS_NSS2( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: STATS_NSS2( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: STATS_NSS2( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: STATS_NSS2( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: STATS_NSS2( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: STATS_NSS2( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`gal_statistics_maximum'", input->type);
+ }
+ *((double *)(out->array)) = n ? sqrt( (s2-s*s/n)/n ) : GAL_BLANK_FLOAT64;
+ return out;
}
+
+
+
+
+/* Return the mean and standard deviation of a dataset in one run in type
+ float64. The output is a two element data structure, with the first
+ value being the mean and the second value the standard deviation. */
gal_data_t *
-gal_statistics_std(gal_data_t *data)
+gal_statistics_mean_std(gal_data_t *input)
{
- return gal_arithmetic(GAL_ARITHMETIC_OP_STDVAL, SIMP_FLAGS, data);
+ size_t dsize=2, n=0;
+ double s=0.0f, s2=0.0f;
+ gal_data_t *out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &dsize,
+ NULL, 1, -1, NULL, NULL, NULL);
+ switch(input->type)
+ {
+ case GAL_DATA_TYPE_UINT8: STATS_NSS2( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: STATS_NSS2( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: STATS_NSS2( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: STATS_NSS2( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: STATS_NSS2( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: STATS_NSS2( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: STATS_NSS2( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: STATS_NSS2( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: STATS_NSS2( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: STATS_NSS2( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`gal_statistics_maximum'", input->type);
+ }
+ ((double *)(out->array))[0] = n ? s/n : GAL_BLANK_FLOAT64;
+ ((double *)(out->array))[1] = n ? sqrt( (s2-s*s/n)/n ) : GAL_BLANK_FLOAT64;
+ return out;
}
+
+
+
+
+/* The input is a sorted array with no blank values, we want the median
+ value to be put inside the already allocated space which is pointed to
+ by `median'. It is in the same type as the input. */
+#define MED_IN_SORTED(IT) { \
+ IT *a=sorted->array; \
+ *(IT *)median = n%2 ? a[n/2] : (a[n/2]+a[n/2-1])/2; \
+ }
+static void
+statistics_median_in_sorted_no_blank(gal_data_t *sorted, void *median)
+{
+ size_t n=sorted->size;
+
+ switch(sorted->type)
+ {
+ case GAL_DATA_TYPE_UINT8: MED_IN_SORTED( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: MED_IN_SORTED( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: MED_IN_SORTED( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: MED_IN_SORTED( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: MED_IN_SORTED( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: MED_IN_SORTED( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: MED_IN_SORTED( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: MED_IN_SORTED( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: MED_IN_SORTED( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: MED_IN_SORTED( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`gal_statistics_median_in_sorted_no_blank'", sorted->type);
+ }
+}
+
+
+
+
+
+/* Return the median value of the dataset in the same type as the input as
+ a one element dataset. If the `inplace' flag is set, the input data
+ structure will be modified: it will have no blank values and will be
+ sorted (increasing). */
gal_data_t *
-gal_statistics_median(gal_data_t *data)
+gal_statistics_median(gal_data_t *input, int inplace)
{
- return gal_arithmetic(GAL_ARITHMETIC_OP_MEDIANVAL, SIMP_FLAGS, data);
+ size_t dsize=1;
+ gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);
+ gal_data_t *out=gal_data_alloc(NULL, input->type, 1, &dsize,
+ NULL, 1, -1, NULL, NULL, NULL);
+
+ /* Write the median. */
+ statistics_median_in_sorted_no_blank(nbs, out->array);
+
+ /* Clean up (if necessary), then return the output */
+ if(nbs!=input) gal_data_free(nbs);
+ return out;
+}
+
+
+
+
+
+/* Return a single element dataset of the same type as input keeping the
+ value that has the given quantile. */
+#define STATS_QUANT(IT) { IT *o=out->array, *a=nbs->array; *o = a[ind]; }
+gal_data_t *
+gal_statistsics_quantile(gal_data_t *input, float quantile, int inplace)
+{
+ size_t dsize=1, ind;
+ gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);
+ gal_data_t *out=gal_data_alloc(NULL, input->type, 1, &dsize,
+ NULL, 1, -1, NULL, NULL, NULL);
+
+ /* A small sanity check. */
+ if(quantile<0 || quantile>1)
+ error(EXIT_FAILURE, 0, "the `quantile' input to "
+ "`gal_statistics_quantile' must be between 0 and 1 (inclusive)");
+
+ /* Find the index of the quantile. */
+ ind=gal_statistics_quantile_index(nbs->size, quantile);
+
+ /* Set the value. */
+ switch(input->type)
+ {
+ case GAL_DATA_TYPE_UINT8: STATS_QUANT( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: STATS_QUANT( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: STATS_QUANT( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: STATS_QUANT( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: STATS_QUANT( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: STATS_QUANT( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: STATS_QUANT( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: STATS_QUANT( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: STATS_QUANT( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: STATS_QUANT( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`gal_statistics_maximum'", input->type);
+ }
+
+ /* Clean up and return. */
+ if(nbs!=input) gal_data_free(nbs);
+ return out;
+}
+
+
+
+
+
+size_t
+gal_statistics_quantile_index(size_t size, float quant)
+{
+ float floatindex;
+
+ if(quant>1.0f)
+ error(EXIT_FAILURE, 0, "the quantile in "
+ "gal_statistics_index_from_quantile should be smaller than 1.0");
+
+ /* Find the index of the quantile. */
+ floatindex=(float)size*quant;
+
+ /*
+ printf("quant: %f, size: %zu, findex: %f\n", quant, size, floatindex);
+ */
+ /* Note that in the conversion from float to size_t, the floor
+ integer value of the float will be used. */
+ if( floatindex - (int)floatindex > 0.5 )
+ return floatindex+1;
+ else
+ return floatindex;
}
@@ -249,6 +631,57 @@ gal_statistics_sort_decreasing(gal_data_t *data)
+/* Return a dataset that has doesn't have blank values and is sorted. If
+ the `inplace' value is set to 1, then the input array will be modified,
+ otherwise, a new array will be allocated with the desired properties. So
+ if it is already sorted and has blank values, the `inplace' variable is
+ irrelevant.*/
+gal_data_t *
+gal_statistics_no_blank_sorted(gal_data_t *input, int inplace)
+{
+ gal_data_t *noblank, *sorted;
+
+ /* Make sure there is no blanks in the array that will be used. After
+ this step, we won't be dealing with `input' any more, but with
+ `noblank'.*/
+ if( gal_blank_present(input) )
+ {
+ if(inplace) /* If we can free the input, then */
+ { /* just remove the blank elements */
+ gal_blank_remove(input); /* in place. */
+ noblank=input;
+ }
+ else
+ {
+ noblank=gal_data_copy(input); /* We aren't allowed to touch */
+ gal_blank_remove(input); /* the input, so make a copy. */
+ }
+ }
+ else noblank=input;
+
+ /* Make sure the array is sorted. After this step, we won't be dealing
+ with `noblank' any more but with `sorted'. */
+ if( gal_statistics_is_sorted(noblank) ) sorted=noblank;
+ else
+ {
+ if(inplace) sorted=noblank;
+ else
+ {
+ if(noblank!=input) /* no-blank has already been allocated. */
+ sorted=noblank;
+ else
+ sorted=gal_data_copy(noblank);
+ }
+ gal_statistics_sort_increasing(sorted);
+ }
+
+ return sorted;
+}
+
+
+
+
+
@@ -297,69 +730,83 @@ gal_statistics_sort_decreasing(gal_data_t *data)
* `onebinstart' A desired value for onebinstart. Note that with this
option, the bins won't start and end exactly on the given range
values, it will be slightly shifted to accommodate this
- request. */
+ request.
+
+ The output is a 1D array (column) of type double, it has to be double to
+ account for small differences on the bin edges.
+*/
gal_data_t *
-gal_statistics_regular_bins(gal_data_t *data, gal_data_t *range,
+gal_statistics_regular_bins(gal_data_t *data, gal_data_t *inrange,
size_t numbins, float onebinstart)
{
size_t i;
- gal_data_t *bins, *tmp;
- float *b, *ra, min, max, hbw, diff, binwidth;
+ gal_data_t *bins, *tmp, *range;
+ double *b, *ra, min, max, hbw, diff, binwidth;
/* Some sanity checks. */
if(numbins==0)
error(EXIT_FAILURE, 0, "`numbins' in `gal_statistics_regular_bins' "
"cannot be given a value of 0");
- if(range && range->type!=GAL_DATA_TYPE_FLOAT32)
- error(EXIT_FAILURE, 0, "gal_statistics_regular_bins currently on works "
- "on ranges of type float32");
- if(data->ndim!=1)
- error(EXIT_FAILURE, 0, "gal_statistics_regular_bins currently on works "
- "in 1D data");
- if(data->type!=GAL_DATA_TYPE_FLOAT32)
- error(EXIT_FAILURE, 0, "gal_statistics_regular_bins currently on works "
- "on float32 type data");
/* Set the minimum and maximum values. */
- if(range && range->size)
+ if(inrange && inrange->size)
{
+ /* Make sure we are dealing with a double type range. */
+ if(inrange->type==GAL_DATA_TYPE_FLOAT64)
+ range=inrange;
+ else
+ range=gal_data_copy_to_new_type(inrange, GAL_DATA_TYPE_FLOAT64);
+
+ /* Set the minimum and maximum of the bins. */
ra=range->array;
if( (range->size)%2 )
error(EXIT_FAILURE, 0, "Quantile ranges are not implemented in "
"`gal_statistics_regular_bins' yet.");
else
{
+ /* If the minimum isn't set (is blank), find it. */
if( isnan(ra[0]) )
{
- tmp=gal_statistics_minimum(data);
- min=*((float *)(tmp->array));
+ tmp=gal_data_copy_to_new_type_free(gal_statistics_minimum(data),
+ GAL_DATA_TYPE_FLOAT64);
+ min=*((double *)(tmp->array));
gal_data_free(tmp);
}
else min=ra[0];
- if( isnan(ra[1]) ) /* When the maximum */
- { /* isn't set, we'll */
- tmp=gal_statistics_maximum(data); /* Add a very small */
- max=*((float *)(tmp->array)) + 1e-5; /* value so all the */
- gal_data_free(tmp); /* points are counted. */
+
+ /* For the maximum, when it isn't set, we'll add a very small
+ value, so all points are included. */
+ if( isnan(ra[1]) )
+ {
+ tmp=gal_data_copy_to_new_type_free(gal_statistics_maximum(data),
+ GAL_DATA_TYPE_FLOAT64);
+ max=*((double *)(tmp->array))+1e-6;
+ gal_data_free(tmp);
}
else max=ra[1];
}
+
+ /* Clean up: if `range' was allocated. */
+ if(range!=inrange) gal_data_free(range);
}
+ /* No range was given, find the minimum and maximum. */
else
{
- tmp=gal_statistics_minimum(data);
- min=*((float *)(tmp->array));
+ tmp=gal_data_copy_to_new_type_free(gal_statistics_minimum(data),
+ GAL_DATA_TYPE_FLOAT64);
+ min=*((double *)(tmp->array));
gal_data_free(tmp);
- tmp=gal_statistics_maximum(data);
- max=*((float *)(tmp->array));
+ tmp=gal_data_copy_to_new_type_free(gal_statistics_maximum(data),
+ GAL_DATA_TYPE_FLOAT64);
+ max=*((double *)(tmp->array)) + 1e-6;
gal_data_free(tmp);
}
/* Allocate the space for the bins. */
- bins=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, 1, &numbins, NULL,
+ bins=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &numbins, NULL,
0, data->minmapsize, "bin_center", data->unit,
"Center value of each bin.");
@@ -402,63 +849,85 @@ gal_statistics_regular_bins(gal_data_t *data, gal_data_t
*range,
/* Make a histogram of all the elements in the given dataset with bin
- values that are defined in the `bins' structure (see
- `gal_statistics_regular_bins'). */
+ values that are defined in the `inbins' structure (see
+ `gal_statistics_regular_bins'). `inbins' is not mandatory, if you pass a
+ NULL pointer, the bins structure will be built within this function
+ based on the `numbins' input. As a result, when you have already defined
+ the bins, `numbins' is not used. */
+
+#define HISTOGRAM_TYPESET(IT) { \
+ IT *a=data->array, *af=a+data->size; \
+ do if( *a>=min && *a<max) ++h[ (size_t)( (*a-min)/binwidth ) ]; \
+ while(++a<af); \
+ }
+
gal_data_t *
-gal_statistics_histogram(gal_data_t *data, gal_data_t *bins,
- int normalize, int maxhistone)
+gal_statistics_histogram(gal_data_t *data, gal_data_t *bins, int normalize,
+ int maxone)
{
- size_t i, *h;
- double ref=NAN;
+ size_t *h;
+ float *f, *ff;
gal_data_t *hist;
- float *f, *ff, min, max, binwidth;
-
-
- /* Some basic sanity checks for now (until it is generalized). */
- if(data->ndim!=1)
- error(EXIT_FAILURE, 0, "gal_statistics_histogram currently on works "
- "in 1D data");
- if(data->type!=GAL_DATA_TYPE_FLOAT32)
- error(EXIT_FAILURE, 0, "gal_statistics_histogram currently on works "
- "on float32 type data");
+ double *d, min, max, ref=NAN, binwidth;
/* Check if the bins are regular or not. For irregular bins, we can
either use the old implementation, or GSL's histogram
functionality. */
+ if(bins==NULL)
+ error(EXIT_FAILURE, 0, "no `bins' in `gal_statistics_histogram");
if(bins->status!=GAL_STATISTICS_BINS_REGULAR)
error(EXIT_FAILURE, 0, "the input bins to `gal_statistics_histogram' "
"are not regular. Currently it is only implemented for regular "
"bins");
- /* Check if normalize and maxhistone are not called together. */
- if(normalize && maxhistone)
- error(EXIT_FAILURE, 0, "only one of `normalize' and `maxhistone' may "
+ /* Check if normalize and `maxone' are not called together. */
+ if(normalize && maxone)
+ error(EXIT_FAILURE, 0, "only one of `normalize' and `maxone' may "
"be given to `gal_statistics_histogram'");
- /* Allocate the histogram (note that we are clearning it. */
+
+ /* Allocate the histogram (note that we are clearning it so all values
+ are zero. */
hist=gal_data_alloc(NULL, GAL_DATA_TYPE_SIZE_T, bins->ndim, bins->dsize,
NULL, 1, data->minmapsize, "hist_number", "counts",
"Number of data points within each bin.");
- /* Set the minimum and maximum values: */
- f=bins->array;
- binwidth=f[1]-f[0];
- max = f[bins->size - 1] + binwidth/2;
- min = f[0] - binwidth/2;
+ /* Set the minimum and maximum range of the histogram from the bins. */
+ d=bins->array;
+ binwidth=d[1]-d[0];
+ max = d[bins->size - 1] + binwidth/2;
+ min = d[0] - binwidth/2;
/* Go through all the elements and find out which bin they belong to. */
h=hist->array;
- f=data->array;
- for(i=0;i<data->size;++i)
- if(f[i]>=min && f[i]<max)
- {
- ++h[ (size_t)((f[i]-min)/binwidth) ];
- /*printf("%-10.3f%zu\n", f[i], (size_t)((f[i]-min)/binwidth) );*/
- }
+ switch(data->type)
+ {
+ case GAL_DATA_TYPE_UINT8: HISTOGRAM_TYPESET(uint8_t); break;
+ case GAL_DATA_TYPE_INT8: HISTOGRAM_TYPESET(int8_t); break;
+ case GAL_DATA_TYPE_UINT16: HISTOGRAM_TYPESET(uint16_t); break;
+ case GAL_DATA_TYPE_INT16: HISTOGRAM_TYPESET(int16_t); break;
+ case GAL_DATA_TYPE_UINT32: HISTOGRAM_TYPESET(uint32_t); break;
+ case GAL_DATA_TYPE_INT32: HISTOGRAM_TYPESET(int32_t); break;
+ case GAL_DATA_TYPE_UINT64: HISTOGRAM_TYPESET(uint64_t); break;
+ case GAL_DATA_TYPE_INT64: HISTOGRAM_TYPESET(int64_t); break;
+ case GAL_DATA_TYPE_FLOAT32: HISTOGRAM_TYPESET(float); break;
+ case GAL_DATA_TYPE_FLOAT64: HISTOGRAM_TYPESET(double); break;
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`gal_statistics_histogram'", data->type);
+ }
+
+
+ /* For a check:
+ {
+ size_t i, *hh=hist->array;
+ for(i=0;i<hist->size;++i) printf("%-10.4f%zu\n", f[i], hh[i]);
+ }
+ */
/* Find the reference to correct the histogram if necessary. */
@@ -476,7 +945,7 @@ gal_statistics_histogram(gal_data_t *data, gal_data_t *bins,
gal_checkset_allocate_copy("Normalized histogram value for this bin",
&hist->comment);
}
- if(maxhistone)
+ if(maxone)
{
/* Calculate the reference. */
ref=-FLT_MAX;
@@ -498,6 +967,7 @@ gal_statistics_histogram(gal_data_t *data, gal_data_t *bins,
if( !isnan(ref) )
{ ff=(f=hist->array)+hist->size; do *f++ /= ref; while(f<ff); }
+
/* Return the histogram. */
return hist;
}
@@ -531,15 +1001,6 @@ gal_statistics_cfp(gal_data_t *data, gal_data_t *bins,
int normalize)
size_t *s, *sf, *hs, sums;
- /* Some basic sanity checks for now (until it is generalized). */
- if(data->ndim!=1)
- error(EXIT_FAILURE, 0, "gal_statistics_cfp currently on works "
- "in 1D data");
- if(data->type!=GAL_DATA_TYPE_FLOAT32)
- error(EXIT_FAILURE, 0, "gal_statistics_cfp currently on works "
- "on float32 type data");
-
-
/* Check if the bins are regular or not. For irregular bins, we can
either use the old implementation, or GSL's histogram
functionality. */
@@ -557,8 +1018,8 @@ gal_statistics_cfp(gal_data_t *data, gal_data_t *bins, int
normalize)
/* If the histogram has float32 type it was given by the user and is
either normalized or its maximum was set to 1. We can only use it if
- it was normalized. If its maximum is 1, then we must ignore it and
- build the histogram again.*/
+ it was normalized. If it isn't normalized, then we must ignore it and
+ build the histogram here.*/
if(hist->type==GAL_DATA_TYPE_FLOAT32)
{
sum=0.0f;
@@ -567,6 +1028,7 @@ gal_statistics_cfp(gal_data_t *data, gal_data_t *bins, int
normalize)
hist=gal_statistics_histogram(data, bins, 0, 0);
}
+
/* Allocate the cumulative frequency plot's necessary space. */
cfp=gal_data_alloc( NULL, hist->type, bins->ndim, bins->dsize,
NULL, 1, data->minmapsize,
@@ -644,203 +1106,214 @@ gal_statistics_cfp(gal_data_t *data, gal_data_t *bins,
int normalize)
/****************************************************************
- ***************** Quantiles ********************
+ ***************** Sigma clip ********************
****************************************************************/
-/* Find the index corresponding to a certain quantile, considering the
- rounding that might be needed. */
-size_t
-gal_statistics_index_from_quantile(size_t size, float quant)
-{
- float floatindex;
-
- if(quant>1.0f)
- error(EXIT_FAILURE, 0, "the quantile in gal_statistics_index_from_quantile
"
- "(statistics.c) Should be smaller");
-
- /* Find the index of the quantile. */
- floatindex=(float)size*quant;
-
- /*
- printf("quant: %f, size: %zu, findex: %f\n", quant, size, floatindex);
- */
- /* Note that in the conversion from float to size_t, the floor
- integer value of the float will be used. */
- if( floatindex - (int)floatindex > 0.5 )
- return floatindex+1;
- else
- return floatindex;
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
+/* Return a data structure with an array of four values: the final number
+ of points used, the median, average and standard deviation. The number
+ of clips is put into the `status' element of the data structure.
+ Inputs:
+ - `multip': multiple of the standard deviation,
+ - `param' must be positive and determines the type of clipping:
+ - param<1.0: interpretted as a tolerance level to stop clipping.
+ - param>=1.0 and an integer: a specific number of times to do the
+ clippping.
+ The way this function works is very simple: first it will sort the input
+ (if it isn't sorted). Afterwards, it will recursively change the starting
+ point of the array and its size, calcluating the basic statistics in each
+ round to define the new starting point and size.
+*/
-/****************************************************************
- ***************** Sigma clip ********************
- ****************************************************************/
-/* This function will repeatedly sigma clip the data and return the
- median. It is assumed that the data have been ordered.
+#define SIGCLIP(IT) { \
+ IT *a = sorted->array, *af = a + sorted->size; \
+ IT *bf = sorted->array, *b = bf + sorted->size - 1; \
+ \
+ /* Remove all out-of-range elements from the start of the array. */ \
+ if(sortstatus==GAL_STATISTICS_SORTED_INCREASING) \
+ do if( *a > (*med - (multip * *std)) ) \
+ { start=a; break; } \
+ while(++a<af); \
+ else \
+ do if( *a < (*med + (multip * *std)) ) \
+ { start=a; break; } \
+ while(++a<af); \
+ \
+ /* Remove all out-of-range elements from the end of the array. */ \
+ if(sortstatus==GAL_STATISTICS_SORTED_INCREASING) \
+ do if( *b < (*med + (multip * *std)) ) \
+ { size=b-a+1; break; } \
+ while(--b>=bf); \
+ else \
+ do if( *b > (*med - (multip * *std)) ) \
+ { size=b-a+1; break; } \
+ while(--b>=bf); \
+ }
- o1_n0: Ordered (1), not ordered (0).
-*/
-int
-gal_statistics_sigma_clip_converge(float *array, int o1_n0, size_t num_elem,
- float sigma_multiple, float accuracy,
- float *outave, float *outmed,
- float *outstd, int print)
+gal_data_t *
+gal_statistics_sigma_clip(gal_data_t *input, float multip, float param,
+ int quiet)
{
- printf("\n ... in gal_statistics_sigma_clip_converge ... \n");
- exit(1);
-#if 0
- size_t counter=0;
- float *start, *oldstart, *dpt;
- float ave=*outave, med=*outmed, std=*outstd;
- float oldstd=0, oldmed=0, oldave=0, *orderedarray;
+ int sortstatus;
+ double *med, *mean, *std;
+ void *start, *sorted_array;
+ double oldmed, oldmean, oldstd;
+ size_t num=0, dsize=4, size, oldsize;
+ size_t maxnum = param>=1.0f ? param : 50; /* for failing to converge */
+ uint8_t bytolerance = param>=1.0f ? 0 : 1;
+ gal_data_t *sorted, *median_i, *median_d, *out, *meanstd, *noblank;
- if(o1_n0==0)
+ /* Some sanity checks. */
+ if( multip<=0 )
+ error(EXIT_FAILURE, 0, "`multip', must be greater than zero in "
+ "`gal_statistics_sigma_clip'. The given value was %g", multip);
+ if( param<=0 )
+ error(EXIT_FAILURE, 0, "`param', must be greater than zero in "
+ "`gal_statistics_sigma_clip'. The given value was %g", param);
+ if( param >= 1.0f && ceil(param) != param )
+ error(EXIT_FAILURE, 0, "when `param' is larger than 1.0, it is "
+ "interpretted as an absolute number of clips in "
+ "`gal_statistics_sigma_clip'. So it must be an integer. However, "
+ "your given value %g", param);
+
+
+ /* If there are blank elements, remove them (from a copied array). NOTE:
+ we don't want to change the input. */
+ if( gal_blank_present(input) )
{
- gal_array_float_copy(array, num_elem, &orderedarray);
- qsort(orderedarray, num_elem, sizeof*orderedarray,
- gal_qsort_float32_increasing);
+ noblank = gal_data_copy(input);
+ gal_blank_remove(noblank);
}
- else orderedarray=array;
+ else noblank=input;
- start=orderedarray;
- while(counter<GAL_STATISTICS_MAX_SIG_CLIP_CONVERGE)
+
+ /* We want to have sorted (increasing) array. */
+ sortstatus=gal_statistics_is_sorted(noblank);
+ switch(sortstatus)
{
- oldstart=start;
+ case GAL_STATISTICS_SORTED_NOT:
+ /* Only copy (or allocate) space if it wasn't already allocated. */
+ sorted = (noblank==input) ? gal_data_copy(input) : noblank;
+ gal_statistics_sort_increasing(sorted);
+ sortstatus=GAL_STATISTICS_SORTED_INCREASING;
+ break;
- med=*(start+num_elem/2);
- gal_statistics_f_ave_std(start, num_elem, &ave, &std, NULL);
+ case GAL_STATISTICS_SORTED_DECREASING:
+ case GAL_STATISTICS_SORTED_INCREASING:
+ sorted=noblank;
+ break;
- if(print)
- printf(" %zu: %f %f %f %zu\n",
- counter+1, med, ave, std, num_elem);
+ default:
+ error(EXIT_FAILURE, 0, "sorted code %d not recognized in "
+ "`gal_statistics_sigma_clip'", gal_statistics_is_sorted(input));
+ }
- /* It might happen that ave and std are NaN. If so, stop the
- process here (the user has not given a mask and some pixels
- have nan values!). */
- if(isnan(ave) || isnan(std))
- return 0;
- /* Normally, oldstd should be larger than std, because the
- possible outliers have been removed. If it is not, it means
- that we have clipped too much! */
- if(counter>0 && (oldstd-std)/std<accuracy)
- {
- *outstd=oldstd; *outave=oldave; *outmed=oldmed;
- return 1;
- }
+ /* Allocate the necessary spaces. */
+ out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, 1, &dsize, NULL, 0,
+ input->minmapsize, NULL, NULL, NULL);
+ median_i=gal_data_alloc(NULL, sorted->type, 1, &dsize, NULL, 0,
+ input->minmapsize, NULL, NULL, NULL);
- for(dpt=start; dpt<start+num_elem; ++dpt)
- if (*dpt>med-sigma_multiple*std)
- {
- start=dpt;
- break;
- }
-
- for(dpt=oldstart+num_elem-1;dpt>start;dpt--)
- if (*dpt<med+sigma_multiple*std)
- {
- num_elem=dpt-start+1;
- break;
- }
-
- oldave=ave;
- oldmed=med;
- oldstd=std;
- ++counter;
- }
-#endif
- return 0;
-}
+ /* Print the comments. */
+ if(!quiet)
+ printf("%-8s %-10s %-15s %-15s %-15s\n",
+ "round", "number", "median", "mean", "STD");
+ /* Do the clipping, but first initialize the values that will be changed
+ during the clipping: the start of the array and the array's size. */
+ size=sorted->size;
+ sorted_array=start=sorted->array;
+ while(num<maxnum)
+ {
+ /* Find the median. */
+ statistics_median_in_sorted_no_blank(sorted, median_i->array);
+ median_d=gal_data_copy_to_new_type(median_i, GAL_DATA_TYPE_FLOAT64);
+
+ /* Find the average and Standard deviation, note that both `start'
+ and `size' will be different in the next round. */
+ sorted->array = start;
+ sorted->size = oldsize = size;
+ meanstd=gal_statistics_mean_std(sorted);
+
+ /* Put the three final values in usable (with a type) pointers. */
+ med = median_d->array;
+ mean = meanstd->array;
+ std = &((double *)(meanstd->array))[1];
+
+ /* If the user wanted to view the steps, show it to them. */
+ if(!quiet)
+ printf("%-8zu %-10zu %-15g %-15g %-15g\n",
+ num+1, size, *med, *mean, *std);
+
+ /* If we are to work by tolerance, then check if we should jump out
+ of the loop. Normally, `oldstd' should be larger than std, because
+ the possible outliers have been removed. If it is not, it means
+ that we have clipped too much and must stop anyway, so we don't
+ need an absolute value on the difference! */
+ if( bytolerance && num>0 && ((oldstd - *std) / *std) < param )
+ break;
+
+ /* Clip all the elements outside of the desired range: since the
+ array is sorted, this means to just change the starting pointer
+ and size of the array. */
+ switch(sorted->type)
+ {
+ case GAL_DATA_TYPE_UINT8: SIGCLIP( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: SIGCLIP( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: SIGCLIP( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: SIGCLIP( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: SIGCLIP( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: SIGCLIP( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: SIGCLIP( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: SIGCLIP( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: SIGCLIP( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: SIGCLIP( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`gal_statistics_sigma_clip'", sorted->type);
+ }
+ /* Set the values from this round in the old elements, so the next
+ round can compare with, and return then if necessary. */
+ oldmed = *med;
+ oldstd = *std;
+ oldmean = *mean;
+ ++num;
-/* This function will do a certain number of sigma clips and return
- the final average, median and std. o1_n0: 1: initially ordered. 2:
- initially not ordered.*/
-int
-gal_statistics_sigma_clip_certain_num(float *array, int o1_n0, size_t num_elem,
- float sigma_multiple, size_t numtimes,
- float *outave, float *outmed,
- float *outstd, int print)
-{
- printf("\n ... in gal_statistics_sigma_clip_certain_num ... \n");
- exit(1);
-#if 0
- size_t counter=0;
- float ave=*outave, med=*outmed, std=*outstd;
- float *start, *oldstart, *dpt, *orderedarray;
+ /* Clean up: */
+ gal_data_free(meanstd);
+ gal_data_free(median_d);
+ }
- if(o1_n0==0)
+ /* If we were in tolerance mode and `num' and `maxnum' are equal (the
+ loop didn't stop by tolerance), so the outputs should be NaN. */
+ out->status=num;
+ if( bytolerance && num==maxnum )
{
- gal_array_float_copy(array, num_elem, &orderedarray);
- qsort(orderedarray, num_elem, sizeof*orderedarray,
- gal_qsort_float32_increasing);
+ ((float *)(out->array))[0] = NAN;
+ ((float *)(out->array))[1] = NAN;
+ ((float *)(out->array))[2] = NAN;
+ ((float *)(out->array))[3] = NAN;
}
- else orderedarray=array;
-
- start=orderedarray;
- for(counter=0;counter<numtimes;++counter)
+ else
{
- oldstart=start;
-
- med=*(start+num_elem/2);
- gal_statistics_f_ave_std(start, num_elem, &ave, &std, NULL);
-
- if(print)
- printf(" %zu: %f %f %f %zu\n",
- counter+1, med, ave, std, num_elem);
-
- /* It might happen that ave and std are nan. If so, stop the
- process here (the user has not given a mask and some pixels
- have nan values!). */
- if(isnan(ave) || isnan(std))
- return 0;
-
-
- for(dpt=start; dpt<start+num_elem; ++dpt)
- if (*dpt>med-sigma_multiple*std)
- {
- start=dpt;
- break;
- }
-
-
- for(dpt=oldstart+num_elem-1;dpt>start;dpt--)
- if (*dpt<med+sigma_multiple*std)
- {
- num_elem=dpt-start+1;
- break;
- }
+ ((float *)(out->array))[0] = size;
+ ((float *)(out->array))[1] = oldmed;
+ ((float *)(out->array))[2] = oldmean;
+ ((float *)(out->array))[3] = oldstd;
}
- if(o1_n0==0)
- free(orderedarray);
-
- *outave=ave;
- *outmed=med;
- *outstd=std;
-#endif
- return 1;
+ /* Clean up and return. */
+ sorted->array=sorted_array;
+ if(sorted!=input) gal_data_free(sorted);
+ return out;
}
@@ -861,6 +1334,7 @@ gal_statistics_sigma_clip_certain_num(float *array, int
o1_n0, size_t num_elem,
+
/****************************************************************/
/************* Identify outliers ****************/
/****************************************************************/
@@ -1089,7 +1563,7 @@ gal_statistics_mode_value_from_sym(float *sorted, size_t
size,
{
float mf=sorted[modeindex];
float af=
- sorted[gal_statistics_index_from_quantile(2*modeindex+1,
+ sorted[gal_statistics_quantile_index(2*modeindex+1,
GAL_STATISTICS_MODE_SYMMETRICITY_LOW_QUANT)];
return sym*(mf-af)+mf;
}
@@ -1123,9 +1597,9 @@ gal_statistics_mode_index_in_sorted(float *sorted, size_t
size, float errorstdm,
if(mp.numcheck>1000)
mp.interval=mp.numcheck/1000;
else mp.interval=1;
- mp.lowi = gal_statistics_index_from_quantile(size,
+ mp.lowi = gal_statistics_quantile_index(size,
GAL_STATISTICS_MODE_LOW_QUANTILE);
- mp.highi = gal_statistics_index_from_quantile(size,
+ mp.highi = gal_statistics_quantile_index(size,
GAL_STATISTICS_MODE_HIGH_QUANTILE);
mp.midi = (((float)mp.highi
+ MODE_GOLDEN_RATIO*(float)mp.lowi)
diff --git a/lib/table.c b/lib/table.c
index 597880c..eb8569f 100644
--- a/lib/table.c
+++ b/lib/table.c
@@ -29,10 +29,12 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <stdlib.h>
#include <string.h>
+#include <gnuastro/git.h>
#include <gnuastro/txt.h>
#include <gnuastro/blank.h>
#include <gnuastro/table.h>
+#include <timing.h>
#include <checkset.h>
@@ -699,16 +701,10 @@ make_list_of_indexs(struct gal_linkedlist_stll *cols,
gal_data_t *allcols,
columns in the input catalog (note that the user is
counting from 1, not 0!) */
if(tlong>numcols)
- {
- if(gal_fits_name_is_fits(filename))
- error(EXIT_FAILURE, 0, "%s (hdu %s): has %zu columns, "
- "but you have asked for column number %zu",
- filename, hdu, numcols, tlong);
- else
- error(EXIT_FAILURE, 0, "%s: has %zu columns, but you "
- "have asked for column number %zu", filename,
- numcols, tlong);
- }
+ error(EXIT_FAILURE, 0, "%s: has %zu columns, but you "
+ "have asked for column number %zu",
+ gal_fits_name_save_as_string(filename, hdu),
+ numcols, tlong);
/* Everything seems to be fine, put this column number in the
output column numbers linked list. Note that internally,
@@ -860,13 +856,48 @@ gal_table_read(char *filename, char *hdu, struct
gal_linkedlist_stll *cols,
/************************************************************************/
/*************** Write a table ***************/
/************************************************************************/
+/* Write the basic information that is necessary by each program into the
+ comments field. Note that the `comments' has to be already sorted in the
+ proper order. */
+void
+gal_table_comments_add_intro(struct gal_linkedlist_stll **comments,
+ char *program_string, time_t *rawtime)
+{
+ char gitdescribe[100], *tmp;
+
+ /* Get the Git description in the running folder. */
+ tmp=gal_git_describe();
+ if(tmp) sprintf(gitdescribe, " from %s,", tmp);
+ else gitdescribe[0]='\0';
+ free(tmp);
+
+
+ /* Git version and time of program's starting, this will be the second
+ line. Note that ctime puts a `\n' at the end of its string, so we'll
+ have to remove that. Also, note that since we are allocating `msg', we
+ are setting the allocate flag of `gal_linkedlist_add_to_stll' to 0. */
+ asprintf(&tmp, "Created%s on %s", gitdescribe, ctime(rawtime));
+ tmp[ strlen(tmp)-1 ]='\0';
+ gal_linkedlist_add_to_stll(comments, tmp, 0);
+
+
+ /* Program name: this will be the top of the list (first line). We will
+ need to set the allocation flag for this one, because program_string
+ is usually statically allocated.*/
+ gal_linkedlist_add_to_stll(comments, program_string, 1);
+}
+
+
+
+
+
/* The input is a linked list of data structures and some comments. The
table will then be written into `filename' with a format that is
specified by `tableformat'. If it already exists, and `dontdelete' has a
value of 1, then it won't be deleted and an error will be printed. */
void
-gal_table_write(gal_data_t *cols, char *comments, int tableformat,
- char *filename, int dontdelete)
+gal_table_write(gal_data_t *cols, struct gal_linkedlist_stll *comments,
+ int tableformat, char *filename, int dontdelete)
{
/* If a filename was given, then the tableformat is relevant and must be
used. When the filename is empty, a text table must be printed on the
@@ -882,3 +913,30 @@ gal_table_write(gal_data_t *cols, char *comments, int
tableformat,
else
gal_txt_write(cols, comments, filename, dontdelete);
}
+
+
+
+
+
+void
+gal_table_write_log(gal_data_t *logll, char *program_string,
+ time_t *rawtime, struct gal_linkedlist_stll *comments,
+ char *filename, int dontdelete, int quiet)
+{
+ char *msg;
+
+ /* Write all the comments into */
+ gal_table_comments_add_intro(&comments, program_string, rawtime);
+
+ /* Write the log file to disk */
+ gal_table_write(logll, comments, GAL_TABLE_FORMAT_TXT, filename,
+ dontdelete);
+
+ /* In verbose mode, print the information. */
+ if(!quiet)
+ {
+ asprintf(&msg, "%s created.", filename);
+ gal_timing_report(NULL, msg, 1);
+ free(msg);
+ }
+}
diff --git a/lib/txt.c b/lib/txt.c
index 684972d..4d7d68b 100644
--- a/lib/txt.c
+++ b/lib/txt.c
@@ -1164,15 +1164,18 @@ txt_print_value(FILE *fp, void *array, int type, size_t
ind, char *fmt)
static FILE *
-txt_open_file_write_info(gal_data_t *datall, char **fmts, char *comment,
+txt_open_file_write_info(gal_data_t *datall, char **fmts,
+ struct gal_linkedlist_stll *comment,
char *filename, int dontdelete)
{
FILE *fp;
gal_data_t *data;
char *tmp, *nstr;
size_t i, j, num=0;
+ struct gal_linkedlist_stll *strt;
int nlen, nw=0, uw=0, tw=0, bw=0;
+
/* Check the file and open it. */
gal_checkset_check_remove_file(filename, dontdelete);
errno=0;
@@ -1183,7 +1186,8 @@ txt_open_file_write_info(gal_data_t *datall, char **fmts,
char *comment,
/* Write the comments if there were any. */
- if(comment) fprintf(fp, "%s\n", comment);
+ for(strt=comment; strt!=NULL; strt=strt->next)
+ fprintf(fp, "# %s\n", strt->v);
/* Get the maximum width for each information field. */
@@ -1253,8 +1257,8 @@ txt_open_file_write_info(gal_data_t *datall, char **fmts,
char *comment,
void
-gal_txt_write(gal_data_t *datall, char *comment, char *filename,
- int dontdelete)
+gal_txt_write(gal_data_t *datall, struct gal_linkedlist_stll *comment,
+ char *filename, int dontdelete)
{
FILE *fp;
char **fmts;
diff --git a/tests/statistics/basicstats.sh b/tests/statistics/basicstats.sh
index e330862..c875ea4 100755
--- a/tests/statistics/basicstats.sh
+++ b/tests/statistics/basicstats.sh
@@ -24,7 +24,7 @@
# file exists (basicchecks.sh is in the source tree).
prog=statistics
execname=../bin/$prog/ast$prog
-img=convolve_spatial_warped_noised.fits
+img=convolve_spatial_scaled_noised.fits
@@ -48,4 +48,4 @@ if [ ! -f $execname ] || [ ! -f $img ]; then exit 77; fi
# Actual test script
# ==================
-$execname $img --checkmode --mirrorquant=0.5
+$execname $img -g9500 -l11000 --numasciibins=65
- [gnuastro-commits] master e90b490 116/125: Internal headers moved to special directory, (continued)
- [gnuastro-commits] master e90b490 116/125: Internal headers moved to special directory, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 29d7e46 091/125: Statistics now successfully uses gal_data_t, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 60c126f 120/125: New implementation of over-segmentation to find clumps, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master c250a6d 123/125: NoiseChisel's clumps grown and objects defined, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master ffaa059 101/125: Tessellation options now common to all programs, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master dec2b22 078/125: ImageCrop works with gal_data_t features, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master f8b1820 061/125: Other configuration files read immediately when seen, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 46b4651 085/125: Image prefix removed from some program names, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master df716bc 097/125: Tesselation with the new data structure: work started, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master db39cee 088/125: Statistics returns basic results, no array lib, new blank lib, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 60cd76e 090/125: Statistics mostly complete,
Mohammad Akhlaghi <=
- [gnuastro-commits] master 15c524d 124/125: New implementation of NoiseChisel complete, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 3b26d5c 114/125: SubtractSky program removed, Statistics estimates Sky, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 99ac7dc 083/125: Gnuastro now uses fixed width (from stdint.h) types, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master eaf19a3 115/125: NoiseChisel's initial detection step implemented, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 3339b92 125/125: New implementation of MakeCatalog now complete, Mohammad Akhlaghi, 2017/04/23