[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 3c6bdfc: Statistics program and library: new 2
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 3c6bdfc: Statistics program and library: new 2D histogram feature added |
|
Date: |
Fri, 14 Aug 2020 00:55:39 -0400 (EDT) |
branch: master
commit 3c6bdfcb0dbef8f92988fc2f87cea8f917d4a65b
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Statistics program and library: new 2D histogram feature added
In many contexts, a 2D histogram is necessary to study the behavoir of two
features in a large distribution and a scatter plot is usually not enough.
With this commit, the 'statistics.h' library now has the new
'gal_statistics_histogram2d' function and it is used in the Statistics
program to get the necessary columns to draw a 2D histogram. A new section
has also been added in the book, describing 2D histograms that also
contains a fully working minimal LaTeX source to draw the result.
---
NEWS | 13 +++
bin/statistics/args.h | 71 ++++++++++++++-
bin/statistics/aststatistics.conf | 1 +
bin/statistics/main.h | 7 +-
bin/statistics/statistics.c | 76 +++++++++++++---
bin/statistics/ui.c | 179 +++++++++++++++++++++++++++++++-------
bin/statistics/ui.h | 5 ++
doc/gnuastro.texi | 112 ++++++++++++++++++++++--
lib/gnuastro/statistics.h | 3 +
lib/statistics.c | 142 +++++++++++++++++++++++++++++-
10 files changed, 554 insertions(+), 55 deletions(-)
diff --git a/NEWS b/NEWS
index 787bded..15517c6 100644
--- a/NEWS
+++ b/NEWS
@@ -40,6 +40,18 @@ See the end of the file for license conditions.
distortion, you can simply run 'astfits --wcsdistortion=SIP' on the
file. The inverse conversion is also supported (from SIP to TPV).
+ Statistics
+ - New feature to create a 2D-histogram using two input columns (useful
+ when you have lots of points that are too dense and may hide important
+ features). This mode can be activated with the new '--histogram2d'
+ option. The binning of the first (X axis) column is specified with the
+ same 1D histogram options. The second column's binning is configured
+ with the following options:
+ --numbins2: Number of bins along the second column.
+ --greaterequal2: Only second column points that are larger than this.
+ --lessthan2; Only second column points that are less than this.
+ --onebinstart2: Make sure one bin starts at the value given here.
+
Table:
- New '--catcolumns' to specify which columns to concatenate (or append)
to the output. You can specify the file name containing the columns to
@@ -67,6 +79,7 @@ See the end of the file for license conditions.
- GAL_INTERPOLATE_NEIGHBORS_FUNC_INVALID: for error-handling/completeness.
- GAL_CONFIG_HAVE_WCSLIB_DIS_H: if the host's WCSLIB supports distortions.
- GAL_CONFIG_HAVE_WCSLIB_MJDREF: if the host's WCSLIB reads MJDREF keyword.
+ - gal_statistics_histogram2d: Generate 2D histogram from two columns.
- GAL_WCS_FLTERROR: Limit to identify a floating point error for WCS.
- gal_wcs_write: Write the given WCS into a FITS extension with no data.
- gal_wcs_clean_errors: clean major WCS components from errors specified
diff --git a/bin/statistics/args.h b/bin/statistics/args.h
index 252d390..710a753 100644
--- a/bin/statistics/args.h
+++ b/bin/statistics/args.h
@@ -36,10 +36,10 @@ struct argp_option program_options[] =
UI_KEY_COLUMN,
"STR",
0,
- "Column name or number if input is a table.",
+ "Column number (counting from 1) or search string.",
GAL_OPTIONS_GROUP_INPUT,
- &p->column,
- GAL_TYPE_STRING,
+ &p->columns,
+ GAL_TYPE_STRLL,
GAL_OPTIONS_RANGE_ANY,
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET
@@ -424,6 +424,19 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
+ "histogram2d",
+ UI_KEY_HISTOGRAM2D,
+ 0,
+ 0,
+ "Save the 2D histogram in output.",
+ UI_GROUP_PARTICULAR_STAT,
+ &p->histogram2d,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
"mirror",
UI_KEY_MIRROR,
"FLT",
@@ -655,6 +668,19 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
+ "numbins2",
+ UI_KEY_NUMBINS2,
+ "INT",
+ 0,
+ "No. of bins in second-dim of 2D histogram.",
+ UI_GROUP_HIST_CFP,
+ &p->numbins2,
+ GAL_TYPE_SIZE_T,
+ GAL_OPTIONS_RANGE_GT_0,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
"numasciibins",
UI_KEY_NUMASCIIBINS,
"INT",
@@ -732,6 +758,45 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET
},
+ {
+ "greaterequal2",
+ UI_KEY_GREATEREQUAL2,
+ "FLT",
+ 0,
+ "Upper range of 2nd dim in 2D histograms.",
+ GAL_OPTIONS_GROUP_INPUT,
+ &p->greaterequal2,
+ GAL_TYPE_FLOAT32,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "lessthan2",
+ UI_KEY_LESSTHAN2,
+ "FLT",
+ 0,
+ "Lower range of 2nd dim in 2D histograms.",
+ GAL_OPTIONS_GROUP_INPUT,
+ &p->lessthan2,
+ GAL_TYPE_FLOAT32,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "onebinstart2",
+ UI_KEY_ONEBINSTART2,
+ "FLT",
+ 0,
+ "Similar to --onebinstart, for 2D histogram",
+ UI_GROUP_HIST_CFP,
+ &p->onebinstart2,
+ GAL_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 4838c98..ee26c09 100644
--- a/bin/statistics/aststatistics.conf
+++ b/bin/statistics/aststatistics.conf
@@ -33,4 +33,5 @@
numasciibins 70
asciiheight 10
numbins 100
+ numbins2 20
mirrordist 1.5
diff --git a/bin/statistics/main.h b/bin/statistics/main.h
index fadcaf2..488ae36 100644
--- a/bin/statistics/main.h
+++ b/bin/statistics/main.h
@@ -58,10 +58,12 @@ struct statisticsparams
gal_list_i32_t *singlevalue; /* Single value calculations. */
gal_list_f64_t *tp_args; /* Arguments for printing. */
char *inputname; /* Input filename. */
- char *column; /* Column name or number if input is table. */
+ gal_list_str_t *columns; /* Column name or number if input is table. */
char *refcol; /* Reference column name or number. */
float greaterequal; /* Only use values >= this value. */
+ float greaterequal2; /* Only use values >= this value (2D hist). */
float lessthan; /* Only use values < this value. */
+ float lessthan2; /* Only use values < this value (2D hist). */
float quantmin; /* Quantile min or range: from Q to 1-Q. */
float quantmax; /* Quantile maximum. */
uint8_t ontile; /* Do single value calculations on tiles. */
@@ -70,6 +72,7 @@ struct statisticsparams
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 histogram2d; /* Save 2D-histogram in output. */
uint8_t cumulative; /* Save cumulative distibution in output. */
double mirror; /* Mirror value for hist and CFP. */
uint8_t sky; /* Find the Sky value over the image. */
@@ -77,11 +80,13 @@ struct statisticsparams
gal_data_t *contour; /* Levels to show contours. */
size_t numbins; /* Number of bins in histogram or CFP. */
+ size_t numbins2; /* No. of second-dim bins in 2D histogram. */
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. */
uint8_t manualbinrange; /* Set bin min/max manually, not from data. */
float onebinstart; /* Shift bins to start at this value. */
+ float onebinstart2; /* Shift bins to start at this value. */
uint8_t maxbinone; /* Set the maximum bin to 1. */
float mirrordist; /* Maximum distance after mirror for mode. */
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index 9e5f6b7..5fbc03a 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -506,7 +506,7 @@ print_ascii_plot(struct statisticsparams *p, gal_data_t
*plot,
/* Data structure that must be fed into 'gal_statistics_regular_bins'.*/
static gal_data_t *
-set_bin_range_params(struct statisticsparams *p)
+set_bin_range_params(struct statisticsparams *p, size_t dim)
{
size_t rsize=2;
gal_data_t *range=NULL;
@@ -516,8 +516,21 @@ set_bin_range_params(struct statisticsparams *p)
/* Allocate the range data structure. */
range=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &rsize, NULL, 0, -1, 1,
NULL, NULL, NULL);
- ((float *)(range->array))[0]=p->greaterequal;
- ((float *)(range->array))[1]=p->lessthan;
+ switch(dim)
+ {
+ case 1:
+ ((float *)(range->array))[0]=p->greaterequal;
+ ((float *)(range->array))[1]=p->lessthan;
+ break;
+ case 2:
+ ((float *)(range->array))[0]=p->greaterequal2;
+ ((float *)(range->array))[1]=p->lessthan2;
+ break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "address the problem. The value %zu for 'dim' isn't "
+ "recogized", __func__, PACKAGE_BUGREPORT, dim);
+ }
}
return range;
}
@@ -532,7 +545,7 @@ ascii_plots(struct statisticsparams *p)
gal_data_t *bins, *hist, *cfp=NULL, *range=NULL;
/* Make the bins and the respective plot. */
- range=set_bin_range_params(p);
+ range=set_bin_range_params(p, 1);
bins=gal_statistics_regular_bins(p->input, range, p->numasciibins, NAN);
hist=gal_statistics_histogram(p->input, bins, 0, 0);
if(p->asciicfp)
@@ -616,7 +629,7 @@ write_output_table(struct statisticsparams *p, gal_data_t
*table,
gal_list_str_add(&comments, tmp, 0);
if(strcmp(fix, "fits")) /* The intro info will be in FITS files anyway.*/
- gal_table_comments_add_intro(&comments, PROGRAM_NAME, &p->rawtime);
+ gal_table_comments_add_intro(&comments, PROGRAM_STRING, &p->rawtime);
/* Write the table. */
@@ -657,7 +670,7 @@ save_hist_and_or_cfp(struct statisticsparams *p)
/* 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). */
- range=set_bin_range_params(p);
+ range=set_bin_range_params(p, 1);
bins=gal_statistics_regular_bins(p->input, range, p->numbins,
p->onebinstart);
hist=gal_statistics_histogram(p->input, bins, p->normalize, p->maxbinone);
@@ -693,11 +706,11 @@ save_hist_and_or_cfp(struct statisticsparams *p)
/* Prepare the contents. */
if(p->histogram && p->cumulative)
- { suf="_hist_cfp"; contents="Histogram and cumulative frequency plot"; }
+ { suf="-hist-cfp"; contents="Histogram and cumulative frequency plot"; }
else if(p->histogram)
- { suf="_hist"; contents="Histogram"; }
+ { suf="-hist"; contents="Histogram"; }
else
- { suf="_cfp"; contents="Cumulative frequency plot"; }
+ { suf="-cfp"; contents="Cumulative frequency plot"; }
/* Set the output file name. */
@@ -711,6 +724,38 @@ save_hist_and_or_cfp(struct statisticsparams *p)
+static void
+histogram_2d(struct statisticsparams *p)
+{
+ char *suf, *contents;
+ gal_data_t *hist2d, *bins;
+ gal_data_t *range1, *range2;
+
+ /* Set the bins for each dimension */
+ range1=set_bin_range_params(p, 1);
+ range2=set_bin_range_params(p, 2);
+ bins=gal_statistics_regular_bins(p->input, range1, p->numbins,
+ p->onebinstart);
+ bins->next=gal_statistics_regular_bins(p->input->next, range2,
+ p->numbins2,
+ p->onebinstart2);
+
+ /* Build the 2D histogram. */
+ hist2d=gal_statistics_histogram2d(p->input, bins);
+
+ /* Write the output. */
+ suf="-hist2d";
+ contents="2D Histogram";
+ write_output_table(p, hist2d, suf, contents);
+
+ /* Clean up. */
+ gal_list_data_free(hist2d);
+}
+
+
+
+
+
void
print_mirror_hist_cfp(struct statisticsparams *p)
{
@@ -785,8 +830,8 @@ print_input_info(struct statisticsparams *p)
printf("Input: %s\n", name);
/* If a table was given, print the column. */
- if(p->column) printf("Column: %s\n",
- p->input->name ? p->input->name : p->column);
+ if(p->columns) printf("Column: %s\n",
+ p->input->name ? p->input->name : p->columns->v);
/* Range. */
str=NULL;
@@ -899,7 +944,7 @@ print_basics(struct statisticsparams *p)
range of the histogram. In that case, we want to print the histogram
information. */
printf("-------");
- range=set_bin_range_params(p);
+ range=set_bin_range_params(p, 1);
p->asciiheight = p->asciiheight ? p->asciiheight : 10;
p->numasciibins = p->numasciibins ? p->numasciibins : 70;
bins=gal_statistics_regular_bins(p->input, range, p->numasciibins, NAN);
@@ -1046,6 +1091,13 @@ statistics(struct statisticsparams *p)
save_hist_and_or_cfp(p);
}
+ /* 2D histogram. */
+ if(p->histogram2d)
+ {
+ print_basic_info=0;
+ histogram_2d(p);
+ }
+
/* Print the sigma-clipped results. */
if( p->sigmaclip )
{
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index 20ccf91..a35e727 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -127,8 +127,11 @@ ui_initialize_options(struct statisticsparams *p,
/* Program-specific initializers */
p->lessthan = NAN;
+ p->lessthan2 = NAN;
p->onebinstart = NAN;
+ p->onebinstart2 = NAN;
p->greaterequal = NAN;
+ p->greaterequal2 = NAN;
p->quantmin = NAN;
p->quantmax = NAN;
p->mirror = NAN;
@@ -397,8 +400,8 @@ ui_read_check_only_options(struct statisticsparams *p)
if( p->ontile || p->sky )
{
/* The tile or sky mode cannot be called with any other modes. */
- if(p->asciihist || p->asciicfp || p->histogram || p->cumulative
- || p->sigmaclip || !isnan(p->mirror) )
+ if( p->asciihist || p->asciicfp || p->histogram || p->histogram2d
+ || p->cumulative || p->sigmaclip || !isnan(p->mirror) )
error(EXIT_FAILURE, 0, "'--ontile' or '--sky' cannot be called with "
"any of the 'particular' calculation options, for example "
"'--histogram'. This is because the latter work over the whole "
@@ -495,10 +498,15 @@ ui_read_check_only_options(struct statisticsparams *p)
/* When binned outputs are requested, make sure that 'numbins' is set. */
- if( (p->histogram || p->cumulative || !isnan(p->mirror)) && p->numbins==0)
+ if( (p->histogram || p->histogram2d || p->cumulative || !isnan(p->mirror))
+ && 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( p->histogram2d && p->numbins2==0 )
+ error(EXIT_FAILURE, 0, "'--numbins2' isn't set. When a 2D histogram "
+ "is requested, the number of bins in the second dimension "
+ "('--numbins2') is also necessary");
/* If an ascii plot is requested, check if the ascii number of bins and
@@ -539,7 +547,7 @@ ui_check_options_and_arguments(struct statisticsparams *p)
/* If its an image, make sure column isn't given (in case the
user confuses an image with a table). */
p->hdu_type=gal_fits_hdu_format(p->inputname, p->cp.hdu);
- if(p->hdu_type==IMAGE_HDU && p->column)
+ if(p->hdu_type==IMAGE_HDU && p->columns)
error(EXIT_FAILURE, 0, "%s (hdu: %s): is a FITS image "
"extension. The '--column' option is only applicable "
"to tables.", p->inputname, p->cp.hdu);
@@ -574,8 +582,7 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
{
size_t one=1;
unsigned char flags=GAL_ARITHMETIC_NUMOK;
- unsigned char flagsor = ( GAL_ARITHMETIC_FREE
- | GAL_ARITHMETIC_INPLACE
+ unsigned char flagsor = ( GAL_ARITHMETIC_INPLACE
| GAL_ARITHMETIC_NUMOK );
gal_data_t *tmp, *cond_g=NULL, *cond_l=NULL, *cond, *blank, *ref;
@@ -650,13 +657,22 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
/* 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.*/
+ 'where' operator. Also reset the blank flags so they are checked again
+ if necessary.*/
gal_arithmetic(GAL_ARITHMETIC_OP_WHERE, 1, flagsor, p->input, cond, blank);
-
-
- /* Reset the blank flags so they are checked again if necessary. */
p->input->flag &= ~GAL_DATA_FLAG_BLANK_CH;
p->input->flag &= ~GAL_DATA_FLAG_HASBLANK;
+ if(p->input->next)
+ {
+ gal_arithmetic(GAL_ARITHMETIC_OP_WHERE, 1, flagsor, p->input->next,
+ cond, blank);
+ p->input->next->flag &= ~GAL_DATA_FLAG_BLANK_CH;
+ p->input->next->flag &= ~GAL_DATA_FLAG_HASBLANK;
+ }
+
+ /* Clean up. */
+ gal_data_free(cond);
+ gal_data_free(blank);
}
@@ -711,25 +727,82 @@ ui_make_sorted_if_necessary(struct statisticsparams *p)
+/* Merge all given columns into one list. This is because people may call
+ for different columns in one call to '--column' ('-c', separated by a
+ comma), or multiple calls to '-c'. */
+gal_list_str_t *
+ui_read_columns_in_one(gal_list_str_t *incolumns)
+{
+ size_t i;
+ gal_data_t *strs;
+ char *c, **strarr;
+ gal_list_str_t *tmp, *final=NULL;
+
+ /* Go over the separate calls to the '-c' option, see the explaination in
+ Table's 'ui_columns_prepare' function for more on every step. */
+ for(tmp=incolumns;tmp!=NULL;tmp=tmp->next)
+ {
+ /* Remove any new-line character. */
+ for(c=tmp->v;*c!='\0';++c)
+ if(*c=='\\' && *(c+1)=='\n') { *c=' '; *(++c)=' '; }
+
+ /* Read the different comma-separated strings into an array (within a
+ 'gal_data_t'). */
+ strs=gal_options_parse_csv_strings_raw(tmp->v, NULL, 0);
+ strarr=strs->array;
+
+ /* Add each array element to the final list of columns. */
+ for(i=0;i<strs->size;++i)
+ gal_list_str_add(&final, strarr[i], 1);
+ }
+
+ /* Reverse the list to be in the same order. */
+ gal_list_str_reverse(&final);
+
+ /* For a check.
+ for(tmp=final;tmp!=NULL;tmp=tmp->next)
+ printf("%s\n", tmp->v);
+ */
+
+ /* Return the final unified list. */
+ return final;
+}
+
+
+
+
+
void
ui_read_columns(struct statisticsparams *p)
{
int toomanycols=0, tformat;
- gal_list_str_t *column=NULL;
gal_data_t *cols, *tmp, *cinfo;
size_t size, ncols, nrows, counter=0;
+ gal_list_str_t *incols, *columnlist=NULL;
gal_list_str_t *lines=gal_options_check_stdin(p->inputname,
p->cp.stdintimeout, "input");
+ /* Merge possibly multiple calls to '-c' (each possibly separated by a
+ coma) into one list. */
+ incols=ui_read_columns_in_one(p->columns);
+
+ /* If any columns are specified, make sure there is a maximum of two
+ columns. */
+ if(gal_list_str_number(incols)>2)
+ error(EXIT_FAILURE, 0, "%zu input columns were given but currently a "
+ "maximum of two columns are supported (two columns only for "
+ "special operations, the majority of operations are on a single "
+ "column)", gal_list_str_number(incols));
+
/* If a reference column is also given, add it to the list of columns to
read. */
if(p->refcol)
- gal_list_str_add(&column, p->refcol, 0);
+ gal_list_str_add(&columnlist, p->refcol, 0);
/* If no column is specified, Statistics will abort and an error will be
printed when the table has more than one column. If there is only one
column, there is no need to specify any, so Statistics will use it. */
- if(p->column==NULL)
+ if(incols==NULL)
{
/* Get the basic table information. */
cinfo=gal_table_info(p->inputname, p->cp.hdu, lines, &ncols, &nrows,
@@ -740,12 +813,12 @@ ui_read_columns(struct statisticsparams *p)
switch(ncols)
{
case 0:
- error(EXIT_FAILURE, 0, "%s contains no usable information",
+ error(EXIT_FAILURE, 0, "%s contains no usable columns",
( p->inputname
? gal_checkset_dataset_name(p->inputname, p->cp.hdu)
: "Standard input" ));
case 1:
- gal_checkset_allocate_copy("1", &p->column);
+ gal_list_str_add(&incols, "1", 1);
break;
default:
error(EXIT_FAILURE, 0, "%s is a table containing more than one "
@@ -765,18 +838,19 @@ ui_read_columns(struct statisticsparams *p)
}
}
- gal_list_str_add(&column, p->column, 0);
+ if(columnlist) columnlist->next=incols;
+ else columnlist=incols;
/* Read the desired column(s). */
- cols=gal_table_read(p->inputname, p->cp.hdu, lines, column, p->cp.searchin,
- p->cp.ignorecase, p->cp.minmapsize, p->cp.quietmmap,
- NULL);
+ cols=gal_table_read(p->inputname, p->cp.hdu, lines, columnlist,
+ p->cp.searchin, p->cp.ignorecase, p->cp.minmapsize,
+ p->cp.quietmmap, NULL);
gal_list_str_free(lines, 1);
- /* If the input was from standard input, we can actually write this into
- it (for future reporting). */
+ /* If the input was from standard input, we'll set the input name to be
+ 'stdin' (for future reporting). */
if(p->inputname==NULL)
- gal_checkset_allocate_copy("statistics", &p->inputname);
+ gal_checkset_allocate_copy("stdin", &p->inputname);
/* Put the columns into the proper gal_data_t. */
size=cols->size;
@@ -806,8 +880,15 @@ ui_read_columns(struct statisticsparams *p)
/* Put the column into the proper pointer. */
switch(++counter)
{
- case 1: p->input=tmp; break;
- case 2: if(p->refcol) p->reference=tmp; else toomanycols=1; break;
+ case 1: p->input=tmp; break;
+ case 2:
+ if(gal_list_str_number(incols)==2)
+ p->input->next=tmp;
+ else
+ if(p->refcol) p->reference=tmp; else toomanycols=1;
+ break;
+ case 3:
+ if(p->refcol) p->reference=tmp; else toomanycols=1; break;
default: toomanycols=1;
}
@@ -822,7 +903,17 @@ ui_read_columns(struct statisticsparams *p)
}
/* Clean up. */
- gal_list_str_free(column, 0);
+ gal_list_str_free(incols, 1);
+ if(columnlist!=incols)
+ {
+ columnlist->next=NULL;
+ gal_list_str_free(columnlist, 0);
+ }
+
+ /* For a check:
+ printf("Number of input columns: %zu\n", gal_list_data_number(p->input));
+ exit(0);
+ */
}
@@ -832,8 +923,11 @@ ui_read_columns(struct statisticsparams *p)
void
ui_preparations(struct statisticsparams *p)
{
- gal_data_t *check;
+ unsigned char flagsor = ( GAL_ARITHMETIC_FREE
+ | GAL_ARITHMETIC_INPLACE
+ | GAL_ARITHMETIC_NUMOK );
int keepinputdir=p->cp.keepinputdir;
+ gal_data_t *check, *flag1, *flag2, *flag;
struct gal_options_common_params *cp=&p->cp;
struct gal_tile_two_layer_params *tl=&cp->tl;
char *checkbasename = p->cp.output ? p->cp.output : p->inputname;
@@ -855,8 +949,14 @@ ui_preparations(struct statisticsparams *p)
}
else
{
+ /* Read the table columns. */
ui_read_columns(p);
p->inputformat=INPUT_FORMAT_TABLE;
+
+ /* Two columns can only be given with 2D histogram mode. */
+ if(p->histogram2d==0 && p->input->next!=NULL)
+ error(EXIT_FAILURE, 0, "two column input is currently only "
+ "supported for 2D histogram mode");
}
/* Read the convolution kernel if necessary. */
@@ -906,17 +1006,34 @@ ui_preparations(struct statisticsparams *p)
/* If we are not to work on tiles, then re-order and change the input. */
if(p->ontile==0 && p->sky==0 && p->contour==NULL)
{
- /* Only keep the elements we want. */
- gal_blank_remove(p->input);
+ /* Only keep the elements we want. Note that if we have more than one
+ column, we need to move the same rows in both (otherwise their
+ widths won't be equal). */
+ if(p->input->next)
+ {
+ flag1=gal_blank_flag(p->input);
+ flag2=gal_blank_flag(p->input->next);
+ flag=gal_arithmetic(GAL_ARITHMETIC_OP_OR, 1, flagsor, flag1, flag2);
+
+ gal_blank_flag_apply(p->input, flag);
+ gal_blank_flag_apply(p->input->next, flag);
+
+ gal_blank_remove(p->input);
+ gal_blank_remove(p->input->next);
+
+ p->input->next->flag &= ~GAL_DATA_FLAG_HASBLANK ;
+ p->input->next->flag |= GAL_DATA_FLAG_BLANK_CH ;
+ }
+ else
+ gal_blank_remove(p->input);
+ p->input->flag &= ~GAL_DATA_FLAG_HASBLANK ;
+ p->input->flag |= GAL_DATA_FLAG_BLANK_CH ;
/* Make sure there actually are any (non-blank) elements left. */
if(p->input->size==0)
error(EXIT_FAILURE, 0, "%s: all elements are blank",
gal_fits_name_save_as_string(p->inputname, cp->hdu));
- p->input->flag &= ~GAL_DATA_FLAG_HASBLANK ;
- p->input->flag |= GAL_DATA_FLAG_BLANK_CH ;
-
/* Make sure there is data remaining: */
if(p->input->size==0)
error(EXIT_FAILURE, 0, "%s: no data, maybe the '--greaterequal' or "
diff --git a/bin/statistics/ui.h b/bin/statistics/ui.h
index 2e0ea6a..f257480 100644
--- a/bin/statistics/ui.h
+++ b/bin/statistics/ui.h
@@ -83,8 +83,10 @@ enum option_keys_enum
UI_KEY_MODESYMVALUE,
UI_KEY_QUANTFUNC,
UI_KEY_ASCIICFP,
+ UI_KEY_HISTOGRAM2D,
UI_KEY_MIRROR,
UI_KEY_NUMBINS,
+ UI_KEY_NUMBINS2,
UI_KEY_NUMASCIIBINS,
UI_KEY_ASCIIHEIGHT,
UI_KEY_LOWERBIN,
@@ -105,6 +107,9 @@ enum option_keys_enum
UI_KEY_SIGCLIPMEDIAN,
UI_KEY_SIGCLIPMEAN,
UI_KEY_SIGCLIPSTD,
+ UI_KEY_GREATEREQUAL2,
+ UI_KEY_LESSTHAN2,
+ UI_KEY_ONEBINSTART2,
};
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 9a694fd..294c3db 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -259,8 +259,8 @@ General program usage tutorial
* Aperture photometry:: Doing photometry on a fixed aperture.
* Matching catalogs:: Easily find corresponding rows from two
catalogs.
* Finding reddest clumps and visual inspection:: Selecting some targets and
inspecting them.
-* Citing and acknowledging Gnuastro:: How to cite and acknowledge Gnuastro in
your papers.
* Writing scripts to automate the steps:: Scripts will greatly help in
re-doing things fast.
+* Citing and acknowledging Gnuastro:: How to cite and acknowledge Gnuastro in
your papers.
Detecting large extended targets
@@ -471,6 +471,7 @@ Data analysis
Statistics
* Histogram and Cumulative Frequency Plot:: Basic definitions.
+* 2D Histograms:: Plotting the distribution of two variables.
* Sigma clipping:: Definition of @mymath{\sigma}-clipping.
* Sky value:: Definition and derivation of the Sky value.
* Invoking aststatistics:: Arguments and options to Statistics.
@@ -1947,8 +1948,8 @@ This will help simulate future situations when you are
processing your own datas
* Aperture photometry:: Doing photometry on a fixed aperture.
* Matching catalogs:: Easily find corresponding rows from two
catalogs.
* Finding reddest clumps and visual inspection:: Selecting some targets and
inspecting them.
-* Citing and acknowledging Gnuastro:: How to cite and acknowledge Gnuastro in
your papers.
* Writing scripts to automate the steps:: Scripts will greatly help in
re-doing things fast.
+* Citing and acknowledging Gnuastro:: How to cite and acknowledge Gnuastro in
your papers.
@end menu
@node Calling Gnuastro's programs, Accessing documentation, General program
usage tutorial, General program usage tutorial
@@ -3414,7 +3415,7 @@ $ astfits matched.fits
-@node Finding reddest clumps and visual inspection, Citing and acknowledging
Gnuastro, Matching catalogs, General program usage tutorial
+@node Finding reddest clumps and visual inspection, Writing scripts to
automate the steps, Matching catalogs, General program usage tutorial
@subsection Finding reddest clumps and visual inspection
@cindex GNU AWK
As a final step, let's go back to the original clumps-based color measurement
we generated in @ref{Working with catalogs estimating colors}.
@@ -3516,7 +3517,7 @@ $ ds9 -mecube seg/xdf-f160w.fits -zscale -zoom to fit \
-regions load all reddest.reg
@end example
-@node Writing scripts to automate the steps, , Citing and acknowledging
Gnuastro, General program usage tutorial
+@node Writing scripts to automate the steps, Citing and acknowledging
Gnuastro, Finding reddest clumps and visual inspection, General program usage
tutorial
@subsection Writing scripts to automate the steps
In the previous sub-sections, we went through a series of steps like
downloading the necessary datasets (in @ref{Setup and data download}),
detecting the objects in the image, and finally selecting a particular subset
of them to inspect visually (in @ref{Finding reddest clumps and visual
inspection}).
@@ -3830,7 +3831,7 @@ if ! [ -f $f160w_flat ]; then
fi
@end example
-@node Citing and acknowledging Gnuastro, Writing scripts to automate the
steps, Finding reddest clumps and visual inspection, General program usage
tutorial
+@node Citing and acknowledging Gnuastro, , Writing scripts to automate the
steps, General program usage tutorial
@subsection Citing and acknowledging Gnuastro
In conclusion, we hope this extended tutorial has been a good starting point
to help in your exciting research.
If this book or any of the programs in Gnuastro have been useful for your
research, please cite the respective papers, and acknowledge the funding
agencies that made all of this possible.
@@ -12701,13 +12702,14 @@ The Statistics program is designed for such
situations.
@menu
* Histogram and Cumulative Frequency Plot:: Basic definitions.
+* 2D Histograms:: Plotting the distribution of two variables.
* Sigma clipping:: Definition of @mymath{\sigma}-clipping.
* Sky value:: Definition and derivation of the Sky value.
* Invoking aststatistics:: Arguments and options to Statistics.
@end menu
-@node Histogram and Cumulative Frequency Plot, Sigma clipping, Statistics,
Statistics
+@node Histogram and Cumulative Frequency Plot, 2D Histograms, Statistics,
Statistics
@subsection Histogram and Cumulative Frequency Plot
@cindex Histogram
@@ -12746,8 +12748,87 @@ 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 @option{--qrange} 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 2D Histograms, Sigma clipping, Histogram and Cumulative Frequency Plot,
Statistics
+@subsection 2D Histograms
+
+In @ref{Histogram and Cumulative Frequency Plot} the concept of histograms
were introduced on a single dataset.
+However, especially when doing high-level science on tables, the distribution
in a 2D space may be of interest (for example a color-magnitude diagram).
+But the number of points may be too large for a simple scatter plot to show
the concentration of the points: they will all fall over each other and just
make a large connected region that will hide potentially interesting behaviors.
+This is where 2D histograms can become very useful.
+The desired 2D region is broken up into 2D bins (boxes) and the number of
points falling within each box is returned.
+Added with a color-bar, you can now clearly see the distribution.
+
+Gnuastro's Statistics program has the @option{--histogram2d} option for this
task.
+Its output will be three columns that have the centers of every box in both
dimensions.
+The first column is the central box coordinates in the first dimension, the
second has values along the second dimension and the third has the number of
input points that fall within each box.
+You can specify the number of bins along each dimension through the
@option{--numbins} (for first input column) an @option{--numbins2} (for second
input column).
+the output file from this command can then be given to any plotting tool to
visualize the distribution.
+
+You can then make high-quality plots within your paper (all in @LaTeX{}), you
can use @url{https://ctan.org/pkg/pgfplots, PGFPlots} with the following
minimal @LaTeX{} source:
+
+@example
+%% Replace 'XXXXXXXXX' with your selected number of bins in the first
+%% dimension.
+%%
+%% Then run these commands to build the plot in a LaTeX command.
+%% mkdir tikz
+%% pdflatex -shell-escape -halt-on-error plot.tex
+\documentclass@{article@}
+
+%% Load PGFPlots and set it to build the figure separately in a 'tikz'
+%% directory (which has to exist before LaTeX is run). This
+%% "externalization" is very useful to include the commands of multiple
+%% plots in the middle of your paper/report, but also have the plots
+%% separately to use in slides or other scenarios.
+\usepackage@{pgfplots@}
+\usetikzlibrary@{external@}
+\tikzexternalize
+\tikzsetexternalprefix@{tikz/@}
+
+%% Start the document
+\begin@{document@}
+
+ You can actually write a full paper here and include many figures!
+ Feel free to change this text.
+
+ %% Define the colormap.
+ \pgfplotsset@{
+ /pgfplots/colormap=@{coldredux@}@{
+ [1cm]
+ rgb255(0cm)=(255,255,255)
+ rgb255(2cm)=(0,192,255)
+ rgb255(4cm)=(0,0,255)
+ rgb255(6cm)=(0,0,0)
+ @}
+ @}
-@node Sigma clipping, Sky value, Histogram and Cumulative Frequency Plot,
Statistics
+ %% Draw the plot.
+ \begin@{tikzpicture@}
+ \small
+ \begin@{axis@}[
+ width=\linewidth,
+ view=@{0@}@{90@},
+ colorbar horizontal,
+ xlabel=X axis,
+ ylabel=Y axis,
+ ylabel shift=-0.1cm,
+ colorbar style=@{at=@{(0,1.01)@}, anchor=south west,
+ xticklabel pos=upper@},
+ ]
+ \addplot3[
+ surf,
+ shader=flat corner,
+ mesh/ordering=rowwise,
+ mesh/cols=XXXXXXXXX, %% <---- Correct this!
+ ] file @{cm-hist2d.txt@};
+
+ \end@{axis@}
+\end@{tikzpicture@}
+
+\end@{document@}
+@end example
+
+@node Sigma clipping, Sky value, 2D Histograms, Statistics
@subsection Sigma clipping
Let's assume that you have pure noise (centered on zero) with a clear
@url{https://en.wikipedia.org/wiki/Normal_distribution,Gaussian distribution},
or see @ref{Photon counting noise}.
@@ -13405,6 +13486,11 @@ By default (when no @option{--output} is specified) a
plain text table will be c
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 @option{--asciihist} option).
+@item --histogram2d
+Save the 2D histogram of two input columns into an output file, see @ref{2D
Histograms}.
+The output will have three columns: the first two are the coordinates of each
box's center in the first and second dimensions/columns.
+The third will be number of input points that fall within that box.
+
@item -C
@itemx --cumulative
Save the cumulative frequency plot of the usable values in the input dataset
into a table, similar to @option{--histogram}.
@@ -13477,6 +13563,18 @@ The @option{--onebinstart} option has a higher
priority than this option.
In other words, @option{--onebinstart} takes effect after the range has been
finalized and the initial bins have been defined, therefore it has the power to
(possibly) shift the bins.
If you want to manually set the range of the bins @emph{and} have one bin on a
special value, it is thus better to avoid @option{--onebinstart}.
+@item --numbins2=INT
+Similar to @option{--numbins}, but for the second column when a 2D histogram
is requested, see @option{--histogram2d}.
+
+@item --greaterequal2=FLT
+Similar to @option{--greaterequal}, but for the second column when a 2D
histogram is requested, see @option{--histogram2d}.
+
+@item --lessthan2=FLT
+Similar to @option{--lessthan}, but for the second column when a 2D histogram
is requested, see @option{--histogram2d}.
+
+@item --onebinstart2=FLT
+Similar to @option{--onebinstart}, but for the second column when a 2D
histogram is requested, see @option{--histogram2d}.
+
@end table
diff --git a/lib/gnuastro/statistics.h b/lib/gnuastro/statistics.h
index 2b4ead5..675bfd0 100644
--- a/lib/gnuastro/statistics.h
+++ b/lib/gnuastro/statistics.h
@@ -163,6 +163,9 @@ gal_statistics_histogram(gal_data_t *data, gal_data_t *bins,
int normalize, int maxhistone);
gal_data_t *
+gal_statistics_histogram2d(gal_data_t *input, gal_data_t *bins);
+
+gal_data_t *
gal_statistics_cfp(gal_data_t *data, gal_data_t *bins, int normalize);
diff --git a/lib/statistics.c b/lib/statistics.c
index c151674..bafab53 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -1569,7 +1569,7 @@ gal_statistics_no_blank_sorted(gal_data_t *input, int
inplace)
* The 'inrange' data structure keeps the desired range along each
dimension of the input data structure, it has to be in float32
- type. Note that if
+ type. Note these points:
- If you want the full range of the dataset (in any dimensions,
then just set 'range' to NULL and the range will be specified
@@ -1852,6 +1852,146 @@ gal_statistics_histogram(gal_data_t *input, gal_data_t
*bins, int normalize,
+/* Build a 2D histogram from the two input columns (a list) and two bins
+ (also a list). */
+#define HISTOGRAM2D_TYPESET(AT, BT) { \
+ BT *b=input->next->array; \
+ AT *a=input->array, *af=a+input->size; \
+ do \
+ { \
+ if(*a>=mina && *a<=maxa && *b>=minb && *b<=maxb) \
+ { \
+ i=(*a-mina)/binwidtha; \
+ j=(*b-minb)/binwidthb; \
+ /* When '*a' is the largest element (within floating */ \
+ /* point errors), 'ii' can be one element larger than */ \
+ /* the number of bins. But since its in the dataset, we */ \
+ /* need to count it. So we'll put it in the last bin. */ \
+ if(i==bsizea) --i; \
+ if(j==bsizeb) --j; \
+ ++h[ i*bsizeb+j ]; \
+ } \
+ ++b; \
+ } \
+ while(++a<af); \
+ }
+
+#define HISTOGRAM2D_TYPESET_A(AT) { \
+ switch(input->next->type) \
+ { \
+ case GAL_TYPE_UINT8: HISTOGRAM2D_TYPESET(AT, uint8_t); break; \
+ case GAL_TYPE_INT8: HISTOGRAM2D_TYPESET(AT, int8_t); break; \
+ case GAL_TYPE_UINT16: HISTOGRAM2D_TYPESET(AT, uint16_t); break; \
+ case GAL_TYPE_INT16: HISTOGRAM2D_TYPESET(AT, int16_t); break; \
+ case GAL_TYPE_UINT32: HISTOGRAM2D_TYPESET(AT, uint32_t); break; \
+ case GAL_TYPE_INT32: HISTOGRAM2D_TYPESET(AT, int32_t); break; \
+ case GAL_TYPE_UINT64: HISTOGRAM2D_TYPESET(AT, uint64_t); break; \
+ case GAL_TYPE_INT64: HISTOGRAM2D_TYPESET(AT, int64_t); break; \
+ case GAL_TYPE_FLOAT32: HISTOGRAM2D_TYPESET(AT, float); break; \
+ case GAL_TYPE_FLOAT64: HISTOGRAM2D_TYPESET(AT, double); break; \
+ default: \
+ error(EXIT_FAILURE, 0, "%s: type code %d not recognized", \
+ __func__, input->type); \
+ } \
+ }
+
+gal_data_t *
+gal_statistics_histogram2d(gal_data_t *input, gal_data_t *bins)
+{
+ double *o1, *o2;
+ gal_data_t *tmp, *out;
+ size_t i, j, *h, bsizea, bsizeb, outsize;
+ double *da, *db, binwidtha, binwidthb, mina, minb, maxa, maxb;
+
+ /* Basic sanity checks */
+ if(input->next==NULL)
+ error(EXIT_FAILURE, 0, "%s: 'input' has to be a list of two datasets",
+ __func__);
+ if(bins->next==NULL)
+ error(EXIT_FAILURE, 0, "%s: 'bins' has to be a list of two datasets",
+ __func__);
+ if(input->next->next)
+ error(EXIT_FAILURE, 0, "%s: 'input' should only contain two datasets, "
+ "not more", __func__);
+ if(bins->next->next)
+ error(EXIT_FAILURE, 0, "%s: 'bins' should only contain two datasets, "
+ "not more", __func__);
+ if(input->size != input->next->size)
+ error(EXIT_FAILURE, 0, "the two input datasets have to have the "
+ "same size");
+ if(bins->status!=GAL_STATISTICS_BINS_REGULAR
+ || bins->next->status!=GAL_STATISTICS_BINS_REGULAR)
+ error(EXIT_FAILURE, 0, "%s: the input bins are not regular. Currently "
+ "it is only implemented for regular bins", __func__);
+
+ /* For easy reading of bin sizes. */
+ da=bins->array;
+ bsizea=bins->size;
+ db=bins->next->array;
+ bsizeb=bins->next->size;
+
+ /* Allocate the output. */
+ outsize=bsizea*bsizeb;
+ out=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &outsize,
+ NULL, 1, input->minmapsize, input->quietmmap,
+ "bin_dim1", input->unit,
+ "Bin centers along first axis.");
+ tmp=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &outsize,
+ NULL, 1, input->minmapsize, input->quietmmap,
+ "bin_dim2", input->next->unit,
+ "Bin centers along second axis.");
+ out->next=tmp;
+ tmp=gal_data_alloc(NULL, GAL_TYPE_SIZE_T, 1, &outsize,
+ NULL, 1, input->minmapsize, input->quietmmap,
+ "hist_number", "counts",
+ "Number of data points within each 2D-bin (box).");
+ out->next->next=tmp;
+
+ /* Fill in the first two output columns and set the histogram pointer. */
+ o1=out->array;
+ o2=out->next->array;
+ h=out->next->next->array;
+ for(i=0;i<bsizea;++i)
+ for(j=0;j<bsizeb;++j)
+ {
+ o1[i*bsizeb+j]=da[i];
+ o2[i*bsizeb+j]=db[j];
+ }
+
+ /* Set the minimum and maximum range of the histogram from the bins. */
+ binwidtha=da[1]-da[0];
+ binwidthb=db[1]-db[0];
+ mina=da[0]-binwidtha/2;
+ minb=db[0]-binwidthb/2;
+ maxa=da[ bins->size - 1 ] + binwidtha/2;
+ maxb=db[ bins->next->size - 1] + binwidthb/2;
+
+ /* Fill the histogram column. */
+ switch(input->type)
+ {
+ case GAL_TYPE_UINT8: HISTOGRAM2D_TYPESET_A(uint8_t); break;
+ case GAL_TYPE_INT8: HISTOGRAM2D_TYPESET_A(int8_t); break;
+ case GAL_TYPE_UINT16: HISTOGRAM2D_TYPESET_A(uint16_t); break;
+ case GAL_TYPE_INT16: HISTOGRAM2D_TYPESET_A(int16_t); break;
+ case GAL_TYPE_UINT32: HISTOGRAM2D_TYPESET_A(uint32_t); break;
+ case GAL_TYPE_INT32: HISTOGRAM2D_TYPESET_A(int32_t); break;
+ case GAL_TYPE_UINT64: HISTOGRAM2D_TYPESET_A(uint64_t); break;
+ case GAL_TYPE_INT64: HISTOGRAM2D_TYPESET_A(int64_t); break;
+ case GAL_TYPE_FLOAT32: HISTOGRAM2D_TYPESET_A(float); break;
+ case GAL_TYPE_FLOAT64: HISTOGRAM2D_TYPESET_A(double); break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
+ __func__, input->type);
+ }
+
+ /* Return the final output */
+ return out;
+}
+
+
+
+
+
/* Make a cumulative frequency plot (CFP) of all the elements in the given
dataset with bin values that are defined in the 'bins' structure (see
'gal_statistics_regular_bins').
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 3c6bdfc: Statistics program and library: new 2D histogram feature added,
Mohammad Akhlaghi <=