[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master dc8615a 2/2: Library (arithmetic): new quantil
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master dc8615a 2/2: Library (arithmetic): new quantile operator for multiple datasets |
|
Date: |
Wed, 8 Jan 2020 20:00:06 -0500 (EST) |
branch: master
commit dc8615aa81cae85fed2933d5bc9c5b0bd8fabe0b
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Library (arithmetic): new quantile operator for multiple datasets
The quantile operator on the pixels of many input datasets is useful in
some data analysis scenarios, with this commit this operator is now
available in the arithmetic library and thus the programs: Arithmetic (for
images/cubes) and Table (for columns).
This was suggested by Alejandro Serrano Borlaff.
---
NEWS | 5 ++++
bin/arithmetic/arithmetic.c | 7 +++--
doc/announce-acknowledge.txt | 1 +
doc/gnuastro.texi | 23 +++++++++++++--
lib/arithmetic.c | 67 ++++++++++++++++++++++++++++++++++++++++++++
lib/gnuastro/arithmetic.h | 1 +
6 files changed, 100 insertions(+), 4 deletions(-)
diff --git a/NEWS b/NEWS
index 6e6c4e2..301134c 100644
--- a/NEWS
+++ b/NEWS
@@ -7,6 +7,9 @@ See the end of the file for license conditions.
** New features
+ Arithmetic:
+ - New `quantile' operator for coadding datasets.
+
CosmicCalculator:
--listlines: list the pre-defined spectral line wavelengths and names
(which you can use with the `--obsline' and `--lineatz' options). This
@@ -19,10 +22,12 @@ See the end of the file for license conditions.
Table:
--equal: Can now work on columns with string type also.
--notequal: Can now work on columns with string type also.
+ - New `quantile' operator in column arithmetic.
Library:
- GAL_SPECLINES_INVALID_MAX: Total number of spectral lines, plus 1.
- gal_txt_trim_space: trim white space before and after a string.
+ - GAL_ARITHMETIC_OP_QUANTILE: operator for `gal_arithmetic'.
** Removed features
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 28526d4..74b7931 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -77,6 +77,9 @@ pop_number_of_operands(struct arithmeticparams *p, int op,
char *token_string,
/* See if this operator needs any parameters. If so, pop them. */
switch(op)
{
+ case GAL_ARITHMETIC_OP_QUANTILE:
+ numparams=1;
+ break;
case GAL_ARITHMETIC_OP_SIGCLIP_STD:
case GAL_ARITHMETIC_OP_SIGCLIP_MEAN:
case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN:
@@ -98,8 +101,8 @@ pop_number_of_operands(struct arithmeticparams *p, int op,
char *token_string,
tmp=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
gal_list_data_add(params, tmp);
- /* A small sanity check (none of the parameters for sigma-clipping
- can be negative.. */
+ /* A small sanity check (none of the parameters for sigma-clipping,
+ or quantile estimation can be negative. */
if( ((float *)(tmp->array))[0]<=0.0 )
error(EXIT_FAILURE, 0, "the %s popped operand of the \"%s\" "
"operator cannot be negative", cstring, token_string);
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index 69a0278..f604d58 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -1,5 +1,6 @@
Alphabetically ordered list to acknowledge in the next release.
+Alejandro Serrano Borlaff
Joseph Putko
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 6eb4a79..1fc5945 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -10108,6 +10108,19 @@ For each pixel, find the median in all given datasets.
The output will have the a single-precision (32-bit) floating point type.
This operator is called similar to the @command{min} operator, please see
there for more.
+@item quantile
+For each pixel, find the quantile from all given datasets.
+The output will have the same numeric data type and size as the input datasets.
+Besides the input datasets, the quantile operator also needs a single
parameter (the requested quantile).
+The parameter should be the first popped operand, with a value between (and
including) 0 and 1.
+The second popped operand must be the number of datasets to use.
+
+In the example below, the first-popped operand (@command{0.7}) is the
quantile, the second-popped operand (@command{3}) is the number of datasets to
pop.
+
+@example
+astarithmetic a.fits b.fits c.fits 3 0.7 quantile
+@end example
+
@item sigclip-number
For each pixel, find the sigma-clipped number (after removing outliers) in all
given datasets.
The output will have the an unsigned 32-bit integer type (see @ref{Numeric
data types}).
@@ -10117,7 +10130,7 @@ This operator is very similar to @command{min}, with
the exception that it expec
The first popped operand is the termination criteria and the second is the
multiple of @mymath{\sigma}.
For example in the command below, the first popped operand (@command{0.2}) is
the sigma clipping termination criteria.
-If the termination criteria is larger than 1 it is interpreted as the number
of clips to do.
+If the termination criteria is larger than, or equal to, 1 it is interpreted
as the number of clips to do.
But if it is between 0 and 1, then it is the tolerance level on the standard
deviation (see @ref{Sigma clipping}).
The second popped operand (@command{5}) is the multiple of sigma to use in
sigma-clipping.
The third popped operand (@command{10}) is number of datasets that will be
used (similar to the first popped operand to @command{min}).
@@ -22226,6 +22239,12 @@ output will be @code{GAL_TYPE_UINT32} and for the
rest, it will be
@code{GAL_TYPE_FLOAT32}.
@end deffn
+@deffn Macro GAL_ARITHMETIC_OP_QUANTILE
+Similar to the operands above (including @code{GAL_ARITHMETIC_MIN}), except
that when @code{gal_arithmetic} is called with these operators, it requires two
arguments.
+The first is the list of datasets like before, and the second is the 1-element
dataset with the quantile value.
+The output type is the same as the inputs.
+@end deffn
+
@deffn Macro GAL_ARITHMETIC_OP_SIGCLIP_STD
@deffnx Macro GAL_ARITHMETIC_OP_SIGCLIP_MEAN
@deffnx Macro GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN
@@ -23521,7 +23540,7 @@ Return the index of the element that has a quantile of
@code{quantile}
assuming the dataset has @code{size} elements.
@end deftypefun
-@deftypefun size_t gal_statistics_quantile (gal_data_t @code{*input}, double
@code{quantile}, int @code{inplace})
+@deftypefun {gal_data_t *} gal_statistics_quantile (gal_data_t @code{*input},
double @code{quantile}, int @code{inplace})
Return a single-element dataset containing the value with in a quantile
@code{quantile} of the non-blank values in @code{input}. The numerical
datatype of the output is the same as @code{input}. See
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index c391f86..9cb21e4 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -932,6 +932,51 @@ struct multioperandparams
+#define MULTIOPERAND_QUANTILE(TYPE) { \
+ size_t n, j; \
+ gal_data_t *quantile; \
+ TYPE *o=p->out->array; \
+ TYPE *pixs=gal_pointer_allocate(p->list->type, p->dnum, 0, \
+ __func__, "pixs"); \
+ gal_data_t *cont=gal_data_alloc(pixs, p->list->type, 1, &p->dnum, \
+ NULL, 0, -1, 1, NULL, NULL, NULL); \
+ \
+ /* Go over all the pixels assigned to this thread. */ \
+ for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
+ { \
+ /* Initialize, `j' is desired pixel's index. */ \
+ n=0; \
+ j=tprm->indexs[tind]; \
+ \
+ /* Read the necessay values from each input. */ \
+ for(i=0;i<p->dnum;++i) pixs[n++]=a[i][j]; \
+ \
+ /* If there are any elements, measure the */ \
+ if(n) \
+ { \
+ /* Calculate the quantile and put it in the output. */ \
+ quantile=gal_statistics_quantile(cont, p->p1, 1); \
+ memcpy(&o[j], quantile->array, \
+ gal_type_sizeof(p->list->type)); \
+ gal_data_free(quantile); \
+ \
+ /* Since we are doing sigma-clipping in place, the size, */ \
+ /* and flags need to be reset. */ \
+ cont->flag=0; \
+ cont->size=cont->dsize[0]=p->dnum; \
+ } \
+ else \
+ o[j]=b; \
+ } \
+ \
+ /* Clean up. */ \
+ gal_data_free(cont); \
+ }
+
+
+
+
+
#define MULTIOPERAND_SIGCLIP(TYPE) { \
size_t n, j; \
gal_data_t *sclip; \
@@ -955,6 +1000,7 @@ struct multioperandparams
/* If there are any elements, measure the */ \
if(n) \
{ \
+ /* Calculate the sigma-clip and write it in. */ \
sclip=gal_statistics_sigma_clip(cont, p->p1, p->p2, 1, 1); \
sarr=sclip->array; \
switch(p->operator) \
@@ -968,6 +1014,7 @@ struct multioperandparams
"valid for sigma-clipping results", __func__, \
p->operator); \
} \
+ gal_data_free(sclip); \
\
/* Since we are doing sigma-clipping in place, the size, */ \
/* and flags need to be reset. */ \
@@ -1036,6 +1083,10 @@ struct multioperandparams
MULTIOPERAND_MEDIAN(TYPE, QSORT_F); \
break; \
\
+ case GAL_ARITHMETIC_OP_QUANTILE: \
+ MULTIOPERAND_QUANTILE(TYPE); \
+ break; \
+ \
case GAL_ARITHMETIC_OP_SIGCLIP_STD: \
case GAL_ARITHMETIC_OP_SIGCLIP_MEAN: \
case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN: \
@@ -1144,6 +1195,17 @@ arithmetic_multioperand(int operator, int flags,
gal_data_t *list,
/* Write them */
if(isnan(p1)) p1=((float *)(tmp->array))[0];
else p2=((float *)(tmp->array))[0];
+
+ /* Operator specific, parameter sanity checks. */
+ switch(operator)
+ {
+ case GAL_ARITHMETIC_OP_QUANTILE:
+ if(p1<0 || p1>1)
+ error(EXIT_FAILURE, 0, "%s: the parameter given to the `quantile' "
+ "operator must be between (and including) 0 and 1. The "
+ "given value is: %g", __func__, p1);
+ break;
+ }
}
@@ -1178,6 +1240,7 @@ arithmetic_multioperand(int operator, int flags,
gal_data_t *list,
case GAL_ARITHMETIC_OP_MEAN: otype=GAL_TYPE_FLOAT32; break;
case GAL_ARITHMETIC_OP_STD: otype=GAL_TYPE_FLOAT32; break;
case GAL_ARITHMETIC_OP_MEDIAN: otype=GAL_TYPE_FLOAT32; break;
+ case GAL_ARITHMETIC_OP_QUANTILE: otype=list->type; break;
case GAL_ARITHMETIC_OP_SIGCLIP_STD: otype=GAL_TYPE_FLOAT32; break;
case GAL_ARITHMETIC_OP_SIGCLIP_MEAN: otype=GAL_TYPE_FLOAT32; break;
case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN: otype=GAL_TYPE_FLOAT32; break;
@@ -1620,6 +1683,8 @@ gal_arithmetic_set_operator(char *string, size_t
*num_operands)
{ op=GAL_ARITHMETIC_OP_STD; *num_operands=-1; }
else if (!strcmp(string, "median"))
{ op=GAL_ARITHMETIC_OP_MEDIAN; *num_operands=-1; }
+ else if (!strcmp(string, "quantile"))
+ { op=GAL_ARITHMETIC_OP_QUANTILE; *num_operands=-1; }
else if (!strcmp(string, "sigclip-number"))
{ op=GAL_ARITHMETIC_OP_SIGCLIP_NUMBER; *num_operands=-1; }
else if (!strcmp(string, "sigclip-mean"))
@@ -1750,6 +1815,7 @@ gal_arithmetic_operator_string(int operator)
case GAL_ARITHMETIC_OP_MEAN: return "mean";
case GAL_ARITHMETIC_OP_STD: return "std";
case GAL_ARITHMETIC_OP_MEDIAN: return "median";
+ case GAL_ARITHMETIC_OP_QUANTILE: return "quantile";
case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER: return "sigclip-number";
case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN: return "sigclip-median";
case GAL_ARITHMETIC_OP_SIGCLIP_MEAN: return "sigclip-mean";
@@ -1860,6 +1926,7 @@ gal_arithmetic(int operator, size_t numthreads, int
flags, ...)
case GAL_ARITHMETIC_OP_MEAN:
case GAL_ARITHMETIC_OP_STD:
case GAL_ARITHMETIC_OP_MEDIAN:
+ case GAL_ARITHMETIC_OP_QUANTILE:
case GAL_ARITHMETIC_OP_SIGCLIP_STD:
case GAL_ARITHMETIC_OP_SIGCLIP_MEAN:
case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN:
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index 7a34d0d..b3c61eb 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -121,6 +121,7 @@ enum gal_arithmetic_operators
GAL_ARITHMETIC_OP_MEAN, /* Mean per pixel of multiple arrays. */
GAL_ARITHMETIC_OP_STD, /* STD per pixel of multiple arrays. */
GAL_ARITHMETIC_OP_MEDIAN, /* Median per pixel of multiple arrays. */
+ GAL_ARITHMETIC_OP_QUANTILE, /* Quantile per pixel of multiple arrays.*/
GAL_ARITHMETIC_OP_SIGCLIP_NUMBER,/* Sigma-clipped number of mult. arrays.*/
GAL_ARITHMETIC_OP_SIGCLIP_MEAN, /* Sigma-clipped mean of multiple arrays.*/
GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN,/* Sigma-clipped median of mult. arrays.*/