[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 29d7e46 091/125: Statistics now successfully u
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 29d7e46 091/125: Statistics now successfully uses gal_data_t |
|
Date: |
Sun, 23 Apr 2017 22:36:45 -0400 (EDT) |
branch: master
commit 29d7e46bf3f0202a58f59e84e8d197bd2672ae2c
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Statistics now successfully uses gal_data_t
The Statistics program now fully implements the previous calculations in a
much more robust and user friendly manner along with many new features. The
main features added with this commit are:
- Using a reference column with the `--refcol' option to select the
elements from the input dataset.
- The `--quantile' option which will print the values at the given
quantiles.
- The `--quantfunc' option which will print the quantile function, or
inverse of the `--quantile'. So you give it a value and it will say what
quantile that value has in the input distribution.
- The `--mode', `--modequant', `--modesym', and `--modesymvalue' options
as described in Appendix C of Akhlaghi and Ichikawa (2015, ApJS 220,
1). The new `--mirror' option is also now available to make the
histogram and cumulative frequency plot of the mirror distribution.
Many functions were also added to the statistics library and also the
Statistics program to make all of these possible.
---
bin/mkprof/ui.c | 2 +-
bin/noisechisel/thresh.c | 9 -
bin/statistics/args.h | 117 +++++-
bin/statistics/main.h | 8 +-
bin/statistics/statistics.c | 331 ++++++++++++---
bin/statistics/ui.c | 221 +++++++---
bin/statistics/ui.h | 9 +-
doc/gnuastro.texi | 217 +++++-----
lib/Makefile.am | 8 +-
lib/data.c | 28 ++
lib/gnuastro/data.h | 3 +
lib/gnuastro/linkedlist.h | 3 +
lib/gnuastro/statistics.h | 81 ++--
lib/linkedlist.c | 22 +
lib/mode.c | 402 ------------------
lib/mode.h | 46 ---
lib/options.c | 2 +-
lib/statistics.c | 963 +++++++++++++++++++++++++++++++-------------
18 files changed, 1442 insertions(+), 1030 deletions(-)
diff --git a/bin/mkprof/ui.c b/bin/mkprof/ui.c
index f763128..09cb801 100644
--- a/bin/mkprof/ui.c
+++ b/bin/mkprof/ui.c
@@ -395,7 +395,7 @@ ui_read_cols(struct mkprofparams *p)
the same as those that were wanted (it might be more). */
while(cols!=NULL)
{
- /* Pop out the top node. */
+ /* Pop out the top column. */
tmp=gal_data_pop_from_ll(&cols);
/* By default check if the column has blank values, but it can be
diff --git a/bin/noisechisel/thresh.c b/bin/noisechisel/thresh.c
index 71836df..b9b1f28 100644
--- a/bin/noisechisel/thresh.c
+++ b/bin/noisechisel/thresh.c
@@ -355,15 +355,6 @@ snthresh(struct noisechiselparams *p, float *sntable,
size_t size,
/* Sort the signal to noise ratios and remove their outliers */
qsort(sntable, size, sizeof *sntable, gal_qsort_float_increasing);
- /* The removal of outliers was useful when the S/N was calculated
- separately on each mesh (thus there were very few
- points). However, now that the S/N is calculated over the full
- image, the number of pseudo-detections and clumps is so high that
- these outliers will not play a significant role in the S/N
- threshold (unless you set unreasonably high quantiles!). So This
- is commented now. */
- /*gal_statistics_remove_outliers_flat_cdf(sntable, &size);*/
-
/* Store the SN value. */
sn=sntable[gal_statistics_index_from_quantile(size, quant)];
diff --git a/bin/statistics/args.h b/bin/statistics/args.h
index 7ae18c3..759230a 100644
--- a/bin/statistics/args.h
+++ b/bin/statistics/args.h
@@ -45,6 +45,19 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
+ "refcol",
+ ARGS_OPTION_KEY_REFCOL,
+ "STR",
+ 0,
+ "Reference column name or number.",
+ GAL_OPTIONS_GROUP_INPUT,
+ &p->refcol,
+ GAL_DATA_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
"greaterequal",
ARGS_OPTION_KEY_GREATEREQUAL,
"FLT",
@@ -99,10 +112,10 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_NUMBER,
0,
0,
- "Print the number of data-points.",
+ "Number (non-blank).",
ARGS_GROUP_IN_ONE_ROW,
&p->toprint,
- GAL_DATA_TYPE_INVALID,
+ GAL_OPTIONS_NO_ARG_TYPE,
GAL_OPTIONS_RANGE_ANY,
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET,
@@ -113,7 +126,7 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_MINIMUM,
0,
0,
- "Print the minimum.",
+ "Minimum.",
ARGS_GROUP_IN_ONE_ROW,
&p->toprint,
GAL_OPTIONS_NO_ARG_TYPE,
@@ -127,7 +140,7 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_MAXIMUM,
0,
0,
- "Print the maximum.",
+ "Maximum.",
ARGS_GROUP_IN_ONE_ROW,
&p->toprint,
GAL_OPTIONS_NO_ARG_TYPE,
@@ -141,7 +154,7 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_SUM,
0,
0,
- "Print the sum of all elements.",
+ "Sum.",
ARGS_GROUP_IN_ONE_ROW,
&p->toprint,
GAL_OPTIONS_NO_ARG_TYPE,
@@ -155,7 +168,7 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_MEAN,
0,
0,
- "Print the mean.",
+ "Mean.",
ARGS_GROUP_IN_ONE_ROW,
&p->toprint,
GAL_OPTIONS_NO_ARG_TYPE,
@@ -169,7 +182,7 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_STD,
0,
0,
- "Print the standad deviation.",
+ "Standad deviation.",
ARGS_GROUP_IN_ONE_ROW,
&p->toprint,
GAL_OPTIONS_NO_ARG_TYPE,
@@ -183,7 +196,7 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_MEDIAN,
0,
0,
- "Print the median.",
+ "Median.",
ARGS_GROUP_IN_ONE_ROW,
&p->toprint,
GAL_OPTIONS_NO_ARG_TYPE,
@@ -193,11 +206,81 @@ struct argp_option program_options[] =
ui_add_to_print_in_row
},
{
+ "quantile",
+ ARGS_OPTION_KEY_QUANTILE,
+ "FLT[,...]",
+ 0,
+ "Quantile (multiple values acceptable).",
+ ARGS_GROUP_IN_ONE_ROW,
+ &p->toprint,
+ GAL_DATA_TYPE_FLOAT32,
+ GAL_OPTIONS_RANGE_GE_0_LE_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_add_to_print_in_row
+ },
+ {
+ "quantfunc",
+ ARGS_OPTION_KEY_QUANTFUNC,
+ "FLT[,...]",
+ 0,
+ "Quantile function (multiple values acceptable).",
+ ARGS_GROUP_IN_ONE_ROW,
+ &p->toprint,
+ GAL_DATA_TYPE_FLOAT32,
+ GAL_OPTIONS_RANGE_GE_0_LE_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_add_to_print_in_row
+ },
+ {
"mode",
ARGS_OPTION_KEY_MODE,
0,
0,
- "Print the mode (for large datasets).",
+ "Mode (Appendix D of arXiv:1505.011664).",
+ ARGS_GROUP_IN_ONE_ROW,
+ &p->toprint,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_add_to_print_in_row
+ },
+ {
+ "modequant",
+ ARGS_OPTION_KEY_MODEQUANT,
+ 0,
+ 0,
+ "Mode quantile (see --mode)",
+ ARGS_GROUP_IN_ONE_ROW,
+ &p->toprint,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_add_to_print_in_row
+ },
+ {
+ "modesym",
+ ARGS_OPTION_KEY_MODESYM,
+ 0,
+ 0,
+ "Mode symmetricity (see --mode).",
+ ARGS_GROUP_IN_ONE_ROW,
+ &p->toprint,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_add_to_print_in_row
+ },
+ {
+ "modesymvalue",
+ ARGS_OPTION_KEY_MODESYMVALUE,
+ 0,
+ 0,
+ "Value at mode symmetricity (see --mode).",
ARGS_GROUP_IN_ONE_ROW,
&p->toprint,
GAL_OPTIONS_NO_ARG_TYPE,
@@ -206,7 +289,6 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET,
ui_add_to_print_in_row
},
-
@@ -282,6 +364,19 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET,
ui_parse_numbers
},
+ {
+ "mirror",
+ ARGS_OPTION_KEY_MIRROR,
+ "FLT",
+ 0,
+ "Save the histogram and CFP of the mirror dist.",
+ ARGS_GROUP_PARTICULAR_STAT,
+ &p->mirror,
+ GAL_DATA_TYPE_FLOAT64,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
@@ -289,7 +384,7 @@ struct argp_option program_options[] =
{
0, 0, 0, 0,
- "Histogram or Cumulative frequency plot settings",
+ "Settings",
ARGS_GROUP_HIST_CFP
},
{
diff --git a/bin/statistics/main.h b/bin/statistics/main.h
index af87ccc..59f6db1 100644
--- a/bin/statistics/main.h
+++ b/bin/statistics/main.h
@@ -44,9 +44,11 @@ struct statisticsparams
{
/* From command-line */
struct gal_options_common_params cp; /* Common parameters. */
- struct gal_linkedlist_ill *toprint; /* Values to print in one row. */
+ struct gal_linkedlist_ill *toprint; /* Options to print in one row. */
+ struct gal_linkedlist_dll *tp_args; /* Arguments for printing. */
char *inputname; /* Input filename. */
char *column; /* Column name or number if input is table. */
+ char *refcol; /* Reference column name or number. */
float greaterequal; /* Only use values >= this value. */
float lessthan; /* Only use values < this value. */
float quantmin; /* Quantile min or range: from Q to 1-Q. */
@@ -58,6 +60,8 @@ struct statisticsparams
uint8_t cumulative; /* Save cumulative distibution in output. */
float sigclipmultip; /* Multiple of sigma in sigma clipping. */
float sigclipparam; /* Tolerance to stop or number of clips. */
+ double mirror; /* Mirror value for hist and CFP. */
+
size_t numbins; /* Number of bins in histogram or CFP. */
size_t numasciibins; /* Number of bins in ASCII plots. */
size_t asciiheight; /* Height of ASCII histogram or CFP plots. */
@@ -66,8 +70,10 @@ struct statisticsparams
uint8_t maxbinone; /* Set the maximum bin to 1. */
/* Internal */
+ int numoutfiles; /* Number of output files made in this run. */
uint8_t needssort; /* If sorting is needed. */
gal_data_t *input; /* Input data structure. */
+ gal_data_t *reference; /* Reference 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). */
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index cc5abc2..5bb9b58 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -49,13 +49,31 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/*******************************************************************/
/************** Print in one row ***************/
/*******************************************************************/
-static void
-ui_print_one_number(gal_data_t *out)
+static gal_data_t *
+statistics_pull_out_element(gal_data_t *input, size_t index)
+{
+ size_t dsize=1;
+ gal_data_t *out=gal_data_alloc(NULL, input->type, 1, &dsize,
+ NULL, 1, -1, NULL, NULL, NULL);
+ gal_data_copy_element_same_type(input, index, out->array);
+ return out;
+}
+
+
+
+
+
+static double
+statistics_read_check_args(struct statisticsparams *p)
{
- char *toprint=gal_data_write_to_string(out->array, out->type, 0);
- printf("%s ", toprint);
- gal_data_free(out);
- free(toprint);
+ double d;
+ if(p->tp_args==NULL)
+ error(EXIT_FAILURE, 0, "a bug! Not enough arguments for the "
+ "requested options to print in one row. Please contact "
+ "us at %s so we can address the problem",
+ PACKAGE_BUGREPORT);
+ gal_linkedlist_pop_from_dll(&p->tp_args, &d);
+ return d;
}
@@ -63,39 +81,130 @@ ui_print_one_number(gal_data_t *out)
static void
-ui_print_one_row(struct statisticsparams *p)
+statistics_print_one_row(struct statisticsparams *p)
{
+ int mustfree;
+ char *toprint;
+ size_t dsize=1;
+ double arg, *d;
+ float mirrdist=1.5;
struct gal_linkedlist_ill *tmp;
+ gal_data_t *out, *tmpv, *num=NULL, *min=NULL, *max=NULL;
+ gal_data_t *sum=NULL, *med=NULL, *meanstd=NULL, *modearr=NULL;
+
- /* Print the numbers. */
+ /* The user can ask for any of the operators more than once, also some
+ operators might return more than one usable value (like mode). So we
+ will calculate the desired values once, and then print them. */
for(tmp=p->toprint; tmp!=NULL; tmp=tmp->next)
switch(tmp->v)
{
+ /* Calculate respective values. Checking with `if(num==NULL)' gives
+ compiler warnings of `this if clause does not guard ...'. So we
+ are using this empty-if and else statement. */
case ARGS_OPTION_KEY_NUMBER:
- ui_print_one_number( gal_statistics_number(p->input) ); break;
+ num = num ? num : gal_statistics_number(p->input); break;
case ARGS_OPTION_KEY_MINIMUM:
- ui_print_one_number( gal_statistics_minimum(p->input) ); break;
+ min = min ? min : gal_statistics_minimum(p->input); break;
case ARGS_OPTION_KEY_MAXIMUM:
- ui_print_one_number( gal_statistics_maximum(p->input) ); break;
+ max = max ? max : gal_statistics_maximum(p->input); break;
case ARGS_OPTION_KEY_SUM:
- ui_print_one_number( gal_statistics_sum(p->input) ); break;
+ sum = sum ? sum : gal_statistics_sum(p->input); break;
+ case ARGS_OPTION_KEY_MEDIAN:
+ med = med ? med : gal_statistics_median(p->sorted, 0); 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;
+ meanstd = meanstd ? meanstd : gal_statistics_mean_std(p->input);
+ break;
case ARGS_OPTION_KEY_MODE:
- error(EXIT_FAILURE, 0, "mode isn't implemented yet!");
+ case ARGS_OPTION_KEY_MODEQUANT:
+ case ARGS_OPTION_KEY_MODESYM:
+ case ARGS_OPTION_KEY_MODESYMVALUE:
+ modearr = ( modearr ? modearr
+ : gal_statistics_mode(p->sorted, mirrdist, 0) );
+ d=modearr->array;
+ if(d[2]<GAL_STATISTICS_MODE_GOOD_SYM) d[0]=d[1]=NAN;
+ break;
+
+ /* Will be calculated as printed. */
+ case ARGS_OPTION_KEY_QUANTILE:
+ case ARGS_OPTION_KEY_QUANTFUNC:
break;
+
+ /* The option isn't recognized. */
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);
+ "`statistics_print_one_row'. Please contact us at %s so we "
+ "can address the problem", tmp->v, PACKAGE_BUGREPORT);
}
+
+ /* Print every requested number. */
+ for(tmp=p->toprint; tmp!=NULL; tmp=tmp->next)
+ {
+ /* By default don't free anything. */
+ mustfree=0;
+
+ /* Get the output. */
+ switch(tmp->v)
+ {
+ /* Previously calculated values. */
+ case ARGS_OPTION_KEY_NUMBER: out=num; break;
+ case ARGS_OPTION_KEY_MINIMUM: out=min; break;
+ case ARGS_OPTION_KEY_MAXIMUM: out=max; break;
+ case ARGS_OPTION_KEY_SUM: out=sum; break;
+ case ARGS_OPTION_KEY_MEDIAN: out=med; break;
+ case ARGS_OPTION_KEY_MEAN:
+ out=statistics_pull_out_element(meanstd, 0); mustfree=1; break;
+ case ARGS_OPTION_KEY_STD:
+ out=statistics_pull_out_element(meanstd, 1); mustfree=1; break;
+ case ARGS_OPTION_KEY_MODE:
+ out=statistics_pull_out_element(modearr, 0); mustfree=1; break;
+ case ARGS_OPTION_KEY_MODEQUANT:
+ out=statistics_pull_out_element(modearr, 1); mustfree=1; break;
+ case ARGS_OPTION_KEY_MODESYM:
+ out=statistics_pull_out_element(modearr, 2); mustfree=1; break;
+ case ARGS_OPTION_KEY_MODESYMVALUE:
+ out=statistics_pull_out_element(modearr, 3); mustfree=1; break;
+
+ /* Not previously calculated. */
+ case ARGS_OPTION_KEY_QUANTILE:
+ arg = statistics_read_check_args(p);
+ out = gal_statistsics_quantile(p->sorted, arg, 0);
+ break;
+
+ case ARGS_OPTION_KEY_QUANTFUNC:
+ arg = statistics_read_check_args(p);
+ tmpv = gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &dsize,
+ NULL, 1, -1, NULL, NULL, NULL);
+ *((double *)(tmpv->array)) = arg;
+ tmpv = gal_data_copy_to_new_type_free(tmpv, p->input->type);
+ out = gal_statistics_quantile_function(p->sorted, tmpv, 0);
+ break;
+ }
+
+ /* Print the number. */
+ toprint=gal_data_write_to_string(out->array, out->type, 0);
+ printf("%s ", toprint);
+ free(toprint);
+
+ /* Clean up (if necessary). */
+ if(mustfree) gal_data_free(out);
+ }
+
+
/* Print a new line. */
printf("\n");
+
+
+ /* Clean any of the allocated arrays. */
+ if(num) gal_data_free(num);
+ if(min) gal_data_free(min);
+ if(max) gal_data_free(max);
+ if(sum) gal_data_free(sum);
+ if(med) gal_data_free(med);
+ if(meanstd) gal_data_free(meanstd);
+ if(modearr) gal_data_free(modearr);
}
@@ -212,12 +321,78 @@ ascii_plots(struct statisticsparams *p)
/*******************************************************************/
/******* Histogram and cumulative frequency tables ***********/
/*******************************************************************/
+void
+write_output_table(struct statisticsparams *p, gal_data_t *table,
+ char *suf, char *contents)
+{
+ char *output;
+ int use_auto_output=0;
+ char *fix, *suffix, *tmp;
+ struct gal_linkedlist_stll *comments=NULL;
+
+
+ /* See if we should use automatic output. */
+ if(p->cp.output)
+ {
+ if(p->numoutfiles>1) use_auto_output=1;
+ }
+ else use_auto_output=1;
+
+
+ /* Set the type suffix */
+ if(use_auto_output)
+ {
+ fix = ( p->cp.output
+ ? gal_fits_name_is_fits(p->cp.output) ? "fits" : "txt"
+ : "txt" );
+ asprintf(&suffix, "%s.%s", suf, fix);
+ }
+
+
+ /* Make the output name. */
+ output= ( use_auto_output
+ ? gal_checkset_automatic_output(&p->cp, p->inputname, suffix)
+ : p->cp.output );
+
+
+ /* 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(table, 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);
+}
+
+
+
+
+
static void
save_hist_and_or_cfp(struct statisticsparams *p)
{
+ char *suf, *contents;
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
@@ -265,47 +440,47 @@ save_hist_and_or_cfp(struct statisticsparams *p)
{ suf="_cfp"; contents="Cumulative frequency plot"; }
- /* 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) );
-
+ /* Set the output file name. */
+ write_output_table(p, bins, suf, contents);
+}
- /* 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);
+void
+print_mirror_hist_cfp(struct statisticsparams *p)
+{
+ size_t dsize=1;
+ gal_data_t *table;
+ double mirror_val;
+ gal_data_t *mirror=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &dsize,
+ NULL, 1, -1, NULL, NULL, NULL);
+ /* Convert the given mirror value into the type of the input dataset. */
+ *((double *)(mirror->array)) = p->mirror;
+ mirror=gal_data_copy_to_new_type_free(mirror, p->input->type);
- /* Let the user know, if we aren't in quiet mode. */
- if(!p->cp.quiet)
- printf("%s created.\n", output);
+ /* Make the table columns. */
+ table=gal_statistics_mode_mirror_plots(p->sorted, mirror, p->numbins, 0,
+ &mirror_val);
+ if(p->mirror!=mirror_val)
+ {
+ fprintf(stderr, "Warning: Mirror value is %f.\n", mirror_val);
+ if(!p->cp.quiet)
+ fprintf(stderr, "\nNote that the mirror distribution is discrete "
+ "and depends on the input data. So the closest point in "
+ "the data to your desired mirror at %f was %f.\n\n",
+ p->mirror, mirror_val);
+ }
- /* Clean up. */
- free(suffix);
- if(output!=p->cp.output) free(output);
- gal_linkedlist_free_stll(comments, 1);
+ /* If the mirror value was out-of-range, then no table will be made. */
+ if(table)
+ write_output_table(p, table, "_mirror_hist_cfp",
+ "Histogram and CFP of mirror distribution");
+ else
+ error(EXIT_FAILURE, 0, "mirror value %g is out of range", p->mirror);
}
@@ -325,7 +500,6 @@ save_hist_and_or_cfp(struct statisticsparams *p)
-
/*******************************************************************/
/************** Basic information ***************/
/*******************************************************************/
@@ -347,7 +521,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->column);
+ if(p->column) printf("Column: %s\n",
+ p->input->name ? p->input->name : p->column);
/* Range. */
str=NULL;
@@ -360,7 +535,11 @@ print_input_info(struct statisticsparams *p)
asprintf(&str, "upto (exclusive) %g", p->lessthan);
if(str)
{
- printf("Range: %s.\n", str);
+ printf("Range: ");
+ if(p->refcol)
+ printf("[on column %s] ",
+ p->reference->name ? p->reference->name : p->refcol);
+ printf("%s.\n", str);
free(str);
}
@@ -384,8 +563,9 @@ void
print_basics(struct statisticsparams *p)
{
char *str;
- double mean, std;
int namewidth=40;
+ float mirrdist=1.5;
+ double mean, std, *d;
gal_data_t *tmp, *bins, *hist;
/* Define the input dataset. */
@@ -415,19 +595,30 @@ print_basics(struct statisticsparams *p)
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.*/
+ /* Mode of the distribution (if it is valid). we want the mode and 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_mode(p->input, mirrdist, 0);
+ d=tmp->array;
+ if(d[2]>GAL_STATISTICS_MODE_GOOD_SYM)
+ { /* Same format as `gal_data_write_to_string' */
+ printf(" %-*s %.10g\n", namewidth, "Mode:", d[0]);
+ printf(" %-*s %.10g\n", namewidth, "Mode quantile:", d[1]);
+ }
+ gal_data_free(tmp);
+
+ /* Find and print the median: */
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);
+ /* Print the mean and standard deviation. Same format as
+ `gal_data_write_to_string' */
+ printf(" %-*s %.10g\n", namewidth, "Mean:", mean);
+ printf(" %-*s %.10g\n", namewidth, "Standard deviation:", std);
/* Ascii histogram. Note that we don't want to force the user to have the
plotting parameters. */
@@ -527,7 +718,6 @@ print_sigma_clip(struct statisticsparams *p)
-
/*******************************************************************/
/************** Main function ***************/
/*******************************************************************/
@@ -540,7 +730,7 @@ statistics(struct statisticsparams *p)
if(p->toprint)
{
print_basic_info=0;
- ui_print_one_row(p);
+ statistics_print_one_row(p);
}
/* Print the ASCII plots if requested. */
@@ -558,12 +748,19 @@ statistics(struct statisticsparams *p)
}
/* Print the sigma-clipped results. */
- if( !isnan(p->sigclipmultip ) )
+ if( !isnan(p->sigclipmultip) )
{
print_basic_info=0;
print_sigma_clip(p);
}
+ /* Make the mirror table. */
+ if( !isnan(p->mirror) )
+ {
+ print_basic_info=0;
+ print_mirror_hist_cfp(p);
+ }
+
/* 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 1a9b7a7..1abe260 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -131,6 +131,7 @@ ui_initialize_options(struct statisticsparams *p,
p->greaterequal = NAN;
p->quantmin = NAN;
p->quantmax = NAN;
+ p->mirror = NAN;
p->sigclipparam = NAN;
p->sigclipmultip = NAN;
@@ -205,8 +206,12 @@ static void *
ui_add_to_print_in_row(struct argp_option *option, char *arg,
char *filename, size_t lineno, void *params)
{
+ size_t i;
+ double *d;
+ gal_data_t *inputs=NULL;
struct statisticsparams *p=(struct statisticsparams *)params;
+ /* In case of printing the option values. */
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 "
@@ -215,19 +220,62 @@ ui_add_to_print_in_row(struct argp_option *option, char
*arg,
"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' )
- error_at_line(EXIT_FAILURE, 0, filename, lineno, "the `--%s' "
- "option takes no arguments. In a configuration file "
- "it can only have the values `1' or `0', indicating "
- "if it should be used or not", option->name);
+ /* Some of these options take values and some don't. */
+ if(option->type==GAL_OPTIONS_NO_ARG_TYPE)
+ {
+ /* 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)
+ {
+ /* Make sure the value is only `0' or `1'. */
+ if( arg[1]!='\0' && *arg!='0' && *arg!='1' )
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "the `--%s' "
+ "option takes no arguments. In a configuration "
+ "file it can only have the values `1' or `0', "
+ "indicating if it should be used or not",
+ option->name);
+
+ /* Only proceed if the (possibly given) argument is 1. */
+ if(arg[0]=='0' && arg[1]=='\0') return NULL;
+ }
+
+ /* Add this option to the print list. */
+ gal_linkedlist_add_to_ill(&p->toprint, option->key);
+ }
+ else
+ {
+ /* Read the string of numbers. */
+ inputs=gal_options_parse_list_of_numbers(arg, filename, lineno);
+ d=inputs->array;
- /* Only proceed if the (possibly given) argument is 1. */
- if(arg && *arg=='0') return NULL;
+ /* Do the appropriate operations with the */
+ switch(option->key)
+ {
+ case ARGS_OPTION_KEY_QUANTILE:
+ case ARGS_OPTION_KEY_QUANTFUNC:
+ /* For the quantile and the quantile function, its possible to
+ give any number of arguments, so add the operation index and
+ the argument once for each given number. */
+ for(i=0;i<inputs->size;++i)
+ {
+ if(option->key==ARGS_OPTION_KEY_QUANTILE && (d[i]<0 || d[i]>1) )
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "values "
+ "to `--quantile' (`-u') must be between 0 "
+ "and 1, you had asked for %g (read from `%s')",
+ d[i], arg);
+ gal_linkedlist_add_to_dll(&p->tp_args, d[i]);
+ gal_linkedlist_add_to_ill(&p->toprint, option->key);
+ }
+ break;
- /* Add this option to the print list. */
- gal_linkedlist_add_to_ill(&p->toprint, option->key);
+ default:
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "a bug! please "
+ "contact us at %s so we can address the problem. "
+ "the option given to `ui_add_to_print_in_row' is "
+ "marked as requiring a value, but is not recognized",
+ PACKAGE_BUGREPORT);
+ }
+ }
return NULL;
}
@@ -411,7 +459,7 @@ ui_read_check_only_options(struct statisticsparams *p)
/* When binned outputs are requested, make sure that `numbins' is set. */
- if( (p->histogram || p->cumulative) && p->numbins==0)
+ if( (p->histogram || 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");
@@ -426,9 +474,10 @@ ui_read_check_only_options(struct statisticsparams *p)
"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. */
+ /* Reverse the list of statistics to print in one row and also the
+ arguments, so it has the same order the user wanted. */
gal_linkedlist_reverse_ill(&p->toprint);
+ gal_linkedlist_reverse_dll(&p->tp_args);
}
@@ -514,10 +563,15 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
{
size_t one=1;
unsigned char flags=GAL_ARITHMETIC_NUMOK;
- gal_data_t *tmp, *cond_g=NULL, *cond_l=NULL, *cond;
unsigned char flagsor = ( GAL_ARITHMETIC_FREE
| GAL_ARITHMETIC_INPLACE
| GAL_ARITHMETIC_NUMOK );
+ gal_data_t *tmp, *cond_g=NULL, *cond_l=NULL, *cond, *blank, *ref;
+
+
+ /* Set the dataset that should be used for the condition. */
+ ref = p->reference ? p->reference : p->input;
+
/* If the user has given a quantile range, then set the `greaterequal'
and `lessthan' values. */
@@ -527,16 +581,17 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
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_statistsics_quantile(ref, 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_statistsics_quantile(ref, 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 */
@@ -545,20 +600,22 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
tmp=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, 1, &one, NULL, 0, -1,
NULL, NULL, NULL);
*((float *)(tmp->array)) = p->greaterequal;
- cond_g=gal_arithmetic(GAL_ARITHMETIC_OP_LT, flags, p->input, tmp);
+ cond_g=gal_arithmetic(GAL_ARITHMETIC_OP_LT, flags, ref, tmp);
gal_data_free(tmp);
}
+
/* Same reasoning as above for `p->greaterthan'. */
if(!isnan(p->lessthan))
{
tmp=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, 1, &one, NULL, 0, -1,
NULL, NULL, NULL);
*((float *)(tmp->array)) = p->lessthan;
- cond_l=gal_arithmetic(GAL_ARITHMETIC_OP_GE, flags, p->input, tmp);
+ cond_l=gal_arithmetic(GAL_ARITHMETIC_OP_GE, flags, ref, tmp);
gal_data_free(tmp);
}
+
/* Now, set the final condition. If both values were specified, then use
the GAL_ARITHMETIC_OP_OR to merge them into one. */
switch( !isnan(p->greaterequal) + !isnan(p->lessthan) )
@@ -572,15 +629,18 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
break;
}
- /* Allocate a blank value to mask all input pixels. */
- tmp=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, 1, &one, NULL,
+
+ /* Allocate a blank value to mask all pixels that don't satisfy the
+ condition. */
+ blank=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, 1, &one, NULL,
0, -1, NULL, NULL, NULL);
- *((float *)(tmp->array)) = NAN;
+ *((float *)(blank->array)) = NAN;
+
/* 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);
+ gal_arithmetic(GAL_ARITHMETIC_OP_WHERE, flagsor, p->input, cond, blank);
}
@@ -599,27 +659,32 @@ ui_make_sorted_if_necessary(struct statisticsparams *p)
for(tmp=p->toprint; tmp!=NULL; tmp=tmp->next)
switch(tmp->v)
{
- case ARGS_OPTION_KEY_MEDIAN:
case ARGS_OPTION_KEY_MODE:
+ case ARGS_OPTION_KEY_MEDIAN:
+ case ARGS_OPTION_KEY_QUANTILE:
+ case ARGS_OPTION_KEY_QUANTFUNC:
is_necessary=1;
break;
}
/* Check in the rest of the outputs. */
- if( is_necessary==0 && !isnan(p->sigclipmultip) )
+ if( is_necessary==0
+ && ( !isnan(p->sigclipmultip) || !isnan(p->mirror) ) )
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
+ if(is_necessary)
{
- p->sorted=gal_data_copy(p->input);
- gal_statistics_sort_increasing(p->sorted);
+ 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);
+ }
}
}
@@ -628,36 +693,81 @@ ui_make_sorted_if_necessary(struct statisticsparams *p)
void
-ui_preparations(struct statisticsparams *p)
+ui_read_columns(struct statisticsparams *p)
{
- gal_data_t *tmp;
- char *errorstring;
- size_t numcolmatches=0;
+ int toomanycols=0;
+ size_t size, counter=0;
+ gal_data_t *cols, *tmp;
struct gal_linkedlist_stll *column=NULL;
- /* Read the input: no matter if its an image or a table column. */
- if(p->isfits && p->hdu_type==IMAGE_HDU)
- p->input=gal_fits_img_read(p->inputname, p->cp.hdu, p->cp.minmapsize);
- else
+ /* Define the columns that we want, note that they should be added to the
+ list in reverse. */
+ if(p->refcol)
+ gal_linkedlist_add_to_stll(&column, p->refcol, 0);
+ gal_linkedlist_add_to_stll(&column, p->column, 0);
+
+ /* Read the desired column(s). */
+ cols=gal_table_read(p->inputname, p->cp.hdu, column, p->cp.searchin,
+ p->cp.ignorecase, p->cp.minmapsize);
+
+ /* Put the columns into the proper gal_data_t. */
+ size=cols->size;
+ while(cols!=NULL)
{
- /* Read the input column. */
- gal_linkedlist_add_to_stll(&column, p->column, 0);
- p->input=gal_table_read(p->inputname, p->cp.hdu, column, p->cp.searchin,
- p->cp.ignorecase, p->cp.minmapsize);
- if(p->input->next)
+ /* Pop out the top column. */
+ tmp=gal_data_pop_from_ll(&cols);
+
+ /* Make sure it has the proper size. */
+ if(tmp->size!=size)
+ error(EXIT_FAILURE, 0, " read column number %zu has a %zu elements, "
+ "while previous column(s) had %zu", counter, tmp->size, size);
+
+ /* Make sure it is a usable datatype. */
+ switch(tmp->type)
+ {
+ case GAL_DATA_TYPE_BIT:
+ case GAL_DATA_TYPE_STRLL:
+ case GAL_DATA_TYPE_STRING:
+ case GAL_DATA_TYPE_COMPLEX32:
+ case GAL_DATA_TYPE_COMPLEX64:
+ error(EXIT_FAILURE, 0, " read column number %zu has a %s type, "
+ "which is not currently supported by %s", counter,
+ gal_data_type_as_string(tmp->type, 1), PROGRAM_NAME);
+ }
+
+ /* Put the column into the proper pointer. */
+ switch(++counter)
{
- for(tmp=p->input;tmp!=NULL;tmp=tmp->next) ++numcolmatches;
- asprintf(&errorstring, "%zu columns were selected with `%s' "
- "(value to `--column' option). In this context, "
- "Statistics can only work on one data-set (column in a "
- "table).", numcolmatches, p->column);
- gal_table_error_col_selection(p->inputname, p->cp.hdu, errorstring);
+ case 1: p->input=tmp; break;
+ case 2: if(p->refcol) p->reference=tmp; else toomanycols=1; break;
+ default: toomanycols=1;
}
- /* Clean up. */
- gal_linkedlist_free_stll(column, 0);
+ /* Print an error if there are too many columns: */
+ if(toomanycols)
+ gal_table_error_col_selection(p->inputname, p->cp.hdu, "too many "
+ "columns were selected by the given "
+ "values to the `--column' and/or "
+ "`--refcol' options. Only one "
+ "is acceptable for each.");
}
+ /* Clean up. */
+ gal_linkedlist_free_stll(column, 0);
+}
+
+
+
+
+
+void
+ui_preparations(struct statisticsparams *p)
+{
+ /* Read the input. */
+ if(p->isfits && p->hdu_type==IMAGE_HDU)
+ p->input=gal_fits_img_read(p->inputname, p->cp.hdu, p->cp.minmapsize);
+ else ui_read_columns(p);
+
/* Set the out-of-range values in the input to blank. */
ui_out_of_range_to_blank(p);
@@ -672,6 +782,10 @@ ui_preparations(struct statisticsparams *p)
/* Make the sorted array if necessary. */
ui_make_sorted_if_necessary(p);
+
+ /* Set the number of output files. */
+ if( !isnan(p->mirror) ) ++p->numoutfiles;
+ if( p->histogram || p->cumulative ) ++p->numoutfiles;
}
@@ -779,4 +893,7 @@ ui_free_report(struct statisticsparams *p)
if(p->sorted!=p->input)
gal_data_free(p->sorted);
gal_data_free(p->input);
+ gal_data_free(p->reference);
+ gal_linkedlist_free_ill(p->toprint);
+ gal_linkedlist_free_dll(p->tp_args);
}
diff --git a/bin/statistics/ui.h b/bin/statistics/ui.h
index 5162519..519d64d 100644
--- a/bin/statistics/ui.h
+++ b/bin/statistics/ui.h
@@ -29,13 +29,14 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Available letters for short options:
- a b c d e f i j k p r u v w x y z
+ a b c d e f i j k p 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_REFCOL = 'r',
ARGS_OPTION_KEY_GREATEREQUAL = 'g',
ARGS_OPTION_KEY_LESSTHAN = 'l',
ARGS_OPTION_KEY_QRANGE = 'Q',
@@ -43,6 +44,7 @@ enum option_keys_enum
ARGS_OPTION_KEY_STD = 't',
ARGS_OPTION_KEY_MEDIAN = 'M',
ARGS_OPTION_KEY_MODE = 'O',
+ ARGS_OPTION_KEY_QUANTILE = 'u',
ARGS_OPTION_KEY_ASCIIHIST = 'A',
ARGS_OPTION_KEY_HISTOGRAM = 'H',
ARGS_OPTION_KEY_CUMULATIVE = 'C',
@@ -55,7 +57,12 @@ enum option_keys_enum
ARGS_OPTION_KEY_MINIMUM,
ARGS_OPTION_KEY_MAXIMUM,
ARGS_OPTION_KEY_SUM,
+ ARGS_OPTION_KEY_MODEQUANT,
+ ARGS_OPTION_KEY_MODESYM,
+ ARGS_OPTION_KEY_MODESYMVALUE,
+ ARGS_OPTION_KEY_QUANTFUNC,
ARGS_OPTION_KEY_ASCIICFP,
+ ARGS_OPTION_KEY_MIRROR,
ARGS_OPTION_KEY_NUMBINS,
ARGS_OPTION_KEY_NUMASCIIBINS,
ARGS_OPTION_KEY_ASCIIHEIGHT,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index f360814..107758a 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -438,7 +438,6 @@ Statistics
* Histogram and Cumulative Frequency Plot:: Basic definitions.
* Sigma clipping:: Definition of @mymath{\sigma}-clipping
-* Mirror distribution:: Used for finding the mode.
* Invoking aststatistics:: Arguments and options to Statistics.
NoiseChisel
@@ -11101,7 +11100,6 @@ Statistics program is designed for such situations.
@menu
* Histogram and Cumulative Frequency Plot:: Basic definitions.
* Sigma clipping:: Definition of @mymath{\sigma}-clipping
-* Mirror distribution:: Used for finding the mode.
* Invoking aststatistics:: Arguments and options to Statistics.
@end menu
@@ -11179,7 +11177,7 @@ 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.
address@hidden Sigma clipping, Mirror distribution, Histogram and Cumulative
Frequency Plot, Statistics
address@hidden Sigma clipping, Invoking aststatistics, Histogram and
Cumulative Frequency Plot, Statistics
@subsection Sigma clipping
Let's assume that you have pure noise (centered on zero) with a clear
@@ -11274,102 +11272,7 @@ ASCII histogram that is printed on the command-line
with
@option{--asciihist} is good enough in most cases.
@end cartouche
address@hidden Mirror distribution, Invoking aststatistics, Sigma clipping,
Statistics
address@hidden Mirror distribution
-
address@hidden Mirror distribution
-The mirror distribution of a data set was defined in Appendix C of
-Akhlaghi and Ichikawa (2015). It is best visualized by mentally
-placing a mirror on the histogram of a distribution at any point
-within the distribution (which we call the mirror point).
-
-Through the @option{--mirrorquant} in Statistics, you can check the mirror
-of a distribution when the mirror is placed on any given quantile. The
-mirror distribution is plotted along with the input distribution both as
-histograms and cumulative frequency plots, see @ref{Histogram and
-Cumulative Frequency Plot}. Unlike the rest of the histograms and
-cumulative frequency plots in Statistics, the text files created with the
address@hidden and @option{--checkmode} will contain 3 columns. The
-first is the horizontal axis similar to all other histograms and cumulative
-frequency plots. The second column shows the respective value for the
-actual data distribution and the third shows the value for the mirror
-distribution.
-
address@hidden Python
-The value for each bin of both histogram is divided by the maximum of
-both. For the cumulative frequency plot, the value in each bin is
-divided by the maximum number of elements. So one of the cumulative
-frequency plots reaches the maximum vertical axis of 1. The outputs
-will have the @file{_mirrorhist.txt} and @file{_mirrorcfp.txt} suffixes
-respectively. You can use a simple Python script like the one below to
-display the histograms and cumulative frequency plots in one plot:
-
address@hidden
-#! /usr/bin/env python3
-
-# Import the necessary modules:
-import sys
-import numpy as np
-import matplotlib.pyplot as plt
-
-# Load the two files:
-a=np.loadtxt(sys.argv[1]+"_mirrorhist.txt")
-b=np.loadtxt(sys.argv[1]+"_mirrorcfp.txt")
-
-# Calculate the bin width:
-w=a[1,0]-a[0,0]
-
-# Plot the two histograms and cumulative frequency plots:
-plt.bar(a[:,0], a[:,1], width=w, color="blue", linewidth=0, alpha=0.6)
-plt.bar(a[:,0], a[:,2], width=w, color="green", linewidth=0, alpha=0.4)
-plt.plot(b[:,0], b[:,1], linewidth=2, color="blue")
-plt.plot(b[:,0], b[:,2], linewidth=2, color="green")
-
-# Write the axis labels:
-plt.ylim([0,np.amax(a[:,1])])
-plt.xlim([np.amin(a[:,0]),np.amax(a[:,0])])
-
-# Save the output to any name from the command-line:
-plt.savefig(sys.argv[1]+"_plot.pdf")
address@hidden example
-
address@hidden PNG
address@hidden Matplotlib
address@hidden
-The output format can be anything that Python's Matplotlib recognizes
-(for example png or jpg are also acceptable). So, if your input data
-file name was @file{input.fits}, and you want to see how its mirror
-distribution would look like if the mirror was placed at the 0.8
-quantile (or the value which is above 80 percent of your data) you can
-run the following sequence of commands to see the combined cumulative
-frequency plot and histogram together in one PDF file. Let's assume you
-have put the Python script above into the file
address@hidden The second command makes the Python script an
-executable file.
-
address@hidden
-$ ls
-input.fits mirrorplot.py
-$ chmod +x mirrorplot.py
-$ aststatistics input.fits --nohist --nocfp --mirrorquant=0.8
-$ ls
-input.fits input_mirrorcfp.txt input_mirrorhist.txt mirrorplot.py
-$ ./mirrorplot.py input
address@hidden example
-
-By default, the range of the mirror distribution is set to the 0.01
-quantile of the image to 2 times the distance of the mode to that
-point. This is done so that the mode (which will have a value of zero
-in this plot) is positioned exactly on 1/3rd point of the x axis
-plot. The number of bins is the same value as the
address@hidden option. Alternatively, through the option
address@hidden, the histogram properties set with the
address@hidden and @option{--histmax} can be used. In these mirror
-plots, the mirror value is going to have the value of zero and one of
-the bins is going to start at zero (if @option{--histrangeformirror}
-is called, the given ranges will also shift accordingly).
-
address@hidden Invoking aststatistics, , Mirror distribution, Statistics
address@hidden Invoking aststatistics, , Sigma clipping, Statistics
@subsection Invoking Statistics
Statistics will print statistical measures of an input dataset (table
@@ -11387,8 +11290,8 @@ One line examples:
$ 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
+$ aststatistics tab.txt -rPHOTO_Z -g3 -cMAG_F160W --median
@end example
@noindent
@@ -11431,10 +11334,12 @@ Unit: Brightness
-------
Number of elements: 9074
Minimum: 9622.35
- Maximum: 10999.8
+ Maximum: 10999.7
+ Mode: 10055.45996
+ Mode quantile: 0.4001983908
Median: 10093.7
- Mean: 10144
- Standard deviation: 221.81
+ Mean: 10143.98257
+ Standard deviation: 221.80834
-------
Histogram:
| **
@@ -11460,10 +11365,20 @@ programs including Statistics, see @ref{Common
options} for those.
@item -c STR/INT
@itemx --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}.
+The input column selector when the input file is a table. See
address@hidden 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 -r STR/INT
address@hidden --refcol=STR/INT
+The reference column selector when the input file is a table. When a
+reference column is given, the range options below will be applied to this
+column and only elements in the input column that have a reference value in
+the correct range will be used. In practice this option allows you to
+select a subset of the input column based on values in another (the
+reference) column. All the statistical calculations will be done on the
+selected input column, not the reference column.
@item -g FLT
@itemx --greaterequal=FLT
@@ -11541,9 +11456,70 @@ Print the standard deviation of all used elements.
@itemx --median
Print the median of all used elements.
address@hidden -u FLT[,FLT[,...]]
address@hidden --quantile=FLT[,FLT[,...]]
+Print the values at the given quantiles of the input dataset. Any number of
+quantiles may be given and one number will be printed for each. Values can
+either be written as a single number or as fractions, but must be between
+zero and one (inclusive). Hence, in effect @command{--quantile=0.25
+--quantile=0.75} is equivalent to @option{--quantile=0.25,3/4}, or
address@hidden/4,3/4}.
+
+The returned value is one of the elements from the dataset. Taking
address@hidden to be your desired quantile, and @mymath{N} to be the total
+number of used (non-blank and within the given range) elements, the
+returned value is at the following position in the sorted array:
address@hidden(q\times{}N}).
+
address@hidden --quantfunc=FLT[,FLT[,...]]
+Print the quantiles of the given values in the dataset. This option is the
+inverse of the @option{--quantile} and operates similarly except that the
+acceptable values are within the range of the dataset, not between 0 and
+1. Formally it is known as the ``Quantile function''.
+
+Since the dataset is not continuous this function will find the nearest
+element of the dataset and use its position to estimate the quantile
+function.
+
@item -O
@itemx --mode
-Print the mode of all used elements.
+Print the mode of all used elements. The mode is found through the mirror
+distribution which is fully described in Appendix C of
address@hidden://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa 2015}. See
+that section for a full description.
+
+This mode calculation algorithm is non-parametric, so when the dataset is
+not large enough (larger than about 1000 elements usually), or doesn't have
+a clear mode it can fail. In such cases, this option will return a value of
address@hidden (for the floating point NaN value).
+
+As described in that paper, the easiest way to assess the quality of this
+mode calculation method is to use it's symmetricity (see @option{--modesym}
+below). A better way would be to use the @option{--mirror} option to
+generate the histogram and cumulative frequency tables for any given mirror
+value (the mode in this case) as a table. If you generate plots like those
+shown in Figure 21 of that paper, then your mode is accurate.
+
address@hidden --modequant
+Print the quantile of the mode. You can get the actual mode value from the
address@hidden described above. In many cases, the absolute value of the
+mode is irrelevant, but its position within the distribution is
+important. In such cases, this option will become handy.
+
address@hidden --modesym
+Print the symmetricity of the calculated mode. See the description of
address@hidden for more. This mode algorithm finds the mode based on how
+symmetric it is, so if the symmetricity returned by this option is too low,
+the mode is not too accurate. See Appendix C of
address@hidden://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa 2015} for a
+full description. In practice, symmetricity values larger than 0.2 are
+mostly good.
+
address@hidden --modesymvalue
+Print the value in the distribution where the mirror and input
+distributions are no longer symmetric, see @option{--mode} and Appendix C
+of @url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa 2015} for
+more.
@end table
@@ -11659,6 +11635,31 @@ 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.
address@hidden --mirror=FLT
+Make a histogram and cumulative frequency plot of the mirror distribution
+for the given dataset when the mirror is located at the value to this
+option. The mirror distribution is fully described in Appendix C of
address@hidden://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa 2015} and
+currently it is only used to calculate the mode (see @option{--mode}).
+
+Just note that the mirror distribtution is a discrete distribution like the
+input, so while you may give any number as the value to this option, the
+actual mirror value is the closest number in the input dataset to this
+value. If the two numbers are different, Statistics will warn you of the
+actual mirror value used.
+
+This option will make a table as output. Depending on your selected name
+for the output, it will be either a FITS table or a plain text table (which
+is the default). It contains three columns: the first is the center of the
+bins, the second is the histogram (with the largest value set to 1) and the
+third is the normalized cumulative frequency plot of the mirror
+distribution. The bins will be positioned such that the mode is on the
+starting interval of one of the bins to make it symmetric around the
+mirror. With this output file and the input histogram (that you can
+generate in another run of Statistics, using the @option{--onebinvalue}),
+it is possible to make plots like Figure 21 of
address@hidden://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa 2015}.
+
@end table
The final set of options in Statistics describe the histogram and
diff --git a/lib/Makefile.am b/lib/Makefile.am
index cb587d6..0d2eefa 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -47,9 +47,9 @@ libgnuastro_la_LIBADD =
$(top_builddir)/bootstrapped/lib/libgnu.la
# Specify the library .c files
-libgnuastro_la_SOURCES = arithmetic.c arithmetic-binary.c \
- arithmetic-onlyint.c blank.c box.c checkset.c data.c fits.c git.c \
- linkedlist.c mesh.c mode.c options.c polygon.c qsort.c spatialconvolve.c \
+libgnuastro_la_SOURCES = arithmetic.c arithmetic-binary.c \
+ arithmetic-onlyint.c blank.c box.c checkset.c data.c fits.c git.c \
+ linkedlist.c mesh.c options.c polygon.c qsort.c spatialconvolve.c \
statistics.c table.c threads.c timing.c txt.c wcs.c
@@ -81,7 +81,7 @@ pkginclude_HEADERS = gnuastro/config.h
$(headersdir)/arithmetic.h \
# distribute them here.
EXTRA_DIST = $(headersdir)/README gnuastro.pc.in arithmetic-binary.h \
arithmetic-onlyint.h arithmetic-other.h config.h.in checkset.h \
- fixedstringmacros.h mode.h neighbors.h options.h timing.h
+ fixedstringmacros.h neighbors.h options.h timing.h
diff --git a/lib/data.c b/lib/data.c
index 322965e..a73058e 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -1504,6 +1504,34 @@ gal_data_to_same_type(gal_data_t *f, gal_data_t *s,
+/* Copy/read the element at `index' of the array in `data' into the space
+ pointed to by `ptr'. */
+#define COPY_ELEM(IT) { IT *o=ptr, *a=input->array; *o = a[index]; }
+void
+gal_data_copy_element_same_type(gal_data_t *input, size_t index, void *ptr)
+{
+ /* Set the value. */
+ switch(input->type)
+ {
+ case GAL_DATA_TYPE_UINT8: COPY_ELEM( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: COPY_ELEM( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: COPY_ELEM( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: COPY_ELEM( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: COPY_ELEM( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: COPY_ELEM( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: COPY_ELEM( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: COPY_ELEM( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: COPY_ELEM( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: COPY_ELEM( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`gal_data_copy_elem'", input->type);
+ }
+}
+
+
+
+
diff --git a/lib/gnuastro/data.h b/lib/gnuastro/data.h
index 7c7fae5..31a614c 100644
--- a/lib/gnuastro/data.h
+++ b/lib/gnuastro/data.h
@@ -263,6 +263,9 @@ void
gal_data_to_same_type(gal_data_t *f, gal_data_t *s, gal_data_t **of,
gal_data_t **os, uint8_t type, int freeinputs);
+void
+gal_data_copy_element_same_type(gal_data_t *input, size_t index, void *ptr);
+
/*************************************************************
diff --git a/lib/gnuastro/linkedlist.h b/lib/gnuastro/linkedlist.h
index e5121ee..296c175 100644
--- a/lib/gnuastro/linkedlist.h
+++ b/lib/gnuastro/linkedlist.h
@@ -104,6 +104,9 @@ gal_linkedlist_dll_to_array(struct gal_linkedlist_dll *list,
double **d, size_t *num);
void
+gal_linkedlist_reverse_dll(struct gal_linkedlist_dll **list);
+
+void
gal_linkedlist_free_dll(struct gal_linkedlist_dll *list);
diff --git a/lib/gnuastro/statistics.h b/lib/gnuastro/statistics.h
index f9fce77..de6a9d0 100644
--- a/lib/gnuastro/statistics.h
+++ b/lib/gnuastro/statistics.h
@@ -47,8 +47,11 @@ __BEGIN_C_DECLS /* From C++ preparations */
-/* Maximum number of tests for sigma-clipping convergence */
-#define GAL_STATISTICS_MAX_SIG_CLIP_CONVERGE 50
+/* Maximum number of tests for sigma-clipping convergence. */
+#define GAL_STATISTICS_SIG_CLIP_MAX_CONVERGE 50
+
+/* Least acceptable mode symmetricity.*/
+#define GAL_STATISTICS_MODE_GOOD_SYM 0.2f
@@ -102,10 +105,34 @@ gal_data_t *
gal_statistics_median(gal_data_t *input, int inplace);
gal_data_t *
-gal_statistsics_quantile(gal_data_t *input, float quantile, int inplace);
+gal_statistsics_quantile(gal_data_t *input, double quantile, int inplace);
+
+size_t
+gal_statistics_quantile_index(size_t size, double quant);
size_t
-gal_statistics_quantile_index(size_t size, float quant);
+gal_statistics_quantile_function_index(gal_data_t *input, gal_data_t *value,
+ int inplace);
+
+gal_data_t *
+gal_statistics_quantile_function(gal_data_t *input, gal_data_t *value,
+ int inplace);
+
+
+
+
+
+/****************************************************************
+ ******** Mode *******
+ ****************************************************************/
+
+gal_data_t *
+gal_statistics_mode(gal_data_t *input, float errorstd, int inplace);
+
+gal_data_t *
+gal_statistics_mode_mirror_plots(gal_data_t *input, gal_data_t *value,
+ size_t numbins, int inplace,
+ double *mirror_val);
@@ -129,12 +156,14 @@ gal_statistics_no_blank_sorted(gal_data_t *input, int
inplace);
+
+
/****************************************************************
******** Histogram and Cumulative Frequency Plot *******
****************************************************************/
gal_data_t *
gal_statistics_regular_bins(gal_data_t *data, gal_data_t *range,
- size_t numbins, float onebinstart);
+ size_t numbins, double onebinstart);
gal_data_t *
gal_statistics_histogram(gal_data_t *data, gal_data_t *bins,
@@ -167,48 +196,6 @@ gal_statistics_remove_outliers_flat_cdf(float *sorted,
size_t *outsize);
-
-/****************************************************************/
-/************* Mode ****************/
-/****************************************************************/
-#define GAL_STATISTICS_MODE_LOW_QUANTILE 0.01f
-#define GAL_STATISTICS_MODE_HIGH_QUANTILE 0.51f
-
-#define GAL_STATISTICS_MODE_SYM_GOOD 0.2f
-#define GAL_STATISTICS_MODE_LOW_QUANT_GOOD 0.02f
-
-#define GAL_STATISTICS_MODE_SYMMETRICITY_LOW_QUANT 0.01f
-
-struct gal_statistics_mode_params
-{
- float *sorted; /* Sorted array to be used. */
- size_t size; /* Number of elements in the sorted array. */
- size_t lowi; /* Lower quantile of interval. */
- size_t midi; /* Middle quantile of interval. */
- size_t midd; /* Middle index of inteval. */
- size_t highi; /* Higher quantile of interval. */
- float tolerance; /* Tolerance level to terminate search. */
- size_t numcheck; /* Number of pixels after mode to check. */
- size_t interval; /* Interval to check pixels. */
- float errorstdm; /* Multiple of standard deviation. */
-};
-
-void
-gal_statistics_mode_mirror_plots(float *sorted, size_t size,
- size_t mirrorindex, float min, float max,
- size_t numbins, char *histsname,
- char *cfpsname, float mirrorplotdist);
-
-float
-gal_statistics_mode_value_from_sym(float *sorted, size_t size,
- size_t modeindex, float sym);
-
-void
-gal_statistics_mode_index_in_sorted(float *sorted, size_t size,
- float errorstdm, size_t *modeindex,
- float *modesym);
-
-
__END_C_DECLS /* From C++ preparations */
#endif /* __GAL_STATISTICS_H__ */
diff --git a/lib/linkedlist.c b/lib/linkedlist.c
index 0e8ca06..b02b873 100644
--- a/lib/linkedlist.c
+++ b/lib/linkedlist.c
@@ -260,6 +260,28 @@ gal_linkedlist_dll_to_array(struct gal_linkedlist_dll
*list,
void
+gal_linkedlist_reverse_dll(struct gal_linkedlist_dll **list)
+{
+ double thisvalue;
+ struct gal_linkedlist_dll *correctorder=NULL;
+
+ /* Only do the reversal if there is more than one element. */
+ if( *list && (*list)->next )
+ {
+ while(*list!=NULL)
+ {
+ gal_linkedlist_pop_from_dll(list, &thisvalue);
+ gal_linkedlist_add_to_dll(&correctorder, thisvalue);
+ }
+ *list=correctorder;
+ }
+}
+
+
+
+
+
+void
gal_linkedlist_free_dll(struct gal_linkedlist_dll *list)
{
struct gal_linkedlist_dll *tmp, *ttmp;
diff --git a/lib/mode.c b/lib/mode.c
deleted file mode 100644
index ac37f4a..0000000
--- a/lib/mode.c
+++ /dev/null
@@ -1,402 +0,0 @@
-/*********************************************************************
-mode -- Find the mode of a distribution.
-This is part of GNU Astronomy Utilities (Gnuastro) package.
-
-Original author:
- Mohammad Akhlaghi <address@hidden>
-Contributing author(s):
-Copyright (C) 2015, Free Software Foundation, Inc.
-
-Gnuastro is free software: you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation, either version 3 of the License, or (at your
-option) any later version.
-
-Gnuastro is distributed in the hope that it will be useful, but
-WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
-**********************************************************************/
-#include <config.h>
-
-#include <math.h>
-#include <stdio.h>
-#include <errno.h>
-#include <error.h>
-#include <float.h>
-#include <stdlib.h>
-
-#include <gnuastro/array.h>
-#include <gnuastro/statistics.h>
-
-#include "mode.h"
-
-
-/****************************************************************
- ***************** Mode plots ********************
- ****************************************************************/
-/*
-This is the python code you can put in `plot.py` so the plot
-command mentioned here works:
-
------------------------------------------------------------------
-#! /usr/bin/env python3
-
-import sys
-import numpy as np
-import matplotlib.pyplot as plt
-
-a=np.loadtxt(sys.argv[1])
-b=np.loadtxt(sys.argv[2])
-w=a[1,0]-a[0,0]
-
-plt.bar(a[:,0], a[:,1], width=w, color="blue", linewidth=0, alpha=0.6)
-plt.bar(a[:,0], a[:,2], width=w, color="green", linewidth=0, alpha=0.4)
-plt.plot(b[:,0], b[:,1], linewidth=2, color="blue")
-plt.plot(b[:,0], b[:,2], linewidth=2, color="green")
-
-plt.ylim([0,np.amax(a[:,1])])
-plt.xlim([np.amin(a[:,0]),np.amax(a[:,0])])
-
-plt.savefig(sys.argv[3])
----------------------------------------------------------
-
-Run it with `./plot.py histsname.txt cfpsname.txt outputpdfname.pdf`
-It will plot the corresponding histograms and cumulative frequency
-plots. If you like to, this call is available as a system() call in
-the functions below.
-*/
-
-/* This is used for the plots, it will allocate an array and put the
- mirrored array values in it. `mi` is the index the mirror is to
- be placed on. */
-void
-makemirrored(float *in, size_t mi, float **outmirror, size_t *outsize)
-{
- size_t i, size;
- float *mirror, zf;
-
- zf=in[mi];
- size=2*mi+1;
-
- errno=0;
- mirror=malloc(size*sizeof *mirror);
- if(mirror==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for mirror array in "
- "makemirrored (mode.c)", size*sizeof *mirror);
-
- for(i=0;i<=mi;++i)
- mirror[i] = in[i];
- for(i=1;i<=mi;++i)
- mirror[mi+i] = 2*zf - mirror[mi-i];
-
- *outmirror=mirror;
- *outsize=size;
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/****************************************************************
- ***************** Mode ********************
- ****************************************************************/
-
-/*Find the index where the two CDFs dirverge beyond CDFSDIVERGEDIFF.
- T he input CDF and a CDF of a mirrored distribution about the
- quantile `mirrorquant` are compared.
-
- The basic idea behind finding the mode is comparing the mirrored CDF
- (about a test for the mode) with the actual CDF for a given
- point. This is the job of this function. It takes the ordered array
- and the quantile that is to be checked, it then finds the maximum
- difference between the mirrored CDF about the given point and the
- input CDF.
-
- `zf` keeps the flux at the mirror (zero) point. `i` is used to
- count the pixels before m. So `m+i` is the index of the mirrored
- distribution and mf=zf+(zf-a[m-i])=2*zf-a[m-i] is the mirrored flux
- at this point. `j` is found such that a[m+j] has the nearest flux to
- `mf`.
-
- The desired difference between the input CDF and the mirrored one
- for each `i` is then simply: `j-i`.
-
- Once `i` is incremented, `mf` will increase, so to find the new `j`
- we don't need to begin looking from `j=0`. Remember that the array
- is sorted, so the desired `j` is definitely larger than the previous
- `j`. So, if we keep the previous `j` in `prevj` then, all we have to
- do is to start incrementing `j` from `prevj`. This will really help
- in speeding up the job :-D. Only for the first element, `prevj=0`.*/
-size_t
-mirrormaxdiff(float *a, size_t size, size_t m,
- size_t numcheck, size_t interval, size_t stdm)
-{
- /* The variables:
- i: Index on mirror distribution.
- j: Index on input distribution.
- prevj: Index of previously checked point in the actual array.
- mf: Flux that is approximately equal in both distributions.*/
- float mf, zf;
- size_t maxdiff=0, errordiff;
- size_t i, j, absdiff, prevj=0;
-
- /* Find the error and mirror flux value */
- zf=a[m];
- errordiff=stdm*sqrt(m);
-
- /*
- printf("###############\n###############\n");
- printf("### Mirror pixel: %zu\n", m);
- printf("###############\n###############\n");
- */
- /* Go over the mirrored points. */
- for(i=1; i<numcheck && i<=m && m+i<size ;i+=interval)
- {
- mf=2*zf-a[m-i];
-
- /* The moment a[m+j]>mf, we have reached the last pixel to
- check. Now we just have to see if a[m+j-1] is closer to mf or
- if a[m+j]. We change `j` accordingly and break out of the `j`
- loop. */
- for(j=prevj;j<size-m;++j)
- if(a[m+j]>mf)
- {
- if( a[m+j]-mf < mf-a[m+j-1] )
- break;
- else
- {
- j--;
- break;
- }
- }
- /*
- printf("i:%-5zu j:%-5zu diff:%-5d maxdiff: %zu\n",
- i, j, (int)j-(int)i, maxdiff);
- */
- /* The index of the actual CDF corresponding the the mirrored
- flux has been found. We want the mirrored distribution to be
- within the actual distribution, not beyond it, so the only
- acceptable results are when i<j. If i>j+errordiff then the
- result is not acceptable! */
- if(i>j+errordiff)
- {
- maxdiff = MODE_MIRROR_IS_ABOVE_RESULT;
- break;
- }
- absdiff = i>j ? i-j : j-i;
- if(absdiff>maxdiff)
- maxdiff=absdiff;
-
- prevj=j;
- }
- return maxdiff;
-}
-
-
-
-
-
-/* Find the mode of a float array of size `size`. I assume that
- mirrormaxdiff() has one minimum (within the statistical errors) in the
- function. To find that minimum, the golden section search algorithm is
- going to used. Read the Wikipedia article for a very nice
- introduction. In summary we will constantly be finding middle points in
- the given interval and thus decreasing the interval until a certain
- tolerance is reached.
-
- If the input interval is on points `a` and `b`, then the middle point
- (lets call it `c`, where c>a and c<b) to test should be positioned such
- that (b-c)/(c-a)=GAL_MODE_GOLDEN_RATIO. Once we open up this relation,
- we can find c using:
-
- c=(b+GAL_MODE_GOLDEN_RATIO*a)/(1+GAL_MODE_GOLDEN_RATIO).
-
- We need a fourth point to be placed between. With this configuration,
- the probing point is located at: */
-size_t
-modegoldenselection(struct gal_statistics_mode_params *mp)
-{
- size_t di, dd;
- /*static int counter=1;*/
- /*------------------------------------------------------------------
- char outname[500], command[1000];
- char histsname[500], cfpsname[500];
- ------------------------------------------------------------------*/
-
- /* Find the probing point in the larger interval. */
- if(mp->highi-mp->midi > mp->midi-mp->lowi)
- di = mp->midi + MODE_TWO_TAKE_GOLDEN_RATIO *(float)(mp->highi-mp->midi);
- else
- di = mp->midi - MODE_TWO_TAKE_GOLDEN_RATIO * (float)(mp->midi-mp->lowi);
-
- /* Since these are all indexs (and positive) we don't need an
- absolute value, highi is also always larger than lowi! In some
- cases, the first (standard) condition might be satisfied, while
- highi-lowi<=2. In such cases, also jump out! */
- if( (mp->highi - mp->lowi) < mp->tolerance*(mp->midi+di)
- || (mp->highi - mp->lowi) <= 3)
- return (mp->highi+mp->lowi)/2;
-
- /* Find the maximum difference for this quantile. */
- dd = mirrormaxdiff(mp->sorted, mp->size, di, mp->numcheck,
- mp->interval, mp->errorstdm);
-
- /*------------------------------------------------------------------
- sprintf(outname, "%dcmp.pdf", counter);
- sprintf(cfpsname, "%dcfps.txt", counter);
- sprintf(histsname, "%dhists.txt", counter);
- gal_mode_make_mirror_plots(mp->sorted, mp->size, di, histsname, cfpsname);
- sprintf(command, "./plot.py %s %s %s", histsname, cfpsname, outname);
- system(command);
- -------------------------------------------------------------------*/
- /*
- printf("%-5zu\t%-5zu(%d)\t%-5zu ----> dq: %-5zu di: %d\n",
- mp->lowi, mp->midi, (int)mp->midd, mp->highi,
- di, (int)dd);
- */
- /* +++++++++++++ The mirrored distribution's cumulative frequency plot
- has be lower than the actual's cfp. If it isn't, `di` will be
- GAL_MODE_MIRROR_IS_ABOVE_RESULT. In this case, the normal golden
- section minimization is not going to give us what we want. So I have
- added this modification to it. In such cases, we want the search to go
- to the lower intervals.*/
- if(dd==MODE_MIRROR_IS_ABOVE_RESULT)
- {
- if(mp->midi < di)
- {
- mp->highi=di;
- return modegoldenselection(mp);
- }
- else
- {
- mp->highi=mp->midi;
- mp->midi=di;
- mp->midd=dd;
- return modegoldenselection(mp);
- }
- }
- /* +++++++++++++ End of my addition to the golden section search. */
-
- /* This is the standard golden section search: */
- if(dd<mp->midd)
- {
- if(mp->highi-mp->midi > mp->midi-mp->lowi)
- {
- mp->lowi = mp->midi;
- mp->midi = di;
- mp->midd = dd;
- return modegoldenselection(mp);
- }
- else
- {
- mp->highi = mp->midi;
- mp->midi = di;
- mp->midd = dd;
- return modegoldenselection(mp);
- }
- }
- else
- {
- if(mp->highi-mp->midi > mp->midi-mp->lowi)
- {
- mp->highi = di;
- return modegoldenselection(mp);
- }
- else
- {
- mp->lowi = di;
- return modegoldenselection(mp);
- }
- }
-}
-
-
-
-
-
-/* Once the mode is found, we need to do a quality control. This
- quality control is the measure of symmetricity. Lets assume the
- mode index is at `m`, the error in `m` can be assumed to be
- sqrt(m). So lets call the first point that the difference between
- the cumulative distribution of the mirror and actual data deviate
- above sqrt(m), is at index `b`. For a scale parameter, lets assume
- that the index of 5% of `m` is `a`. We could have taken the
- distribution minimum, but the scatter in that can be too high! Now
- the symmetricity of the mode can be quantified as: (b-m)/(m-a). For
- a completly symmetric mode, this should be 1. Note that the search
- for `b` only goes to the 95% of the distribution. */
-void
-modesymmetricity(float *a, size_t size, size_t mi, float errorstdm,
- float *sym)
-{
- float af, bf, mf, fi;
- size_t i, j, bi=0, topi, errdiff, prevj=0;
-
- mf=a[mi];
- errdiff=errorstdm*sqrt(mi);
- topi = 2*mi>size-1 ? size-1 : 2*mi;
- 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
- find the index where the difference between the two cumulative
- frequency plots exceeds that of the error in the mirror index. */
- for(i=1; i<topi-mi ;i+=1)
- {
- fi=2*mf-a[mi-i];
-
- for(j=prevj;j<size-mi;++j)
- if(a[mi+j]>fi)
- {
- if( a[mi+j]-fi < fi-a[mi+j-1] )
- break;
- else
- {
- j--;
- break;
- }
- }
-
- if(i>j+errdiff || j>i+errdiff)
- {
- bi=mi+i;
- break;
- }
- prevj=j;
- }
-
- /* bi==0 shows that no point with a larger difference could be
- found. So bi should be set to the end of the search region. */
- if(bi==0) bi=topi;
-
- bf=a[bi];
- /*
- printf("%f, %f, %f\n", af, mf, bf);
- */
-
- *sym=(bf-mf)/(mf-af);
-
- /* This is mainly used for plotting, which subtracts `mf`.
- printf("SymmetricFlux: %f\n", a[bi]-mf);
- printf("symmetricity: %f\n", sym);
- */
-}
diff --git a/lib/mode.h b/lib/mode.h
deleted file mode 100644
index 3838f4c..0000000
--- a/lib/mode.h
+++ /dev/null
@@ -1,46 +0,0 @@
-/*********************************************************************
-mode -- Find the mode of a distribution.
-This is part of GNU Astronomy Utilities (Gnuastro) package.
-
-Original author:
- Mohammad Akhlaghi <address@hidden>
-Contributing author(s):
-Copyright (C) 2015, Free Software Foundation, Inc.
-
-Gnuastro is free software: you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation, either version 3 of the License, or (at your
-option) any later version.
-
-Gnuastro is distributed in the hope that it will be useful, but
-WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
-**********************************************************************/
-#ifndef __MODE_H__
-#define __MODE_H__
-
-#include <gnuastro/statistics.h>
-
-#define MODE_GOLDEN_RATIO 1.618034f
-#define MODE_TWO_TAKE_GOLDEN_RATIO 0.38197f
-#define MODE_MIRROR_IS_ABOVE_RESULT (size_t)(-1)
-
-size_t
-mirrormaxdiff(float *a, size_t size, size_t m,
- size_t numcheck, size_t interval, size_t stdm);
-
-void
-modesymmetricity(float *a, size_t size, size_t mi, float errorstdm,
- float *sym);
-
-size_t
-modegoldenselection(struct gal_statistics_mode_params *mp);
-
-void
-makemirrored(float *in, size_t mi, float **outmirror, size_t *outsize);
-
-#endif
diff --git a/lib/options.c b/lib/options.c
index b66681f..0cd0c7c 100644
--- a/lib/options.c
+++ b/lib/options.c
@@ -270,7 +270,7 @@ gal_options_parse_list_of_numbers(char *string, char
*filename, size_t lineno)
out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &num, NULL, 0,
minmapsize, NULL, NULL, NULL);
for(tdll=list;tdll!=NULL;tdll=tdll->next)
- ((double *)out->array)[--i]=tdll->v;
+ ((double *)(out->array))[--i]=tdll->v;
/* Clean up and return. */
diff --git a/lib/statistics.c b/lib/statistics.c
index 7604418..ce64b1e 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -38,8 +38,6 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <checkset.h>
-#include "mode.h"
-
@@ -396,7 +394,6 @@ gal_statistics_median(gal_data_t *input, int inplace)
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);
@@ -408,46 +405,136 @@ gal_statistics_median(gal_data_t *input, int inplace)
+/* For a given size, return the index (starting from zero) that is at the
+ given quantile. */
+size_t
+gal_statistics_quantile_index(size_t size, double quant)
+{
+ double floatindex;
+
+ if(quant<0.0f || quant>1.0f)
+ error(EXIT_FAILURE, 0, "the quantile in `gal_statistics_quantile_index' "
+ "should be between 0.0 and 1.0 (inclusive). You have asked for "
+ "%g", quant);
+
+ /* Find the index of the quantile. */
+ floatindex=(double)(size-1)*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 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)
+gal_statistsics_quantile(gal_data_t *input, double quantile, int inplace)
{
- size_t dsize=1, ind;
+ size_t dsize=1, index;
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);
+ index=gal_statistics_quantile_index(nbs->size, quantile);
+
+ /* Read the value at this index into the output. */
+ gal_data_copy_element_same_type(nbs, index, out->array);
+
+ /* Clean up and return. */
+ if(nbs!=input) gal_data_free(nbs);
+ return out;
+}
+
+
+
+
+
+/* Return the index of the (first) point in the sorted dataset that has the
+ closest value to `value' (which has to be the same type as the `input'
+ dataset). */
+#define STATS_QFUNC(IT) { \
+ IT *r, *a=nbs->array, *af=a+nbs->size, v=*((IT *)(value->array)); \
+ \
+ /* For a reference. Since we are comparing with the previous. */ \
+ /* element, we need to start with the second element.*/ \
+ r=a++; \
+ \
+ /* Increasing array: */ \
+ if( *a < *(a+1) ) \
+ do if(*a>v) { if( v - *(a-1) < *a - v ) --a; break; } while(++a<af); \
+ \
+ /* Decreasing array. */ \
+ else \
+ do if(*a<v) { if( *(a-1) - v < v - *a ) --a; break; } while(++a<af); \
+ \
+ /* Set the difference. */ \
+ if(a<af) index=a-r; \
+ }
+size_t
+gal_statistics_quantile_function_index(gal_data_t *input, gal_data_t *value,
+ int inplace)
+{
+ size_t index=-1;
+ gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);
- /* Set the value. */
+ /* A sanity check. */
+ if(input->type!=value->type)
+ error(EXIT_FAILURE, 0, "the types of the input dataset and value to "
+ "`gal_statistics_quantile_function' have to be the same");
+
+ /* Find the result: */
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;
+ case GAL_DATA_TYPE_UINT8: STATS_QFUNC( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: STATS_QFUNC( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: STATS_QFUNC( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: STATS_QFUNC( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: STATS_QFUNC( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: STATS_QFUNC( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: STATS_QFUNC( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: STATS_QFUNC( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: STATS_QFUNC( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: STATS_QFUNC( double ); break;
default:
error(EXIT_FAILURE, 0, "type code %d not recognized in "
- "`gal_statistics_maximum'", input->type);
+ "`gal_statistics_quantile_function'", input->type);
}
/* Clean up and return. */
if(nbs!=input) gal_data_free(nbs);
+ return index;
+}
+
+
+
+
+
+/* Return the quantile function of the given value as float64. */
+gal_data_t *
+gal_statistics_quantile_function(gal_data_t *input, gal_data_t *value,
+ int inplace)
+{
+ double *d;
+ size_t dsize=1;
+ gal_data_t *out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &dsize,
+ NULL, 1, -1, NULL, NULL, NULL);
+ size_t ind=gal_statistics_quantile_function_index(input, value, inplace);
+
+ /* Note that counting of the index starts from 0, so for the quantile we
+ should divided by (size - 1). */
+ d=out->array;
+ d[0] = ( ind==-1 ? NAN : ((double)ind) / ((double)(input->size - 1)) );
return out;
}
@@ -455,33 +542,581 @@ gal_statistsics_quantile(gal_data_t *input, float
quantile, int inplace)
-size_t
-gal_statistics_quantile_index(size_t size, float quant)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*********************************************************************/
+/***************** Mode ***********************/
+/*********************************************************************/
+/* Main structure to keep mode parameters. */
+struct statistics_mode_params
{
- float floatindex;
+ gal_data_t *data; /* Sorted input dataset with no blank values. */
+ size_t lowi; /* Lower quantile of interval. */
+ size_t midi; /* Index of the mid-interval point. */
+ size_t midd; /* Maximum CDF distance at the middle point. */
+ size_t highi; /* Higher quantile of interval. */
+ float tolerance; /* Tolerance level to terminate search. */
+ size_t numcheck; /* Number of pixels after mode to check. */
+ size_t interval; /* Interval to check pixels. */
+ float mirrordist; /* Distance after mirror to check ( x STD). */
+};
+
- 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;
+
+
+/* Macros for the mode finding algorithm. */
+#define MODE_MIN_Q 0.01f /* Mode search lower interval quantile. */
+#define MODE_MAX_Q 0.90f /* Mode search higher interval quantile. */
+#define MODE_GOOD_LQ 0.02f /* Least acceptable mode quantile. */
+#define MODE_SYM_LOW_Q 0.01f /* Lower quantile to get symmetricity. */
+#define MODE_GOLDEN_RATIO 1.618034f /* Golden ratio: (1+sqrt(5))/2. */
+#define MODE_TWO_TAKE_GR 0.38197f /* 2 - Golden ratio. */
+#define MODE_MIRROR_ABOVE (size_t)(-1) /* Mirror is above the result. */
+
+
+
+
+/*
+ Given a mirror point (`m'), return the maximum distance between the
+ mirror distribution and the original distribution.
+
+ The basic idea behind finding the mode is comparing the mirrored CDF
+ (where the mirror is a test for the mode) with the original CDF for a
+ given point. The job of this function is to return the maximum distance,
+ given a mirror point. It takes the index of the mirror that is to be
+ checked, it then finds the maximum difference between the mirrored CDF
+ about the given point and the input CDF.
+
+ `zf` keeps the value at the mirror (zero) point. `i` is used to count
+ the pixels after the mirror in the mirror distribution. So `m+i` is the
+ index of the mirrored distribution and mf=zf+(zf-a[m-i])=2*zf-a[m-i] is
+ the mirrored flux at this point. Having found `mf', we find the `j` such
+ that a[m+j] has the nearest flux to `mf`.
+
+ The desired difference between the input CDF and the mirrored one
+ for each `i` is then simply: `j-i`.
+
+ Once `i` is incremented, `mf` will increase, so to find the new `j` we
+ don't need to begin looking from `j=0`. Remember that the array is
+ sorted, so the desired `j` is definitely larger than the previous
+ `j`. So, if we keep the previous `j` in `prevj` then, all we have to do
+ is to start incrementing `j` from `prevj`. This will really help in
+ speeding up the job :-D. Only for the first element, `prevj=0`. */
+#define MIRR_MAX_DIFF(IT) { \
+ IT *a=p->data->array, zf=a[m], mf=2*zf-a[m-i]; \
+ \
+ /* When a[m+j]>mf, we have reached the last pixel to check. Now, */ \
+ /* we just have to see which one of a[m+j-1] or a[m+j] is closer */ \
+ /* to `mf'. We then change `j` accordingly and break out of the */ \
+ /* `j' loop. */ \
+ for(j=prevj;j<size-m;++j) \
+ if(a[m+j]>mf) \
+ { \
+ if( a[m+j]-mf < mf-a[m+j-1] ) \
+ break; \
+ else \
+ { \
+ j--; \
+ break; \
+ } \
+ } \
+ }
+
+static size_t
+mode_mirror_max_index_diff(struct statistics_mode_params *p, size_t m)
+{
+ /* The variables:
+ i: Index on mirror distribution.
+ j: Index on input distribution.
+ prevj: Index of previously checked point in the actual array.
+ mf: (in macro) Value that is approximately equal in both
+ distributions. */
+ size_t i, j, absdiff, prevj=0, size=p->data->size;
+ size_t maxdiff=0, errordiff=p->mirrordist*sqrt(m);
/*
- printf("quant: %f, size: %zu, findex: %f\n", quant, size, floatindex);
+ printf("###############\n###############\n");
+ printf("### Mirror pixel: %zu\n", m);
+ printf("###############\n###############\n");
*/
- /* 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;
+ /* Go over the mirrored points. */
+ for(i=1; i<p->numcheck && i<=m && m+i<size ;i+=p->interval)
+ {
+ /* Find `j': the index of the closest point in the original
+ distribution that has a value similar to the mirror
+ distribution. */
+ switch(p->data->type)
+ {
+ case GAL_DATA_TYPE_UINT8: MIRR_MAX_DIFF( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: MIRR_MAX_DIFF( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: MIRR_MAX_DIFF( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: MIRR_MAX_DIFF( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: MIRR_MAX_DIFF( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: MIRR_MAX_DIFF( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: MIRR_MAX_DIFF( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: MIRR_MAX_DIFF( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: MIRR_MAX_DIFF( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: MIRR_MAX_DIFF( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`mode_mirror_max_diff'", p->data->type);
+ }
+
+ /*
+ printf("i:%-5zu j:%-5zu diff:%-5d maxdiff: %zu\n",
+ i, j, (int)j-(int)i, maxdiff);
+ */
+ /* The index of the actual CDF corresponding the the mirrored flux
+ has been found. We want the mirrored distribution to be within the
+ actual distribution, not beyond it, so the only acceptable results
+ are when i<j. But we also have noise, so we can't simply use that
+ as the criterion, small `j's with `i>j' are acceptable. So, only
+ when `i>j+errordiff' the result is not acceptable! */
+ if(i>j+errordiff)
+ {
+ maxdiff = MODE_MIRROR_ABOVE;
+ break;
+ }
+ absdiff = i>j ? i-j : j-i;
+ if(absdiff>maxdiff) maxdiff=absdiff;
+
+ prevj=j;
+ }
+
+ /* Return the maximum difference */
+ return maxdiff;
+}
+
+
+
+
+
+/* Find the mode through the Golden-section search. It is assumed that
+ `mode_mirror_max_index_diff' has one minimum (within the statistical
+ errors) in the function. To find that minimum, the golden section search
+ algorithm is going to used. Read the Wikipedia article for a very nice
+ introduction.
+
+ In summary we will constantly be finding middle points in the given
+ interval and thus decreasing the interval until a certain tolerance is
+ reached.
+
+ If the input interval is on points `a' and `b', then the middle point
+ (lets call it `c', where c>a and c<b) to test should be positioned such
+ that (b-c)/(c-a)=MODE_GOLDEN_RATIO. Once we open up this relation, we
+ can find c using:
+
+ c = ( b + MODE_GOLDEN_RATIO * a ) / ( 1 + MODE_GOLDEN_RATIO )
+
+ We need a fourth point to be placed between. With this configuration,
+ the probing point is located at: */
+static size_t
+mode_golden_section(struct statistics_mode_params *p)
+{
+ size_t di, dd;
+
+ /* Find the probing point in the larger interval. */
+ if(p->highi-p->midi > p->midi-p->lowi)
+ di = p->midi + MODE_TWO_TAKE_GR * (float)(p->highi-p->midi);
else
- return floatindex;
+ di = p->midi - MODE_TWO_TAKE_GR * (float)(p->midi-p->lowi);
+
+ /* Since these are all indexs (and positive) we don't need an absolute
+ value, highi is also always larger than lowi! In some cases, the first
+ (standard) condition might be satisfied, while highi-lowi<=2. In such
+ cases, also jump out! */
+ if( (p->highi - p->lowi) < p->tolerance*(p->midi+di)
+ || (p->highi - p->lowi) <= 3)
+ return (p->highi+p->lowi)/2;
+
+ /* Find the maximum difference for this mirror point. */
+ dd = mode_mirror_max_index_diff(p, di);
+
+ /*------------------------------------------------------------------
+ {
+ static int counter=1;
+ char outname[500], command[1000];
+ char histsname[500], cfpsname[500];
+ sprintf(outname, "%dcmp.pdf", counter);
+ sprintf(cfpsname, "%dcfps.txt", counter);
+ sprintf(histsname, "%dhists.txt", counter);
+ gal_mode_make_mirror_plots(p->sorted, p->size, di, histsname, cfpsname);
+ sprintf(command, "./plot.py %s %s %s", histsname, cfpsname, outname);
+ system(command);
+ }
+ -------------------------------------------------------------------*/
+
+ /*
+ printf("%-5zu\t%-5zu(%d)\t%-5zu ----> dq: %-5zu di: %d\n",
+ p->lowi, p->midi, (int)p->midd, p->highi,
+ di, (int)dd);
+ */
+
+ /* +++++++++++++ Start of addition to the golden section search.
+
+ The mirrored distribution's cumulative frequency plot has be lower
+ than the actual's cfp. If it isn't, `di` will be MODE_MIRROR_ABOVE. In
+ this case, the normal golden section minimization is not going to give
+ us what we want. So we have this modification. In such cases, we want
+ the search to go to the lower interval. */
+ if(dd==MODE_MIRROR_ABOVE)
+ {
+ if( p->midi < di )
+ {
+ p->highi=di;
+ return mode_golden_section(p);
+ }
+ else
+ {
+ p->highi=p->midi;
+ p->midi=di;
+ p->midd=dd;
+ return mode_golden_section(p);
+ }
+ }
+ /* End of addition to the golden section search. +++++++++++++*/
+
+ /* This is the standard golden section search: */
+ if(dd<p->midd)
+ {
+ if(p->highi-p->midi > p->midi-p->lowi)
+ {
+ p->lowi = p->midi;
+ p->midi = di;
+ p->midd = dd;
+ return mode_golden_section(p);
+ }
+ else
+ {
+ p->highi = p->midi;
+ p->midi = di;
+ p->midd = dd;
+ return mode_golden_section(p);
+ }
+ }
+ else
+ {
+ if(p->highi-p->midi > p->midi-p->lowi)
+ {
+ p->highi = di;
+ return mode_golden_section(p);
+ }
+ else
+ {
+ p->lowi = di;
+ return mode_golden_section(p);
+ }
+ }
}
+/* Once the mode is found, we need to do a quality control. This quality
+ control is the measure of its symmetricity. Lets assume the mode index
+ is at `m', since an index is just a count, from the Poisson
+ distribution, the error in `m' is sqrt(m).
+
+ Now, let's take `b' to be the first point that the difference between
+ the cumulative distribution of the mirror and actual data deviate more
+ than sqrt(m). For a scale parameter, lets assume that the index of 5% of
+ `m` is `a`. We could have taken the distribution minimum, but the
+ scatter in the minimum can be too high!
+
+ Now, the "symmetricity" of the mode can be defined as: (b-m)/(m-a). For
+ a completly symmetric mode, this should be 1. Note that the search for
+ `b` only goes to the 95% of the distribution. */
+#define MODE_SYM(IT) { \
+ IT *a=p->data->array, af, bf, mf, fi; \
+ \
+ /* Set the values at the mirror and at `a' (see above). */ \
+ mf=a[m]; \
+ af=a[ gal_statistics_quantile_index(2*m+1, MODE_SYM_LOW_Q) ]; \
+ \
+ /* This loop is very similar to that of */ \
+ /* `mode_mirror_max_index_diff'. It will find the index where the */\
+ /* difference between the two cumulative frequency plots exceeds */ \
+ /* that of the error in the mirror index.*/ \
+ for(i=1; i<topi-m ;i+=1) \
+ { \
+ fi=2*mf-a[m-i]; \
+ \
+ for(j=prevj;j<size-m;++j) \
+ if(a[m+j]>fi) \
+ { \
+ if( a[m+j]-fi < fi-a[m+j-1] ) \
+ break; \
+ else \
+ { \
+ j--; \
+ break; \
+ } \
+ } \
+ \
+ if(i>j+errdiff || j>i+errdiff) \
+ { \
+ bi=m+i; \
+ break; \
+ } \
+ prevj=j; \
+ } \
+ \
+ /* bi==0 shows that no point with a larger difference could be */ \
+ /* found. So bi should be set to the end of the search region. */ \
+ if(bi==0) bi=topi; \
+ \
+ bf = *(IT *)b_val = a[bi]; \
+ /* printf("%f, %f, %f\n", af, mf, bf); */ \
+ \
+ return (bf-mf)/(mf-af); \
+ }
+static double
+mode_symmetricity(struct statistics_mode_params *p, size_t m, void *b_val)
+{
+ size_t i, j, bi=0, topi, errdiff, prevj=0, size=p->data->size;
+
+ /* Set the basic constants. */
+ topi = 2*m>size-1 ? size-1 : 2*m;
+ errdiff = p->mirrordist * sqrt(m);
+
+ /* Do the process. */
+ switch(p->data->type)
+ {
+ case GAL_DATA_TYPE_UINT8: MODE_SYM( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: MODE_SYM( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: MODE_SYM( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: MODE_SYM( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: MODE_SYM( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: MODE_SYM( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: MODE_SYM( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: MODE_SYM( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: MODE_SYM( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: MODE_SYM( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`mode_symmetricity'", p->data->type);
+ }
+
+ /* Control shouldn't reach here! */
+ error(EXIT_FAILURE, 0, "a bug! please contact us at %s so we can address "
+ "the problem. For some reason control has reached the end of "
+ "`mode_symmetricity', this should not have happened",
+ PACKAGE_BUGREPORT);
+ return NAN;
+}
+
+
+
+
+
+/* Return the mode and related parameters in a float64 `gal_data_t' with
+ the following elements in its array, the array:
+
+ array[0]: mode
+ array[1]: mode quantile.
+ array[2]: symmetricity.
+ array[3]: value at the end of symmetricity.
+
+ The inputs are:
+
+ - `input' is the input dataset, it doesn't have to be sorted and can
+ have blank values.
+
+ - `mirrordist' is the maximum distance after the mirror point to check
+ as a multiple of sigma.
+
+ - `inplace' is either 0 or 1. If it is 1 and the input array has blank
+ values and is not sorted, then the removal of blank values and
+ sorting will occur in-place (input will be modified): all blank
+ elements in the input array will be removed and it will be sorted. */
+gal_data_t *
+gal_statistics_mode(gal_data_t *input, float mirrordist, int inplace)
+{
+ double *oa;
+ size_t modeindex;
+ size_t dsize=4, mdsize=1;
+ struct statistics_mode_params p;
+ gal_data_t *mode=gal_data_alloc(NULL, input->type, 1, &mdsize,
+ NULL, 1, -1, NULL, NULL, NULL);
+ gal_data_t *b_val=gal_data_alloc(NULL, input->type, 1, &mdsize,
+ NULL, 1, -1, NULL, NULL, NULL);
+ gal_data_t *out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &dsize,
+ NULL, 1, -1, NULL, NULL, NULL);
+
+ /* Make sure the input doesn't have blank values and is sorted. */
+ p.data=gal_statistics_no_blank_sorted(input, inplace);
+
+
+ /* Basic constants. */
+ p.tolerance = 0.01;
+ p.mirrordist = mirrordist;
+ p.numcheck = p.data->size/2;
+
+
+ /* Fill in the interval: Checking every single element is over-kill, so
+ if the dataset is large enough, we'll set an interval to only check
+ elements at an interval (so only 1000 elements are checked). */
+ p.interval = p.numcheck>1000 ? p.numcheck/1000 : 1;
+
+
+ /* Set the lower and higher acceptable indexes for the mode based on
+ quantiles. */
+ p.lowi = gal_statistics_quantile_index(p.data->size, MODE_MIN_Q);
+ p.highi = gal_statistics_quantile_index(p.data->size, MODE_MAX_Q);
+
+
+ /* Having set the low and higher interval values, we will set the first
+ middle point and also the maximum distance on that point. This is
+ necessary to start the iteration. */
+ p.midi = ( ( (float)p.highi + MODE_GOLDEN_RATIO * (float)p.lowi )
+ / ( 1 + MODE_GOLDEN_RATIO ) );
+ p.midd = mode_mirror_max_index_diff(&p, p.midi);
+
+
+ /* Do the golden-section search iteration, read the mode value from the
+ input array and save it as the first element of the output dataset's
+ array, then free the `mode' structure. */
+ oa=out->array;
+ modeindex = mode_golden_section(&p);
+ gal_data_copy_element_same_type(p.data, modeindex, mode->array);
+ mode=gal_data_copy_to_new_type_free(mode, GAL_DATA_TYPE_FLOAT64);
+ gal_data_copy_element_same_type(mode, 0, &oa[0]);
+
+
+ /* Put in the rest of the values of the output structure. */
+ oa[1] = ((double)modeindex) / ((double)(p.data->size-1));
+ oa[2] = mode_symmetricity(&p, modeindex, b_val->array);
+
+
+ /* Put the end of the symmetricity in. */
+ b_val=gal_data_copy_to_new_type_free(mode, GAL_DATA_TYPE_FLOAT64);
+ gal_data_copy_element_same_type(mode, 0, &oa[3]);
+
+
+ /* For a check:
+ printf("mode: %g\nquantile: %g\nsymmetricity: %g\nsym value: %g\n",
+ oa[0], oa[1], oa[2], oa[3]);
+ */
+
+ /* Clean up (if necessary), then return the output */
+ if(p.data!=input) gal_data_free(p.data);
+ gal_data_free(mode);
+ return out;
+}
+
+
+
+
+
+/* Make the mirror array. */
+#define STATS_MKMIRROR(IT) { \
+ IT *a=noblank_sorted->array, *m=mirror->array; \
+ IT zf=a[index]; \
+ *mirror_val=zf; \
+ for(i=0;i<=index;++i) m[i] = a[i]; \
+ for(i=1;i<=index;++i) m[index+i] = 2 * zf - m[index - i]; \
+ }
+static gal_data_t *
+statistics_make_mirror(gal_data_t *noblank_sorted, size_t index,
+ double *mirror_val)
+{
+ size_t i, dsize = 2*index+1;
+ gal_data_t *mirror=gal_data_alloc(NULL, noblank_sorted->type, 1, &dsize,
+ NULL, 1, -1, NULL, NULL, NULL);
+
+ /* Make sure the index is less than or equal to the number of
+ elements. */
+ if( index >= noblank_sorted->size )
+ error(EXIT_FAILURE, 0, "the index value to `statistics_make_mirror' "
+ "must be less than or equal to the number of elements in the "
+ "input, but it isn't: index: %zu, size of input: %zu", index,
+ noblank_sorted->size);
+
+ /* Fill in the mirror array. */
+ switch(noblank_sorted->type)
+ {
+ case GAL_DATA_TYPE_UINT8: STATS_MKMIRROR( uint8_t ); break;
+ case GAL_DATA_TYPE_INT8: STATS_MKMIRROR( int8_t ); break;
+ case GAL_DATA_TYPE_UINT16: STATS_MKMIRROR( uint16_t ); break;
+ case GAL_DATA_TYPE_INT16: STATS_MKMIRROR( int16_t ); break;
+ case GAL_DATA_TYPE_UINT32: STATS_MKMIRROR( uint32_t ); break;
+ case GAL_DATA_TYPE_INT32: STATS_MKMIRROR( int32_t ); break;
+ case GAL_DATA_TYPE_UINT64: STATS_MKMIRROR( uint64_t ); break;
+ case GAL_DATA_TYPE_INT64: STATS_MKMIRROR( int64_t ); break;
+ case GAL_DATA_TYPE_FLOAT32: STATS_MKMIRROR( float ); break;
+ case GAL_DATA_TYPE_FLOAT64: STATS_MKMIRROR( double ); break;
+ }
+
+ /* Return the mirrored distribution. */
+ return mirror;
+}
+
+
+
+
+
+/* Make a mirrored histogram and cumulative frequency plot with the mirror
+ distribution of the input with a value at `value'.
+
+ The output is a linked list of data structures: the first is the bins
+ with on bin at the mirror point, the second is the histogram with a
+ maximum of one and the third is the cumulative frequency plot. */
+gal_data_t *
+gal_statistics_mode_mirror_plots(gal_data_t *input, gal_data_t *value,
+ size_t numbins, int inplace,
+ double *mirror_val)
+{
+ gal_data_t *mirror, *bins, *hist, *cfp;
+ gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);
+ size_t ind=gal_statistics_quantile_function_index(nbs, value, inplace);
+
+
+ /* If the given mirror was outside the range of the input, then index
+ will be 0 (below the range) or -1 (above the range), in that case, we
+ should return NULL. */
+ if(ind==-1 || ind==0)
+ return NULL;
+
+
+ /* Make the mirror array. */
+ mirror=statistics_make_mirror(nbs, ind, mirror_val);
+
+
+ /* Set the bins for histogram and cdf. */
+ bins=gal_statistics_regular_bins(mirror, NULL, numbins, *mirror_val);
+
+
+ /* Make the histogram: set it's maximum value to 1 for a nice comparison
+ with the CDF. */
+ hist=gal_statistics_histogram(mirror, bins, 0, 1);
+
+
+ /* Make the cumulative frequency plot. */
+ cfp=gal_statistics_cfp(mirror, bins, 1);
+
+
+ /* Set the pointers to make a table and return. */
+ bins->next=hist;
+ hist->next=cfp;
+ return bins;
+}
+
+
@@ -737,7 +1372,7 @@ gal_statistics_no_blank_sorted(gal_data_t *input, int
inplace)
*/
gal_data_t *
gal_statistics_regular_bins(gal_data_t *data, gal_data_t *inrange,
- size_t numbins, float onebinstart)
+ size_t numbins, double onebinstart)
{
size_t i;
gal_data_t *bins, *tmp, *range;
@@ -825,18 +1460,19 @@ gal_statistics_regular_bins(gal_data_t *data, gal_data_t
*inrange,
if( (b[i]-hbw) < onebinstart && (b[i+1]-hbw) > onebinstart) break;
if( i != numbins-1 )
{
- diff=onebinstart-b[i];
+ diff = onebinstart - (b[i]-hbw);
for(i=0;i<numbins;++i)
- b[i]+=diff;
+ b[i] += diff;
}
}
- /* For a check
+ /* For a check:
printf("min: %g\n", min);
printf("max: %g\n", max);
+ printf("onebinstart: %.10f\n", onebinstart);
printf("binwidth: %g\n", binwidth);
for(i=0;i<numbins;++i)
- printf("%zu: %.4f\n", i, b[i]);
+ printf("%zu: %.4f\t(%f, %f)\n", i, b[i], b[i]-hbw, b[i]+hbw);
*/
/* Set the status of the bins to regular and return. */
@@ -942,7 +1578,7 @@ gal_statistics_histogram(gal_data_t *data, gal_data_t
*bins, int normalize,
free(hist->name); free(hist->unit); free(hist->comment);
gal_checkset_allocate_copy("hist_normalized", &hist->name);
gal_checkset_allocate_copy("frac", &hist->unit);
- gal_checkset_allocate_copy("Normalized histogram value for this bin",
+ gal_checkset_allocate_copy("Normalized histogram value for this bin.",
&hist->comment);
}
if(maxone)
@@ -958,7 +1594,7 @@ gal_statistics_histogram(gal_data_t *data, gal_data_t
*bins, int normalize,
gal_checkset_allocate_copy("hist_maxone", &hist->name);
gal_checkset_allocate_copy("frac", &hist->unit);
gal_checkset_allocate_copy("Fractional histogram value for this bin "
- "when maximum bin value is 1.0",
+ "when maximum bin value is 1.0.",
&hist->comment);
}
@@ -1106,7 +1742,7 @@ gal_statistics_cfp(gal_data_t *data, gal_data_t *bins,
int normalize)
/****************************************************************
- ***************** Sigma clip ********************
+ ***************** Outliers ********************
****************************************************************/
/* Return a data structure with an array of four values: the final number
of points used, the median, average and standard deviation. The number
@@ -1163,9 +1799,9 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
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;
+ size_t maxnum = param>=1.0f ? param : GAL_STATISTICS_SIG_CLIP_MAX_CONVERGE;
/* Some sanity checks. */
if( multip<=0 )
@@ -1320,24 +1956,6 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/****************************************************************/
-/************* Identify outliers ****************/
-/****************************************************************/
/* Using the cumulative distribution function this funciton will
remove outliers from a dataset. */
void
@@ -1397,218 +2015,3 @@ gal_statistics_remove_outliers_flat_cdf(float *sorted,
size_t *outsize)
free(slopes);
#endif
}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/****************************************************************/
-/************* Mode calculation ****************/
-/****************************************************************/
-void
-gal_statistics_mode_mirror_plots(float *sorted, size_t size, size_t
mirrorindex,
- float min, float max, size_t numbins,
- char *histsname, char *cfpsname,
- float mirrorplotdist)
-{
- printf("\n... in gal_statistics_mode_mirror_plots ...\n");
- exit(1);
-#if 0
- FILE *fp;
- size_t i, msize;
- float *out, maxhist=-FLT_MAX, maxcfp, d;
- int normhist=0, maxhistone=0, normcfp=0;
- float *bins, *mirror, *actual, mf, onebinvalue=0.0f;
-
-
- /* Find the index of the mirror and make the mirror array: */
- mf=sorted[mirrorindex];
- gal_array_float_copy(sorted, size, &actual);
- makemirrored(sorted, mirrorindex, &mirror, &msize);
-
-
- /* Set the best range if asked for, such that the mirror is on the
- 1/3 of the image scale. */
- if(mirrorplotdist!=0.0f)
- {
- min=actual[gal_statistics_index_from_quantile(size, 0.001f)];
- max=mf+mirrorplotdist*(mf-min);
- }
-
-
- /* set the mirror pixel to have the value of zero.*/
- min-=mf;
- max-=mf;
- gal_array_fsum_const(actual, size, -1.0f*mf);
- gal_array_fsum_const(mirror, msize, -1.0f*mf);
-
-
- /* Allocate space for the 3 column array keeping the histograms:*/
- errno=0;
- out=malloc(numbins*3*sizeof *out);
- if(out==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for out in "
- "gal_mode_make_mirror_plots (mode.c)", numbins*3*sizeof *out);
-
-
- /* Define the bin sides: */
- gal_statistics_set_bins(actual, size, numbins, min,
- max, onebinvalue, 0, &bins);
-
-
- /* Find the histogram of the actual data and put it in out. Note
- that maxhistone=0, because here we want to use one value for both
- histograms so they are comparable. */
- gal_statistics_histogram(actual, size, bins, numbins,
- normhist, maxhistone);
- for(i=0;i<numbins;++i)
- if(bins[i*2+1]>maxhist) maxhist=bins[i*2+1];
- for(i=0;i<numbins;++i)
- { out[i*3]=bins[i*2]; out[i*3+1]=bins[i*2+1]/maxhist; bins[i*2+1]=0.0f;}
- bins[i*2+1]=0.0f; /* bins[] actually has numbins+1 elements. */
- d=(out[3]-out[0])/2;
-
-
- /* Find the histogram of the mirrired distribution and put it in
- out: */
- gal_statistics_histogram(mirror, msize, bins, numbins, normhist,
- maxhistone);
- for(i=0;i<numbins;++i)
- { out[i*3+2]=bins[i*2+1]/maxhist; bins[i*2+1]=0.0f;}
- bins[i*2+1]=0.0f; /* bins[] actually has numbins+1 elements. */
-
-
- /* Print out the histogram: */
- errno=0;
- fp=fopen(histsname, "w");
- if(fp==NULL)
- error(EXIT_FAILURE, errno, "could not open file %s", histsname);
- fprintf(fp, "# Histogram of actual and mirrored distributions.\n");
- fprintf(fp, "# Column 0: Value in the middle of this bin.\n");
- fprintf(fp, "# Column 1: Input data.\n");
- fprintf(fp, "# Column 2: Mirror distribution.\n");
- for(i=0;i<numbins;++i)
- fprintf(fp, "%-25.6f%-25.6f%-25.6f\n", out[i*3]+d,
- out[i*3+1], out[i*3+2]);
- fclose(fp);
-
-
-
-
- /* Find the cumulative frequency plots of the two distributions: */
- gal_statistics_cumulative_fp(actual, size, bins, numbins, normcfp);
- for(i=0;i<numbins;++i)
- { out[i*3+1]=bins[i*2+1]; bins[i*2+1]=0.0f; }
- bins[i*2+1]=0.0f; /* bins[] actually has numbins+1 elements. */
- gal_statistics_cumulative_fp(mirror, msize, bins, numbins, normcfp);
- for(i=0;i<numbins;++i)
- { out[i*3+2]=bins[i*2+1]; bins[i*2+1]=0.0f;}
- bins[i*2+1]=0.0f; /* bins[] actually has numbins+1 elements. */
-
-
- /* Since the two cumultiave frequency plots have to be on scale, see
- which one is larger and divided both CFPs by the size of the
- larger one. Then print the CFPs. */
- if(size>msize) maxcfp=size;
- else maxcfp=msize;
- errno=0;
- fp=fopen(cfpsname, "w");
- if(fp==NULL)
- error(EXIT_FAILURE, errno, "could not open file %s", cfpsname);
- fprintf(fp, "# Cumulative frequency plot (average index in bin) of\n"
- "# Actual and mirrored distributions.\n");
- fprintf(fp, "# Column 0: Value in the middle of this bin.\n");
- fprintf(fp, "# Column 1: Actual data.\n");
- fprintf(fp, "# Column 2: Mirror distribution.\n");
- for(i=0;i<numbins;++i)
- fprintf(fp, "%-25.6f%-25.6f%-25.6f\n", out[i*3],
- out[i*3+1]/maxcfp, out[i*3+2]/maxcfp);
- fclose(fp);
-
-
- /* Clean up. */
- free(out);
- free(bins);
- free(mirror);
- free(actual);
-#endif
-}
-
-
-
-
-
-/* It happens that you have the symmetricity and you want the flux
- value at that point, this function will do that job. In practice,
- it just finds bf from the equation to calculate symmetricity in
- modesymmetricity. */
-float
-gal_statistics_mode_value_from_sym(float *sorted, size_t size,
- size_t modeindex, float sym)
-{
- float mf=sorted[modeindex];
- float af=
- sorted[gal_statistics_quantile_index(2*modeindex+1,
- GAL_STATISTICS_MODE_SYMMETRICITY_LOW_QUANT)];
- return sym*(mf-af)+mf;
-}
-
-
-
-
-
-/* Find the quantile of the mode of a sorted distribution. The return
- value is either 0 (not accurate) or 1 (accurate). Accuracy is
- defined based on the difference between the maximum and minimum
- maxdiffs that were found during the golden section search.
-
- A good mode will have:
-
- modequant=(float)(modeindex)/(float)size;
- modesym>GAL_MODE_SYM_GOOD && modequant>GAL_MODE_LOW_QUANT_GOOD
-*/
-void
-gal_statistics_mode_index_in_sorted(float *sorted, size_t size, float
errorstdm,
- size_t *modeindex, float *modesym)
-{
- struct gal_statistics_mode_params mp;
-
- /* Initialize the gal_mode_params structure: */
- mp.size=size;
- mp.sorted=sorted;
- mp.tolerance=0.01;
- mp.numcheck=size/2;
- mp.errorstdm=errorstdm;
- if(mp.numcheck>1000)
- mp.interval=mp.numcheck/1000;
- else mp.interval=1;
- mp.lowi = gal_statistics_quantile_index(size,
- GAL_STATISTICS_MODE_LOW_QUANTILE);
- mp.highi = gal_statistics_quantile_index(size,
- GAL_STATISTICS_MODE_HIGH_QUANTILE);
- mp.midi = (((float)mp.highi
- + MODE_GOLDEN_RATIO*(float)mp.lowi)
- /(1+MODE_GOLDEN_RATIO));
- mp.midd = mirrormaxdiff(mp.sorted, mp.size, mp.midi, mp.numcheck,
- mp.interval, mp.errorstdm);
-
- /* Do the golden section search and find the resulting
- symmetricity. */
- *modeindex=modegoldenselection(&mp);
- modesymmetricity(sorted, size, *modeindex, errorstdm, modesym);
-}
- [gnuastro-commits] master 38f401b 121/125: Clump S/N threshold is now calculated on Sky clumps, (continued)
- [gnuastro-commits] master 38f401b 121/125: Clump S/N threshold is now calculated on Sky clumps, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 1e961f5 067/125: Starting to use gal_data_t in MakeProfiles, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master e58ed29 108/125: One value per tile in multi-channel tessellation, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 0e344fe 086/125: Warp now uses gal_data_t and all its fruits, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 8348f1e 119/125: NoiseChisel pseudo-detections identified on threads, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master fca5e80 117/125: First re-implementation of NoiseChisel's complete detection, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master a1ff97b 111/125: Permutations and neighbor finding algorithms, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 0bc96c9 122/125: Clumps and their S/N tables over detections, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 80cdd4f 118/125: Correction of convolution over channel borders, Mohammad Akhlaghi, 2017/04/23
- [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 <=
- [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, 2017/04/23
- [gnuastro-commits] master 15c524d 124/125: New implementation of NoiseChisel complete, Mohammad Akhlaghi, 2017/04/23