[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 97cb586b: Library (fit.h): least squares fitti
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 97cb586b: Library (fit.h): least squares fitting, accessible via Statistics |
|
Date: |
Wed, 28 Sep 2022 14:11:09 -0400 (EDT) |
branch: master
commit 97cb586b0c1e0533ce716aed10e6f73a4b6f7a12
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Library (fit.h): least squares fitting, accessible via Statistics
Until now, Gnuastro didn't have methods for least squares fitting (which is
usually necessary in astronomical data analysis). Fortunately the GNU
Scientific Library (GSL) has some powerful tools for doing this, so it was
just necessary to add a high-level interface to those functions.
With this commit, a high-level interface to some of GSL's linear least
squares fitting functions have been added in the new 'fit.h' library, and
an even higher-level interface has been written in the Statistics program
(calling Gnuastro's 'fit.h' library).
This task was proposed by Sepideh Eskandarlou.
---
NEWS | 44 ++
bin/statistics/args.h | 115 +++-
bin/statistics/aststatistics.conf | 5 +
bin/statistics/main.h | 22 +-
bin/statistics/statistics.c | 619 ++++++++++++++++++-
bin/statistics/ui.c | 474 ++++++++++-----
bin/statistics/ui.h | 11 +-
bin/table/arithmetic.c | 34 +-
doc/gnuastro.texi | 844 ++++++++++++++++++++++----
lib/Makefile.am | 2 +
lib/arithmetic.c | 77 +--
lib/blank.c | 10 +-
lib/fit.c | 667 ++++++++++++++++++++
lib/fits.c | 10 +-
lib/gnuastro/fit.h | 126 ++++
lib/gnuastro/jpeg.h | 2 +-
lib/options.c | 21 +-
lib/txt.c | 15 +-
tests/Makefile.am | 9 +-
tests/statistics/fitting-data.txt | 63 ++
tests/statistics/fitting-polynomial-robust.sh | 59 ++
21 files changed, 2840 insertions(+), 389 deletions(-)
diff --git a/NEWS b/NEWS
index 3def0eb2..bb466aba 100644
--- a/NEWS
+++ b/NEWS
@@ -77,6 +77,33 @@ See the end of the file for license conditions.
surface brightness image using the new 'counts-to-sb' operator of
Arithmetic.
+ Statistics:
+ - Linear and Polynomial least squares fitting are now available and very
+ easy to call on the command-line. They are wrappers over the respective
+ least squares fitting functions of the GNU Scientific Library. The
+ interface is pretty simple, like the example below:
+ aststatistics table.fits -cX,Y,Yerr --fit=linear-weighted
+ It is also possible to estimate values and errors of the fitted model
+ on a new X axis. A complete example has been added to the newly added
+ "Least squares fitting" section of the book (under the Statistics
+ program documentation). Please see that tutorial to easily use this
+ feature. The following new options have been added to the Statistics
+ program for this purpose:
+ --fit: the model to use. Currently the following models are supported:
+ linear
+ linear-weighted
+ linear-no-constant
+ linear-no-constant-weighted
+ polynomial
+ polynomial-weighted
+ polynomial-robust
+ --fitweight: nature of the "weight" column (default: standard dev).
+ --fitmaxpower: maximum power of X in polynomial models.
+ --fitrobust: weight function to use in the "robust" polynomial model.
+ --fitestimate: File name, or number to estimate the fit on.
+ --fitestimatehdu: HDU containing table in file given to '--fitestimate'.
+ --fitestimatecol: Column containing X axis values for '--fitestimate'.
+
astscript-fits-view:
--ds9colorbarmulti: show a separate color-bar for each image in DS9. By
default this script will show a single color-bar for all the images to
@@ -138,6 +165,16 @@ See the end of the file for license conditions.
- gal_dimension_collapse_sigclip_number: collapse with num. after sig-clip.
- gal_eps_shape_id_to_name: return the name of a shape from its ID.
- gal_eps_shape_name_to_id: return the ID of a shape from its name.
+ - gal_fit_name_to_id: Convert string name to ID of fitting model.
+ - gal_fit_name_from_id: Convert ID of fitting model to string name.
+ - gal_fit_name_robust_to_id: Convert name of robust weights to ID.
+ - gal_fit_name_robust_from_id: Convert ID of robust weights to name.
+ - gal_fit_1d_linear: linear fit of input columns.
+ - gal_fit_1d_linear_no_constant: linear fit with no constant.
+ - gal_fit_1d_linear_estimate: estimate a linear fit on a new X column.
+ - gal_fit_1d_polynomial: polynomial fit of input columns.
+ - gal_fit_1d_polynomial_robust: robust polynomial fit of input columns.
+ - gal_fit_1d_polynomial_estimate: estimate a polynomial fit on new X column.
- gal_fits_unique_keyvalues: extract all unique values to a certain
keyword in many files.
- gal_fits_with_keyvalue: select FITS image with a certain key value.
@@ -154,6 +191,13 @@ See the end of the file for license conditions.
** Removed features
+ Statistics:
+ --refcol has been removed because it breaks the modularity principle
+ (given that it is the job of Gnuastro's Table program to limit rows
+ from a larger table based on many different criteria). The output of
+ Table can be directly piped to Statistics to achieve the same (and
+ much more feature-rich effect).
+
** Changed features
Book:
diff --git a/bin/statistics/args.h b/bin/statistics/args.h
index 8ad9b2e7..9a22e980 100644
--- a/bin/statistics/args.h
+++ b/bin/statistics/args.h
@@ -44,19 +44,6 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET
},
- {
- "refcol",
- UI_KEY_REFCOL,
- "STR",
- 0,
- "Reference column name or number.",
- GAL_OPTIONS_GROUP_INPUT,
- &p->refcol,
- GAL_TYPE_STRING,
- GAL_OPTIONS_RANGE_ANY,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
{
"greaterequal",
UI_KEY_GREATEREQUAL,
@@ -815,6 +802,108 @@ struct argp_option program_options[] =
},
+ {
+ 0, 0, 0, 0,
+ "Fitting (1 dimensional)",
+ UI_GROUP_FIT
+ },
+ {
+ "fit",
+ UI_KEY_FIT,
+ "STR",
+ 0,
+ "Type of fitting: 'linear', 'linear-weighted', "
+ "'linear-no-constant', 'linear-no-constant-weighted', "
+ "'polynomial', 'polynomial-weighted', "
+ "'polynomial-robust'.",
+ UI_GROUP_FIT,
+ &p->fitname,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "fitweight",
+ UI_KEY_FITWEIGHT,
+ "STR",
+ 0,
+ "Input weight: 'std', 'var' or 'inv-variance'.",
+ UI_GROUP_FIT,
+ &p->fitweight,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "fitmaxpower",
+ UI_KEY_FITMAXPOWER,
+ "INT",
+ 0,
+ "Maximum power in polynomial to fit.",
+ UI_GROUP_FIT,
+ &p->fitmaxpower,
+ GAL_TYPE_SIZE_T,
+ GAL_OPTIONS_RANGE_GE_0,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "fitrobust",
+ UI_KEY_FITROBUST,
+ "STR",
+ 0,
+ "The robust function name to use.",
+ UI_GROUP_FIT,
+ &p->fitrobustname,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "fitestimate",
+ UI_KEY_FITESTIMATE,
+ "STR/FLT",
+ 0,
+ "Estimate fitted func. (number, file or 'self').",
+ UI_GROUP_FIT,
+ &p->fitestimate,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "fitestimatehdu",
+ UI_KEY_FITESTIMATEHDU,
+ "STR/INT",
+ 0,
+ "HDU containing the --fitestimate values.",
+ UI_GROUP_FIT,
+ &p->fitestimatehdu,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "fitestimatecol",
+ UI_KEY_FITESTIMATECOL,
+ "STR/INT",
+ 0,
+ "Column containing the --fitestimate values.",
+ UI_GROUP_FIT,
+ &p->fitestimatecol,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+
+
+
{0}
};
diff --git a/bin/statistics/aststatistics.conf
b/bin/statistics/aststatistics.conf
index ffe7ccec..7b5bdcc9 100644
--- a/bin/statistics/aststatistics.conf
+++ b/bin/statistics/aststatistics.conf
@@ -35,3 +35,8 @@
numbins 100
numbins2 100
mirrordist 1.5
+
+# Fitting
+ fitweight std
+ fitestimatehdu 1
+ fitrobust bisquare
diff --git a/bin/statistics/main.h b/bin/statistics/main.h
index 4781ea9a..ab591f51 100644
--- a/bin/statistics/main.h
+++ b/bin/statistics/main.h
@@ -46,6 +46,15 @@ enum statistics_input_format
INPUT_FORMAT_IMAGE,
};
+enum statistics_fit_weight_types
+ {
+ STATISTICS_FIT_WHT_INVALID,
+
+ STATISTICS_FIT_WHT_STD,
+ STATISTICS_FIT_WHT_VAR,
+ STATISTICS_FIT_WHT_INVVAR,
+ };
+
@@ -59,7 +68,6 @@ struct statisticsparams
gal_list_f64_t *tp_args; /* Arguments for printing. */
char *inputname; /* Input filename. */
gal_list_str_t *columns; /* Column name or number if input is table. */
- char *refcol; /* Reference column name or number. */
float greaterequal; /* Only use values >= this value. */
float greaterequal2; /* Only use values >= this value (2D hist). */
float lessthan; /* Only use values < this value. */
@@ -68,6 +76,13 @@ struct statisticsparams
float quantmax; /* Quantile maximum. */
uint8_t ontile; /* Do single value calculations on tiles. */
uint8_t interpolate; /* Use interpolation to fill blank tiles. */
+ char *fitname; /* Name of fitting function to use. */
+ char *fitweight; /* Input weight is 'std' or 'invvar'. */
+ char *fitestimate; /* Values to estimate the fit. */
+ char *fitestimatecol; /* Column for values to estimate the fit. */
+ char *fitestimatehdu; /* HDU for values to estimate the fit. */
+ size_t fitmaxpower; /* Maximum power of polynomial fit. */
+ char *fitrobustname; /* Name of robust function. */
uint8_t asciihist; /* Print an ASCII histogram. */
uint8_t asciicfp; /* Print an ASCII cumulative frequency plot.*/
@@ -107,7 +122,6 @@ struct statisticsparams
uint8_t needssort; /* If sorting is needed. */
gal_data_t *input; /* Input data structure. */
gal_data_t *sorted; /* Sorted input data structure. */
- gal_data_t *reference; /* Reference data structure. */
int isfits; /* Input is a FITS file. */
int hdu_type; /* Type of HDU (image or table). */
gal_data_t *kernel; /* Kernel for convolution of input for Sky. */
@@ -116,6 +130,10 @@ struct statisticsparams
gal_data_t *std_t; /* Sky standard deviation on each tile. */
char *checkskyname; /* Name of file for Sky calculation steps. */
time_t rawtime; /* Starting time of the program. */
+ uint8_t fitid; /* ID of desired fit. */
+ uint8_t fitrobustid; /* ID of robust fit type. */
+ gal_data_t *fitestval; /* Values to estimate over fit. */
+ int fitwhtid; /* Code for the nature of the weight column.*/
};
#endif
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index 978fc18c..348c1da3 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -30,6 +30,9 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <string.h>
#include <stdlib.h>
+#include <gsl/gsl_fit.h>
+
+#include <gnuastro/fit.h>
#include <gnuastro/wcs.h>
#include <gnuastro/fits.h>
#include <gnuastro/tile.h>
@@ -980,11 +983,7 @@ print_input_info(struct statisticsparams *p)
}
if(str)
{
- printf("Range: ");
- if(p->refcol)
- printf("[on column %s] ",
- p->reference->name ? p->reference->name : p->refcol);
- printf("%s.\n", str);
+ printf("Range: %s.\n", str);
free(str);
}
@@ -1076,8 +1075,8 @@ print_basics(struct statisticsparams *p)
p->numasciibins = p->numasciibins ? p->numasciibins : 70;
bins=gal_statistics_regular_bins(p->input, range, p->numasciibins, NAN);
hist=gal_statistics_histogram(p->input, bins, 0, 0);
- if(p->refcol==NULL) printf("\nHistogram:\n");
- print_ascii_plot(p, hist, bins, 1, p->refcol ? 1 : 0);
+ printf("\nHistogram:\n");
+ print_ascii_plot(p, hist, bins, 1, 0);
gal_data_free(bins);
gal_data_free(hist);
gal_data_free(range);
@@ -1175,6 +1174,605 @@ print_sigma_clip(struct statisticsparams *p)
+
+
+/*******************************************************************/
+/************** Fitting ***************/
+/*******************************************************************/
+static struct gal_fits_list_key_t *
+statistics_fit_params_to_keys(struct statisticsparams *p, gal_data_t *fit,
+ char *whtnat, double *redchisq)
+{
+ size_t i, j;
+ char *kname, *kcomm;
+ struct gal_fits_list_key_t *out=NULL;
+ double *c=fit->array, *cov=fit->next?fit->next->array:NULL;
+
+ /* Set the title and basic info (independent of the type of fit). */
+ gal_fits_key_list_title_add(&out, "Fit results", 0);
+ gal_fits_key_list_add(&out, GAL_TYPE_STRING, "FITTYPE", 0,
+ gal_fit_name_from_id(p->fitid), 0,
+ "Functional form of the fitting.", 0, NULL, 0);
+ if( p->fitid==GAL_FIT_POLYNOMIAL
+ || p->fitid==GAL_FIT_POLYNOMIAL_ROBUST
+ || p->fitid==GAL_FIT_POLYNOMIAL_WEIGHTED )
+ gal_fits_key_list_add(&out, GAL_TYPE_SIZE_T, "FITMAXP", 0,
+ &p->fitmaxpower, 0,
+ "Maximum power of polynomial.", 0, NULL, 0);
+ if( p->fitid==GAL_FIT_POLYNOMIAL_ROBUST )
+ gal_fits_key_list_add(&out, GAL_TYPE_STRING, "FITRTYP", 0,
+ p->fitrobustname, 0,
+ "Function for removing outliers", 0, NULL, 0);
+ gal_fits_key_list_add(&out, GAL_TYPE_STRING, "FITIN", 0,
+ p->inputname, 0,"Name of file with input columns.",
+ 0, NULL, 0);
+ if(p->isfits)
+ gal_fits_key_list_add(&out, GAL_TYPE_STRING, "FITINHDU", 0,
+ p->cp.hdu, 0,"Name or Number of HDU with input "
+ "columns.", 0, NULL, 0);
+ gal_fits_key_list_add(&out, GAL_TYPE_STRING, "FITXCOL", 0,
+ p->columns->v, 0,"Name or Number of independent "
+ "(X) column.", 0, NULL, 0);
+ gal_fits_key_list_add(&out, GAL_TYPE_STRING, "FITYCOL", 0,
+ p->columns->next->v, 0,"Name or Number of "
+ "measured (Y) column.", 0, NULL, 0);
+ if(p->columns->next->next)
+ {
+ gal_fits_key_list_add(&out, GAL_TYPE_STRING, "FITWCOL", 0,
+ p->columns->next->next->v, 0,
+ "Name or Number of weight column.", 0, NULL, 0);
+ gal_fits_key_list_add(&out, GAL_TYPE_STRING, "FITWNAT", 0,
+ whtnat, 0, "Nature of weight column.",
+ 0, NULL, 0);
+ }
+ if(p->fitid==GAL_FIT_POLYNOMIAL_ROBUST)
+ gal_fits_key_list_add(&out, GAL_TYPE_STRING, "FITROBST", 0,
+ p->fitrobustname, 0,
+ "Robust fitting (rejecting outliers) function.",
+ 0, NULL, 0);
+ gal_fits_key_list_add(&out, GAL_TYPE_FLOAT64, "FRDCHISQ", 0,
+ redchisq, 0, "Reduced chi^2 of fit.", 0, NULL, 0);
+
+ /* Add the Fitting results. */
+ switch(p->fitid)
+ {
+ /* Linear with no constant */
+ case GAL_FIT_LINEAR_NO_CONSTANT:
+ case GAL_FIT_LINEAR_NO_CONSTANT_WEIGHTED:
+ gal_fits_key_list_add(&out, GAL_TYPE_FLOAT64, "FITC1", 0,
+ c, 0, "C1: Multiple of X in linear fit "
+ "(y=C1*x).", 0, NULL, 0);
+ gal_fits_key_list_add(&out, GAL_TYPE_FLOAT64, "FCOV11", 0,
+ c+1, 0, "Variance of C1 (only element of "
+ "cov. matrix).", 0, NULL, 0);
+ break;
+
+ /* Basic linear. */
+ case GAL_FIT_LINEAR:
+ case GAL_FIT_LINEAR_WEIGHTED:
+ gal_fits_key_list_add(&out, GAL_TYPE_FLOAT64, "FITC0", 0,
+ c, 0, "C0: Constant in linear fit "
+ "(y=C0+C1*x).", 0, NULL, 0);
+ gal_fits_key_list_add(&out, GAL_TYPE_FLOAT64, "FITC1", 0,
+ c+1, 0, "C1: Multiple of X in linear fit "
+ "(y=C0+C1*x).", 0, NULL, 0);
+ gal_fits_key_list_add(&out, GAL_TYPE_FLOAT64, "FCOV11", 0,
+ c+2, 0, "Covariance matrix element (1,1).",
+ 0, NULL, 0);
+ gal_fits_key_list_add(&out, GAL_TYPE_FLOAT64, "FCOV12", 0,
+ c+3, 0, "Covariance matrix element "
+ "(1,2)=(2,1).", 0, NULL, 0);
+ gal_fits_key_list_add(&out, GAL_TYPE_FLOAT64, "FCOV22", 0,
+ c+4, 0, "Covariance matrix element (2,2).",
+ 0, NULL, 0);
+ break;
+
+ /* Polynomial fit. */
+ case GAL_FIT_POLYNOMIAL:
+ case GAL_FIT_POLYNOMIAL_ROBUST:
+ case GAL_FIT_POLYNOMIAL_WEIGHTED:
+ for(i=0;i<fit->size;++i)
+ {
+ if( asprintf(&kname, "FITC%zu", i)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf in FITCxx name",
+ __func__);
+ if( asprintf(&kcomm, "C%zu: multiple of x^%zu in polynomial",
+ i, i)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf in FITCxx comment",
+ __func__);
+ gal_fits_key_list_add(&out, GAL_TYPE_FLOAT64, kname, 1,
+ c+i, 0, kcomm, 1, NULL, 0);
+ }
+ for(i=0;i<fit->size;++i)
+ for(j=0;j<fit->size;++j)
+ {
+ if( asprintf(&kname, "FCOV%zu%zu", i+1, j+1)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf in FCOVxx name",
+ __func__);
+ if( asprintf(&kcomm, "Covariance matrix element (%zu,%zu).",
+ i+1, j+1)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf in FCOVxx comment",
+ __func__);
+ gal_fits_key_list_add(&out, GAL_TYPE_FLOAT64, kname, 1,
+ cov+(i*fit->size+j), 0, kcomm, 1,
+ NULL, 0);
+ }
+ break;
+
+ /* Unrecognized FIT ID. */
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to "
+ "fix the problem. The code '%d' isn't recognized for 'fitid'",
+ __func__, PACKAGE_BUGREPORT, p->fitid);
+ }
+
+ /* Reverse the last-in-first-out list to be in the same logical order we
+ inserted the items here. */
+ gal_fits_key_list_reverse(&out);
+
+ /* Return the key list. */
+ return out;
+}
+
+
+
+
+
+static void
+statistics_fit_estimate(struct statisticsparams *p, gal_data_t *fit,
+ char *whtnat, double *redchisq)
+{
+ gal_data_t *est=NULL;
+ double *x, *y, *yerr;
+ struct gal_fits_list_key_t *keys=NULL;
+
+ /* If the input had no metadata, add them. */
+ if(p->fitestval->name==NULL)
+ gal_checkset_allocate_copy("X-INPUT", &p->fitestval->name);
+ if(p->fitestval->comment==NULL)
+ gal_checkset_allocate_copy("Requested values to estimate fit.",
+ &p->fitestval->comment);
+
+ /* Estimations are done on a per-row level. */
+ switch(p->fitid)
+ {
+ case GAL_FIT_LINEAR:
+ case GAL_FIT_LINEAR_WEIGHTED:
+ case GAL_FIT_LINEAR_NO_CONSTANT:
+ case GAL_FIT_LINEAR_NO_CONSTANT_WEIGHTED:
+ est=gal_fit_1d_linear_estimate(fit, p->fitestval);
+ break;
+
+ case GAL_FIT_POLYNOMIAL:
+ case GAL_FIT_POLYNOMIAL_ROBUST:
+ case GAL_FIT_POLYNOMIAL_WEIGHTED:
+ est=gal_fit_1d_polynomial_estimate(fit, p->fitestval);
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to "
+ "fix the problem. The code '%d' isn't recognized for 'fitid'",
+ __func__, PACKAGE_BUGREPORT, p->fitid);
+ }
+
+ /* Set the estimated columns to be after the input's columns. */
+ p->fitestval->next=est;
+
+ /* Non-quiet title. */
+ if(p->cp.quiet==0)
+ printf("\nRequested estimation:\n");
+
+ /* If only one value was estimated and no column was given, then the
+ user wanted the value on the command-line. */
+ if(p->fitestimatecol)
+ {
+ if(p->cp.output)
+ {
+ if(p->cp.quiet==0)
+ printf(" Written to: %s\n", p->cp.output);
+ }
+ keys=statistics_fit_params_to_keys(p, fit, whtnat, redchisq);
+ gal_table_write(p->fitestval, &keys, NULL, p->cp.tableformat,
+ p->cp.output, "FIT_ESTIMATE", 0);
+ }
+
+ /* Print estimated value on the commandline. */
+ else
+ {
+ x=p->fitestval->array;
+ y=p->fitestval->next->array;
+ yerr=p->fitestval->next->next->array;
+ if(p->cp.quiet)
+ printf("%f %f %f\n", x[0], y[0], yerr[0]);
+ else
+ printf(" X: %f (given on command-line)\n"
+ " Y: %f\n"
+ " Y_error: %f\n", x[0], y[0], yerr[0]);
+ }
+
+ /* Clean up. */
+ gal_list_data_free(p->fitestval);
+ p->fitestval=NULL; /* Was freed with 'out', avoid generic freeing later. */
+}
+
+
+
+
+
+static char *
+statistics_fit_whtnat(struct statisticsparams *p)
+{
+ switch(p->fitwhtid)
+ {
+ case STATISTICS_FIT_WHT_STD: return "Standard deviation";
+ case STATISTICS_FIT_WHT_VAR: return "Variance";
+ case STATISTICS_FIT_WHT_INVVAR: return "Inverse variance";
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' "
+ "to find and fix the problem. The value '%d' isn't a "
+ "recognized weight type identifier", __func__,
+ PACKAGE_BUGREPORT, p->fitwhtid);
+ }
+
+ /* Control should not reach here. */
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' "
+ "to find and fix the problem. Control should not reach the "
+ "end of this function", __func__, PACKAGE_BUGREPORT);
+ return NULL;
+}
+
+
+
+
+
+static char *
+statistics_fit_print_intro(struct statisticsparams *p, char **whtnat)
+{
+ char *colspace;
+ char *filename, *intro, *wcolstr=NULL;
+ gal_list_str_t *xn=p->columns, *yn=xn->next?xn->next:NULL,
+ *wn=yn->next?yn->next:NULL;
+
+ /* Set the full file name (for easy reading later!).*/
+ filename=gal_fits_name_save_as_string(p->inputname, p->cp.hdu);
+
+ /* Prepare string for nature of weight (if any weight was given!). */
+ *whtnat = p->input->next->next ? statistics_fit_whtnat(p) : NULL;
+
+ /* Set the Weight column string(s). */
+ if( p->fitid==GAL_FIT_LINEAR_WEIGHTED
+ || p->fitid==GAL_FIT_POLYNOMIAL_WEIGHTED
+ || p->fitid==GAL_FIT_LINEAR_NO_CONSTANT_WEIGHTED )
+ {
+ colspace=" ";
+ if( asprintf(&wcolstr, " Weight column: %s "
+ "[%s of Y in each row]\n", wn->v, *whtnat)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else colspace=" ";
+
+ /* Put everything into one string. */
+ if( asprintf(&intro,
+ "%s\n"
+ "-------\n"
+ "Fitting results (remove extra info with '--quiet' "
+ "or '-q)\n"
+ " Input file: %s with %zu non-blank rows.\n"
+ " X%scolumn: %s\n"
+ " Y%scolumn: %s\n"
+ "%s",
+ PROGRAM_STRING, filename, p->input->size,
+ colspace, xn->v, colspace, yn->v,
+ wcolstr ? wcolstr : "")<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+
+ /* Clean up and return. */
+ free(filename);
+ free(wcolstr);
+ return intro;
+}
+
+
+
+
+
+static int
+statistics_fit_linear(struct statisticsparams *p)
+{
+ double *f, redchisq;
+ gal_data_t *fit=NULL;
+ char *intro, *funcvals, *whtnat=NULL;
+
+ /* These have to be defined after each other. */
+ gal_data_t *x=p->input;
+ gal_data_t *y=x->next?x->next:NULL, *w=y->next?y->next:NULL;
+
+ /* If the required number of columns aren't given, */
+ switch(p->fitid)
+ {
+
+ /* Linear fit. */
+ case GAL_FIT_LINEAR:
+ if(gal_list_data_number(p->input)!=2) return 2;
+ fit=gal_fit_1d_linear(x, y, NULL);
+ break;
+
+ /* Linear (weighted). */
+ case GAL_FIT_LINEAR_WEIGHTED:
+ if(gal_list_data_number(p->input)!=3) return 3;
+ fit=gal_fit_1d_linear(x, y, w);
+ break;
+
+ /* Linear (no constant) */
+ case GAL_FIT_LINEAR_NO_CONSTANT:
+ if(gal_list_data_number(p->input)!=2) return 2;
+ fit=gal_fit_1d_linear_no_constant(x, y, NULL);
+ break;
+
+ /* Linear (no constant, weighted) */
+ case GAL_FIT_LINEAR_NO_CONSTANT_WEIGHTED:
+ if(gal_list_data_number(p->input)!=3) return 3;
+ fit=gal_fit_1d_linear_no_constant(x, y, w);
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to "
+ "fix the problem. '%d' isn't recognized as a fitting ID",
+ __func__, PACKAGE_BUGREPORT, p->fitid);
+ }
+
+ /* Print the output. */
+ f=fit->array;
+ if(p->cp.quiet)
+ {
+ if(p->fitestval) whtnat=statistics_fit_whtnat(p);
+ else
+ {
+ /* The covariance matrix will only be present if a constant was
+ present. After it, just print the residual (chi-squared or sum
+ of squares of residuals). */
+ if( p->fitid==GAL_FIT_LINEAR
+ || p->fitid==GAL_FIT_LINEAR_WEIGHTED )
+ printf("%+-.10f %+-.10f\n" /* The Two coefficients. */
+ "%+-20.10f %+-20.10f\n"/* First row of cov. matrix. */
+ "%+-20.10f %+-20.10f\n"/* Second row of cov. matrix. */
+ "%+-.10f\n", /* Residual. */
+ f[0], f[1], f[2], f[3], f[3], f[4], f[5]);
+ else
+ printf("%+-.10f\n%+-.10f\n%+-.10f\n", f[0], f[1], f[2]);
+ }
+ }
+ else
+ {
+ /* With or without constant. */
+ if( p->fitid==GAL_FIT_LINEAR
+ || p->fitid==GAL_FIT_LINEAR_WEIGHTED )
+ {
+ redchisq=f[5];
+ if( asprintf(&funcvals,
+ "Fit function: Y = c0 + (c1 * X)\n"
+ " c0: %+-.10f\n"
+ " c1: %+-.10f\n\n"
+ "Covariance matrix (off-diagonal are identical "
+ "same):\n"
+ " %+-20.10f %+-20.10f\n"
+ " %+-20.10f %+-20.10f\n", f[0], f[1], f[2], f[3],
+ f[3], f[4])<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ redchisq=f[2];
+ if( asprintf(&funcvals,
+ "Fit function: Y = c1 * X\n"
+ " c1: %+-.10f\n\n"
+ "Variance of 'c1':\n"
+ " %+-.10f\n", f[0], f[1])<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+
+ /* Statistics version, and input filenames. */
+ intro=statistics_fit_print_intro(p, &whtnat);
+
+ /* Final printed report.*/
+ printf("%s\n%s\nReduced chi^2 of fit:\n %+f\n", intro,
+ funcvals, redchisq);
+
+ /* Clean up. */
+ free(intro);
+ free(funcvals);
+ }
+
+ /* Estimate values (if requested), note that it involves writing the
+ fitted parameters in the header. */
+ if(p->fitestval)
+ statistics_fit_estimate(p, fit, whtnat, fit->size==6?f+5:f+2);
+
+ /* Clean up. */
+ gal_data_free(fit);
+ return 0; /* Means everything is good! */
+}
+
+
+
+
+
+static void
+statistics_fit_polynomial_print(struct statisticsparams *p, gal_data_t *fit,
+ double redchisq, char **whtnat)
+{
+ size_t i, j;
+ char *intro;
+ size_t nconst=p->fitmaxpower+1;
+ double *farr=fit->array, *carr=fit->next->array;
+
+ /* Print fitted constants */
+ if(p->cp.quiet)
+ {
+ if(p->fitestval==NULL)
+ {
+ for(i=0;i<nconst;++i) printf("%+-.10f ", farr[i]);
+ printf("\b\n");
+ }
+ }
+ else
+ {
+ /* Statistics version, and input filenames. */
+ intro=statistics_fit_print_intro(p, whtnat);
+
+ /* Final printed report.*/
+ printf("%s\n"
+ "Fit function: Y = c0 + (c1 * X^1) + (c2 * X^2) "
+ "+ ... (cN * X^N)\n", intro);
+
+ /* Notice for the (possible) robust function and maximum power of
+ polynomial. */
+ if(p->fitid==GAL_FIT_POLYNOMIAL_ROBUST)
+ printf(" Robust function: %s\n", p->fitrobustname);
+ printf(" N: %zu\n", p->fitmaxpower);
+
+ /* Print the fitted values. */
+ for(i=0;i<nconst;++i)
+ printf(" c%zu: %s%+-.10f\n", i, i<10?" ":"", farr[i]);
+
+ /* Print the information on the covariance matrix. */
+ printf("\nCovariance matrix:\n");
+
+ /* Clean up. */
+ free(intro);
+ }
+
+ /* Print the covariance matrix (when no estimation is requested, this is
+ the same for quiet or non-quiet mode). But when estimation is
+ requested, they will be written in the FITS keywords of the output so
+ there is no more need to have them on the command-line.*/
+ if(p->cp.quiet==0 || p->fitestval==NULL)
+ {
+ for(i=0;i<nconst;++i)
+ {
+ if(p->cp.quiet==0) printf(" ");
+ for(j=0;j<nconst;++j)
+ printf("%+-20.10f ", carr[i*nconst+j]);
+ printf("\b\n");
+ }
+
+ /* Print the chi^2. */
+ if(p->cp.quiet==0)
+ printf("\nReduced chi^2 of fit:\n");
+ printf("%s%+-.10f\n", p->cp.quiet?"":" ", redchisq);
+ }
+}
+
+
+
+
+
+static int
+statistics_fit_polynomial(struct statisticsparams *p)
+{
+ char *whtnat;
+ gal_data_t *fit=NULL;
+ double redchisq=NAN; /* Important to initialize. */
+ gal_data_t *x=p->input, *y=x->next?x->next:NULL, *w=y->next?y->next:NULL;
+
+ /* Sanity check and call the fitting functions. */
+ switch(p->fitid)
+ {
+ case GAL_FIT_POLYNOMIAL:
+ if(gal_list_data_number(p->input)!=2) return 2; /* Sanity check. */
+ fit=gal_fit_1d_polynomial(x, y, NULL, p->fitmaxpower, &redchisq);
+ break;
+
+ case GAL_FIT_POLYNOMIAL_ROBUST:
+ if(gal_list_data_number(p->input)!=2) return 2; /* Sanity check. */
+ fit=gal_fit_1d_polynomial_robust(x, y, p->fitmaxpower,
+ p->fitrobustid, &redchisq);
+ break;
+
+ case GAL_FIT_POLYNOMIAL_WEIGHTED:
+ if(gal_list_data_number(p->input)!=3) return 3; /* Sanity check. */
+ fit=gal_fit_1d_polynomial(x, y, w, p->fitmaxpower, &redchisq);
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to "
+ "fix the problem. '%d' isn't recognized as a fitting ID",
+ __func__, PACKAGE_BUGREPORT, p->fitid);
+ }
+
+ /* Print the output. */
+ statistics_fit_polynomial_print(p, fit, redchisq, &whtnat);
+
+ /* Estimate values (if requested), note that it involves writing the
+ fitted parameters in the header. */
+ if(p->fitestval)
+ statistics_fit_estimate(p, fit, whtnat, &redchisq);
+
+ /* Clean up. */
+ gal_data_free(fit);
+ gal_list_data_free(p->fitestval); p->fitestval=NULL;
+ return 0; /* Means everything is good! */
+}
+
+
+
+
+
+static void
+statistics_fit(struct statisticsparams *p)
+{
+ size_t neededcols=0;
+
+ /* Make sure that at least two columns are provided. */
+ if(p->input->next==NULL)
+ error(EXIT_FAILURE, 0, "at least two columns are necessary for "
+ "the fitting operations");
+
+ /* Do the fitting. */
+ switch(p->fitid)
+ {
+ case GAL_FIT_LINEAR:
+ case GAL_FIT_LINEAR_WEIGHTED:
+ case GAL_FIT_LINEAR_NO_CONSTANT:
+ case GAL_FIT_LINEAR_NO_CONSTANT_WEIGHTED:
+ neededcols=statistics_fit_linear(p);
+ break;
+
+ case GAL_FIT_POLYNOMIAL:
+ case GAL_FIT_POLYNOMIAL_ROBUST:
+ case GAL_FIT_POLYNOMIAL_WEIGHTED:
+ neededcols=statistics_fit_polynomial(p);
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to "
+ "fix the problem. '%s' is not a recognized as a fit type",
+ __func__, PACKAGE_BUGREPORT, p->fitname);
+ }
+
+ /* If the number of columns is not sufficient, then 'neededcols' will be
+ non-zero. In this case, we should abort and inform the user. */
+ if(neededcols)
+ error(EXIT_FAILURE, 0, "'%s' fitting requires %zu columns as input, "
+ "but %zu columns have been given", p->fitname, neededcols,
+ gal_list_data_number(p->input));
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
/*******************************************************************/
/************** Main function ***************/
/*******************************************************************/
@@ -1239,6 +1837,13 @@ statistics(struct statisticsparams *p)
print_mirror_hist_cfp(p);
}
+ /* Fitting. */
+ if( p->fitname )
+ {
+ print_basic_info=0;
+ statistics_fit(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 e2ee5347..d39a8ac6 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -28,6 +28,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <stdio.h>
#include <string.h>
+#include <gnuastro/fit.h>
#include <gnuastro/txt.h>
#include <gnuastro/wcs.h>
#include <gnuastro/fits.h>
@@ -73,11 +74,12 @@ const char
doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" will do statistical "
"analysis on the input dataset (table column or image). All blank "
"pixels or pixels outside of the given range are ignored. You can "
- "either directly ask for certain statistics in one line/row as shown "
- "below with the same order as requested, or get tables of different "
- "statistical measures like the histogram, cumulative frequency style "
- "and etc. If no particular statistic is requested, some basic "
- "information about the dataset is printed on the command-line.\n"
+ "either directly ask for certain statistics in one line/row as "
+ "shown below with the same order as requested, or get tables of "
+ "different statistical measures like the histogram, cumulative "
+ "frequency style and etc. If no particular statistic is "
+ "requested, some basic information about the dataset is printed "
+ "on the command-line.\n"
GAL_STRINGS_MORE_HELP_INFO
/* After the list of options: */
"\v"
@@ -140,6 +142,7 @@ ui_initialize_options(struct statisticsparams *p,
p->meanmedqdiff = NAN;
p->sclipparams[0] = NAN;
p->sclipparams[1] = NAN;
+ p->fitmaxpower = GAL_BLANK_SIZE_T;
/* Set the mandatory common options. */
for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
@@ -178,10 +181,10 @@ parse_opt(int key, char *arg, struct argp_state *state)
check if the first character of arg is the equal sign, then the
user is warned and the program is stopped: */
if(arg && arg[0]=='=')
- argp_error(state, "incorrect use of the equal sign ('='). For short "
- "options, '=' should not be used and for long options, "
- "there should be no space between the option, equal sign "
- "and value");
+ argp_error(state, "incorrect use of the equal sign ('='). For "
+ "short options, '=' should not be used and for long "
+ "options, there should be no space between the "
+ "option, equal sign and value");
/* Set the key to this option. */
switch(key)
@@ -192,7 +195,8 @@ parse_opt(int key, char *arg, struct argp_state *state)
'arg' will be an empty string! We don't want to account for such
cases (and give a clear error that no input has been given). */
if(p->inputname)
- argp_error(state, "only one argument (input file) should be given");
+ argp_error(state, "only one argument (input file) should be "
+ "given");
else
if(arg[0]!='\0') p->inputname=arg;
break;
@@ -220,19 +224,19 @@ ui_add_to_single_value(struct argp_option *option, char
*arg,
/* 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 "
- "with the '--printparams' ('-P'), or writing into configuration "
- "files due to lack of time when implementing these features. "
- "You can put them into configuration files manually. Please get "
- "in touch with us at '%s', so we can implement it",
- PACKAGE_BUGREPORT);
+ error(EXIT_FAILURE, 0, "currently the options to be printed in one "
+ "row (like '--number', '--mean', and etc) do not support "
+ "printing with the '--printparams' ('-P'), or writing into "
+ "configuration files due to lack of time when implementing "
+ "these features. You can put them into configuration files "
+ "manually. Please get in touch with us at '%s', so we can "
+ "implement it", PACKAGE_BUGREPORT);
/* 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 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'. */
@@ -271,8 +275,8 @@ ui_add_to_single_value(struct argp_option *option, char
*arg,
if(option->key==UI_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);
+ "and 1, you had asked for %g (read from "
+ "'%s')", d[i], arg);
gal_list_f64_add(&p->tp_args, d[i]);
gal_list_i32_add(&p->singlevalue, option->key);
}
@@ -282,8 +286,8 @@ ui_add_to_single_value(struct argp_option *option, char
*arg,
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);
+ "marked as requiring a value, but is not "
+ "recognized", PACKAGE_BUGREPORT);
}
}
@@ -403,12 +407,13 @@ ui_read_check_only_options(struct statisticsparams *p)
{
/* The tile or sky mode cannot be called with any other modes. */
if( p->asciihist || p->asciicfp || p->histogram || p->histogram2d
- || p->cumulative || p->sigmaclip || !isnan(p->mirror) )
- error(EXIT_FAILURE, 0, "'--ontile' or '--sky' cannot be called with "
- "any of the 'particular' calculation options, for example "
- "'--histogram'. This is because the latter work over the whole "
- "dataset and element positions are changed, but in the former "
- "positions are significant");
+ || p->cumulative || p->sigmaclip || p->fitname
+ || !isnan(p->mirror) )
+ error(EXIT_FAILURE, 0, "'--ontile' or '--sky' cannot be called "
+ "with any of the 'particular' calculation options, for "
+ "example '--histogram'. This is because the latter work "
+ "over the whole dataset and element positions are changed, "
+ "but in the former positions are significant");
/* Make sure the tessellation defining options are given. */
if( tl->tilesize==NULL || tl->numchannels==NULL
@@ -432,14 +437,15 @@ ui_read_check_only_options(struct statisticsparams *p)
/* Make sure a reasonable value is given to '--meanmedqdiff'. */
if(p->meanmedqdiff>0.5)
- error(EXIT_FAILURE, 0, "%g is not acceptable for '--meanmedqdiff'! "
- "This option is the quantile-difference between the quantile "
- "of the mean and 0.5 (the quantile of the median). Since the "
- "range of quantile is from 0.0 to 1.0, the maximum difference "
- "can therefore be 0.5. For more on this option, please see "
- "the \"Quantifying signal in a tile\" section of the Gnuastro "
- "book (with this command: 'info gnuastro \"Quantifying "
- "signal in a tile\"'", p->meanmedqdiff);
+ error(EXIT_FAILURE, 0, "%g is not acceptable for "
+ "'--meanmedqdiff'! This option is the quantile-difference "
+ "between the quantile of the mean and 0.5 (the quantile of "
+ "the median). Since the range of quantile is from 0.0 to "
+ "1.0, the maximum difference can therefore be 0.5. For "
+ "more on this option, please see the \"Quantifying signal "
+ "in a tile\" section of the Gnuastro book (with this "
+ "command: 'info gnuastro \"Quantifying signal in a "
+ "tile\"'", p->meanmedqdiff);
if(p->meanmedqdiff>0.3 && p->cp.quiet==0)
error(EXIT_SUCCESS, 0, "WARNING: %g is very high for "
"'--meanmedqdiff'! This option is the quantile-difference "
@@ -453,19 +459,20 @@ ui_read_check_only_options(struct statisticsparams *p)
/* If a kernel name has been given, we need the HDU. */
if(p->kernelname && gal_fits_file_recognized(p->kernelname)
&& p->khdu==NULL )
- error(EXIT_FAILURE, 0, "no HDU specified for the kernel image. When "
- "A HDU is necessary for FITS files. You can use the '--khdu' "
- "('-u') option and give it the HDU number (starting from "
- "zero), extension name, or anything acceptable by CFITSIO");
+ error(EXIT_FAILURE, 0, "no HDU specified for the kernel image. "
+ "When A HDU is necessary for FITS files. You can use the "
+ "'--khdu' ('-u') option and give it the HDU number "
+ "(starting from zero), extension name, or anything "
+ "acceptable by CFITSIO");
}
/* Sigma-clipping needs 'sclipparams'. */
if(p->sigmaclip && isnan(p->sclipparams[0]))
- error(EXIT_FAILURE, 0, "'--sclipparams' is necessary with '--sigmaclip'. "
- "'--sclipparams' takes two values (separated by a comma) for "
- "defining the sigma-clip: the multiple of sigma, and tolerance "
- "(<1) or number of clips (>1).");
+ error(EXIT_FAILURE, 0, "'--sclipparams' is necessary with "
+ "'--sigmaclip'. '--sclipparams' takes two values (separated "
+ "by a comma) for defining the sigma-clip: the multiple of "
+ "sigma, and tolerance (<1) or number of clips (>1).");
/* If any of the mode measurements are requested, then 'mirrordist' is
@@ -479,8 +486,8 @@ ui_read_check_only_options(struct statisticsparams *p)
case UI_KEY_MODESYMVALUE:
if( isnan(p->mirrordist) )
error(EXIT_FAILURE, 0, "'--mirrordist' is required for the "
- "mode-related single measurements ('--mode', '--modequant', "
- "'--modesym', and '--modesymvalue')");
+ "mode-related single measurements ('--mode', "
+ "'--modequant', '--modesym', and '--modesymvalue')");
break;
case UI_KEY_SIGCLIPSTD:
case UI_KEY_SIGCLIPMEAN:
@@ -489,9 +496,9 @@ ui_read_check_only_options(struct statisticsparams *p)
if( isnan(p->sclipparams[0]) )
error(EXIT_FAILURE, 0, "'--sclipparams' is necessary with "
"sigma-clipping measurements.\n\n"
- "'--sclipparams' takes two values (separated by a comma) for "
- "defining the sigma-clip: the multiple of sigma, and tolerance
"
- "(<1) or number of clips (>1).");
+ "'--sclipparams' takes two values (separated by a comma) "
+ "for defining the sigma-clip: the multiple of sigma, "
+ "and tolerance (<1) or number of clips (>1).");
break;
}
@@ -514,7 +521,10 @@ ui_read_check_only_options(struct statisticsparams *p)
/* When binned outputs are requested, make sure that 'numbins' is set. */
- if( (p->histogram || p->histogram2d || p->cumulative || !isnan(p->mirror))
+ if( (p->histogram
+ || p->histogram2d
+ || p->cumulative
+ || !isnan(p->mirror))
&& p->numbins==0)
error(EXIT_FAILURE, 0, "'--numbins' isn't set. When the histogram or "
"cumulative frequency plots are requested, the number of bins "
@@ -522,12 +532,13 @@ ui_read_check_only_options(struct statisticsparams *p)
if( p->histogram2d )
{
if( p->numbins2==0 )
- error(EXIT_FAILURE, 0, "'--numbins2' isn't set. When a 2D histogram "
- "is requested, the number of bins in the second dimension "
- "('--numbins2') is also necessary");
- if( strcmp(p->histogram2d,"table") && strcmp(p->histogram2d,"image") )
- error(EXIT_FAILURE, 0, "the value to '--histogram2d' can either be "
- "'table' or 'image'");
+ error(EXIT_FAILURE, 0, "'--numbins2' isn't set. When a 2D "
+ "histogram is requested, the number of bins in the "
+ "second dimension ('--numbins2') is also necessary");
+ if( strcmp(p->histogram2d,"table")
+ && strcmp(p->histogram2d,"image") )
+ error(EXIT_FAILURE, 0, "the value to '--histogram2d' can "
+ "either be 'table' or 'image'");
}
/* If an ascii plot is requested, check if the ascii number of bins and
@@ -535,9 +546,52 @@ ui_read_check_only_options(struct statisticsparams *p)
if( (p->asciihist || p->asciicfp)
&& (p->numasciibins==0 || p->asciiheight==0) )
error(EXIT_FAILURE, 0, "when an ascii plot is requested, "
- "'--numasciibins' and '--asciiheight' are mandatory, but atleast "
- "one of these has not been given");
+ "'--numasciibins' and '--asciiheight' are mandatory, but "
+ "at least one of these has not been given");
+
+ /* Find the fitting type code from the input string. */
+ if(p->fitname)
+ {
+ /* Interpret the name. */
+ p->fitid=gal_fit_name_to_id(p->fitname);
+
+ /* Wrong name? */
+ if(p->fitid==GAL_FIT_INVALID)
+ error(EXIT_FAILURE, 0, "'%s' is not a recognized name for "
+ "the type of fitting, please see the description of "
+ "'--fit' in the manual (you can run `info %s`",
+ p->fitname, PROGRAM_EXEC);
+
+ /* Options for polynomial fits. */
+ if( p->fitid==GAL_FIT_POLYNOMIAL
+ || p->fitid==GAL_FIT_POLYNOMIAL_ROBUST
+ || p->fitid==GAL_FIT_POLYNOMIAL_WEIGHTED )
+ {
+
+ /* '--fitmaxpower' is mandatory. */
+ if(p->fitmaxpower==GAL_BLANK_SIZE_T)
+ error(EXIT_FAILURE, 0, "'--fitmaxpower' is necessary for "
+ "polynomial fitting. This is the maximum power of X "
+ "in the fitted polynomial");
+ /* For the robust types, '--fitrobustname' is mandatory. */
+ if( p->fitid==GAL_FIT_POLYNOMIAL_ROBUST )
+ {
+ /* Make sure '--fitrobust' is given at all. */
+ if(p->fitrobustname==NULL)
+ error(EXIT_FAILURE, 0, "'--fitrobust' is mandatory for "
+ "robust fittings");
+
+ /* Find the ID of the fit. */
+ p->fitrobustid=gal_fit_name_robust_to_id(p->fitrobustname);
+ if(p->fitrobustid==GAL_FIT_ROBUST_INVALID)
+ error(EXIT_FAILURE, 0, "'%s' is not a recognized robust "
+ "function in the polynomial fittings, please see "
+ "the manual for the acceptable names",
+ p->fitrobustname);
+ }
+ }
+ }
/* Reverse the list of statistics to print in one row and also the
arguments, so it has the same order the user wanted. */
@@ -609,7 +663,7 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
/* Set the dataset that should be used for the condition. */
- ref = p->reference ? p->reference : p->input;
+ ref = p->input;
/* If the user has given a quantile range, then set the 'greaterequal'
@@ -647,8 +701,10 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
tmp=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL, 0, -1, 1,
NULL, NULL, NULL);
*((float *)(tmp->array)) = p->greaterequal2;
- tmp2=gal_arithmetic(GAL_ARITHMETIC_OP_LT, 1, flags, p->input->next, tmp);
- cond_g=gal_arithmetic(GAL_ARITHMETIC_OP_OR, 1, flagsor, cond_g, tmp2);
+ tmp2=gal_arithmetic(GAL_ARITHMETIC_OP_LT, 1, flags, p->input->next,
+ tmp);
+ cond_g=gal_arithmetic(GAL_ARITHMETIC_OP_OR, 1, flagsor, cond_g,
+ tmp2);
gal_data_free(tmp);
gal_data_free(tmp2);
}
@@ -668,8 +724,10 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
tmp=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL, 0, -1, 1,
NULL, NULL, NULL);
*((float *)(tmp->array)) = p->lessthan2;
- tmp2=gal_arithmetic(GAL_ARITHMETIC_OP_GE, 1, flags, p->input->next, tmp);
- cond_l=gal_arithmetic(GAL_ARITHMETIC_OP_OR, 1, flagsor, cond_l, tmp2);
+ tmp2=gal_arithmetic(GAL_ARITHMETIC_OP_GE, 1, flags, p->input->next,
+ tmp);
+ cond_l=gal_arithmetic(GAL_ARITHMETIC_OP_OR, 1, flagsor, cond_l,
+ tmp2);
gal_data_free(tmp);
gal_data_free(tmp2);
}
@@ -684,7 +742,8 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
cond = isnan(p->greaterequal) ? cond_l : cond_g;
break;
case 2:
- cond = gal_arithmetic(GAL_ARITHMETIC_OP_OR, 1, flagsor, cond_l, cond_g);
+ cond = gal_arithmetic(GAL_ARITHMETIC_OP_OR, 1, flagsor, cond_l,
+ cond_g);
gal_data_free(cond_g);
break;
}
@@ -701,7 +760,8 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
blank value will be used in the proper type of the input in the
'where' operator. Also reset the blank flags so they are checked again
if necessary.*/
- gal_arithmetic(GAL_ARITHMETIC_OP_WHERE, 1, flagsor, p->input, cond, blank);
+ gal_arithmetic(GAL_ARITHMETIC_OP_WHERE, 1, flagsor, p->input,
+ cond, blank);
p->input->flag &= ~GAL_DATA_FLAG_BLANK_CH;
p->input->flag &= ~GAL_DATA_FLAG_HASBLANK;
if(p->input->next)
@@ -773,8 +833,8 @@ ui_make_sorted_if_necessary(struct statisticsparams *p)
/* Merge all given columns into one list. This is because people may call
for different columns in one call to '--column' ('-c', separated by a
comma), or multiple calls to '-c'. */
-gal_list_str_t *
-ui_read_columns_in_one(gal_list_str_t *incolumns)
+void
+ui_read_columns_in_one(struct statisticsparams *p)
{
size_t i;
gal_data_t *strs;
@@ -783,7 +843,7 @@ ui_read_columns_in_one(gal_list_str_t *incolumns)
/* Go over the separate calls to the '-c' option, see the explaination in
Table's 'ui_columns_prepare' function for more on every step. */
- for(tmp=incolumns;tmp!=NULL;tmp=tmp->next)
+ for(tmp=p->columns;tmp!=NULL;tmp=tmp->next)
{
/* Remove any new-line character. */
for(c=tmp->v;*c!='\0';++c)
@@ -810,8 +870,9 @@ ui_read_columns_in_one(gal_list_str_t *incolumns)
printf("%s\n", tmp->v);
*/
- /* Return the final unified list. */
- return final;
+ /* Free the input list and replace it with 'final'. */
+ gal_list_str_free(p->columns, 1);
+ p->columns=final;
}
@@ -821,35 +882,29 @@ ui_read_columns_in_one(gal_list_str_t *incolumns)
void
ui_read_columns(struct statisticsparams *p)
{
- int toomanycols=0, tformat;
+ int tformat;
gal_data_t *cols, *tmp, *cinfo;
- size_t size, ncols, nrows, counter=0;
- gal_list_str_t *incols, *columnlist=NULL;
+ size_t ncols, nrows, counter=0;
gal_list_str_t *lines=gal_options_check_stdin(p->inputname,
p->cp.stdintimeout,
"input");
/* Merge possibly multiple calls to '-c' (each possibly separated by a
coma) into one list. */
- incols=ui_read_columns_in_one(p->columns);
+ ui_read_columns_in_one(p);
- /* If any columns are specified, make sure there is a maximum of two
- columns. */
- if(gal_list_str_number(incols)>2)
+ /* If any columns are specified, and fitting hasn't been requested, make
+ sure there is a maximum of two columns. */
+ if(p->fitname==NULL && gal_list_str_number(p->columns)>2)
error(EXIT_FAILURE, 0, "%zu input columns were given but currently a "
"maximum of two columns are supported (two columns only for "
"special operations, the majority of operations are on a single "
- "column)", gal_list_str_number(incols));
-
- /* If a reference column is also given, add it to the list of columns to
- read. */
- if(p->refcol)
- gal_list_str_add(&columnlist, p->refcol, 0);
+ "column)", gal_list_str_number(p->columns));
/* If no column is specified, Statistics will abort and an error will be
printed when the table has more than one column. If there is only one
column, there is no need to specify any, so Statistics will use it. */
- if(incols==NULL)
+ if(p->columns==NULL)
{
/* Get the basic table information. */
cinfo=gal_table_info(p->inputname, p->cp.hdu, lines, &ncols, &nrows,
@@ -865,7 +920,7 @@ ui_read_columns(struct statisticsparams *p)
? gal_checkset_dataset_name(p->inputname, p->cp.hdu)
: "Standard input" ));
case 1:
- gal_list_str_add(&incols, "1", 1);
+ gal_list_str_add(&p->columns, "1", 1);
break;
default:
error(EXIT_FAILURE, 0, "%s is a table containing more than one "
@@ -873,11 +928,11 @@ ui_read_columns(struct statisticsparams *p)
"specified.\n\n"
"Please use the '--column' ('-c') option to specify a "
"column. You can either give it the column number "
- "(couting from 1), or a match/search in its meta-data (e.g., "
- "column names).\n\n"
+ "(couting from 1), or a match/search in its meta-data "
+ "(e.g., column names).\n\n"
"For more information, please run the following command "
- "(press the 'SPACE' key to go down and 'q' to return to the "
- "command-line):\n\n"
+ "(press the 'SPACE' key to go down and 'q' to return to "
+ "the command-line):\n\n"
" $ info gnuastro \"Selecting table columns\"\n",
( p->inputname
? gal_checkset_dataset_name(p->inputname, p->cp.hdu)
@@ -885,33 +940,31 @@ ui_read_columns(struct statisticsparams *p)
}
}
- if(columnlist) columnlist->next=incols;
- else columnlist=incols;
/* Read the desired column(s). */
- cols=gal_table_read(p->inputname, p->cp.hdu, lines, columnlist,
+ cols=gal_table_read(p->inputname, p->cp.hdu, lines, p->columns,
p->cp.searchin, p->cp.ignorecase, p->cp.numthreads,
p->cp.minmapsize, p->cp.quietmmap, NULL);
gal_list_str_free(lines, 1);
+ p->input=cols;
/* If the input was from standard input, we'll set the input name to be
'stdin' (for future reporting). */
if(p->inputname==NULL)
gal_checkset_allocate_copy("stdin", &p->inputname);
- /* Put the columns into the proper gal_data_t. */
- size=cols->size;
- while(cols!=NULL)
+ /* Print an error if there are too many columns: */
+ if(gal_list_str_number(p->columns)!=gal_list_data_number(cols))
+ gal_tableintern_error_col_selection(p->inputname, p->cp.hdu, "too "
+ "many columns were selected by the given values to the "
+ "'--column' (or '-c'). In other words, the number of "
+ "columns that were read, are more than the number of "
+ "columns given to '--column' (this can happen due to "
+ "columns with same name for example)");
+
+ /* Make sure all columns have a usable datatype. */
+ for(tmp=cols; tmp!=NULL; tmp=tmp->next)
{
- /* Pop out the top column. */
- tmp=gal_list_data_pop(&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_TYPE_BIT:
@@ -923,44 +976,141 @@ ui_read_columns(struct statisticsparams *p)
"which is not currently supported by %s", counter,
gal_type_name(tmp->type, 1), PROGRAM_NAME);
}
+ }
+
+ /* For a check:
+ printf("Number of input columns: %zu\n", gal_list_data_number(p->input));
+ exit(0);
+ */
+}
+
+
+
- /* Put the column into the proper pointer. */
- switch(++counter)
+
+void
+ui_preparations_fitestimate(struct statisticsparams *p)
+{
+ size_t one=1;
+ double dbl, *dptr=&dbl;
+ gal_list_str_t *fecols=NULL;
+
+ if( gal_type_from_string((void **)(&dptr), p->fitestimate,
+ GAL_TYPE_FLOAT64) )
+ {
+ /* If the value was "self", then we should put the input file name,
+ its HDU and the first column in required values so they are fully
+ read. We aren't using the read 'p->input' because rows with a
+ blank value have been removed there. */
+ if( !strcmp(p->fitestimate, "self") )
{
- case 1: p->input=tmp; break;
- case 2:
- if(gal_list_str_number(incols)==2)
- p->input->next=tmp;
- else
- if(p->refcol) p->reference=tmp; else toomanycols=1;
- break;
- case 3:
- if(p->refcol) p->reference=tmp; else toomanycols=1; break;
- default: toomanycols=1;
+ free(p->fitestimate);
+ free(p->fitestimatehdu);
+ free(p->fitestimatecol);
+ gal_checkset_allocate_copy(p->inputname, &p->fitestimate);
+ gal_checkset_allocate_copy(p->cp.hdu, &p->fitestimatehdu);
+ gal_checkset_allocate_copy(p->columns->v, &p->fitestimatecol);
}
- /* Print an error if there are too many columns: */
- if(toomanycols)
- gal_tableintern_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.");
+ /* Make sure a HDU is specified. We need to do this here (not
+ in 'ui_read_check_only_options') because only here we know
+ that the user specified a file not a value. */
+ if( gal_fits_name_is_fits(p->fitestimate) && p->fitestimatehdu==NULL )
+ error(EXIT_FAILURE, 0, "no HDU specified for '%s' (given to "
+ "'--fitestimate'). Please use the '--fitestimatehdu' "
+ "option to specify the HDU", p->fitestimate);
+
+ /* Make sure a column is specified. We need to do this here (not
+ in 'ui_read_check_only_options') because only here we know
+ that the user specified a file not a value. */
+ if( p->fitestimatecol==NULL )
+ error(EXIT_FAILURE, 0, "no column specified for '%s' (given to "
+ "'--fitestimate'). Please use the '--fitestimatecol' "
+ "option to specify the column", p->fitestimate);
+
+ /* Given string couldn't be read as a number, so try reading it
+ as a table. */
+ gal_list_str_add(&fecols, p->fitestimatecol, 1);
+ p->fitestval=gal_table_read(p->fitestimate, p->fitestimatehdu,
+ NULL, fecols,
+ p->cp.searchin, p->cp.ignorecase, 1,
+ p->cp.minmapsize, p->cp.quietmmap,
+ NULL);
+
+ /* If more than one column matched, inform the user. */
+ if(p->fitestval->next)
+ gal_tableintern_error_col_selection(p->fitestimate,
+ p->fitestimatehdu, "More than one column matched "
+ "the value given to '--fitestimatecol'.");
}
-
- /* Clean up. */
- gal_list_str_free(incols, 1);
- if(columnlist!=incols)
+ else
{
- columnlist->next=NULL;
- gal_list_str_free(columnlist, 0);
+ /* Given string could be read as a double (in variable 'd'). */
+ p->fitestval=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &one,
+ NULL, 0, -1, 1, NULL, NULL, NULL);
+ ((double *)(p->fitestval->array))[0]=dbl;
}
- /* For a check:
- printf("Number of input columns: %zu\n", gal_list_data_number(p->input));
- exit(0);
- */
+ /* Make sure 'fitestval' has a double type. */
+ p->fitestval=gal_data_copy_to_new_type_free(p->fitestval,
+ GAL_TYPE_FLOAT64);
+}
+
+
+
+
+
+void
+ui_preparations_fit_weight(struct statisticsparams *p)
+{
+ double *d, *df;
+ gal_data_t *wht;
+
+ /* This is only necessary for fitting models that require a weight. */
+ switch(p->fitid)
+ {
+ case GAL_FIT_LINEAR_WEIGHTED:
+ case GAL_FIT_POLYNOMIAL_WEIGHTED:
+ case GAL_FIT_LINEAR_NO_CONSTANT_WEIGHTED:
+
+ /* Basic sanity check first. */
+ if( gal_list_data_number(p->input)<3)
+ error(EXIT_FAILURE, 0, "no weight column specified! A "
+ "weight-based fit needs a third input column");
+ if(p->fitweight==0)
+ error(EXIT_FAILURE, 0, "the nature of the input weights for "
+ "fitting haven't been specified. Please use the "
+ "'--fitweight' option to specify this. It can take "
+ "values of 'std' (if the weight column is the standard "
+ "deviation), 'var' (for variance) and 'inv-var' "
+ "(inverse-variance) which is direct");
+
+ /* Convert the third dataset to double-precision. */
+ wht=gal_data_copy_to_new_type_free(p->input->next->next,
+ GAL_TYPE_FLOAT64);
+ p->input->next->next=wht;
+
+ /* Based on the input nature, convert it to the inverse of the
+ variance. */
+ d=wht->array;
+ if( !strcmp(p->fitweight, "std") ) /* Standard deviation. */
+ {
+ p->fitwhtid=STATISTICS_FIT_WHT_STD;
+ df=d+wht->size; do *d=1/(*d * *d); while(++d<df);
+ }
+ else if ( !strcmp(p->fitweight, "var") ) /* Variance. */
+ {
+ p->fitwhtid=STATISTICS_FIT_WHT_VAR;
+ df=d+wht->size; do *d = 1 / *d; while(++d<df);
+ }
+ else if ( !strcmp(p->fitweight, "inv-var") ) /* Inverse variance. */
+ p->fitwhtid=STATISTICS_FIT_WHT_INVVAR;
+ else /* Not recognized: error! */
+ error(EXIT_FAILURE, 0, "'%s' is not a recognized weight-type! "
+ "Please use either 'std' (standard deviation), 'var' "
+ "(variance) or 'inv-var' (inverse variance)", p->fitweight);
+ break;
+ }
}
@@ -970,11 +1120,10 @@ ui_read_columns(struct statisticsparams *p)
void
ui_preparations(struct statisticsparams *p)
{
+ gal_data_t *check;
int keepinputdir=p->cp.keepinputdir;
- gal_data_t *check, *flag1, *flag2, *flag;
struct gal_options_common_params *cp=&p->cp;
struct gal_tile_two_layer_params *tl=&cp->tl;
- unsigned char flagsor = GAL_ARITHMETIC_FLAGS_BASIC;
char *checkbasename = p->cp.output ? p->cp.output : p->inputname;
/* Change 'keepinputdir' based on if an output name was given. */
@@ -1000,9 +1149,9 @@ ui_preparations(struct statisticsparams *p)
p->inputformat=INPUT_FORMAT_TABLE;
/* Two columns can only be given with 2D histogram mode. */
- if(p->histogram2d==0 && p->input->next!=NULL)
- error(EXIT_FAILURE, 0, "two column input is currently only "
- "supported for 2D histogram mode");
+ if(p->histogram2d==0 && p->fitname==NULL && p->input->next!=NULL)
+ error(EXIT_FAILURE, 0, "multi-column input is currently only "
+ "supported for 2D histogram or fitting modes");
}
/* Read the convolution kernel if necessary. */
@@ -1025,11 +1174,13 @@ ui_preparations(struct statisticsparams *p)
/* Make the tile check image if requested. */
if(tl->checktiles)
{
- tl->tilecheckname=gal_checkset_automatic_output(cp, checkbasename,
+ tl->tilecheckname=gal_checkset_automatic_output(cp,
+ checkbasename,
"_tiled.fits");
check=gal_tile_block_check_tiles(tl->tiles);
if(p->inputformat==INPUT_FORMAT_IMAGE)
- gal_fits_img_write(check, tl->tilecheckname, NULL, PROGRAM_NAME);
+ gal_fits_img_write(check, tl->tilecheckname, NULL,
+ PROGRAM_NAME);
else
{
gal_checkset_writable_remove(tl->tilecheckname, 0,
@@ -1055,26 +1206,8 @@ ui_preparations(struct statisticsparams *p)
/* Only keep the elements we want. Note that if we have more than one
column, we need to move the same rows in both (otherwise their
widths won't be equal). */
- if(p->input->next)
- {
- flag1=gal_blank_flag(p->input);
- flag2=gal_blank_flag(p->input->next);
- flag=gal_arithmetic(GAL_ARITHMETIC_OP_OR, 1, flagsor, flag1, flag2);
-
- gal_blank_flag_apply(p->input, flag);
- gal_blank_flag_apply(p->input->next, flag);
-
- gal_blank_remove(p->input);
- gal_blank_remove(p->input->next);
-
- gal_data_free(flag);
- p->input->next->flag &= ~GAL_DATA_FLAG_HASBLANK ;
- p->input->next->flag |= GAL_DATA_FLAG_BLANK_CH ;
- }
- else
- gal_blank_remove(p->input);
- p->input->flag &= ~GAL_DATA_FLAG_HASBLANK ;
- p->input->flag |= GAL_DATA_FLAG_BLANK_CH ;
+ if(p->input->next) gal_blank_remove_rows(p->input, NULL);
+ else gal_blank_remove(p->input);
/* Make sure there actually are any (non-blank) elements left. */
if(p->input->size==0)
@@ -1083,8 +1216,8 @@ ui_preparations(struct statisticsparams *p)
/* Make sure there is data remaining: */
if(p->input->size==0)
- error(EXIT_FAILURE, 0, "%s: no data, maybe the '--greaterequal' or "
- "'--lessthan' options need to be adjusted",
+ error(EXIT_FAILURE, 0, "%s: no data, maybe the '--greaterequal' "
+ "or '--lessthan' options need to be adjusted",
gal_fits_name_save_as_string(p->inputname, cp->hdu) );
/* Make the sorted array if necessary. */
@@ -1097,6 +1230,15 @@ ui_preparations(struct statisticsparams *p)
/* Reset 'keepinputdir' to what it originally was. */
p->cp.keepinputdir=keepinputdir;
+
+ /* Make sure the output doesn't already exist. */
+ gal_checkset_writable_remove(p->cp.output, p->cp.keep,
+ p->cp.dontdelete);
+
+ /* Set the fit-estimate column, and prepare the weight based on the
+ user's specification.*/
+ ui_preparations_fit_weight(p);
+ if(p->fitestimate) ui_preparations_fitestimate(p);
}
@@ -1122,7 +1264,8 @@ ui_preparations(struct statisticsparams *p)
/**************************************************************/
void
-ui_read_check_inputs_setup(int argc, char *argv[], struct statisticsparams *p)
+ui_read_check_inputs_setup(int argc, char *argv[],
+ struct statisticsparams *p)
{
struct gal_options_common_params *cp=&p->cp;
@@ -1209,7 +1352,6 @@ ui_free_report(struct statisticsparams *p)
/* Free the allocated arrays: */
free(p->cp.hdu);
free(p->cp.output);
- gal_data_free(p->reference);
gal_list_data_free(p->input);
gal_list_f64_free(p->tp_args);
gal_list_i32_free(p->singlevalue);
diff --git a/bin/statistics/ui.h b/bin/statistics/ui.h
index 1dcaeaa4..3e699e0b 100644
--- a/bin/statistics/ui.h
+++ b/bin/statistics/ui.h
@@ -37,6 +37,7 @@ enum program_args_groups
UI_GROUP_PARTICULAR_STAT,
UI_GROUP_SKY,
UI_GROUP_HIST_CFP,
+ UI_GROUP_FIT,
};
@@ -45,14 +46,13 @@ enum program_args_groups
/* Available letters for short options:
- a b e f j p v w x z
+ a b e j r v w x z
B G J L W X Y
*/
enum option_keys_enum
{
/* With short-option version. */
UI_KEY_COLUMN = 'c',
- UI_KEY_REFCOL = 'r',
UI_KEY_GREATEREQUAL = 'g',
UI_KEY_LESSTHAN = 'l',
UI_KEY_QRANGE = 'Q',
@@ -71,6 +71,8 @@ enum option_keys_enum
UI_KEY_SKY = 'y',
UI_KEY_KERNEL = 'k',
UI_KEY_CONTOUR = 'R',
+ UI_KEY_FIT = 'f',
+ UI_KEY_FITMAXPOWER = 'p',
/* Only with long version (start with a value 1000, the rest will be set
automatically). */
@@ -111,6 +113,11 @@ enum option_keys_enum
UI_KEY_GREATEREQUAL2,
UI_KEY_LESSTHAN2,
UI_KEY_ONEBINSTART2,
+ UI_KEY_FITWEIGHT,
+ UI_KEY_FITESTIMATE,
+ UI_KEY_FITESTIMATEHDU,
+ UI_KEY_FITESTIMATECOL,
+ UI_KEY_FITROBUST,
};
diff --git a/bin/table/arithmetic.c b/bin/table/arithmetic.c
index 901c774d..0bbec8a7 100644
--- a/bin/table/arithmetic.c
+++ b/bin/table/arithmetic.c
@@ -289,17 +289,15 @@ arithmetic_init(struct tableparams *p, struct
arithmetic_token **arith,
size_t one=1;
uint8_t ntype;
gal_data_t *col;
- char *delimiter=" \t";
struct arithmetic_token *atmp, *node=NULL;
- char *token=NULL, *lasttoken=NULL, *saveptr;
struct gal_options_common_params *cp=&p->cp;
+ char *token=NULL, *saveptr, *delimiter=" \t";
/* Parse all the given tokens. */
token=strtok_r(expression, delimiter, &saveptr);
while(token!=NULL)
{
/* Allocate and initialize this arithmetic token. */
- lasttoken=token;
node=arithmetic_add_new_to_end(arith);
/* See if the token is an operator, if not check other cases. */
@@ -351,13 +349,6 @@ arithmetic_init(struct tableparams *p, struct
arithmetic_token **arith,
/* Go to the next token. */
token=strtok_r(NULL, delimiter, &saveptr);
}
-
- /* A small sanity check: the last added token must be an operator (unless
- its a 'loadcol'). */
- if(node==NULL
- || (node->loadcol==NULL && node->operator==GAL_ARITHMETIC_OP_INVALID))
- error(EXIT_FAILURE, 0, "last token in arithmetic column ('%s') "
- "is not a recognized operator", lasttoken);
}
@@ -418,7 +409,6 @@ arithmetic_indexs_final(struct tableparams *p)
tmp->start=startind;
tmp->numsimple=numcols;
}
-
}
}
@@ -957,7 +947,9 @@ gal_data_t *
arithmetic_read_at_usage(struct tableparams *p,
struct arithmetic_token *token)
{
+ char *c;
size_t i;
+ int hasnewline=0;
/* If the number is blank, then we have a name and should return the
first column in the table that has the given name. */
@@ -965,14 +957,26 @@ arithmetic_read_at_usage(struct tableparams *p,
{
/* Search all the existing columns in the table. */
for(i=0; i<p->numcolarray; ++i)
- if( strcasecmp(p->colarray[i]->name, token->id_at_usage)==0 )
+ if( p->colarray[i]->name
+ && strcasecmp(p->colarray[i]->name, token->id_at_usage)==0 )
return p->colarray[i];
/* If control reaches here, then the requested name didn't exist in
- the table. */
+ the table. This may happen because the user broke an arithmetic
+ command into muliple lines and forgot to put a back-slash. If so,
+ the string will have a new-line within it.*/
+ for(c=token->id_at_usage;*c!='\0';++c)
+ if(*c=='\n') hasnewline=1;
+
+ /* print the error message. */
error(EXIT_FAILURE, 0, "'%s' doesn't correspond to the name of "
- "any column in the table until this column arithmetic",
- token->id_at_usage);
+ "any column in the table until this column arithmetic%s",
+ token->id_at_usage, hasnewline?". Because there is a "
+ "new-line in the name, it is highly probable that it "
+ "wasn't actually a name. This can happen when you break "
+ "an arithmetic command into multiple lines, but have "
+ "forgot to put a back-slash ('\\') at the end of a "
+ "broken line":"");
}
/* We already have the desired column number. */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 7e1598e7..7ca0f6a8 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -228,7 +228,7 @@ To go to the sections or subsections, you have to click on
the menu entries that
* Data containers:: Tools to operate on extensions and tables.
* Data manipulation:: Tools for basic image manipulation.
* Data analysis:: Analyze images.
-* Modeling and fittings:: Make and fit models.
+* Data modeling:: Modeling observed data.
* High-level calculations:: Physical calculations.
* Installed scripts:: Installed scripts that operate like programs.
* Makefile extensions:: Use Gnuastro's features as GNU Make functions.
@@ -507,8 +507,8 @@ Arithmetic operators
* Basic mathematical operators:: For example +, -, /, log, pow, and etc.
* Trigonometric and hyperbolic operators:: sin, cos, atan, asinh, and etc
-* Constants::
-* Unit conversion operators:: magnitudes to counts, or parsecs to AUs, and
etc.
+* Constants:: Physical and Mathematical constants.
+* Unit conversion operators:: Various unit conversions necessary.
* Statistical operators:: Statistics of a single dataset (for example
mean).
* Stacking operators:: Coadding or combining multiple datasets into
one.
* Filtering operators:: Smoothing a dataset through mixing pixel with
neighbors.
@@ -570,6 +570,7 @@ Statistics
* Histogram and Cumulative Frequency Plot:: Basic definitions.
* 2D Histograms:: Plotting the distribution of two variables.
* Sigma clipping:: Definition of @mymath{\sigma}-clipping.
+* Least squares fitting::
* Sky value:: Definition and derivation of the Sky value.
* Invoking aststatistics:: Arguments and options to Statistics.
@@ -584,6 +585,15 @@ Sky value
* Sky value misconceptions:: Wrong methods to estimate the Sky value.
* Quantifying signal in a tile:: Method to estimate the presence of signal.
+Invoking Statistics
+
+* Input to Statistics:: How to specify the inputs to Statistics.
+* Single value measurements:: Can be used together (like --mean, or
--maximum).
+* Generating histograms and cumulative frequency plots:: Histogram and CFP
tables.
+* Fitting options:: Least squares fitting.
+* Contour options:: Table of contours.
+* Statistics on tiles:: Possible to do single-valued measurements on
tiles.
+
NoiseChisel
* NoiseChisel changes after publication:: Updates since published papers.
@@ -636,7 +646,7 @@ Match
* Matching algorithms:: Different ways to find the match
* Invoking astmatch:: Inputs, outputs and options of Match
-Modeling and fitting
+Data modeling
* MakeProfiles:: Making mock galaxies and stars.
* MakeNoise:: Make (add) noise to an image.
@@ -771,6 +781,7 @@ Gnuastro library
* Permutations:: Re-order (or permute) the values in a dataset.
* Matching:: Matching catalogs based on position.
* Statistical operations:: Functions for basic statistics.
+* Fitting functions:: Fit independent and measured variables.
* Binary datasets:: Datasets that can only have values of 0 or 1.
* Labeled datasets:: Working with Segmented/labeled datasets.
* Convolution functions:: Library functions to do convolution.
@@ -16337,7 +16348,7 @@ These operators take a single operand.
@node Constants, Unit conversion operators, Trigonometric and hyperbolic
operators, Arithmetic operators
@subsubsection Constants
-@cindex @mymath{\pi}
+@cindex Pi
During your analysis it is often necessary to have certain constants like the
number @mymath{\pi}.
The ``operators'' in this section do not actually take any operand, they just
replace the desired constant into the stack.
So in effect, these are actually operands.
@@ -16345,14 +16356,14 @@ But since their value is not inserted by the user, we
have placed them in the li
@table @code
@item e
-@cindex @mymath{e} (base of natural logarithm)
+@cindex e (base of natural logarithm)
@cindex Euler's number (@mymath{e})
@cindex Base of natural logarithm (@mymath{e})
Euler’s number, or the base of the natural logarithm (no units).
See @url{https://en.wikipedia.org/wiki/E_(mathematical_constant), Wikipedia}.
@item pi
-@cindex @mymath{\pi}
+@cindex Pi
Ratio of circle’s circumference to its diameter (no units).
See @url{https://en.wikipedia.org/wiki/Pi, Wikipedia}.
@@ -17520,6 +17531,9 @@ Without this increment, when the column values are the
same (happens a lot, for
But this feature has a side-effect: that if the order of calling the
@code{mknoise-*} operators changes, the seeds used for each operator will
change@footnote{We have defined @url{https://savannah.gnu.org/task/?15971, Task
15971} in Gnuastro's project management system to address this.
If you need this feature please send us an email at
@code{bug-gnuastro@@gnu.org} (to motivate us in its implementation).}.
+In case each data element should have an independent sigma, the first popped
operand can be a dataset of the same size as the second.
+In this case, for each element, a different noise measure (for example sigma
in @code{mknoise-sigma}) will be used.
+
@item mknoise-poisson
@cindex Poisson noise
Add Poisson noise to each element of the input dataset (see @ref{Photon
counting noise}).
@@ -19601,7 +19615,7 @@ Alternatively, a value of @code{0} will keep output
pixels that are even infinit
-@node Data analysis, Modeling and fittings, Data manipulation, Top
+@node Data analysis, Data modeling, Data manipulation, Top
@chapter Data analysis
Astronomical datasets (images or tables) contain very valuable information,
the tools in this section can help in analyzing, extracting, and quantifying
that information.
@@ -19632,6 +19646,7 @@ The Statistics program is designed for such situations.
* Histogram and Cumulative Frequency Plot:: Basic definitions.
* 2D Histograms:: Plotting the distribution of two variables.
* Sigma clipping:: Definition of @mymath{\sigma}-clipping.
+* Least squares fitting::
* Sky value:: Definition and derivation of the Sky value.
* Invoking aststatistics:: Arguments and options to Statistics.
@end menu
@@ -19931,7 +19946,7 @@ $ pdflatex cmd-report.tex
The improved quality, blending in with the text, vector-graphics resolution
and other features make this plot pleasing to the eye, and let your readers
focus on the main point of your scientific argument.
PGFPlots can also built the PDF of the plot separately from the rest of the
paper/report, see @ref{2D histogram as a table for plotting} for the necessary
changes in the preamble.
-@node Sigma clipping, Sky value, 2D Histograms, Statistics
+@node Sigma clipping, Least squares fitting, 2D Histograms, Statistics
@subsection Sigma clipping
Let's assume that you have pure noise (centered on zero) with a clear
@url{https://en.wikipedia.org/wiki/Normal_distribution,Gaussian distribution},
or see @ref{Photon counting noise}.
@@ -20002,9 +20017,243 @@ The ASCII histogram that is printed on the
command-line with @option{--asciihist
@end cartouche
+@node Least squares fitting, Sky value, Sigma clipping, Statistics
+@subsection Least squares fitting
+
+@cindex Radial profile
+@cindex Least squares fitting
+@cindex Fitting (least squares)
+@cindex Star formation main sequence
+After completing a good observation, doing robust data reduction and
finalizing the measurements, it is commonly necessary to parametrize the
derived correlations.
+For example, you have derived the radial profile of the PSF of your image (see
@ref{Building the extended PSF}).
+You now want to parametrize the radial profile to estimate the slope.
+Alternatively, you may have found the star formation rate and stellar mass of
your sample of galaxies.
+Now, you want to derive the star formation main sequence as a parametric
relation between the two.
+The fitting functions below can be used for such purposes.
+
+@cindex GSL
+@cindex GNU Scientific Library
+Gnuastro's least squares fitting features are just wrappers over the least
squares fitting methods of the
@url{https://www.gnu.org/software/gsl/doc/html/lls.html, linear} and
@url{https://www.gnu.org/software/gsl/doc/html/nls.html, nonlinear}
least-squares fitting functions of the GNU Scientific Library (GSL).
+For the low-level details and equations of the methods, please see the GSL
documentation.
+The names have been preserved here in Gnuastro to simplify the connection with
GSL and follow the details in the detailed documentation there.
+
+GSL is a very low-level library, designed for maximal portability to many
scenarios, and power.
+Therefore calling GSL's functions directly for a fast operation requires a
good knowledge of the C programming language and many lines of code.
+As a low-level library, GSL is designed to be the back-end of higher-level
programs (like Gnuastro).
+Through the Statistics program, in Gnuastro we provide a high-level interface
to access to GSL's very powerful least squares fitting engine to read/write
from/to standard data formats in astronomy.
+A fully working example is shown below.
+
+@cindex Gaussian noise
+@cindex Noise (Gaussian)
+To activate fitting in Statistics, simply give your desired fitting method to
the @option{--fit} option (for the full list of acceptable methods, see
@ref{Fitting options}).
+For example, with the command below, we'll build a fake measurement table
(including noise) from the polynomial @mymath{y=1.23-4.56x+7.89x^2}.
+To understand how this equation translates to the command below (part before
@code{set-y}), see @ref{Reverse polish notation} and @ref{Column arithmetic}.
+We will set the X axis to have values from 0.1 to 2, with steps of 0.01 and
let's assume a random Gaussian noise to each @mymath{y} measurement:
@mymath{\sigma_y=0.1y}.
+To make the random number generation exactly reproducible, we are also setting
the seed (see @ref{Generating random numbers}, which also uses GSL as a
backend).
+To learn more about the @code{mknoise-sigma} operator, see the Arithmetic
program's @ref{Random number generators}.
+
+@example
+$ export GSL_RNG_SEED=1664015492
+$ seq 0.1 0.01 2 \
+ | asttable --output=noisy.fits --envseed -c1 \
+ -c'arith 1.23 -4.56 $1 x + 7.89 $1 x $1 x + set-y \
+ 0.1 y x set-yerr \
+ y yerr mknoise-sigma yerr' \
+ --colmetadata=1,X --colmetadata=2,Y \
+ --colmetadata=3,Yerr
+@end example
+
+@noindent
+Let's have a look at the output plot with TOPCAT using the command below.
+
+@example
+$ astscript-fits-view noisy.fits
+@end example
+
+@noindent
+To see the error-bars, after opening the scatter plot, go into the ``Form''
tab for that plot.
+Click on the button with a green ``+'' sign followed by ``Forms'' and select
``XYError''.
+On the side-menu, in front of ``Y Positive Error'', select the @code{Yerr}
column of the input table.
+
+As you see, the error bars do indeed increase for higher X axis values.
+Since we have error bars in this example (as in any measurement), we can use
weighted fitting.
+Also, this isn't a linear relation, so we'll use a polynomial to second order
(a maximum power of 2 in the form of @mymath{Y=c_0+c_1X+c_2X^2}):
+
+@example
+$ aststatistics noisy.fits -cX,Y,Yerr --fit=polynomial-weighted \
+ --fitmaxpower=2
+Statistics (GNU Astronomy Utilities) @value{VERSION}
+-------
+Fitting results (remove extra info with '--quiet' or '-q)
+ Input file: noisy.fits (hdu: 1) with 191 non-blank rows.
+ X column: X
+ Y column: Y
+ Weight column: Yerr [Standard deviation of Y in each row]
+
+Fit function: Y = c0 + (c1 * X^1) + (c2 * X^2) + ... (cN * X^N)
+ N: 2
+ c0: +1.2286211608
+ c1: -4.5127796636
+ c2: +7.8435883943
+
+Covariance matrix:
+ +0.0010496001 -0.0039928488 +0.0028367390
+ -0.0039928488 +0.0175244127 -0.0138030778
+ +0.0028367390 -0.0138030778 +0.0128129806
+
+Reduced chi^2 of fit:
+ +0.9740670090
+@end example
+
+As you see from the elaborate message, the weighted polynomial fitting has
found return the @mymath{c_0}, @mymath{c_1} and @mymath{c_2} of
@mymath{Y=c_0+c_1X+c_2X^2} that best represents the data we inserted.
+Our input values were @mymath{c_0=1.23}, @mymath{c_1=-4.56} and
@mymath{c_2=7.89}, and the fitted values are @mymath{c_0\approx1.2286},
@mymath{c_1\approx-4.5128} and @mymath{c_2\approx7.8436} (which is
statistically a very good fit! given that we knew the original values
a-priori!).
+The covariance matrix is also calculated, it is necessary to calculate error
bars on the estimations and contains a lot of information (e.g., possible
correlations between parameters).
+Finally, the reduced @mymath{\chi^2} (or @mymath{\chi_{red}^2}) of the fit is
also printed (which was the measure to minimize).
+A @mymath{\chi_{red}^2\approx1} shows a good fit.
+This is good for real-world scenarios when you don't know the original values
a-priori.
+For more on interpretting @mymath{\chi_{red}^2\approx1}, see
@url{https://arxiv.org/abs/1012.3754, Andrae et al (2010)}.
+
+The comparison of fitted and input values look pretty good, but nothing beats
visual inspection!
+To see how this looks compared to the data, let's open the table again:
+
+@example
+$ astscript-fits-view noisy.fits
+@end example
+
+Repeat the steps below to show the scatter plot and error-bars.
+Then, go to the ``Layers'' menu and select ``Add Function Control''.
+Use the results above to fill the box infront of ``Function Expression'':
@code{1.2286+(-4.5128*x)+(7.8436*x*x)}.
+You will see that the second order polynomial falls very nicely over the
points@footnote{After plotting, you will notice that the legend made the plot
too thin.
+Fortunately you have a lot of empty space within the plot.
+To bring the legend in, click on the ``Legend'' item on the bottom-left menu,
in the ``Location'' tab, click on ``Internal'' and hold and move it to the
top-left in the box below.
+To make the functional fit more clear, you can click on the ``Function'' item
of the bottom-left menu.
+In the ``Style'' tab, change the color and thickness.}.
+But this fit is not perfect: it also has errors (inherited from the
measurement errors).
+We need the covariance matrix to estimate the errors on each point, and that
can be complex to do by hand.
+
+Fortunately GSL has the tools to easily estimate the function at any point and
also calculate its corresponding error.
+To access this feature within Gnuastro's Statistics program, you should use
the @option{--fitestimate} option.
+You can either give an independent table file name (with
@option{--fitestimatehdu} and @option{--fitestimatecol} to specify the HDU and
column in that file), or just @code{self} so it uses the same X axis column
that was used in this fit.
+Let's use the easier case:
+
+@example
+$ aststatistics noisy.fits -cX,Y,Yerr --fit=polynomial-weighted \
+ --fitmaxpower=2 --fitestimate=self --output=est.fits
+@end example
+You will see a new line printed in the output, saying that the estimation was
writte in @file{est.fits}.
+You can now inspect the two tables with TOPCAT again with the command below.
+After TOPCAT opens, plot both scatter plots:
-@node Sky value, Invoking aststatistics, Sigma clipping, Statistics
+@example
+$ astscript-fits-view noisy.fits est.fits
+@end example
+
+It is clear that they fall nicely ontop of each other.
+The @file{est.fits} table also has a third column with error bars.
+You can follow the same steps before and draw the error bars to see how they
compare with the scatter of the measured data.
+They are much smaller than the error in each point because we had a very good
sampling of the function in our noisy data.
+
+Another useful point with the estimated output file is that it contains all
the fitting outputs as keywords in the header:
+
+@example
+$ astfits est.fits -h1
+...
+ / Fit results
+FITTYPE = 'polynomial-weighted' / Functional form of the fitting.
+FITMAXP = 2 / Maximum power of polynomial.
+FITIN = 'noisy.fits' / Name of file with input columns.
+FITINHDU= '1 ' / Name or Number of HDU with input cols.
+FITXCOL = 'X ' / Name or Number of independent (X) col.
+FITYCOL = 'Y ' / Name or Number of measured (Y) column.
+FITWCOL = 'Yerr ' / Name or Number of weight column.
+FITWNAT = 'Standard deviation' / Nature of weight column.
+FRDCHISQ= 0.974067008958516 / Reduced chi^2 of fit.
+FITC0 = 1.22862116084727 / C0: multiple of x^0 in polynomial
+FITC1 = -4.51277966356177 / C1: multiple of x^1 in polynomial
+FITC2 = 7.84358839431161 / C2: multiple of x^2 in polynomial
+FCOV11 = 0.00104960011629718 / Covariance matrix element (1,1).
+FCOV12 = -0.00399284880859776 / Covariance matrix element (1,2).
+FCOV13 = 0.00283673901863388 / Covariance matrix element (1,3).
+FCOV21 = -0.00399284880859776 / Covariance matrix element (2,1).
+FCOV22 = 0.0175244126670659 / Covariance matrix element (2,2).
+FCOV23 = -0.0138030778380786 / Covariance matrix element (2,3).
+FCOV31 = 0.00283673901863388 / Covariance matrix element (3,1).
+FCOV32 = -0.0138030778380786 / Covariance matrix element (3,2).
+FCOV33 = 0.0128129806394559 / Covariance matrix element (3,3).
+...
+@end example
+
+In scenarios were you don't want the estimation, but only the fitted
parameters, all that verbose, human-friendly text or FITS keywords can be an
annoying extra step.
+For such cases, you should use the @option{--quiet} option like below.
+It will print the parameters, rows of the covariance matrix and
@mymath{\chi_{red}^2} on separate lines with nothing extra.
+This allows you to parse the values in any way that you would like.
+
+@example
+$ aststatistics noisy.fits -cX,Y,Yerr --fit=polynomial-weighted \
+ --fitmaxpower=2 --quiet
++1.2286211608 -4.5127796636 +7.8435883943
++0.0010496001 -0.0039928488 +0.0028367390
+-0.0039928488 +0.0175244127 -0.0138030778
++0.0028367390 -0.0138030778 +0.0128129806
++0.9740670090
+@end example
+
+As a final example, because real data usually have outliers, let's look at the
``robust'' polynomial fit which has special features to remove outliers.
+First, we need to add some outliers to the table.
+To do this, we'll make a plain-text table with @command{echo}, and use Table's
@option{--catrowfile} to concatenate (or append) those two rows to the original
table.
+Finally, we'll run the same fitting step above:
+
+@example
+$ echo "0.6 20 0.01" > outliers.txt
+$ echo "0.8 20 0.01" >> outliers.txt
+
+$ asttable noisy.fits --catrowfile=outliers.txt \
+ --output=with-outlier.fits
+
+$ aststatistics with-outlier.fits -cX,Y,Yerr --fit=polynomial-weighted \
+ --fitmaxpower=2 --fitestimate=self \
+ --output=est-out.fits --quiet
+-13.6446036899 +66.8463258547 -30.8746303591
++0.0007889160 -0.0027706310 +0.0022208939
+-0.0027706310 +0.0113922468 -0.0100306732
++0.0022208939 -0.0100306732 +0.0094087226
++4501.8356719150
+@end example
+
+We see that the coefficient values have changed significantly and that
@mymath{\chi_{red}^2} has increased to @mymath{4501}!
+Recall that a good fit should have @mymath{\chi_{red}^2\approx1}.
+These numbers clearly show that the fit was bad, but again, nothing beats a
visual inspection.
+To visually see the effect of those outliers, let's plot them with the command
below.
+You see that those two points have clearly caused a turn in the fitted result
which is terrible.
+
+@example
+$ astscript-fits-view with-outlier.fits est-out.fits
+@end example
+
+For such cases, GSL has
@url{https://www.gnu.org/software/gsl/doc/html/lls.html#robust-linear-regression,
Robust linear regression}.
+In Gnuastro's Statistics, you can access it with
@option{--fit=polynomial-robust}, like the example below.
+Just note that the robust method doesn't take an error column (because it
estimates the errors internally while rejecting outliers, based on the method):
+
+@example
+$ aststatistics with-outlier.fits -cX,Y --fit=polynomial-robust \
+ --fitmaxpower=2 --fitestimate=self \
+ --output=est-out.fits --quiet
++1.2042269119 -4.4779253576 +7.8498615369
++0.0237588085 -0.0468869673 +0.0193819707
+-0.0468869673 +0.1126156328 -0.0507842538
++0.0193819707 -0.0507842538 +0.0241623033
++0.3428106540
+
+$ astscript-fits-view with-outlier.fits est-out.fits
+@end example
+
+It is clear that the coefficients are very similar to the no-outlier scenario
above and if you run the second command to view the scatter plots on TOPCAT,
you also see that the fit nicely follows the curve and is not affected by those
two points.
+GSL provides many methods to reject outliers.
+For their full list, see the description of @option{--fitrobust} in
@ref{Fitting options}.
+For a description of the outlier rejection methods, see the
@url{https://www.gnu.org/software/gsl/doc/html/lls.html#c.gsl_multifit_robust_workspace,
GSL manual}.
+
+@node Sky value, Invoking aststatistics, Least squares fitting, Statistics
@subsection Sky value
@cindex Sky
@@ -20292,6 +20541,13 @@ $ aststatistics image.fits --sky --kernel=kernel.fits
## column value between 26 and 27:
$ aststatistics cat.fits -cMAG_F160W -g26 -l27 --sigmaclip=3,0.2
+## Find the polynomial (to third order) that best fits the X and Y
+## columns of 'table.fits'. Robust fitting will be used to reject
+## outliers. Also, estimate the fitted polynomial on the same input
+## column (with errors).
+$ aststatistics table.fits --fit=polynomial-robust --fitmaxpower=3 \
+ -cX,Y --fitestimate=self --output=estimated.fits
+
## Print the median value of all records in column MAG_F160W that
## have a value larger than 3 in column PHOTO_Z:
$ aststatistics tab.txt -rPHOTO_Z -g3 -cMAG_F160W --median
@@ -20315,51 +20571,6 @@ Initially, the full dataset will be read, but it is
possible to select a specifi
You can either directly specify a minimum and maximum value for the range of
data elements to use (with @option{--greaterequal} or @option{--lessthan}), or
specify the range using quantiles (with @option{--qrange}).
If a range is specified, all pixels outside of it are ignored before any
processing.
-The following set of options are for specifying the input/outputs of
Statistics.
-There are many other input/output options that are common to all Gnuastro
programs including Statistics, see @ref{Input output options} for those.
-
-@table @option
-
-@item -c STR/INT
-@itemx --column=STR/INT
-The column to use when the input file is a table with more than one column.
-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}.
-
-@item -r STR/INT
-@itemx --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
-Limit the range of inputs into those with values greater and equal to what is
given to this option.
-None of the values below this value will be used in any of the processing
steps below.
-
-@item -l FLT
-@itemx --lessthan=FLT
-Limit the range of inputs into those with values less-than what is given to
this option.
-None of the values greater or equal to this value will be used in any of the
processing steps below.
-
-@item -Q FLT[,FLT]
-@itemx --qrange=FLT[,FLT]
-Specify the range of usable inputs using the quantile.
-This option can take one or two quantiles to specify the range.
-When only one number is input (let's call it @mymath{Q}), the range will be
those values in the quantile range @mymath{Q} to @mymath{1-Q}.
-So when only one value is given, it must be less than 0.5.
-When two values are given, the first is used as the lower quantile range and
the second is used as the larger quantile range.
-
-@cindex Quantile
-The quantile of a given element in a dataset is defined by the fraction of its
index to the total number of values in the sorted input array.
-So the smallest and largest values in the dataset have a quantile of 0.0 and
1.0.
-The quantile is a very useful non-parametric (making no assumptions about the
input) relative measure to specify a range.
-It can best be understood in terms of the cumulative frequency plot, see
@ref{Histogram and Cumulative Frequency Plot}.
-The quantile of each horizontal axis value in the cumulative frequency plot is
the vertical axis value associate with it.
-
-@end table
-
@cindex ASCII plot
When no operation is requested, Statistics will print some general basic
properties of the input dataset on the command-line like the example below (ran
on one of the output images of @command{make check}@footnote{You can try it by
running the command in the @file{tests} directory, open the image with a FITS
viewer and have a look at it to get a sense of how these statistics relate to
the input image/dataset.}).
This default behavior is designed to help give you a general feeling of how
the data are distributed and help in narrowing down your analysis.
@@ -20411,6 +20622,59 @@ This feature makes these options very useful in
scripts, or to redirect into pro
These are some of the most basic measures, Gnuastro is still under heavy
development and this list will grow.
If you want another statistical parameter, please contact us and we will do
out best to add it to this list, see @ref{Suggest new feature}.
+@menu
+* Input to Statistics:: How to specify the inputs to Statistics.
+* Single value measurements:: Can be used together (like --mean, or
--maximum).
+* Generating histograms and cumulative frequency plots:: Histogram and CFP
tables.
+* Fitting options:: Least squares fitting.
+* Contour options:: Table of contours.
+* Statistics on tiles:: Possible to do single-valued measurements on
tiles.
+@end menu
+
+@node Input to Statistics, Single value measurements, Invoking aststatistics,
Invoking aststatistics
+@subsubsection Input to Statistics
+
+The following set of options are for specifying the input/outputs of
Statistics.
+There are many other input/output options that are common to all Gnuastro
programs including Statistics, see @ref{Input output options} for those.
+
+@table @option
+
+@item -c STR/INT
+@itemx --column=STR/INT
+The column to use when the input file is a table with more than one column.
+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}.
+
+@item -g FLT
+@itemx --greaterequal=FLT
+Limit the range of inputs into those with values greater and equal to what is
given to this option.
+None of the values below this value will be used in any of the processing
steps below.
+
+@item -l FLT
+@itemx --lessthan=FLT
+Limit the range of inputs into those with values less-than what is given to
this option.
+None of the values greater or equal to this value will be used in any of the
processing steps below.
+
+@item -Q FLT[,FLT]
+@itemx --qrange=FLT[,FLT]
+Specify the range of usable inputs using the quantile.
+This option can take one or two quantiles to specify the range.
+When only one number is input (let's call it @mymath{Q}), the range will be
those values in the quantile range @mymath{Q} to @mymath{1-Q}.
+So when only one value is given, it must be less than 0.5.
+When two values are given, the first is used as the lower quantile range and
the second is used as the larger quantile range.
+
+@cindex Quantile
+The quantile of a given element in a dataset is defined by the fraction of its
index to the total number of values in the sorted input array.
+So the smallest and largest values in the dataset have a quantile of 0.0 and
1.0.
+The quantile is a very useful non-parametric (making no assumptions about the
input) relative measure to specify a range.
+It can best be understood in terms of the cumulative frequency plot, see
@ref{Histogram and Cumulative Frequency Plot}.
+The quantile of each horizontal axis value in the cumulative frequency plot is
the vertical axis value associate with it.
+
+@end table
+
+@node Single value measurements, Generating histograms and cumulative
frequency plots, Input to Statistics, Invoking aststatistics
+@subsubsection Single value measurements
+
@table @option
@item -n
@@ -20565,6 +20829,9 @@ Standard deviation after applying
@mymath{\sigma}-clipping (see @ref{Sigma clipp
@end table
+@node Generating histograms and cumulative frequency plots, Fitting options,
Single value measurements, Invoking aststatistics
+@subsubsection Generating histograms and cumulative freq.
+
The list of options below are for those statistical operations that output
more than one value.
So while they can be called together in one run, their outputs will be
distinct (each one's output will usually be printed in more than one line).
@@ -20730,23 +20997,181 @@ Similar to @option{--onebinstart}, but for the
second column when a 2D histogram
@end table
+@node Fitting options, Contour options, Generating histograms and cumulative
frequency plots, Invoking aststatistics
+@subsubsection Fitting options
-All the options described until now were from the first class of operations
discussed above: those that treat the whole dataset as one.
-However, it often happens that the relative position of the dataset elements
over the dataset is significant.
-For example you do not want one median value for the whole input image, you
want to know how the median changes over the image.
-For such operations, the input has to be tessellated (see @ref{Tessellation}).
-Thus this class of options cannot currently be called along with the options
above in one run of Statistics.
+With the options below, you can customize the least squares fitting features
of Statistics.
+For a tutorial of the usage of least squares fitting in Statistics, please see
@ref{Least squares fitting}.
+Here, we will just reivew the details of each option.
+
+To activate least squares fitting in Statistics, it is necessary to use the
@option{--fit} option to specify the type of fit you want to do.
+See the description of @option{--fit} for the various available fitting models.
+The fitting models that account for weights require three input columns, while
the non-weighted ones only take two input columns.
+Here is a summary of the input columns:
+
+@enumerate
+@item
+The first input column is assumed to be the independent variable (on the
horizontal axis of a plot, or @mymath{X} in the equations of each fit).
+@item
+The second input column is assumed to be the measured value (on the vertical
axis of a plot, or @mymath{Y} in the equation above).
+@item
+The third input column is only for fittings with a weight.
+It is assumed to be the ``weight'' of the measurement column.
+The nature of the ``weight'' can be set with the @option{--fitweight} option,
for example if you have the standard deviation of the error in @mymath{Y}, you
can use @option{--fitweight=std} (which is the default, so unless the default
value has been changed, you will not need to set this).
+@end enumerate
+
+If three columns are given to a model without weight, or two columns are given
to a model that requires weights, Statistics will abort and inform you.
+Below you can see an example of fitting with the same linear model, once
weighted and once without weights.
+
+@example
+$ aststatistics table.fits --column=X,Y --fit=linear
+$ aststatistics table.fits --column=X,Y,Yerr --fit=linear-weighted
+@end example
+
+The output of the fitting can be in three modes listed below.
+For a complete example, see the tutorial in @ref{Least squares fitting}).
+@table @asis
+@item Human friendly format
+By default (for example the commands above) the output is an elaborate
description of the model parameters.
+For example @mymath{c_0} and @mymath{c_1} in the linear model
(@mymath{Y=c_0+c_1X}).
+Their covariance matrix and the reduced @mymath{\chi^2} of the fit are also
printed on the output.
+@item Raw numbers
+If you don't need the human friendly components of the output (which are
annoying when you want to parse the outputs in some scenarios), you can use
@option{--quiet} option.
+Only the raw output numbers will be printed.
+@item Estimate on a custom X column
+Through the @option{--fitestimate} option, you can specify an independent
table column to estimate the fit (it can also take a single value).
+See the description of this option for more.
+@end table
@table @option
+@item -f STR
+@itemx --fit=STR
+The name of the fitting method to use.
+They are based on the @url{https://www.gnu.org/software/gsl/doc/html/lls.html,
linear} and @url{https://www.gnu.org/software/gsl/doc/html/nls.html, nonlinear}
least-squares fitting functions of the GNU Scientific Library (GSL).
+@table @code
+@item linear
+@mymath{Y=c_0+c_1X}
+@item linear-weighted
+@mymath{Y=c_0+c_1X}; accounting for ``weights'' in @mymath{Y}.
+@item linear-no-constant
+@mymath{Y=c_1X}.
+@item linear-no-constant-weighted
+@mymath{Y=c_1X}; accounting for ``weights'' in @mymath{Y}.
+@item polynomial
+@mymath{Y=c_0+c_1X+c_2X^2+\cdots+c_nX^n}; the maximum required power
(@mymath{n}) is specifed by @option{--fitmaxpower}.
+@item polynomial-weighted
+@mymath{Y=c_0+c_1X+c_2X^2+\cdots+c_nX^n}; accounting for ``weights'' in
@mymath{Y}.
+The maximum required power (@mymath{n}) is specifed by @option{--fitmaxpower}.
+@item polynomial-robust
+@cindex Robust polynomial fit
+@cindex Polynomial fit (robust)
+@mymath{Y=c_0+c_1X+c_2X^2+\cdots+c_nX^n}; rejects outliers.
+The function to use for outlier removal can be specified with the
@option{--fitrobust} option described below.
+This model doesn't take weights since they are calculated internally based on
the outlier removal function (requires two input columns).
+The maximum required power (@mymath{n}) is specifed by @option{--fitmaxpower}.
+
+For a comprehensive reivew of ``robust'' fitting and the available functions,
please see the
@url{https://www.gnu.org/software/gsl/doc/html/lls.html#robust-linear-regression,
Robust linear regression} section of the GNU Scientific Library.
+@end table
-@item -t
-@itemx --ontile
-Do the respective single-valued calculation over one tile of the input
dataset, not the whole dataset.
-This option must be called with at least one of the single valued options
discussed above (for example @option{--mean} or @option{--quantile}).
-The output will be a file in the same format as the input.
-If the @option{--oneelempertile} option is called, then one element/pixel will
be used for each tile (see @ref{Processing options}).
-Otherwise, the output will have the same size as the input, but each element
will have the value corresponding to that tile's value.
-If multiple single valued operations are called, then for each operation there
will be one extension in the output FITS file.
+@item --fitweight=STR
+The nature of the ``weight'' column (when a weight is necessary for the model).
+It can take one of the following values:
+@table @code
+@item std
+Standard deviation of each @mymath{Y} axis measurement: this is the usual
``error'' associated with a measurement (for example in @ref{MakeCatalog}) and
is the default value to this option.
+@item var
+Variance of each @mymath{Y} axis measurement.
+Assuming a Gaussian distribution with standard deviation @mymath{\sigma}, the
variance is @mymath{\sigma^2}.
+@item inv-var
+Inverse variance of each @mymath{Y} axis measurement.
+Assuming a Gaussian distribution with standard deviation @mymath{\sigma}, the
variance is @mymath{1/\sigma^2}.
+@end table
+
+@item --fitmaxpower=INT
+The maximum power (an integer) in a polynomial (@mymath{n} in
@mymath{Y=c_0+c_1X+c_2X^2+\cdots+c_nX^n}).
+This is only relevant when one of the polynomial models is given to
@option{--fit}.
+The fit will return @mymath{n+1} coefficients.
+
+@item --fitrobust=STR
+The function for rejecting outliers in the @code{polynomial-robust} fitting
model.
+For a comprehensive reivew of ``robust'' fitting and the available functions,
please see the
@url{https://www.gnu.org/software/gsl/doc/html/lls.html#robust-linear-regression,
Robust linear regression} section of the GNU Scientific Library.
+This function can take the following values:
+@table @code
+@item bisquare
+@cindex Tukey’s biweight (bisquare) function
+@cindex Biweight function of Tukey
+@cindex Bisquare function of Tukey
+Tukey’s biweight (bisquare) function, this is the default function.
+According to the GSL manual, this is a good general purpose weight function.
+@item cauchy
+@cindex Cauchy's function (robust weight)
+@cindex Lorentzian function (robust weight)
+Cauchy’s function (also known as the Lorentzian function).
+It doesn't guarantee a unique solution, so it should be used with care.
+@item fair
+@cindex Fair function (robust weight)
+The fair function.
+It guarantees a unique solution and has continuous derivatives to three orders.
+@item huber
+@cindex Huber function (robust weight)
+Huber's @mymath{\rho} function.
+This is also a good general purpose weight function for rejecting outliers,
but can cause difficulty in some special scenarios.
+@item ols
+Ordinary Least Squares (OLS) solution with a constant weight of unity.
+@item welsch
+@cindex Welsch function (robust weight)
+Welsch function which is useful when the residuals follow an exponential
distribution.
+@end table
+
+@item --fitestimate=STR/FLT
+Estimate the fitted function at a single point or a complete column of points.
+The input @mymath{X} axis positions to estimate the function can be specified
in the following ways:
+@itemize
+@item
+A real number: the fitted function will be estimated at that @mymath{X}
position and the corresponding @mymath{Y} and its error will be printed to
standard output.
+@item
+@code{self}: in this mode, the same X axis column that was used in the fit
will be used for estimating the fitted function.
+This can be useful to visually/easily check the fit, see @ref{Least squares
fitting}.
+@item
+A file name: If the value is none of the above, Statistics expects it to be a
file name containing a table.
+If the file is a FITS file, the HDU containing the table should be specified
with the @option{--fitestimatehdu} option.
+The column of the table to use for the @mymath{X} axis points should be
specified with the @option{--fitestimatecol} option.
+@end itemize
+The output in this mode can be customized in the following ways:
+@itemize
+@item
+If a single floating point value is given @option{--fitestimate}, the fitted
function will be estimated on that point and printed to standard output.
+@item
+When nothing is given to @option{--output}, the independent column and the
estimated values and errors are printed on the standard output.
+@item
+If a file name is given to @option{--output}, the estimated table above is
saved in that file.
+It can have any of the formats in @ref{Recognized table formats}.
+As a FITS file, all the fitted numbers (coefficients, covariance matrix and
reduced @mymath{\chi^2} are kept as FITS keywords in the same HDU of the
estimated table.
+For a complete example, see @ref{Least squares fitting}.
+@item
+When @option{--quiet} is given with @option{--fitestimate}, the fitted
parameters are no longer printed on the standard output (they will be available
as FITS keywords in the file given to @option{--output}).
+@end itemize
+
+@item --fitestimatehdu=STR/INT
+HDU name or counter (counting from zero) that contains the table to be used
for the estimating the fitted function over many points through
@option{--fitestimate}.
+For more on selecting a HDU, see the description of @option{--hdu} in
@ref{Input output options}.
+
+@item --fitestimatecol=STR/INT
+Column name or counter (counting from one) that contains the table to be used
for the estimating the fitted function over many points through
@option{--fitestimate}.
+See @ref{Selecting table columns}.
+@end table
+
+
+
+
+
+@node Contour options, Statistics on tiles, Fitting options, Invoking
aststatistics
+@subsubsection Contour options
+
+Contours are useful to highlight the 2D shape of a certain flux level over an
image.
+To derive contours in Statistics, you can use the option below:
+
+@table @option
@item -R FLT[,FLT[,FLT...]]
@itemx --contour=FLT[,FLT[,FLT...]]
@@ -20760,6 +21185,27 @@ If any other format can be useful for your work please
let us know so we can add
If the image has World Coordinate System information, the written coordinates
will be in RA and Dec, otherwise, they will be in pixel coordinates.
Note that currently, this is a very crude/simple implementation, please let us
know if you find problematic situations so we can fix it.
+@end table
+
+@node Statistics on tiles, , Contour options, Invoking aststatistics
+@subsubsection Statistics on tiles
+
+All the options described until now were from the first class of operations
discussed above: those that treat the whole dataset as one.
+However, it often happens that the relative position of the dataset elements
over the dataset is significant.
+For example you do not want one median value for the whole input image, you
want to know how the median changes over the image.
+For such operations, the input has to be tessellated (see @ref{Tessellation}).
+Thus this class of options cannot currently be called along with the options
above in one run of Statistics.
+
+@table @option
+
+@item -t
+@itemx --ontile
+Do the respective single-valued calculation over one tile of the input
dataset, not the whole dataset.
+This option must be called with at least one of the single valued options
discussed above (for example @option{--mean} or @option{--quantile}).
+The output will be a file in the same format as the input.
+If the @option{--oneelempertile} option is called, then one element/pixel will
be used for each tile (see @ref{Processing options}).
+Otherwise, the output will have the same size as the input, but each element
will have the value corresponding to that tile's value.
+If multiple single valued operations are called, then for each operation there
will be one extension in the output FITS file.
@item -y
@itemx --sky
@@ -24256,8 +24702,8 @@ The last three are the three Euler angles in units of
degrees in the ZXZ order a
-@node Modeling and fittings, High-level calculations, Data analysis, Top
-@chapter Modeling and fitting
+@node Data modeling, High-level calculations, Data analysis, Top
+@chapter Data modeling
@cindex Fitting
@cindex Modeling
@@ -24272,7 +24718,7 @@ The tools in this chapter create model galaxies and
will provide 2D fittings to
-@node MakeProfiles, MakeNoise, Modeling and fittings, Modeling and fittings
+@node MakeProfiles, MakeNoise, Data modeling, Data modeling
@section MakeProfiles
@cindex Checking detection algorithms
@@ -25333,7 +25779,7 @@ If an individual image was created, this column will
have a value of @code{1}, o
-@node MakeNoise, , MakeProfiles, Modeling and fittings
+@node MakeNoise, , MakeProfiles, Data modeling
@section MakeNoise
@cindex Noise
@@ -25668,7 +26114,7 @@ This option will be most useful if the input images
were of integer types.
-@node High-level calculations, Installed scripts, Modeling and fittings, Top
+@node High-level calculations, Installed scripts, Data modeling, Top
@chapter High-level calculations
After the reduction of raw data (for example with the programs in @ref{Data
manipulation}) you will have reduced images/data ready for processing/analyzing
(for example with the programs in @ref{Data analysis}).
@@ -28833,6 +29279,7 @@ If you use the Info version of this manual (see
@ref{Info}), you do not have to
* Permutations:: Re-order (or permute) the values in a dataset.
* Matching:: Matching catalogs based on position.
* Statistical operations:: Functions for basic statistics.
+* Fitting functions:: Fit independent and measured variables.
* Binary datasets:: Datasets that can only have values of 0 or 1.
* Labeled datasets:: Working with Segmented/labeled datasets.
* Convolution functions:: Library functions to do convolution.
@@ -30204,71 +30651,52 @@ Free all the individual datasets within the elements
of @code{dataptr}, then fre
@node Copying datasets, , Arrays of datasets, Library data container
@subsubsection Copying datasets
-The functions in this section describes Gnuastro's facilities to copy a
-given dataset into another. The new dataset can have a different type
-(including a string), it can be already allocated (in which case only the
-values will be written into it). In all these cases, if the input dataset
-is a tile or a list, only the data within the given tile, or the given node
-in a list, are copied. If the input is a list, the @code{next} pointer will
-also be copied to the output, see @ref{List of gal_data_t}.
+The functions in this section describes Gnuastro's facilities to copy a given
dataset into another.
+The new dataset can have a different type (including a string), it can be
already allocated (in which case only the values will be written into it).
+In all these cases, if the input dataset is a tile or a list, only the data
within the given tile, or the given node in a list, are copied.
+If the input is a list, the @code{next} pointer will also be copied to the
output, see @ref{List of gal_data_t}.
-In many of the functions here, it is possible to copy the dataset to a new
-numeric data type (see @ref{Numeric data types}. In such cases, Gnuastro's
-library is going to use the native conversion by C. So if you are
-converting to a smaller type, it is up to you to make sure that the values
-fit into the output type.
+In many of the functions here, it is possible to copy the dataset to a new
numeric data type (see @ref{Numeric data types}.
+In such cases, Gnuastro's library is going to use the native conversion by C.
+So if you are converting to a smaller type, it is up to you to make sure that
the values fit into the output type.
@deftypefun {gal_data_t *} gal_data_copy (gal_data_t @code{*in})
-Return a new dataset that is a copy of @code{in}, all of @code{in}'s
-meta-data will also copied into the output, except for @code{block}. If the
-dataset is a tile/list, only the given tile/node will be copied, the
-@code{next} pointer will also be copied however.
+Return a new dataset that is a copy of @code{in}, all of @code{in}'s meta-data
will also copied into the output, except for @code{block}.
+If the dataset is a tile/list, only the given tile/node will be copied, the
@code{next} pointer will also be copied however.
@end deftypefun
@deftypefun {gal_data_t *} gal_data_copy_to_new_type (gal_data_t @code{*in},
uint8_t @code{newtype})
-Return a copy of the dataset @code{in}, converted to @code{newtype}, see
-@ref{Library data types} for Gnuastro library's type identifiers. The
-returned dataset will have all meta-data except their type and @code{block}
-equal to the input's metadata. If the dataset is a tile/list, only the
-given tile/node will be copied, the @code{next} pointer will also be copied
-however.
+Return a copy of the dataset @code{in}, converted to @code{newtype}, see
@ref{Library data types} for Gnuastro library's type identifiers.
+The returned dataset will have all meta-data except their type and
@code{block} equal to the input's metadata.
+If the dataset is a tile/list, only the given tile/node will be copied, the
@code{next} pointer will also be copied however.
@end deftypefun
@deftypefun {gal_data_t *} gal_data_copy_to_new_type_free (gal_data_t
@code{*in}, uint8_t @code{newtype})
-Return a copy of the dataset @code{in} that is converted to @code{newtype}
-and free the input dataset. See @ref{Library data types} for Gnuastro
-library's type identifiers. The returned dataset will have all meta-data,
-except their type, equal to the input's metadata (including
-@code{next}). Note that if the input is a tile within a larger block, it
-will not be freed. This function is similar to
-@code{gal_data_copy_to_new_type}, except that it will free the input
-dataset.
+Return a copy of the dataset @code{in} that is converted to @code{newtype} and
free the input dataset.
+See @ref{Library data types} for Gnuastro library's type identifiers.
+The returned dataset will have all meta-data, except their type, equal to the
input's metadata (including @code{next}).
+Note that if the input is a tile within a larger block, it will not be freed.
+This function is similar to @code{gal_data_copy_to_new_type}, except that it
will free the input dataset.
@end deftypefun
@deftypefun {void} gal_data_copy_to_allocated (gal_data_t @code{*in},
gal_data_t @code{*out})
-Copy the contents of the array in @code{in} into the already allocated
-array in @code{out}. The types of the input and output may be different,
-type conversion will be done internally. When @code{in->size != out->size}
-this function will behave as follows:
+Copy the contents of the array in @code{in} into the already allocated array
in @code{out}.
+The types of the input and output may be different, type conversion will be
done internally.
+When @code{in->size != out->size} this function will behave as follows:
@table @code
@item out->size < in->size
-This function will not re-allocate the necessary space, it will abort with an
-error, so please check before calling this function.
+This function will not re-allocate the necessary space, it will abort with an
error, so please check before calling this function.
@item out->size > in->size
-This function will write the values in @code{out->size} and
-@code{out->dsize} from the same values of @code{in}. So if you want to use
-a pre-allocated space/dataset multiple times with varying input sizes, be
-sure to reset @code{out->size} before every call to this function.
+This function will write the values in @code{out->size} and @code{out->dsize}
from the same values of @code{in}.
+So if you want to use a pre-allocated space/dataset multiple times with
varying input sizes, be sure to reset @code{out->size} before every call to
this function.
@end table
@end deftypefun
@deftypefun {gal_data_t *} gal_data_copy_string_to_number (char @code{*string})
-Read @code{string} into the smallest type that can store the value (see
-@ref{Numeric data types}). This function is just a wrapper for the
-@code{gal_type_string_to_number}, but will put the value into a
-single-element dataset.
+Read @code{string} into the smallest type that can store the value (see
@ref{Numeric data types}).
+This function is just a wrapper for the @code{gal_type_string_to_number}, but
will put the value into a single-element dataset.
@end deftypefun
@@ -34755,7 +35183,7 @@ If internal allocation is necessary and the space is
larger than @code{minmapsiz
@end deftypefun
-@node Statistical operations, Binary datasets, Matching, Gnuastro library
+@node Statistical operations, Fitting functions, Matching, Gnuastro library
@subsection Statistical operations (@file{statistics.h})
After reading a dataset into memory from a file or fully simulating it with
@@ -35187,7 +35615,185 @@ input dataset, so the input may be altered after this
function.
@end deftypefun
-@node Binary datasets, Labeled datasets, Statistical operations, Gnuastro
library
+
+
+
+
+
+
+
+
+@node Fitting functions, Binary datasets, Statistical operations, Gnuastro
library
+@subsection Fitting functions (@file{fit.h})
+
+@cindex Fitting
+@cindex Least squares fitting
+
+After doing a measurement, it is usually necessary to parametrize the relation
that has been found.
+The functions in this section are wrappers over the GNU Scientific Library
(GSL) @url{https://www.gnu.org/software/gsl/doc/html/lls.html, Linear
Least-Squares Fitting}, to make them easily accessible using Gnuastro's
@ref{Generic data container}.
+The respective GSL function is mentioned under each function.
+
+@deffn {Global integer} GAL_FIT_INVALID
+@deffnx {Global integer} GAL_FIT_LINEAR
+@deffnx {Global integer} GAL_FIT_LINEAR_WEIGHTED
+@deffnx {Global integer} GAL_FIT_LINEAR_NO_CONSTANT
+@deffnx {Global integer} GAL_FIT_LINEAR_NO_CONSTANT_WEIGHTED
+@deffnx {Global integer} GAL_FIT_POLYNOMIAL
+@deffnx {Global integer} GAL_FIT_POLYNOMIAL_WEIGHTED
+@deffnx {Global integer} GAL_FIT_POLYNOMIAL_NUMBER
+Identifiers for the various types of fitting functions.
+These can be used by the callers of these functions to select between various
fitting types.
+They can easily be converted to, and from, fixed human-readable strings using
the @code{gal_fit_name_*} functions below.
+The last one @code{GAL_FIT_ROBUST_NUMBER} is the total number of available
fitting methods (can be used to add more macros in the calling program and to
avoid overlaps with existing codes).
+@end deffn
+
+@deffn {Global integer} GAL_FIT_ROBUST_INVALID
+@deffnx {Global integer} GAL_FIT_ROBUST_DEFAULT
+@deffnx {Global integer} GAL_FIT_ROBUST_BISQUARE
+@deffnx {Global integer} GAL_FIT_ROBUST_CAUCHY
+@deffnx {Global integer} GAL_FIT_ROBUST_FAIR
+@deffnx {Global integer} GAL_FIT_ROBUST_HUBER
+@deffnx {Global integer} GAL_FIT_ROBUST_OLS
+@deffnx {Global integer} GAL_FIT_ROBUST_WELSCH
+@deffnx {Global integer} GAL_FIT_ROBUST_NUMBER
+Identifiers for the various types of robust polynomial fitting functions.
+For a description of each, see
@url{https://www.gnu.org/s/gsl/doc/html/lls.html#c.gsl_multifit_robust_alloc}.
+The last one @code{GAL_FIT_ROBUST_NUMBER} is the total number of available
functions (can be used to add more macros in the calling program and to avoid
overlaps with existing codes).
+@end deffn
+
+@deftypefun uint8_t gal_fit_name_to_id (char @code{*name})
+Return the internal code of a standard human-readable name for the various
fitting functions.
+If the name is not recognized, the returned value will be
@code{GAL_FIT_INVALID}.
+@end deftypefun
+
+@deftypefun {char *} gal_fit_name_from_id (uint8_t @code{fitid})
+Return a standard human-readable name for the fitting function identified with
the @code{fitid} (read as ``fitting ID'').
+If the fitting ID couldn't be recognized, a NULL pointer is returned.
+@end deftypefun
+
+@deftypefun uint8_t gal_fit_name_robust_to_id (char @code{*name})
+Return the internal code of a standard human-readable name for the various
robust fitting types.
+If the name is not recognized, the returned value will be
@code{GAL_FIT_INVALID}.
+@end deftypefun
+
+@deftypefun {char *} gal_fit_name_robust_from_id (uint8_t @code{robustid})
+Return a standard human-readable name for the input robust fitting type.
+If the fitting ID couldn't be recognized, a NULL pointer is returned.
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_fit_1d_linear (gal_data_t @code{*xin},
gal_data_t @code{*yin}, gal_data_t @code{*ywht})
+@cindex Weight (in fitting)
+Preform a 1D linear regression fit with a constant
term@footnote{@url{https://www.gnu.org/s/gsl/doc/html/lls.html#linear-regression-with-a-constant-term}}
in the form of @mymath{Y=c_0+c_1X}.
+The input @code{xin} contains the independent variable values and @code{yin}
contains the measured variable values for each independent variable.
+When @code{ywht!=NULL}, it is assumed to contain the ``weight'' of each Y
measurement (if you don't have weights on your measured values, simply set this
to @code{NULL}).
+The weight of each measurement is the inverse of its variance.
+For a Gaussian error distribution with standard deviation @mymath{\sigma}, the
weight is therefore @mymath{1/\sigma^2}.
+
+If any of the values in any of the inputs is blank (NaN in floating point),
the final fitted parameters will all be NaN.
+To remove rows with a NaN/blank, you can use @code{gal_blank_remove_rows}
(which will remove all rows with a blank values in any of the columns with a
single call).
+
+@cindex Chi-squared
+@cindex Covariance matrix
+@cindex Matrix (covariance)
+@cindex Variance-covariance matrix
+The output is a single dataset with a @code{GAL_TYPE_FLOAT64} type with 6
elements:
+@enumerate
+@item
+@mymath{c_0}: the constant in @mymath{Y=c_0+c_1X}.
+@item
+@mymath{c_1}: the multiple in @mymath{Y=c_0+c_1X}.
+@item
+First element of variance-covariance matrix.
+@item
+Second and third (which are equal) elements of the variance-covariance matrix.
+@item
+Fourth element of the variance-covariance matrix.
+@item
+The reduced @mymath{\chi^2} of the fit.
+@end enumerate
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_fit_1d_linear_no_constant (gal_data_t
@code{*xin}, gal_data_t @code{*yin}, gal_data_t @code{*ywht})
+@cindex Weight (in fitting)
+Preform a 1D linear regression fit @emph{without} a constant
term@footnote{@url{https://www.gnu.org/s/gsl/doc/html/lls.html#linear-regression-without-a-constant-term}},
formally: @mymath{Y=c_1X}.
+The input @code{xin} contains the independent variable values and @code{yin}
contains the measured variable values for each independent variable.
+When @code{ywht!=NULL}, it is assumed to contain the ``weight'' of each Y
measurement (if you don't have weights on your measured values, simply set this
to @code{NULL}).
+The weight of each measurement is the inverse of its variance.
+For a Gaussian error distribution with standard deviation @mymath{\sigma}, the
weight is therefore @mymath{1/\sigma^2}.
+
+If any of the values in any of the inputs is blank (NaN in floating point),
the final fitted parameters will all be NaN.
+To remove rows with a NaN/blank, you can use @code{gal_blank_remove_rows}
(which will remove all rows with a blank values in any of the columns with a
single call).
+
+The output is a single dataset with a @code{GAL_TYPE_FLOAT64} type with 3
elements:
+@enumerate
+@item
+@mymath{c_1}: the multiple in @mymath{Y=c_0+c_1X}.
+@item
+Variance of @mymath{c_1}.
+@item
+The reduced @mymath{\chi^2} of the fit.
+@end enumerate
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_fit_1d_linear_estimate (gal_data_t @code{*fit},
gal_data_t @code{*xin})
+Given a linear least squares fit output (@code{fit}), estimate the fit on an
arbitrary number of independent variable (horizontal axis, or X, in an X-Y
plot) within @code{xin}.
+@code{fit} is assumed to be the output of either @code{gal_fit_1d_linear} or
@code{gal_fit_1d_linear_no_constant}.
+In case you haven't used those functions to obtain the constants and
covariance matrix elements, see the description of those functions for the
expected format of @code{fit}.
+
+This function returns two columns (as a @ref{List of gal_data_t}): The top
node of the list is the estimated values at the input X-axis positions, and the
next node is the erros in the estimation.
+Naturally, both have the same number of elements as @code{xin}.
+Being a list, helps in easily printing the output columns to a table (see
@ref{Table input output}).
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_fit_1d_polynomial (gal_data_t @code{*xin},
gal_data_t @code{*yin}, gal_data_t @code{*ywht}, size_t @code{maxpower}, double
@code{*redchisq})
+@cindex Polynomial fit
+@cindex Fitting (polynomial)
+Preform a 1D polynomial fit, formally:
@mymath{Y=c+0+c_1X+c_2X^2+\cdots+c_nX^n} (using GSL's multi-parameter
regression@footnote{@url{https://www.gnu.org/s/gsl/doc/html/lls.html#multi-parameter-regression}}).
+The largest power of @mymath{X} is determined with the @code{maxpower}
argument (which is @mymath{n} in the equation above).
+The reduced @mymath{\chi^2} of the fit is written in the space that
@code{*redchisq} points to.
+
+The input @code{xin} contains the independent variable values and the input
@code{yin} contains the measured variable values for each independent variable.
+When @code{ywht!=NULL}, it is assumed to contain the ``weight'' of each Y
measurement (if you don't have weights on your measured values, simply set this
to @code{NULL}).
+The weight of each measurement is the inverse of its variance.
+For a Gaussian error distribution with standard deviation @mymath{\sigma}, the
weight is therefore @mymath{1/\sigma^2}.
+
+If any of the values in any of the inputs is blank (NaN in floating point),
the final fitted parameters will all be NaN.
+To remove rows with a NaN/blank, you can use @code{gal_blank_remove_rows}
(which will remove all rows with a blank values in any of the columns with a
single call).
+
+The output of this function is a list of two datasets, linked as a list (as a
@ref{List of gal_data_t}).
+Both have a @code{GAL_TYPE_FLOAT64} type, and are described below (in order).
+@enumerate
+@item
+A one dimensional and contains @mymath{n+1} elements (for the @mymath{n+1}
constants that have been found @mymath{(c_0, c_1, c_2, \cdots, c_n)}.
+@item
+A two dimensional variance-covariance matrix with @mymath{(n+1)\times(n+1)}
elements.
+@end enumerate
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_fit_1d_polynomial_robust (gal_data_t
@code{*xin}, gal_data_t @code{*yin}, size_t @code{maxpower}, uint8_t
@code{robustid}, double @code{*redchisq})
+@cindex Robust Polynomial fit
+Preform a 1D robust polynomial fit, formally:
@mymath{Y=c+0+c_1X+c_2X^2+\cdots+c_nX^n} (using GSL's robust linear
regression@footnote{@url{https://www.gnu.org/software/gsl/doc/html/lls.html#robust-linear-regression}}).
+See the description there for the details.
+
+The inputs and outputs of this function are almost identical to
@code{gal_fit_1d_polynomial}, with the difference that you need to specify the
function to reject outliers through the @code{robustid} input argument.
+You can pass any of the @code{GAL_FIT_ROBUST_*} codes defined at the top of
this section to this (the names are identical to the names in GSL).
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_fit_1d_polynomial_estimate (gal_data_t
@code{*fit}, gal_data_t @code{*xin})
+Given a 1D polynomial fit output (@code{fit}), estimate the fit on an
arbitrary number of independent variable (horizontal axis, or X, in an X-Y
plot) within @code{xin}.
+@code{fit} is assumed to be the output of @code{gal_fit_1d_polynomial}.
+In case you haven't used this function to obtain the constants and covariance
matrix, see the description of that function for the expected format of
@code{fit}.
+
+This function returns two columns (as a @ref{List of gal_data_t}): The top
node of the list is the estimated values at the input X-axis positions, and the
next node is the erros in the estimation.
+Naturally, both have the same number of elements as @code{xin}.
+Being a list, helps in easily printing the output columns to a table (see
@ref{Table input output}).
+@end deftypefun
+
+
+
+
+
+@node Binary datasets, Labeled datasets, Fitting functions, Gnuastro library
@subsection Binary datasets (@file{binary.h})
@cindex Thresholding
diff --git a/lib/Makefile.am b/lib/Makefile.am
index 94513bb3..440b2f88 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -116,6 +116,7 @@ libgnuastro_la_SOURCES = \
data.c \
ds9.c \
eps.c \
+ fit.c \
fits.c \
git.c \
interpolate.c \
@@ -166,6 +167,7 @@ pkginclude_HEADERS = gnuastro/config.h \
$(headersdir)/dimension.h \
$(headersdir)/ds9.h \
$(headersdir)/eps.h \
+ $(headersdir)/fit.h \
$(headersdir)/fits.h \
$(headersdir)/git.h \
$(headersdir)/interpolate.h \
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 2bec4b75..6abac87d 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -811,24 +811,26 @@ static gal_data_t *
arithmetic_mknoise(int operator, int flags, gal_data_t *in,
gal_data_t *arg)
{
+ size_t i;
gsl_rng *rng;
const char *rng_name;
- double *d, *df, arg_v;
unsigned long rng_seed;
gal_data_t *out, *targ;
+ double *d, *aarr, arg_v;
/* The dataset may be empty. In this case, the output should also be
empty (we can have tables and images with 0 rows or pixels!). */
if(in->size==0 || in->array==NULL) return in;
/* Sanity checks. */
- if(arg->size!=1)
+ if(arg->size!=1 && arg->size!=in->size)
error(EXIT_FAILURE, 0, "the first popped operand to the '%s' "
- "operator should be a single number (specifying the fixed "
- "sigma, or background level for Poisson noise), but it "
- "has %zu elements, in %zu dimension(s)",
- gal_arithmetic_operator_string(operator), arg->size,
- arg->ndim);
+ "operator should either be a single number or have the "
+ "same number of elements as the input data. Recall that "
+ "it specifies the fixed sigma, or background level for "
+ "Poisson noise). However, it has %zu elements, in %zu "
+ "dimension(s)", gal_arithmetic_operator_string(operator),
+ arg->size, arg->ndim);
if(in->type==GAL_TYPE_STRING)
error(EXIT_FAILURE, 0, "the input dataset to the '%s' operator "
"should have a numerical data type (integer or float), but "
@@ -845,46 +847,49 @@ arithmetic_mknoise(int operator, int flags, gal_data_t
*in,
{ gal_data_free(in); in=NULL; }
}
targ=gal_data_copy_to_new_type(arg, GAL_TYPE_FLOAT64);
- arg_v=((double *)(targ->array))[0];
- gal_data_free(targ);
-
- /* Make sure the noise identifier is positive. */
- if(arg_v<0)
- error(EXIT_FAILURE, 0, "the noise identifier (sigma for "
- "'mknoise-sigma', background for 'mknoise-poisson', or "
- "range for 'mknoise-uniform') must be positive (it is %g)",
- arg_v);
+ aarr=targ->array;
/* Initialize the GSL random number generator. */
rng=arithmetic_gsl_initialize(flags, &rng_name, &rng_seed,
gal_arithmetic_operator_string(operator));
/* Add the noise. */
- df=(d=out->array)+out->size;
- switch(operator)
+ d=out->array;
+ for(i=0;i<out->size;++i)
{
- case GAL_ARITHMETIC_OP_MKNOISE_SIGMA:
- do *d += gsl_ran_gaussian(rng, arg_v); while(++d<df);
- break;
- case GAL_ARITHMETIC_OP_MKNOISE_POISSON:
- do
- *d += arg_v + gsl_ran_gaussian(rng, sqrt( arg_v + *d ));
- while(++d<df);
- break;
- case GAL_ARITHMETIC_OP_MKNOISE_UNIFORM:
- do
- *d += ( (gsl_rng_uniform(rng)*arg_v) - (arg_v/2) );
- while(++d<df);
- break;
- default:
- error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
- "fix the problem. The operator code %d isn't recognized "
- "in this function", __func__, PACKAGE_BUGREPORT, operator);
+ /* Set the argument value. */
+ arg_v = arg->size==1 ? aarr[0] : aarr[i];
+
+ /* Make sure the noise identifier is positive. */
+ if(arg_v<0)
+ error(EXIT_FAILURE, 0, "the noise identifier (sigma for "
+ "'mknoise-sigma', background for 'mknoise-poisson', or "
+ "range for 'mknoise-uniform') must be positive (it is %g)",
+ arg_v);
+
+ /* Do the necessary operation. */
+ switch(operator)
+ {
+ case GAL_ARITHMETIC_OP_MKNOISE_SIGMA:
+ d[i] += gsl_ran_gaussian(rng, arg_v);
+ break;
+ case GAL_ARITHMETIC_OP_MKNOISE_POISSON:
+ d[i] += arg_v + gsl_ran_gaussian(rng, sqrt( arg_v + *d ));
+ break;
+ case GAL_ARITHMETIC_OP_MKNOISE_UNIFORM:
+ d[i] += ( (gsl_rng_uniform(rng)*arg_v) - (arg_v/2) );
+ break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The operator code %d isn't recognized "
+ "in this function", __func__, PACKAGE_BUGREPORT, operator);
+ }
}
/* Clean up and return */
- gsl_rng_free(rng);
if(flags & GAL_ARITHMETIC_FLAG_FREE) gal_data_free(arg);
+ gal_data_free(targ);
+ gsl_rng_free(rng);
return out;
}
diff --git a/lib/blank.c b/lib/blank.c
index 8444d10e..4de4de1a 100644
--- a/lib/blank.c
+++ b/lib/blank.c
@@ -790,8 +790,8 @@ gal_blank_remove(gal_data_t *input)
/* This function currently assumes a contiguous patch of memory. */
if(input->block && input->ndim!=1 )
- error(EXIT_FAILURE, 0, "%s: tiles in datasets with more dimensions than "
- "1 are not yet supported. Your input has %zu dimensions",
+ error(EXIT_FAILURE, 0, "%s: tiles in datasets with more dimensions "
+ "than 1 are not yet supported. Your input has %zu dimensions",
__func__, input->ndim);
/* If the dataset doesn't have blank values, then just get the size. */
@@ -829,7 +829,7 @@ gal_blank_remove(gal_data_t *input)
input->dsize[0]=input->size=num;
/* Set the flags to mark that there is no blanks. */
- input->flag |= GAL_DATA_FLAG_BLANK_CH;
+ input->flag |= GAL_DATA_FLAG_BLANK_CH;
input->flag &= ~GAL_DATA_FLAG_HASBLANK;
}
@@ -913,8 +913,8 @@ gal_blank_remove_rows(gal_data_t *columns, gal_list_sizet_t
*column_indexs)
/* If the given index is larger than the number of elements in the
input list, then print an error. */
if(tmp==NULL)
- error(EXIT_FAILURE, 0, "%s: input list has %zu elements, but the "
- "column %zu (counting from zero) has been requested",
+ error(EXIT_FAILURE, 0, "%s: input list has %zu elements, but "
+ "the column %zu (counting from zero) has been requested",
__func__, gal_list_data_number(columns), tcol->v);
/* Build the flag of blank elements for this column. */
diff --git a/lib/fit.c b/lib/fit.c
new file mode 100644
index 00000000..af889b6a
--- /dev/null
+++ b/lib/fit.c
@@ -0,0 +1,667 @@
+/*********************************************************************
+Functions for parametric fitting.
+This is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+ Mohammad Akhlaghi <mohammad@akhlaghi.org>
+Contributing author(s):
+Copyright (C) 2022 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 <stdio.h>
+#include <errno.h>
+#include <error.h>
+#include <string.h>
+#include <stdlib.h>
+
+#include <gsl/gsl_fit.h>
+#include <gsl/gsl_multifit.h>
+
+#include <gnuastro/fit.h>
+#include <gnuastro/blank.h>
+#include <gnuastro/pointer.h>
+
+#include <gnuastro-internal/checkset.h>
+
+
+
+
+
+/**********************************************************************/
+/**************** Identifiers ****************/
+/**********************************************************************/
+
+/* Read the desired parameters. */
+uint8_t
+gal_fit_name_to_id(char *name)
+{
+ if( !strcmp(name, "linear") )
+ return GAL_FIT_LINEAR;
+ else if( !strcmp(name, "linear-weighted") )
+ return GAL_FIT_LINEAR_WEIGHTED;
+ else if( !strcmp(name, "linear-no-constant") )
+ return GAL_FIT_LINEAR_NO_CONSTANT;
+ else if( !strcmp(name, "linear-no-constant-weighted") )
+ return GAL_FIT_LINEAR_NO_CONSTANT_WEIGHTED;
+ else if( !strcmp(name, "polynomial-weighted") )
+ return GAL_FIT_POLYNOMIAL_WEIGHTED;
+ else if( !strcmp(name, "polynomial") )
+ return GAL_FIT_POLYNOMIAL;
+ else if( !strcmp(name, "polynomial-robust") )
+ return GAL_FIT_POLYNOMIAL_ROBUST;
+ else return GAL_FIT_INVALID;
+
+ /* If control reaches here, there was a bug! */
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to "
+ "find a fix it. Control should not have reached here",
+ __func__, PACKAGE_BUGREPORT);
+ return GAL_FIT_INVALID;
+}
+
+
+
+
+
+char *
+gal_fit_name_from_id(uint8_t fitid)
+{
+ /* Prepare the temporary array. */
+ switch(fitid)
+ {
+ case GAL_FIT_LINEAR: return "linear";
+ case GAL_FIT_LINEAR_WEIGHTED: return "linear-weighted";
+ case GAL_FIT_LINEAR_NO_CONSTANT: return "linear-no-constant";
+ case GAL_FIT_POLYNOMIAL: return "polynomial";
+ case GAL_FIT_POLYNOMIAL_WEIGHTED: return "polynomial-weighted";
+ case GAL_FIT_POLYNOMIAL_ROBUST: return "polynomial-robust";
+ case GAL_FIT_LINEAR_NO_CONSTANT_WEIGHTED:
+ return "linear-no-constant-weighted";
+ default: return NULL;
+ }
+
+ /* If control reaches here, there was a bug! */
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to "
+ "find a fix it. Control should not have reached here",
+ __func__, PACKAGE_BUGREPORT);
+ return NULL;
+}
+
+
+
+
+
+int
+gal_fit_name_robust_to_id(char *name)
+{
+ /* In case 'name' is NULL, then return the invalid type. */
+ if(name==NULL) return GAL_FIT_ROBUST_INVALID;
+
+ /* Match the name. */
+ if( !strcmp(name, "bisquare") ) return GAL_FIT_ROBUST_BISQUARE;
+ else if( !strcmp(name, "cauchy") ) return GAL_FIT_ROBUST_CAUCHY;
+ else if( !strcmp(name, "fair") ) return GAL_FIT_ROBUST_FAIR;
+ else if( !strcmp(name, "huber") ) return GAL_FIT_ROBUST_HUBER;
+ else if( !strcmp(name, "ols") ) return GAL_FIT_ROBUST_OLS;
+ else if( !strcmp(name, "welsch") ) return GAL_FIT_ROBUST_WELSCH;
+ else return GAL_FIT_ROBUST_INVALID;
+
+ /* If control reaches here, there was a bug! */
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to "
+ "find a fix it. Control should not have reached this point",
+ __func__, PACKAGE_BUGREPORT);
+ return GAL_FIT_ROBUST_INVALID;
+}
+
+
+
+
+
+char *
+gal_fit_name_robust_from_id(uint8_t robustid)
+{
+ switch(robustid)
+ {
+ case GAL_FIT_ROBUST_BISQUARE: return "bisquare";
+ case GAL_FIT_ROBUST_CAUCHY: return "cauchy";
+ case GAL_FIT_ROBUST_FAIR: return "fair";
+ case GAL_FIT_ROBUST_HUBER: return "huber";
+ case GAL_FIT_ROBUST_OLS: return "ols";
+ case GAL_FIT_ROBUST_WELSCH: return "welsch";
+ default: return NULL;
+ }
+
+ /* If control reaches here, there was a bug! */
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to "
+ "find a fix it. Control should not have reached this point",
+ __func__, PACKAGE_BUGREPORT);
+ return NULL;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/**********************************************************************/
+/**************** Common to all ****************/
+/**********************************************************************/
+static gal_data_t *
+fit_1d_sanity_check(gal_data_t *in, gal_data_t *ref, const char *func)
+{
+ gal_data_t *out;
+
+ /* Make sure the input is 1-dimensional. */
+ if(in->ndim!=1)
+ error(EXIT_FAILURE, 0, "%s: inputs must have one dimension", func);
+
+ /* Make sure the input has the same size as the reference. */
+ if(in->size != ref->size)
+ error(EXIT_FAILURE, 0, "%s: all inputs must have the same size",
+ func);
+
+ /* Make sure output has a double type. */
+ out = ( in->type==GAL_TYPE_FLOAT64
+ ? in
+ : gal_data_copy_to_new_type(in, GAL_TYPE_FLOAT64) );
+
+ /* If there are blank values, print a warning, then return. */
+ if(gal_blank_present(out, 1))
+ error(EXIT_SUCCESS, 0, "%s: at least one of the input columns "
+ "have a blank value; the fit will become NaN. Within the "
+ "Gnuastro, you can use 'gal_blank_remove_rows' to remove "
+ "all rows that have at least one blank value in any column",
+ func);
+ return out;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/**********************************************************************/
+/**************** Linear fit ****************/
+/**********************************************************************/
+static gal_data_t *
+fit_1d_linear_base(gal_data_t *xin, gal_data_t *yin,
+ gal_data_t *ywht, int fitid)
+{
+ double *o, nparam=NAN;
+ size_t osize, chisqind=GAL_BLANK_SIZE_T;
+ gal_data_t *x=NULL, *y=NULL, *w=NULL, *out;
+
+ /* Basic sanity checks. */
+ x=fit_1d_sanity_check(xin, xin, __func__);
+ y=fit_1d_sanity_check(yin, xin, __func__);
+ if(ywht) w=fit_1d_sanity_check(ywht, xin, __func__);
+
+ /* Allocate the output dataset. */
+ osize = ( fitid==GAL_FIT_LINEAR || fitid==GAL_FIT_LINEAR_WEIGHTED
+ ? 6 : 3 );
+ out=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &osize, NULL, 0,
+ -1, 1, NULL, NULL, NULL);
+
+ /* For a check.
+ {
+ size_t i;
+ double *xa=x->array, *ya=y->array;
+ for(i=0;i<x->size;++i)
+ printf("%-15f %-15f\n", xa[i], ya[i]);
+ } */
+
+ /* Do the fitting. */
+ o=out->array;
+ switch(fitid)
+ {
+ case GAL_FIT_LINEAR:
+ nparam=2;
+ chisqind=5;
+ gsl_fit_linear(x->array, 1, y->array, 1, x->size, o, o+1,
+ o+2, o+3, o+4, o+5);
+ break;
+ case GAL_FIT_LINEAR_WEIGHTED:
+ nparam=2;
+ chisqind=5;
+ gsl_fit_wlinear(x->array, 1, w->array, 1, y->array, 1, x->size,
+ o, o+1, o+2, o+3, o+4, o+5);
+ break;
+ case GAL_FIT_LINEAR_NO_CONSTANT:
+ nparam=1;
+ chisqind=2;
+ gsl_fit_mul(x->array, 1, y->array, 1, x->size, o, o+1, o+2);
+ break;
+ case GAL_FIT_LINEAR_NO_CONSTANT_WEIGHTED:
+ nparam=1;
+ chisqind=2;
+ gsl_fit_wmul(x->array, 1, w->array, 1, y->array, 1, x->size,
+ o, o+1, o+2);
+ break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' "
+ "to fix the problem. The fitting id '%d' isn't recognized",
+ __func__, PACKAGE_BUGREPORT, fitid);
+ }
+
+ /* For a check.
+ {
+ printf("c0: %f\nc1: %f\n"
+ "cov00: %f\ncov01: %f\ncov11: %f\nsumsq: %f\n",
+ o[0], o[1], o[2], o[3], o[4], o[5]);
+ } */
+
+ /* Calculate the reduced chi^2: As mentioned in [1], in case we have the
+ chi^2, then it is simply the chi^2 divided by the degrees of
+ freedom. But GSL only returns the chi^2 for weighted fits. Therefore,
+ according to [1], we can also use the residual sum of squares instead.
+
+ The number of degrees of freedom is defined by the number of
+ observations subtracted from the number of fitted parameters.
+
+ [1] https://en.wikipedia.org/wiki/Reduced_chi-squared_statistic */
+ o[chisqind] /= (x->size - nparam);
+
+ /* Clean up and return. */
+ if(x!=xin) gal_data_free(x);
+ if(y!=yin) gal_data_free(y);
+ if(ywht && w!=ywht) gal_data_free(w);
+ return out;
+}
+
+
+
+
+
+gal_data_t *
+gal_fit_1d_linear(gal_data_t *xin, gal_data_t *yin, gal_data_t *ywht)
+{
+ return fit_1d_linear_base(xin, yin, ywht,
+ ( ywht
+ ? GAL_FIT_LINEAR_WEIGHTED
+ : GAL_FIT_LINEAR));
+}
+
+
+
+
+
+gal_data_t *
+gal_fit_1d_linear_no_constant(gal_data_t *xin, gal_data_t *yin,
+ gal_data_t *ywht)
+{
+ return fit_1d_linear_base(xin, yin, ywht,
+ ( ywht
+ ? GAL_FIT_LINEAR_NO_CONSTANT_WEIGHTED
+ : GAL_FIT_LINEAR_NO_CONSTANT) );
+}
+
+
+
+
+
+gal_data_t *
+fit_1d_estimate_prepare(gal_data_t *xin, gal_data_t *fit, gal_data_t **xd,
+ const char *func)
+{
+ gal_data_t *out=NULL;
+
+ /* The Fit arrays should be double precision. */
+ if(fit->type!=GAL_TYPE_FLOAT64
+ || (fit->next && fit->next->type!=GAL_TYPE_FLOAT64) )
+ error(EXIT_FAILURE, 0, "%s: the 'fit' argument should only "
+ "contain double precision floating point types", func);
+ if(fit->ndim!=1 || (fit->next && fit->next->ndim!=2) )
+ error(EXIT_FAILURE, 0, "%s: the 'fit' argument should only "
+ "contain single-dimensional outputs", func);
+ if(fit->next && (fit->next->dsize[0]!=fit->next->dsize[1]))
+ error(EXIT_FAILURE, 0, "%s: the secont dataset of the 'fit' "
+ "argument should be square (same size in both "
+ "dimensions)", func);
+
+ /* Make sure the input X values are in double precision. */
+ *xd = ( xin->type==GAL_TYPE_FLOAT64
+ ? xin
+ : gal_data_copy_to_new_type(xin, GAL_TYPE_FLOAT64) );
+
+ /* Allocate the output datasets. */
+ gal_list_data_add_alloc(&out, NULL, GAL_TYPE_FLOAT64, 1, xin->dsize,
+ NULL, 1, xin->minmapsize, xin->quietmmap,
+ "Y-ESTIMATED", xin->unit,
+ "Estimated value after fitting.");
+ gal_list_data_add_alloc(&out, NULL, GAL_TYPE_FLOAT64, 1, xin->dsize,
+ NULL, 1, xin->minmapsize, xin->quietmmap,
+ "Y-ESTIMATED-ERR", xin->unit,
+ "Estimated error on value after fitting.");
+ gal_list_data_reverse(&out);
+
+ /* Return the output. */
+ return out;
+}
+
+
+
+
+
+gal_data_t *
+gal_fit_1d_linear_estimate(gal_data_t *fit, gal_data_t *xin)
+{
+ size_t i;
+ gal_data_t *out=NULL, *xd;
+ double *x, *y, *yerr, *f=fit->array;
+
+ /* Do the basic preparations. */
+ out=fit_1d_estimate_prepare(xin, fit, &xd, __func__);
+
+ /* Set the pointers. */
+ x = xd->array;
+ y = out->array;
+ yerr = out->next->array;
+
+ /* Estimate the values. */
+ switch(fit->size)
+ {
+ case 6: /* Linear with constant. */
+ for(i=0;i<out->size;++i)
+ gsl_fit_linear_est(x[i], f[0], f[1], f[2], f[3], f[4],
+ y+i, yerr+i);
+ break;
+
+ case 3: /* Linear WITHOUT constant. */
+ for(i=0;i<out->size;++i)
+ gsl_fit_mul_est(x[i], f[0], f[1], y+i, yerr+i);
+ break;
+
+ default: /* Un-recognized situation! */
+ error(EXIT_FAILURE, 0, "%s: the 'fit' argument should "
+ "either have 6 or 3 elements (be an output of "
+ "'gal_fit_1d_linear' or 'gal_fit_1d_linear_no_constant'"
+ "respectively), but it has %zu elements", __func__,
+ fit->size);
+ }
+
+ /* Clean up. */
+ if(xd!=xin) gal_data_free(xd);
+ return out;
+}
+
+
+
+
+
+static void
+fit_1d_polynomial_prepare(gal_data_t *xin, gal_data_t *yin,
+ gal_data_t *ywht, int nconst,
+ gsl_matrix **x, gsl_vector **c,
+ gsl_matrix **cov, gsl_vector *y,
+ gsl_vector *w)
+{
+ size_t i, j;
+ double *xo, *xi;
+
+ /* Use GSL's own matrix allocation functions for the structures that need
+ allocation and we can't use the same allocated space of the inputs. */
+ *c = gsl_vector_alloc(nconst);
+ *cov = gsl_matrix_alloc(nconst, nconst);
+ *x = gsl_matrix_alloc(xin->size, nconst);
+
+ /* Fill in the X matrix. */
+ xi=xin->array;
+ xo=(*x)->data;
+ for(i=0;i<xin->size;++i)
+ {
+ /* The first column (constant) doesn't depend on X. So we'll give it
+ a value of 1.0. */
+ xo[ i*nconst ] = 1.0f;
+
+ /* Column i is the multiplication of column i-1 with the input
+ horizontal value. This will make it a polynomial. */
+ for(j=1;j<nconst;++j)
+ xo[ i*nconst + j ] = xo[ i*nconst + j-1 ] * xi[i];
+ }
+
+ /* For a check.
+ {
+ size_t checki=5;
+ printf("Row %zu: ", checki);
+ for(j=0;j<maxpower;++j)
+ printf("%.3f ", xo[ checki*maxpower + j ]);
+ printf("\n");
+ exit(0);
+ } */
+
+ /* Set the pointers of the 'y' and 'w' GSL vectors. */
+ y->data=yin->array;
+ if(ywht) w->data=ywht->array;
+}
+
+
+
+
+
+gal_data_t *
+gal_fit_1d_polynomial_base(gal_data_t *xin, gal_data_t *yin,
+ gal_data_t *ywht, size_t maxpower,
+ uint8_t robustid, double *redchisq)
+{
+ /* Low-level variable. */
+ size_t nconst=maxpower+1;
+
+ /* Other variables */
+ gsl_vector *c=NULL;
+ double chisq=NAN, sse=NAN;
+ gsl_matrix *x=NULL, *cov=NULL;
+ size_t covsize[2]={nconst, nconst};
+ gsl_multifit_linear_workspace *work_n;
+ gsl_multifit_robust_workspace *work_r;
+ const gsl_multifit_robust_type *rtype=NULL;
+ gal_data_t *xdata, *ydata, *wdata, *tmp, *out=NULL;
+
+ /* For the 'y' and 'w' GSL vectors, we don't actually need to allocate
+ any space, we can just use the allocated space within the
+ 'gal_data_t'. We can't set the pointers now because we aren't sure
+ they have 'double' type yet. */
+ gsl_vector yvec={yin->size, 1, NULL, NULL, 0}; /* Both have same size, */
+ gsl_vector wvec={yin->size, 1, NULL, NULL, 0}; /* but ywht may be NULL!*/
+ gsl_vector *y=&yvec, *w=&wvec; /* These have to be after the two above.*/
+
+ /* Basic sanity checks. */
+ xdata = fit_1d_sanity_check(xin, xin, __func__);
+ ydata = fit_1d_sanity_check(yin, xin, __func__);
+ wdata = ywht ? fit_1d_sanity_check(ywht, xin, __func__) : NULL;
+
+ /* Fill all the GSL structures. */
+ fit_1d_polynomial_prepare(xdata, ydata, wdata, nconst,
+ &x, &c, &cov, y, w);
+
+ /* Do the fit (depending on if it is robust or not. */
+ if(robustid==GAL_FIT_ROBUST_INVALID)
+ {
+ work_n = gsl_multifit_linear_alloc(xin->size, nconst);
+ if(ywht) gsl_multifit_wlinear(x, w, y, c, cov, &chisq, work_n);
+ else gsl_multifit_linear( x, y, c, cov, &sse, work_n);
+ gsl_multifit_linear_free(work_n);
+ }
+ else
+ {
+ /* Select the robust function type. */
+ switch(robustid)
+ {
+ case GAL_FIT_ROBUST_BISQUARE: rtype=gsl_multifit_robust_bisquare;
+ break;
+ case GAL_FIT_ROBUST_CAUCHY: rtype=gsl_multifit_robust_cauchy;
+ break;
+ case GAL_FIT_ROBUST_FAIR: rtype=gsl_multifit_robust_fair;
+ break;
+ case GAL_FIT_ROBUST_HUBER: rtype=gsl_multifit_robust_huber;
+ break;
+ case GAL_FIT_ROBUST_OLS: rtype=gsl_multifit_robust_ols;
+ break;
+ case GAL_FIT_ROBUST_WELSCH: rtype=gsl_multifit_robust_welsch;
+ break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at "
+ "'%s' to fix the problem. the 'robustid' value '%d' "
+ "isn't recognize", __func__, PACKAGE_BUGREPORT, robustid);
+ }
+
+ /* Initialize the worker and do the fit (depending on if a weight
+ image was provided). */
+ work_r=gsl_multifit_robust_alloc(rtype, x->size1, x->size2);
+ gsl_multifit_robust(x, y, c, cov, work_r);
+
+ /* Get the residual sum of squares and free the worker. */
+ sse=gsl_multifit_robust_statistics(work_r).sse;
+ gsl_multifit_robust_free(work_r);
+ }
+
+ /* For a check:
+ {
+ size_t i;
+ double *ca=c->data;
+ for(i=0;i<=nconst;++i) { printf("%f ", ca[i]); } printf("\n");
+ } */
+
+ /* Allocate the output dataset containing the fit results as first
+ 'gal_data_t'. */
+ tmp=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &nconst, NULL, 0,
+ xin->minmapsize, xin->quietmmap, NULL, NULL,
+ NULL);
+ memcpy(tmp->array, c->data, nconst*sizeof c->data);
+ out=tmp;
+
+ /* Allocate the second element of the output (the covariance matrix). */
+ tmp=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 2, covsize, NULL, 0,
+ xin->minmapsize, xin->quietmmap, NULL, NULL,
+ NULL);
+ memcpy(tmp->array, cov->data, nconst*nconst*sizeof cov->data);
+ out->next=tmp;
+
+ /* Calculate the reduced chi^2, see the description of same step in
+ 'fit_1d_linear_base'. */
+ *redchisq = (isnan(chisq) ? sse : chisq) / (xdata->size-nconst);
+
+ /* Clean up and return. */
+ gsl_matrix_free(x);
+ gsl_vector_free(c);
+ gsl_matrix_free(cov);
+ if(xdata!=xin) gal_data_free(xdata);
+ if(ydata!=yin) gal_data_free(ydata);
+ if(ywht && wdata!=ywht) gal_data_free(wdata);
+ return out;
+}
+
+
+
+
+
+gal_data_t *
+gal_fit_1d_polynomial(gal_data_t *xin, gal_data_t *yin,
+ gal_data_t *ywht, size_t maxpower,
+ double *redchisq)
+{
+ return gal_fit_1d_polynomial_base(xin, yin, ywht, maxpower,
+ GAL_FIT_ROBUST_INVALID,
+ redchisq);
+}
+
+
+
+
+
+gal_data_t *
+gal_fit_1d_polynomial_robust(gal_data_t *xin, gal_data_t *yin,
+ size_t maxpower, uint8_t robustid,
+ double *redchisq)
+{
+ /* Robust fitting doesn't use weights (the functions are effectively the
+ weight). */
+ return gal_fit_1d_polynomial_base(xin, yin, NULL, maxpower,
+ robustid, redchisq);
+}
+
+
+
+
+
+/* Estimate values from a polynomial fit.*/
+gal_data_t *
+gal_fit_1d_polynomial_estimate(gal_data_t *fit, gal_data_t *xin)
+{
+ size_t i, j;
+ size_t nconst=fit->size;
+ gal_data_t *xd, *out=NULL;
+ double *y, *xi, *xo, *yerr;
+
+ /* We don't need to allocate space for the GSL vectors and matrices, we
+ can just use the allocated space within the 'gal_data_t'. We can't set
+ the pointers now because we aren't sure they have 'double' type
+ yet. */
+ gsl_vector xvec={nconst, 1, NULL, NULL, 0};
+ gsl_vector cvec={nconst, 1, NULL, NULL, 0};
+ gsl_matrix cmat={nconst, nconst, nconst, NULL, NULL, 0};
+
+ /* Do the basic preparations. */
+ out=fit_1d_estimate_prepare(xin, fit, &xd, __func__);
+
+ /* Set the pointers. */
+ xo = xvec.data = gal_pointer_allocate(GAL_TYPE_FLOAT64, nconst,
+ 0, __func__, "xvec.data");
+ xi = xd->array;
+ cvec.data = fit->array;
+ y = out->array;
+ yerr = out->next->array;
+ cmat.data = fit->next->array;
+
+ /* Do the estimation. */
+ for(i=0;i<xd->size;++i)
+ {
+ xo[0]=1.0f; for(j=1;j<nconst;++j) xo[j] = xo[j-1] * xi[i];
+ gsl_multifit_linear_est(&xvec, &cvec, &cmat, y+i, yerr+i);
+ }
+
+ /* Clean up and return. */
+ if(xd!=xin) gal_data_free(xd);
+ free(xvec.data);
+ return out;
+}
diff --git a/lib/fits.c b/lib/fits.c
index b0a0e4a1..3901eb0b 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -2086,7 +2086,8 @@ gal_fits_key_write_in_ptr(gal_fits_list_key_t **keylist,
fitsfile *fptr)
gal_fits_key_write_in_ptr_nan_check(tmp);
/* Write/Update the keyword value. */
- if( fits_update_key(fptr, gal_fits_type_to_datatype(tmp->type),
+ if( fits_update_key(fptr,
+ gal_fits_type_to_datatype(tmp->type),
tmp->keyname, tmp->value, tmp->comment,
&status) )
gal_fits_io_error(status, NULL);
@@ -2100,14 +2101,15 @@ gal_fits_key_write_in_ptr(gal_fits_list_key_t
**keylist, fitsfile *fptr)
/* Write the units if it was given. */
if( tmp->unit
- && fits_write_key_unit(fptr, tmp->keyname, tmp->unit, &status) )
+ && fits_write_key_unit(fptr, tmp->keyname, tmp->unit,
+ &status) )
gal_fits_io_error(status, NULL);
/* Free the value pointer if desired: */
- if(tmp->kfree) free(tmp->keyname);
+ if(tmp->ufree) free(tmp->unit);
if(tmp->vfree) free(tmp->value);
+ if(tmp->kfree) free(tmp->keyname);
if(tmp->cfree) free(tmp->comment);
- if(tmp->ufree) free(tmp->unit);
}
/* Keep the pointer to the next keyword and free the allocated
diff --git a/lib/gnuastro/fit.h b/lib/gnuastro/fit.h
new file mode 100644
index 00000000..eea4fa96
--- /dev/null
+++ b/lib/gnuastro/fit.h
@@ -0,0 +1,126 @@
+/*********************************************************************
+fit -- functions for doing fitting to input datasets.
+This is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+ Mohammad Akhlaghi <mohammad@akhlaghi.org>
+Contributing author(s):
+Copyright (C) 2022 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 __GAL_FIT_H__
+#define __GAL_FIT_H__
+
+
+/* Include other headers if necessary here. Note that other header files
+ must be included before the C++ preparations below */
+#include <gnuastro/data.h>
+
+
+
+/* C++ Preparations */
+#undef __BEGIN_C_DECLS
+#undef __END_C_DECLS
+#ifdef __cplusplus
+# define __BEGIN_C_DECLS extern "C" {
+# define __END_C_DECLS }
+#else
+# define __BEGIN_C_DECLS /* empty */
+# define __END_C_DECLS /* empty */
+#endif
+/* End of C++ preparations */
+
+
+
+
+
+/* Actual header contants (the above were for the Pre-processor). */
+__BEGIN_C_DECLS /* From C++ preparations */
+
+
+
+
+/* Definitions. */
+enum gal_fit_types
+{
+ GAL_FIT_INVALID, /* Invalid (=0 by C standard). */
+ GAL_FIT_LINEAR,
+ GAL_FIT_LINEAR_WEIGHTED,
+ GAL_FIT_LINEAR_NO_CONSTANT,
+ GAL_FIT_LINEAR_NO_CONSTANT_WEIGHTED,
+ GAL_FIT_POLYNOMIAL,
+ GAL_FIT_POLYNOMIAL_ROBUST,
+ GAL_FIT_POLYNOMIAL_WEIGHTED,
+
+ /* This will be the total number of shapes (good for scripts). */
+ GAL_FIT_NUMBER
+};
+
+enum gal_fit_robust_types
+{
+ GAL_FIT_ROBUST_INVALID, /* Invalid (=0 by C standard). */
+ GAL_FIT_ROBUST_BISQUARE,
+ GAL_FIT_ROBUST_CAUCHY,
+ GAL_FIT_ROBUST_FAIR,
+ GAL_FIT_ROBUST_HUBER,
+ GAL_FIT_ROBUST_OLS,
+ GAL_FIT_ROBUST_WELSCH,
+
+ /* This will be the total number of shapes (good for scripts). */
+ GAL_FIT_ROBUST_NUMBER
+};
+
+
+
+/* Functions */
+uint8_t
+gal_fit_name_to_id(char *name);
+
+char *
+gal_fit_name_from_id(uint8_t fitid);
+
+int
+gal_fit_name_robust_to_id(char *name);
+
+char *
+gal_fit_name_robust_from_id(uint8_t robustid);
+
+gal_data_t *
+gal_fit_1d_linear(gal_data_t *xin, gal_data_t *yin,
+ gal_data_t *ywht);
+
+gal_data_t *
+gal_fit_1d_linear_no_constant(gal_data_t *xin, gal_data_t *yin,
+ gal_data_t *ywht);
+
+gal_data_t *
+gal_fit_1d_linear_estimate(gal_data_t *fit, gal_data_t *xin);
+
+gal_data_t *
+gal_fit_1d_polynomial(gal_data_t *xin, gal_data_t *yin,
+ gal_data_t *ywht, size_t maxpower,
+ double *redchisq);
+
+gal_data_t *
+gal_fit_1d_polynomial_robust(gal_data_t *xin, gal_data_t *yin,
+ size_t maxpower, uint8_t robustid,
+ double *redchisq);
+
+gal_data_t *
+gal_fit_1d_polynomial_estimate(gal_data_t *fit, gal_data_t *xin);
+
+__END_C_DECLS /* From C++ preparations */
+
+#endif /* __GAL_FIT_H__ */
diff --git a/lib/gnuastro/jpeg.h b/lib/gnuastro/jpeg.h
index 8f42bcd9..3bbe02f4 100644
--- a/lib/gnuastro/jpeg.h
+++ b/lib/gnuastro/jpeg.h
@@ -72,4 +72,4 @@ gal_jpeg_write(gal_data_t *in, char *filename, uint8_t
quality,
__END_C_DECLS /* From C++ preparations */
-#endif /* __GAL_TIFF_H__ */
+#endif /* __GAL_JPEG_H__ */
diff --git a/lib/options.c b/lib/options.c
index da70c522..8d725e71 100644
--- a/lib/options.c
+++ b/lib/options.c
@@ -1034,7 +1034,8 @@ gal_options_parse_list_of_strings(char *string, char
*filename, size_t lineno)
'gal_data_t' contains the array of given strings. You can read the
number of inputs from its 'size' element. */
gal_data_t *
-gal_options_parse_csv_strings_raw(char *string, char *filename, size_t lineno)
+gal_options_parse_csv_strings_raw(char *string, char *filename,
+ size_t lineno)
{
size_t i, num;
gal_data_t *out;
@@ -1156,10 +1157,10 @@ gal_options_parse_csv_strings(struct argp_option
*option, char *arg,
for(i=0;i<values->size;++i)
{
if( nc > GAL_OPTIONS_STATIC_MEM_FOR_VALUES-100 )
- error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s so we "
- "can address the problem. The number of necessary "
- "characters in the statically allocated string has become "
- "too close to %d", __func__, PACKAGE_BUGREPORT,
+ error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s "
+ "so we can address the problem. The number of necessary "
+ "characters in the statically allocated string has "
+ "become too close to %d", __func__, PACKAGE_BUGREPORT,
GAL_OPTIONS_STATIC_MEM_FOR_VALUES);
nc += sprintf(sstr+nc, "%s,", strarr[i]);
}
@@ -2005,8 +2006,9 @@ options_sanity_check(struct argp_option *option, char
*arg,
static void
-gal_options_read_check(struct argp_option *option, char *arg, char *filename,
- size_t lineno, struct gal_options_common_params *cp)
+gal_options_read_check(struct argp_option *option, char *arg,
+ char *filename, size_t lineno,
+ struct gal_options_common_params *cp)
{
void *topass;
@@ -2059,8 +2061,9 @@ gal_options_read_check(struct argp_option *option, char
*arg, char *filename,
/* Read the string argument into the value. */
if( gal_type_from_string(&option->value, arg, option->type) )
- /* Fortunately 'error_at_line' will behave like 'error' when the
- filename is NULL (the option was read from a command-line). */
+ /* Fortunately 'error_at_line' will behave like 'error' when
+ the filename is NULL (the option was read from a
+ command-line). */
error_at_line(EXIT_FAILURE, 0, filename, lineno,
"'%s' (value to option '--%s') couldn't be read "
"into the proper numerical type. Common causes "
diff --git a/lib/txt.c b/lib/txt.c
index 4922d56d..92329ac4 100644
--- a/lib/txt.c
+++ b/lib/txt.c
@@ -1692,9 +1692,9 @@ gal_txt_write(gal_data_t *input, struct
gal_fits_list_key_t **keylist,
/* Check if the dimensionality and size is the same for all the
elements. */
if( input!=data && gal_dimension_is_different(input, data) )
- error(EXIT_FAILURE, 0, "%s: the input list of datasets must have the "
- "same sizes (dimensions and length along each dimension)",
- __func__);
+ error(EXIT_FAILURE, 0, "%s: the input list of datasets must "
+ "have the same sizes (dimensions and length along each "
+ "dimension)", __func__);
}
@@ -1709,9 +1709,10 @@ gal_txt_write(gal_data_t *input, struct
gal_fits_list_key_t **keylist,
{
/* Make sure the file doesn't already exist. */
if( gal_checkset_check_file_return(filename) )
- error(EXIT_FAILURE, 0, "%s: %s already exists. For safety, this "
- "function will not over-write an existing file. Please delete "
- "it before calling this function", __func__, filename);
+ error(EXIT_FAILURE, 0, "%s: %s already exists. For safety, "
+ "this function will not over-write an existing file. "
+ "Please delete it before calling this function",
+ __func__, filename);
/* Open the output file. */
errno=0;
@@ -1755,7 +1756,7 @@ gal_txt_write(gal_data_t *input, struct
gal_fits_list_key_t **keylist,
j=0;
for(data=input;data!=NULL;data=data->next) /* Column. */
txt_print_value(fp, data->array, data->type, i,
- fmts[j++ * FMTS_COLS]);
+ fmts[j++ * FMTS_COLS]);
fprintf(fp, "\n");
}
break;
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 4d9b902a..33d12a7d 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -181,12 +181,15 @@ if COND_SEGMENT
segment/segment-3d.sh: noisechisel/noisechisel-3d.sh.log
endif
if COND_STATISTICS
- MAYBE_STATISTICS_TESTS = statistics/basicstats.sh statistics/estimate_sky.sh
\
- statistics/from-stdin.sh
+ MAYBE_STATISTICS_TESTS = statistics/basicstats.sh \
+ statistics/from-stdin.sh \
+ statistics/estimate_sky.sh \
+ statistics/fitting-polynomial-robust.sh
+ statistics/from-stdin.sh: prepconf.sh.log
statistics/basicstats.sh: mknoise/addnoise.sh.log
statistics/estimate_sky.sh: mknoise/addnoise.sh.log
- statistics/from-stdin.sh: prepconf.sh.log
+ statistics/fitting-polynomial-robust.sh: prepconf.sh.log
endif
if COND_TABLE
MAYBE_TABLE_TESTS = table/txt-to-fits-binary.sh \
diff --git a/tests/statistics/fitting-data.txt
b/tests/statistics/fitting-data.txt
new file mode 100644
index 00000000..f9039553
--- /dev/null
+++ b/tests/statistics/fitting-data.txt
@@ -0,0 +1,63 @@
+# Table to test the fitting features of Statistics
+#
+# This table was generated with the following commands using Gnuastro 0.19.
+#
+# $ export GSL_RNG_SEED=1664015492
+# $ seq 0.1 0.05 2 \
+# | asttable -ofitting-data.txt --envseed -c1 \
+# -c'arith 1.23 -4.56 $1 x + 7.89 $1 x $1 x + set-y \
+# 0.1 y x set-yerr \
+# y yerr mknoise-sigma yerr' \
+# --colmetadata=1,X \
+# --colmetadata=2,Y \
+# --colmetadata=3,Yerr
+#
+# Copyright (C) 2022 Free Software Foundation, Inc.
+#
+# Copying and distribution of this file, with or without modification, are
+# permitted in any medium without royalty provided the copyright notice and
+# this notice are preserved. This file is offered as-is, without any
+# warranty.
+#
+# Column 1: X [ ,f64,]
+# Column 2: Y [arith_unit_7,f64,] Column from arithmetic operation 7
+# Column 3: Yerr [ ,f64,]
+0.1 0.89231850942267 0.085290003616959
+0.15 0.5266417844425 0.072352503543384
+0.2 0.74688387286754 0.063360003461838
+0.25 0.68313579973585 0.058312503372319
+0.3 0.56406673731523 0.057210003274828
+0.35 0.72296542514828 0.060052503169365
+0.4 0.63978968990778 0.06684000305593
+0.45 0.68345010851872 0.077572502934523
+0.5 0.85462157058084 0.092250002805144
+0.55 1.0907144431374 0.11087250266779
+0.6 1.4194495635087 0.13344000252247
+0.65 1.5886836818975 0.15995250236917
+0.7 1.5345607654082 0.19041000220791
+0.75 2.5068820594536 0.22481250203867
+0.8 2.6457411697165 0.26316000186145
+0.85 3.6429002647991 0.30545250167627
+0.9 3.4988944766376 0.35169000148311
+0.95 3.9566235305386 0.40187250128198
+1 4.7519722911804 0.45600000107288
+1.05 5.2957101184984 0.51407250085581
+1.1 6.3431257005222 0.57609000063077
+1.15 5.8520058385084 0.64205250039775
+1.2 6.612813901809 0.71196000015676
+1.25 7.8955679134514 0.7858124999078
+1.3 8.8245735526149 0.86360999965087
+1.35 7.7953468595072 0.94535249938596
+1.4 9.8088627943 1.0310399991131
+1.45 11.561043909455 1.1206724988322
+1.5 13.267813490096 1.2142499985434
+1.55 14.12894084285 1.3117724982466
+1.6 14.259646409046 1.4132399979419
+1.65 15.735055395656 1.5186524976291
+1.7 15.535047056382 1.6280099973084
+1.75 19.426912050868 1.7413124969797
+1.8 17.000366817577 1.8585599966431
+1.85 18.847181551801 1.9797524962984
+1.9 18.270437654511 2.1048899959458
+1.95 24.684339280065 2.2339724955853
+2 24.522383853646 2.3669999952167
diff --git a/tests/statistics/fitting-polynomial-robust.sh
b/tests/statistics/fitting-polynomial-robust.sh
new file mode 100755
index 00000000..df795c10
--- /dev/null
+++ b/tests/statistics/fitting-polynomial-robust.sh
@@ -0,0 +1,59 @@
+# Do a robust polynomial fit to the given data
+#
+# See the Tests subsection of the manual for a complete explanation
+# (in the Installing gnuastro section).
+#
+# Original author:
+# Mohammad Akhlaghi <mohammad@akhlaghi.org>
+# Contributing author(s):
+# Copyright (C) 2022 Free Software Foundation, Inc.
+#
+# Copying and distribution of this file, with or without modification,
+# are permitted in any medium without royalty provided the copyright
+# notice and this notice are preserved. This file is offered as-is,
+# without any warranty.
+
+
+
+
+
+# Preliminaries
+# =============
+#
+# Set the variables (The executable is in the build tree). Do the
+# basic checks to see if the executable is made or if the defaults
+# file exists (basicchecks.sh is in the source tree).
+prog=statistics
+execname=../bin/$prog/ast$prog
+table=$topsrc/tests/statistics/fitting-data.txt
+
+
+
+
+
+# Skip?
+# =====
+#
+# If the dependencies of the test don't exist, then skip it. There are two
+# types of dependencies:
+#
+# - The executable was not made (for example due to a configure option),
+#
+# - The input data was not made (for example the test that created the
+# data file failed).
+if [ ! -f $execname ]; then echo "$execname not created."; exit 77; fi
+if [ ! -f $table ]; then echo "$table does not exist."; exit 77; fi
+
+
+
+
+
+# Actual test script
+# ==================
+#
+# 'check_with_program' can be something like Valgrind or an empty
+# string. Such programs will execute the command if present and help in
+# debugging when the developer doesn't have access to the user's system.
+$check_with_program $execname $table -cX,Y --fit=polynomial-robust \
+ --fitmaxpower=2 --fitestimate=self \
+ --output=fitting-estimate.txt
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 97cb586b: Library (fit.h): least squares fitting, accessible via Statistics,
Mohammad Akhlaghi <=