[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 7bb809c2: Library (arithmetic.h): loading exte
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 7bb809c2: Library (arithmetic.h): loading external columns and random from hist |
|
Date: |
Mon, 27 Jun 2022 08:17:39 -0400 (EDT) |
branch: master
commit 7bb809c2ac872af3aee31106be4946654f7ebc3f
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Library (arithmetic.h): loading external columns and random from hist
Until now, Gnuastro would only use the Gaussian, Poisson and Uniform
distributions of the GNU Scientific Library (GSL) for random number
generation. However, in some scenarios it is necessary to extract random
values from a custom distribution and there was no way to do it easily.
One of the hinderances was this: it was not possible to easily import
external columns into the arithmetic commands of the Arithmetic or Table
programs. In Arithmetic, all input files should be images and in Table,
they columns have to be in one of the given files, and have to have the
same number of rows. However, inserting custom column was necessary to load
custom-length columns from a table that are only used by the operator and
don't depend on the size of the main input to be operated on (to define the
custom distribution in the example above).
With this commit, a new 'load-col-' operator has been defined in the
Arithmetic library (that is available to both the Arithmetic and Table
programs). This operator can be used to load any column into the processing
stack, and is used by the random histogram distribution to load the
histogram.
Having that operator in place, it was possible to define two new arithmetic
operators in Gnuastro: 'random-from-hist-raw' and 'random-from-hist'. The
first one simply returns a random selection of bin centers that have a
distribution similar to the given histogram. The second uses a uniform
distribution within each bin to random select values within each bin (so
they aren't descrete any more, but still flat if the bins are wide).
Adding the 'random-from-hist' operator was suggested by Raul Infante-Sainz.
---
NEWS | 22 ++-
bin/arithmetic/arithmetic.c | 23 ++-
bin/arithmetic/ui.c | 47 +++---
bin/table/arithmetic.c | 34 +++-
bin/table/main.h | 1 +
doc/gnuastro.texi | 283 ++++++++++++++++++++++++++++++++--
lib/arithmetic.c | 368 ++++++++++++++++++++++++++++++++++++++++----
lib/gnuastro/arithmetic.h | 15 ++
8 files changed, 709 insertions(+), 84 deletions(-)
diff --git a/NEWS b/NEWS
index 97437a58..98e4e3f7 100644
--- a/NEWS
+++ b/NEWS
@@ -17,7 +17,7 @@ See the end of the file for license conditions.
Eskandarlou.
Arithmetic:
- - New operators:
+ - New operators (also available in Table's column arithmetic):
- stitch: connect any number of input images along the given dimension
(for example to stitch the images of separate amplifiers of a CCD
into one image).
@@ -27,6 +27,21 @@ See the end of the file for license conditions.
Zeropoint (the inverse operator 'counts-to-jy' already existed).
- mag-to-jy: convert AB magnitudes to Janskys.
- jy-to-mag: convert Janskys to AB magnitudes.
+ - load-col-%-from-%-hdu-%: Load a certain column into the operands
+ stack. The name or number of the column is the first '%'. The file
+ name hosting that column is the second '%', and HDU (if the file is
+ a FITS file), is identified by third '%'. This can be used to load
+ columns outside the main table in Table's column arithmetic or by
+ operators that need columns as input parameters (like the new
+ 'random_from_hist' operator).
+ - random-from-hist-raw: return a set of random points based on a
+ custom distribution defined by its histogram. See the description
+ and complete examples in the book for more.
+ - random-from-hist: Similar to 'random-from-hist-raw', but within each
+ bin, the final returned value is itself the result of a uniform
+ sampling (to make the final distribution less descrete). See the
+ example in the book for selecting random star magnitudes based on
+ Gaia's observed magnitude distribution.
Crop:
--oneelemstdout: when a crop has a single pixel and this option is
@@ -77,6 +92,11 @@ See the end of the file for license conditions.
- GAL_ARITHMETIC_OP_UNIQUE: to find unique elements in dataset.
- GAL_ARITHMETIC_OP_NOBLANK: to remove blank elements from dataset.
- GAL_ARITHMETIC_OP_STITCH: to stitch multiple datasets.
+ - GAL_ARITHMETIC_OP_RANDOM_FROM_HIST: the new 'random-from-hist' operator.
+ - GAL_ARITHMETIC_OP_RANDOM_FROM_HIST_RAW: the new 'random-from-hist-raw' op.
+ - gal_arithmetic_load_col: Low-level function for parsing the string of
+ the newly added 'load-col' operator, mentioned in the Arithmetic
+ section above.
** Removed features
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index b4003836..3f306615 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -1339,10 +1339,10 @@ arithmetic_operator_run(struct arithmeticparams *p, int
operator,
break;
default:
- error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
- "the problem. '%zu' is not recognized as an operand "
- "counter (with '%s')", __func__, PACKAGE_BUGREPORT,
- num_operands, operator_string);
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s "
+ "to fix the problem. '%zu' is not recognized as an "
+ "operand counter (with '%s')", __func__,
+ PACKAGE_BUGREPORT, num_operands, operator_string);
}
/* Run the arithmetic operation. Note that 'gal_arithmetic'
@@ -1455,11 +1455,12 @@ arithmetic_set_name_used_later(void *in, char *name)
void
reversepolish(struct arithmeticparams *p)
{
- gal_data_t *data;
size_t num_operands=0;
gal_list_str_t *token;
+ gal_data_t *data, *col;
char *hdu, *filename, *printnum;
int operator=GAL_ARITHMETIC_OP_INVALID;
+ struct gal_options_common_params *cp=&p->cp;
/* Prepare the processing: */
p->popcounter=0;
@@ -1488,6 +1489,11 @@ reversepolish(struct arithmeticparams *p)
else if( !strncmp(token->v, GAL_ARITHMETIC_SET_PREFIX,
GAL_ARITHMETIC_SET_PREFIX_LENGTH) )
gal_arithmetic_set_name(&p->setprm, token->v);
+ else if( (col=gal_arithmetic_load_col(token->v, cp->searchin,
+ cp->ignorecase,
+ cp->minmapsize,
+ cp->quietmmap))!=NULL )
+ operands_add(p, NULL, col);
else if( gal_array_file_recognized(token->v)
|| gal_arithmetic_set_is_name(p->setprm.named, token->v) )
operands_add(p, token->v, NULL);
@@ -1500,6 +1506,7 @@ reversepolish(struct arithmeticparams *p)
data->minmapsize=p->cp.minmapsize;
operands_add(p, NULL, data);
}
+
/* Last option is an operator: the program will abort if the token
isn't an operator. */
else
@@ -1543,8 +1550,10 @@ reversepolish(struct arithmeticparams *p)
p->cp.minmapsize,
p->cp.quietmmap);
data=p->operands->data;
- data->ndim=gal_dimension_remove_extra(data->ndim, data->dsize, NULL);
- if(!p->cp.quiet) printf(" - %s (hdu %s) is read.\n", filename, hdu);
+ data->ndim=gal_dimension_remove_extra(data->ndim, data->dsize,
+ NULL);
+ if(!p->cp.quiet)
+ printf(" - %s (hdu %s) is read.\n", filename, hdu);
}
else
error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s to fix "
diff --git a/bin/arithmetic/ui.c b/bin/arithmetic/ui.c
index 2aab601f..cd453887 100644
--- a/bin/arithmetic/ui.c
+++ b/bin/arithmetic/ui.c
@@ -34,6 +34,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/table.h>
#include <gnuastro/array.h>
#include <gnuastro/threads.h>
+#include <gnuastro/arithmetic.h>
#include <gnuastro-internal/timing.h>
#include <gnuastro-internal/options.h>
@@ -142,8 +143,6 @@ ui_initialize_options(struct arithmeticparams *p,
break;
case GAL_OPTIONS_KEY_TYPE:
- case GAL_OPTIONS_KEY_SEARCHIN:
- case GAL_OPTIONS_KEY_IGNORECASE:
case GAL_OPTIONS_KEY_STDINTIMEOUT:
cp->coptions[i].flags=OPTION_HIDDEN;
break;
@@ -290,11 +289,20 @@ ui_check_options_and_arguments(struct arithmeticparams *p)
p->outnamerequested = cp->output ? 1 : 0;
for(token=p->tokens; token!=NULL; token=token->next)
{
- /* Strings given to the 'tofile' operator are also considered as
- outputs and we should delete them before starting the parse. */
+ /* Strings given to the 'tofile' or 'tofilefree' operators are also
+ considered as outputs and we should delete them before starting
+ the parse. */
if( strncmp(OPERATOR_PREFIX_TOFILE, token->v,
- OPERATOR_PREFIX_LENGTH_TOFILE) )
+ OPERATOR_PREFIX_LENGTH_TOFILE)==0
+ || strncmp(OPERATOR_PREFIX_TOFILEFREE, token->v,
+ OPERATOR_PREFIX_LENGTH_TOFILEFREE)==0 )
+ {
+ filename=&token->v[ OPERATOR_PREFIX_LENGTH_TOFILE ];
+ gal_checkset_writable_remove(filename, cp->keep, cp->dontdelete);
+ }
+ /* This may be a simple filename. */
+ else
{
/* This token is a file, count how many mult-extension files we
have and use the first to set the output filename (if it has
@@ -303,10 +311,13 @@ ui_check_options_and_arguments(struct arithmeticparams *p)
|| gal_fits_file_recognized(token->v) )
{
/* Increment the counter for FITS files (if they are
- input). Recall that the 'tofile' operator can also have
- '.fits' suffixes (they are the names of the output
- files). */
- if( gal_array_name_recognized_multiext(token->v) )
+ input). Recall that the 'load-col' operator can also have
+ '.fits' suffixes, but they aren't relevant here because
+ they are inputs, and their HDUs are given within the
+ operator string. */
+ if( strncmp(GAL_ARITHMETIC_OPSTR_LOADCOL_PREFIX, token->v,
+ GAL_ARITHMETIC_OPSTR_LOADCOL_PREFIX_LEN)
+ && gal_array_name_recognized_multiext(token->v) )
++nummultiext;
/* If no output name is given, we need to extract the output
@@ -315,19 +326,12 @@ ui_check_options_and_arguments(struct arithmeticparams *p)
basename=token->v;
}
- /* This token is a number. Check if a negative dash was present that
- has been temporarily replaced with 'NEG_DASH_REPLACE' before
- option parsing. */
+ /* This token is a number. Check if a negative dash was present
+ that has been temporarily replaced with 'NEG_DASH_REPLACE'
+ before option parsing. */
else if(token->v[0]==NEG_DASH_REPLACE && isdigit(token->v[1]) )
token->v[0]='-';
}
-
- /* We are on the 'tofile' operator. */
- else
- {
- filename=&token->v[ OPERATOR_PREFIX_LENGTH_TOFILE ];
- gal_checkset_writable_remove(filename, cp->keep, cp->dontdelete);
- }
}
/* In case no output name has been given (can happen with operators like
@@ -342,8 +346,9 @@ ui_check_options_and_arguments(struct arithmeticparams *p)
"_arith.fits");
else
{
- gal_checkset_allocate_copy("arithmetic.fits", &p->cp.output);
- gal_checkset_writable_remove(p->cp.output, cp->keep, cp->dontdelete);
+ gal_checkset_allocate_copy("arithmetic.fits", &cp->output);
+ gal_checkset_writable_remove(cp->output, cp->keep,
+ cp->dontdelete);
}
}
diff --git a/bin/table/arithmetic.c b/bin/table/arithmetic.c
index f5c38ed3..bac5f669 100644
--- a/bin/table/arithmetic.c
+++ b/bin/table/arithmetic.c
@@ -58,6 +58,7 @@ arithmetic_add_new_to_end(struct arithmetic_token **list)
/* Initialize its elements. */
node->next=NULL;
+ node->loadcol=NULL;
node->name_def=NULL;
node->name_use=NULL;
node->constant=NULL;
@@ -246,7 +247,7 @@ arithmetic_init_addcol(char *token, struct arithmetic_token
*node,
if(num==GAL_BLANK_SIZE_T) /* ID is a string. */
{
for(i=0;i<numcols;++i)
- if( !strcasecmp(id, colinfo[i].name) )
+ if( colinfo[i].name && !strcasecmp(id, colinfo[i].name) )
{ inmaininput=1; break; }
}
else /* ID is a number. */
@@ -264,7 +265,7 @@ arithmetic_init_addcol(char *token, struct arithmetic_token
*node,
}
else
{
- gal_checkset_allocate_copy(id, &node->id_at_usage); /* Future usage. */
+ gal_checkset_allocate_copy(id, &node->id_at_usage);/* Future usage. */
node->num_at_usage=num; /* Avoid converting to a number again! */
node->index=GAL_BLANK_SIZE_T; /* Crash if not used properly! */
}
@@ -283,9 +284,11 @@ arithmetic_init(struct tableparams *p, struct
arithmetic_token **arith,
void *num;
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;
/* Parse all the given tokens. */
token=strtok_r(expression, delimiter, &saveptr);
@@ -313,6 +316,13 @@ arithmetic_init(struct tableparams *p, struct
arithmetic_token **arith,
gal_checkset_allocate_copy(token, &node->name_def);
}
+ /* This is a 'load-col' operator. */
+ else if( (col=gal_arithmetic_load_col(token, cp->searchin,
+ cp->ignorecase,
+ cp->minmapsize,
+ cp->quietmmap))!=NULL )
+ node->loadcol=col;
+
/* Non-pre-defined situation.*/
else
{
@@ -336,8 +346,10 @@ arithmetic_init(struct tableparams *p, struct
arithmetic_token **arith,
token=strtok_r(NULL, delimiter, &saveptr);
}
- /* A small sanity check: the last added token must be an operator. */
- if( node==NULL || node->operator==GAL_ARITHMETIC_OP_INVALID )
+ /* 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);
}
@@ -872,9 +884,9 @@ arithmetic_operator_run(struct tableparams *p,
break;
default:
- error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
- "the problem. '%zu' is not recognized as an operand "
- "counter (with '%s')", __func__, PACKAGE_BUGREPORT,
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s "
+ "to fix the problem. '%zu' is not recognized as an "
+ "operand counter (with '%s')", __func__, PACKAGE_BUGREPORT,
token->num_operands,
arithmetic_operator_name(token->operator));
}
@@ -997,12 +1009,20 @@ arithmetic_reverse_polish(struct tableparams *p,
if(token->operator!=GAL_ARITHMETIC_OP_INVALID)
arithmetic_operator_run(p, token, &setprm, &stack);
+
/* We are on a named variable. */
else if( token->name_use )
gal_list_data_add(&stack,
gal_arithmetic_set_copy_named(&setprm,
token->name_use));
+ /* We are on a column loaded from another file. */
+ else if( token->loadcol )
+ {
+ gal_list_data_add(&stack, token->loadcol);
+ token->loadcol=NULL;
+ }
+
/* Constant number: just put it ontop of the stack. */
else if(token->constant)
{
diff --git a/bin/table/main.h b/bin/table/main.h
index cf91975d..73519fe0 100644
--- a/bin/table/main.h
+++ b/bin/table/main.h
@@ -69,6 +69,7 @@ struct arithmetic_token
char *id_at_usage; /* OPERAND: col identifier at usage time. */
size_t num_at_usage; /* OPERAND: col number at usage. */
gal_data_t *constant; /* OPERAND: a constant/single number. */
+ gal_data_t *loadcol; /* OPERAND: column from another file. */
char *name_def; /* Name given to the 'set-' operator. */
char *name_use; /* If this a usage of a name. */
struct arithmetic_token *next; /* Pointer to next token. */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 5ab12c00..68a7e053 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -505,8 +505,9 @@ Arithmetic operators
* Mathematical morphology operators:: Work on binary images, for example
erode.
* Bitwise operators:: Work on bits within one pixel.
* Numerical type conversion operators:: Convert the numeric datatype of a
dataset.
-* Adding noise operators:: Add noise to a dataset.
+* Random number generators:: Random numbers can be used to add noise for
example.
* Elliptical shape operators:: Operations that are focused on an ellipse.
+* Loading external columns:: Read a column from a table into the stack.
* Building new dataset:: How to construct an empty dataset from scratch.
* Operand storage in memory or a file:: Tools for complex operations in one
command.
@@ -13291,6 +13292,15 @@ asttable f160w-cat.fits --output=both.fits \
@noindent
For a more complete example, see @ref{Working with catalogs estimating colors}.
+@cartouche
+@noindent
+@strong{Loading external columns with Arithmetic:} an alternative way to load
external columns into your output is to use column arithmetic (@ref{Column
arithmetic})
+In particular the @option{load-col-} operator described in @ref{Loading
external columns}.
+But this operator will load only one column per file/HDU everytime it is
called.
+So if you have many columns to insert, it is much faster to use
@option{--catcolumnfile}.
+Because @option{--catcolumnfile} will load all the columns in one opening of
the file, and possibly even read them all into memory in parallel!
+@end cartouche
+
@item -u STR/INT
@itemx --catcolumnhdu=STR/INT
The HDU/extension of the FITS file(s) that should be concatenated, or
appended, by column with @option{--catcolumnfile}.
@@ -14994,8 +15004,9 @@ Reading NaN as a floating point number in Gnuastro
isn't case-sensitive.
* Mathematical morphology operators:: Work on binary images, for example
erode.
* Bitwise operators:: Work on bits within one pixel.
* Numerical type conversion operators:: Convert the numeric datatype of a
dataset.
-* Adding noise operators:: Add noise to a dataset.
+* Random number generators:: Random numbers can be used to add noise for
example.
* Elliptical shape operators:: Operations that are focused on an ellipse.
+* Loading external columns:: Read a column from a table into the stack.
* Building new dataset:: How to construct an empty dataset from scratch.
* Operand storage in memory or a file:: Tools for complex operations in one
command.
@end menu
@@ -16123,7 +16134,7 @@ For example (assuming numbers can be written as bit
strings on the command-line)
Note that the bitwise operators only work on integer type datasets/numbers.
@end table
-@node Numerical type conversion operators, Adding noise operators, Bitwise
operators, Arithmetic operators
+@node Numerical type conversion operators, Random number generators, Bitwise
operators, Arithmetic operators
@subsubsection Numerical type conversion operators
With the operators below you can convert the numerical data type of your
input, see @ref{Numeric data types}.
@@ -16178,8 +16189,8 @@ Convert the type of the popped operand to 64-bit
(double precision) floating poi
The internal conversion of C will be used.
@end table
-@node Adding noise operators, Elliptical shape operators, Numerical type
conversion operators, Arithmetic operators
-@subsubsection Adding noise operators
+@node Random number generators, Elliptical shape operators, Numerical type
conversion operators, Arithmetic operators
+@subsubsection Random number generators
When you simulate data (for example see @ref{Sufi simulates a detection}),
everything is ideal and there is no noise!
The final step of the process is to add simulated noise to the data.
@@ -16264,9 +16275,194 @@ echo 12 | asttable -c'arith $1 4 mknoise-uniform'
@end example
Similar to the example in @code{mknoise-sigma}, you can pipe the output of
@command{echo} to @command{awk} before passing it to @command{asttable} to
generate a full column of uniformly selected values within the same interval.
+
+@item random-from-hist-raw
+Generate random values from a custom distribution (defined by a histogram).
+The output will have a double-precision floating point type (see @ref{Numeric
data types}).
+This operator takes three operands:
+@itemize
+@item
+The first popped operand (nearest to the operator) is the histogram values.
+The histogram is a 1-dimensional dataset (a table column) and contains the
probability of obtaining a certain interval of values.
+The histogram doesn't have to be normalized: the GNU Scientific Library (or
GSL, which is used by Gnuastro for this operator), will normalize it internally.
+The value of each bin (whose probability is given in the histogram) is given
in the second popped operand.
+Therefore these two operands have to have the same number of rows.
+@item
+The second popped operand is the bin value (mostly the bin center, but it can
be anything).
+The probability of each bin is defined in the histogram operand (first popped
operand).
+The bins can have any width (don't have to be evenly spaced), and any order.
+Just make sure that the same row in the bins column corresponds to the same
row in the histogram: the number of rows in the bins and histogram must be
equal.
+@item
+The third popped operand is the dataset that the random values should be
written over.
+Effectively only its size will be used by this operator (all values will be
over-written as a double-precision floating point number).
+@end itemize
+
+The first two operands have to be single-dimensional (a table column) and have
the same number of rows, but the last popped operand can have any number of
dimensions.
+You can use the @code{load-col-} operator to load the two bins and histogram
columns from an external file (see @ref{Loading external columns}).
+
+For example, in the command below, we first construct a fake histogram to
represent a @mymath{y=x^2} distribution with AWK.
+We aim to distribute random values from this distribution in a
@mymath{100\times100} image.
+Therefore, we use the @command{makenew} operator to construct an empty image
of that size, use the @command{load-col-} operator to load the histogram
columns into Arithmetic and put the output in @file{random.fits}.
+Finally we visually inspect @file{random.fits} with DS9 and also have a look
at its pixel distribution with @command{aststatistics}.
+
+@example
+$ echo "" | awk '@{for(i=1;i<5;++i) print i, i*i@}' \
+ > histogram.txt
+
+$ cat histogram.txt
+1 1
+2 4
+3 9
+4 16
+
+$ astarithmetic 100 100 2 makenew \
+ load-col-1-from-histogram.txt \
+ load-col-2-from-histogram.txt \
+ random-from-hist-raw \
+ --output=random.fits
+
+$ astscript-fits-view random.fits
+
+$ aststatistics random.fits --asciihist --numasciibins=50
+ | *
+ | *
+ | *
+ | *
+ | * *
+ | * *
+ | * *
+ | * * *
+ | * * *
+ |* * * *
+ |* * * *
+ |--------------------------------------------------
+@end example
+
+As you see, the 10000 pixels in the image only have values 1, 2, 3 or 4 (which
were the values in the bins column of @file{histogram.txt}), and the number of
times each of these values occurs follows the @mymath{y=x^2} distribution.
+
+Generally, any value given in the bins column will be used for the final
output values.
+For example in the command below, in the step that we generated the histogram,
we are adding the bins by 20 (while keeping the same probability distribution
of @mymath{y=x^2}).
+If you re-run the Arithmetic command above after this, you will notice that
the pixels values are now one of the following 21, 22, 23 or 24 (instead of 1,
2, 3, or 4).
+But the shape of the histogram of the resulting random distribution will be
unchanged.
+
+@example
+$ echo "" | awk '@{for(i=1;i<5;++i) print 20+i, i*i@}' \
+ > histogram.txt
+@end example
+
+If you don't want the outputs to have exactly the value of the bin identifier,
but be a randomly selected value from a uniform distribution within the bin,
you should use @command{random-from-hist} (see below).
+
+As mentioned above, the output will have a double-precision floating point
type (see @ref{Numeric data types}).
+Therefore, by default each element of the output will consume 8 bytes
(64-bits) of storage.
+This is usually far more than the statistical error/precision of your data
(and just results in wasted storage in your file system, or wasted RAM when a
program that uses the data is being run, and a slower running time of the
program).
+
+It is therefore recommended to use a type-conversion operator after this
operator to put the output in the smallest type that can be used to safely
store your data without wasting storage, RAM or time.
+For the list of type conversion operators, see @ref{Numerical type conversion
operators}.
+Recall that you already know the values returned by this operator (they are
one of the values in the bins column).
+
+For example, in the example above, the whole image only have values 1, 2, 3 or
4.
+Since they are always positive and are below 255, we can safely place them in
an unsigned 8-bit integer (see @ref{Numeric data types}) with the command below
(note the @code{uint8} after the operator name, and that we are using a
different name for the output).
+After builing the image, we are listing the two outputs, while printing their
size with @command{ls -l}:
+
+@example
+$ astarithmetic 100 100 2 makenew \
+ load-col-1-from-histogram.txt \
+ load-col-2-from-histogram.txt \
+ random-from-hist-raw uint8 \
+ --output=random-u8.fits
+
+$ ls -lh random.fits random-u8.fits
+-rw-r--r-- 1 name name 85K Jan 01 13:40 random.fits
+-rw-r--r-- 1 name name 17K Jan 01 13:45 random-u8.fits
+@end example
+
+As you see, when using the proper data type, we can shrink the size of the
file significantly without loosing any information (from 85 kilobytes to 17
kilobytes).
+This difference can be felt much better for larger (real-world) datasets, so
be sure to consider the output data type when using this operator.
+
+
+@item random-from-hist
+Similar to @code{random-from-hist-raw}, but don't return the exact bin value,
instead return a random value from a uniform distribution within each bin.
+Therefore the following limitations have to be taken into account (compared to
@code{random-from-hist-raw}:
+@itemize
+@item
+The bins should be represented by their center (the number in the bin column).
+@item
+The bins have to be in descending order.
+@item
+The bin widths (distance from one bin to another) have to be fixed.
+@end itemize
+
+For example, in the commands below, let's replace @code{random-from-hist-raw}
with @code{random-from-hist} in the example of the description of
@code{random-from-hist-raw}.
+Note how we are manually converting the output of this operator into
single-precision floating point (32-bit, since 64-bit precision is
statistically meaningless in this scenario and we don't want to waste storage,
memory and running time):
+@example
+$ echo "" | awk '@{for(i=1;i<5;++i) print i, i*i@}' \
+ > histogram.txt
+
+$ astarithmetic 100 100 2 makenew \
+ load-col-1-from-histogram.txt \
+ load-col-2-from-histogram.txt \
+ random-from-hist float32 \
+ --output=random.fits
+
+$ aststatistics random.fits --asciihist --numasciibins=50
+ | *
+ | *** ********
+ | ************
+ | *************
+ | * * *************
+ | * ***********************
+ | *************************
+ | *************************
+ | *************************************
+ |********* * **************************************
+ |**************************************************
+ |--------------------------------------------------
+@end example
+
+You can see that the pixels of @file{histogram.fits} are no longer just 1, 2,
3 or 4.
+Instead, the values within each bin are selected from a uniform distribution
covering that bin.
+This creates the step-like feature in the histogram of the output.
+
+Of course, this extra uniform random number generation can make your program
slower so be sure to check if it is worth it.
+In particular, one way to avoid this (and use @command{random-from-hist-raw}
with a more contiguous-looking output distribution) is to simply use a
higher-resolution histogram (assuming it is possible: you have a sufficient
number of data points, or its an analytical expression that you can sample at
smaller bin sizes).
+
+To better demonstrate this operator and its practical usage in everyday
research, let's look at another example:
+Assume you want to get 100 random star magnitudes that follow the real-world
Gaia Data release 3 magnitude distribution within a radius of 2 degrees around
the (RA,Dec) coordinate of (1.23,4.56).
+Let's further assume that you want to distribute them uniformly over an image
of size 1000 by 1000 pixels.
+So your desired output table should have three columns, the first two are
pixel positions of each star, and the third is the magnitude.
+
+First, we need to query the Gaia database and ask for all the magnitudes in
this region of the sky.
+We know that Gaia is not complete for stars fainter than the 20th magnitude,
so we'll use the @option{--range} option and only ask for those stars that are
brighter than magnitude 20.
+
+@example
+$ astquery gaia --dataset=dr3 --center=1.23,3.45 --radius=2 \
+ --column=phot_g_mean_mag --output=gaia.fits \
+ --range=phot_g_mean_mag,-inf,20
+@end example
+
+We now have more than 25000 magnitudes in @file{gaia.fits}!
+To get a more accurate random sampling of our stars, let's construct a
histogram with 500 bins, and generate our three desired randomly selected
columns:
+
+@example
+$ aststatistics gaia.fits --histogram --numbins=500 \
+ --output=gaia-hist.fits
+
+$ asttable gaia-hist.fits -i
+
+$ echo 1000 \
+ | awk '@{for(i=0;i<100;++i) print $1/2@}' \
+ | asttable -c'arith $1 500 mknoise-uniform' \
+ -c'arith $1 500 mknoise-uniform' \
+ -c'arith $1 \
+ load-col-1-from-gaia-hist.fits-hdu-1 \
+ load-col-2-from-gaia-hist.fits-hdu-1 \
+ random-from-hist float32'
+@end example
+
+These colums can easily be placed in the format for @ref{MakeProfiles} to be
inserted into an image automatically.
@end table
-@node Elliptical shape operators, Building new dataset, Adding noise
operators, Arithmetic operators
+@node Elliptical shape operators, Loading external columns, Random number
generators, Arithmetic operators
@subsubsection Elliptical shape operators
The operators here describe certain functions that will be necessary when
dealing with objects that have a certain elliptical shape.
@@ -16314,7 +16510,45 @@ $ asttable catalog.fits \
@end example
@end table
-@node Building new dataset, Operand storage in memory or a file, Elliptical
shape operators, Arithmetic operators
+@node Loading external columns, Building new dataset, Elliptical shape
operators, Arithmetic operators
+@subsubsection Loading external columns
+
+In the Arithmetic program, you can always load new dataset by simply giving
their name.
+However, they can only be images, not a column.
+In the Table program, you can load columns in @ref{Column arithmetic}, but it
has to be columns within the same table (and thus the same number of rows).
+However, in some situations, it is necessary to use certain columns of a table
in the Arithmetic program, or columns of different rows (from the main input)
in Table.
+
+@table @command
+@item load-col-%-from-%
+@itemx load-col-%-from-%-hdu-%
+Load the requested column (first @command{%}) from the requested file (second
@command{%}).
+If the file is a FITS file, it is also necessary to specify a HDU using the
second form (where the HDU identifier is the third @command{%}.
+For example @command{load-col-MAG-from-catalog.fits-hdu-1} will load the
@code{MAG} column from HDU 1 of @code{catalog.fits}.
+
+For example, let's assume you have the following two tables, and you would
like to add the first column of the first with the second:
+
+@example
+$ asttable tab-1.fits
+1 43.23
+2 21.91
+3 71.28
+4 18.10
+
+$ cat tab-2.txt
+5
+6
+7
+8
+
+$ asttable tab-1.txt -c'arith $1 load-col-1-from-tab-2.txt +'
+6
+8
+10
+12
+@end example
+@end table
+
+@node Building new dataset, Operand storage in memory or a file, Loading
external columns, Arithmetic operators
@subsubsection Building new dataset
With the operator here, you can create a new dataset from scratch to start
certain operations without any input data.
@@ -16331,7 +16565,7 @@ $ astarithmetic 100 200 2 makenew
@end example
@noindent
-To further extend the example, you can use any of the noise-making operators
to add noise to this new dataset (see @ref{Adding noise operators}), like the
command below:
+To further extend the example, you can use any of the noise-making operators
to add noise to this new dataset (see @ref{Random number generators}), like the
command below:
@example
$ astarithmetic 100 200 2 makenew 5 mknoise-sigma
@@ -28223,12 +28457,9 @@ tiles that are created from this pointer.
@node Dataset allocation, Arrays of datasets, Generic data container, Library
data container
@subsubsection Dataset allocation
-
Gnuastro's main data container was defined in @ref{Generic data container}.
-The functions listed in this section describe the most basic operations on
-@code{gal_data_t}: those related to allocation and freeing. These functions
-are declared in @file{gnuastro/data.h} which is also visible from the
-function names (see @ref{Gnuastro library}).
+The functions listed in this section describe the most basic operations on
@code{gal_data_t}: those related to allocation and freeing.
+These functions are declared in @file{gnuastro/data.h} which is also visible
from the function names (see @ref{Gnuastro library}).
@deftypefun {gal_data_t *} gal_data_alloc (void @code{*array}, uint8_t
@code{type}, size_t @code{ndim}, size_t @code{*dsize}, struct wcsprm
@code{*wcs}, int @code{clear}, size_t @code{minmapsize}, int @code{quietmmap},
char @code{*name}, char @code{*unit}, char @code{*comment})
@@ -29814,7 +30045,7 @@ suffix doesn't have to start with `@key{.}': this
function will return
@code{1} (one) for both @code{fits} and @code{.fits}.
@end deftypefun
-@deftypefun in gal_fits_file_recognized (char @code{*name})
+@deftypefun int gal_fits_file_recognized (char @code{*name})
Return @code{1} if the given file name (possibly including its contents) is a
FITS file.
This is necessary in when the contents of a FITS file do follow the FITS
standard, but it the file doesn't have a Gnuastro-recognized FITS suffix.
Therefore, it will first call @code{gal_fits_name_is_fits}, if the result is
negative, then this function will attempt to open the file with CFITSIO and if
it works, it will close it again and return 1.
@@ -31367,6 +31598,12 @@ For more on random number generation in Gnuastro, see
@ref{Generating random num
By default these operators will print the random number generator function and
seed (in case the user wants to reproduce the result later), but this can be
disabled by activating the bit-flag @code{GAL_ARITHMETIC_FLAG_QUIET} described
above.
@end deffn
+@deffn Macro GAL_ARITHMETIC_OP_RANDOM_FROM_HIST
+@deffnx Macro GAL_ARITHMETIC_OP_RANDOM_FROM_HIST_RAW
+Select random values from a custom distribution (defined by a histogram).
+For more, see the description of the respective operators in @ref{Generating
random numbers}.
+@end deffn
+
@deffn Macro GAL_ARITHMETIC_OP_SIZE
Size operator that will return a single value for datasets of any kind. When
@code{gal_arithmetic} is called with this operator, it requires two arguments.
The first is the dataset, and the second is a single integer value.
@@ -31436,6 +31673,16 @@ The number of dimensions is derived from the number of
nodes in the list and the
Just note that the list should be in the reverse of the desired dimensions.
@end deffn
+@deffn Macro GAL_ARITHMETIC_OPSTR_LOADCOL_HDU
+@deffnx Macro GAL_ARITHMETIC_OPSTR_LOADCOL_FILE
+@deffnx Macro GAL_ARITHMETIC_OPSTR_LOADCOL_PREFIX
+@deffnx Macro GAL_ARITHMETIC_OPSTR_LOADCOL_HDU_LEN
+@deffnx Macro GAL_ARITHMETIC_OPSTR_LOADCOL_FILE_LEN
+@deffnx Macro GAL_ARITHMETIC_OPSTR_LOADCOL_PREFIX_LEN
+Constant components of the @command{load-col-} operator (see @ref{Loading
external columns}).
+These are just fixed strings (and their lengths) that are placed in between
the various components of that operator to allow choosing a certain column of a
certain HDU of a certain file.
+@end deffn
+
@deftypefun {gal_data_t *} gal_arithmetic (int @code{operator}, size_t
@code{numthreads}, int @code{flags}, ...)
Apply the requested arithmetic operator on the operand(s).
The @emph{operator} is identified through the macros above (that start with
@code{GAL_ARITHMETIC_OP_}).
@@ -31509,6 +31756,14 @@ Return the human-readable standard string that
corresponds to the given operator
For example when the input is @code{GAL_ARITHMETIC_OP_PLUS} or
@code{GAL_ARITHMETIC_OP_MEAN}, the strings @code{+} or @code{mean} will be
returned.
@end deftypefun
+@deftypefun {gal_data_t *} gal_arithmetic_load_col (char @code{*str}, int
@code{searchin}, int @code{ignorecase}, size_t @code{minmapsize}, int
@code{quietmmap})
+Return the column that corresponds to the identifier in the input string
(@code{str}).
+@code{str} is expected to be in the format of the @command{load-col-} operator
(see @ref{Loading external columns}).
+This function will extract the column identifier, the file name and the HDU
(if necessary) from the string, read the requested column in memory and return
it.
+
+See @ref{Table input output} for the macros that can be given to
@code{searchin} and @code{ignorecase} and @ref{Generic data container} for the
definitions of @code{minmapsize} and @code{quietmmap}.
+@end deftypefun
+
@node Tessellation library, Bounding box, Arithmetic on datasets, Gnuastro
library
@subsection Tessellation library (@file{tile.h})
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 9199894b..e74b05ea 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -759,6 +759,51 @@ arithmetic_to_oned(int operator, int flags, gal_data_t
*input)
/*************** Adding noise **************/
/***********************************************************************/
+static gsl_rng *
+arithmetic_gsl_initialize(int flags, const char **rng_name,
+ unsigned long *rng_seed, char *operator_str)
+{
+ gsl_rng *rng;
+
+ /* Column counter in case '--envseed' is given and we have multiple
+ columns. */
+ static unsigned long colcounter=0;
+
+ /* Setup the random number generator. For 'envseed', we want to pass a
+ boolean value: either 0 or 1. However, when we say 'flags &
+ GAL_ARITHMETIC_ENVSEED', the returned value is the integer positioning
+ of the envseed bit (for example if its on the fourth bit, the value
+ will be 8). This can cause problems if it is on the 8th bit (or any
+ multiple of 8). So to avoid issues with the bit-positioning of the
+ 'ENVSEED', we will return the conditional to see if the result of the
+ bit-wise '&' is larger than 0 or not (binary). */
+ rng=gal_checkset_gsl_rng( (flags & GAL_ARITHMETIC_FLAG_ENVSEED)>0,
+ rng_name, rng_seed);
+
+ /* If '--envseed' was called, we need to add the column counter to the
+ requested seed (so its not the same for all columns. */
+ if(flags & GAL_ARITHMETIC_FLAG_ENVSEED)
+ {
+ *rng_seed += colcounter++;
+ gsl_rng_set(rng, *rng_seed);
+ }
+
+ /* Print the basic RNG information if requested. */
+ if( (flags & GAL_ARITHMETIC_FLAG_QUIET)==0 )
+ {
+ printf(" - Parameters used for '%s':\n", operator_str);
+ printf(" - Random number generator name: %s\n", *rng_name);
+ printf(" - Random number generator seed: %lu\n", *rng_seed);
+ }
+
+ /* Return the GSL random number generator */
+ return rng;
+}
+
+
+
+
+
/* The size operator. Reports the size along a given dimension. */
static gal_data_t *
arithmetic_mknoise(int operator, int flags, gal_data_t *in,
@@ -770,10 +815,6 @@ arithmetic_mknoise(int operator, int flags, gal_data_t *in,
unsigned long rng_seed;
gal_data_t *out, *targ;
- /* Column counter in case '--envseed' is given and we have multiple
- columns. */
- static unsigned long colcounter=0;
-
/* 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;
@@ -812,33 +853,9 @@ arithmetic_mknoise(int operator, int flags, gal_data_t *in,
"range for 'mknoise-uniform') must be positive (it is %g)",
arg_v);
- /* Setup the random number generator. For 'envseed', we want to pass a
- boolean value: either 0 or 1. However, when we say 'flags &
- GAL_ARITHMETIC_ENVSEED', the returned value is the integer positioning
- of the envseed bit (for example if its on the fourth bit, the value
- will be 8). This can cause problems if it is on the 8th bit (or any
- multiple of 8). So to avoid issues with the bit-positioning of the
- 'ENVSEED', we will return the conditional to see if the result of the
- bit-wise '&' is larger than 0 or not (binary). */
- rng=gal_checkset_gsl_rng( (flags & GAL_ARITHMETIC_FLAG_ENVSEED)>0,
- &rng_name, &rng_seed);
-
- /* If '--envseed' was called, we need to add the column counter to the
- requested seed. */
- if(flags & GAL_ARITHMETIC_FLAG_ENVSEED)
- {
- rng_seed += colcounter++;
- gsl_rng_set(rng, rng_seed);
- }
-
- /* Print the basic RNG information if requested. */
- if( (flags & GAL_ARITHMETIC_FLAG_QUIET)==0 )
- {
- printf(" - Parameters used for '%s':\n",
- gal_arithmetic_operator_string(operator));
- printf(" - Random number generator name: %s\n", rng_name);
- printf(" - Random number generator seed: %lu\n", rng_seed);
- }
+ /* 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;
@@ -864,6 +881,7 @@ arithmetic_mknoise(int operator, int flags, gal_data_t *in,
}
/* Clean up and return */
+ gsl_rng_free(rng);
if(flags & GAL_ARITHMETIC_FLAG_FREE) gal_data_free(arg);
return out;
}
@@ -872,6 +890,149 @@ arithmetic_mknoise(int operator, int flags, gal_data_t
*in,
+/* Sanity checks for 'random_from_hist' */
+static void
+arithmetic_random_from_hist_sanity(gal_data_t **inhist, gal_data_t **inbinc,
+ gal_data_t *in, int operator)
+{
+ size_t i;
+ double *d, diff, binwidth, binwidth_e;
+ gal_data_t *hist=*inhist, *binc=*inbinc;
+
+ /* The 'hist' and 'bincenter' arrays should be 1-dimensional and both
+ should have the same size. */
+ if(hist->ndim!=1)
+ error(EXIT_FAILURE, 0, "the first popped (nearest) operand to the "
+ "'%s' operator should only have a single dimension. However, "
+ "it has %zu dimensions",
+ gal_arithmetic_operator_string(operator), hist->ndim);
+ if(binc->ndim!=1)
+ error(EXIT_FAILURE, 0, "the second popped (second nearest) operand "
+ "to the '%s' operator should only have a single "
+ "dimension. However, it has %zu dimensions",
+ gal_arithmetic_operator_string(operator), binc->ndim);
+ if(hist->size!=binc->size)
+ error(EXIT_FAILURE, 0, "the first and second popped operands to the "
+ "'%s' operator must have the same size! However, "
+ "they have %zu and %zu elements respectively",
+ gal_arithmetic_operator_string(operator), hist->size,
+ binc->size);
+
+ /* The input type shouldn't be a string. */
+ 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 it has a string type",
+ gal_arithmetic_operator_string(operator));
+
+ /* The histogram should be in float64 type, we'll also set the bin
+ centers to float64 to enable easy sanity checks. */
+ *inhist=gal_data_copy_to_new_type_free(hist, GAL_TYPE_FLOAT64);
+ binc=*inbinc=gal_data_copy_to_new_type_free(binc, GAL_TYPE_FLOAT64);
+
+ /* For the 'random-from-hist' operator, we will need to assume that the
+ bins are equally spaced and that they are ascending.*/
+ if(operator==GAL_ARITHMETIC_OP_RANDOM_FROM_HIST)
+ {
+ d=binc->array;
+ binwidth=d[1]-d[0];
+ if(binwidth<0)
+ error(EXIT_FAILURE, 0, "the bins given to the '%s' operator "
+ "should be in ascending order, but the second row "
+ "(%g) is smaller than the first (%g)",
+ gal_arithmetic_operator_string(operator), d[1], d[0]);
+
+ /* Make sure the bins are in ascending order and have a fixed bin
+ width (within floating point errors: hence where 1e-6 comes
+ from). */
+ binwidth_e=binwidth*1e-6;
+ for(i=0;i<binc->size-1;++i)
+ {
+ diff=d[i+1]-d[i];
+ if(diff>binwidth+binwidth_e || diff<binwidth-binwidth_e)
+ error(EXIT_FAILURE, 0, "the bins given to the '%s' operator "
+ "should be in ascending order, but row %zu (with value "
+ "%g) is larger than the next row's value (%g)",
+ gal_arithmetic_operator_string(operator), i+1, d[i],
+ d[i+1]);
+ }
+ }
+}
+
+
+
+
+
+/* Build a custom noise distribution. */
+static gal_data_t *
+arithmetic_random_from_hist(int operator, int flags, gal_data_t *in,
+ gal_data_t *binc, gal_data_t *hist)
+{
+ size_t rind;
+ const char *rng_name;
+ gal_data_t *out=NULL;
+ unsigned long rng_seed;
+ gsl_rng *rng=NULL, *rngu=NULL;
+ gsl_ran_discrete_t *disc=NULL;
+ double *b, *d, *df, binwidth, halfbinwidth;
+
+ /* 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;
+
+ /* Basic sanity checks */
+ arithmetic_random_from_hist_sanity(&hist, &binc, in,
+ GAL_ARITHMETIC_OP_RANDOM_FROM_HIST);
+
+ /* Allocate the output dataset (based on the size and dimensions of the
+ input), then free the input. */
+ out=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, in->ndim, in->dsize, in->wcs,
+ 0, in->minmapsize, in->quietmmap, "RANDOM_FROM_HIST",
+ hist->unit, "Random values from custom distribution");
+
+ /* Initialize the GSL random number generator. */
+ rng=arithmetic_gsl_initialize(flags, &rng_name, &rng_seed,
+ gal_arithmetic_operator_string(operator));
+ if(operator==GAL_ARITHMETIC_OP_RANDOM_FROM_HIST)
+ rngu=arithmetic_gsl_initialize(flags, &rng_name, &rng_seed,
+ "uniform component of 'random-"
+ "from-hist'");
+
+ /* Pre-process the descrete random number generator. */
+ disc=gsl_ran_discrete_preproc(hist->size, hist->array);
+
+ /* Apply the random number generator over the output. Note that
+ 'gsl_ran_discrete' returns the index to the bin centers table that we
+ can return. */
+ b=binc->array;
+ binwidth=b[1]-b[0];
+ halfbinwidth=binwidth/2;
+ df=(d=out->array)+out->size;
+ do
+ {
+ rind=gsl_ran_discrete(rng, disc);
+ *d = ( operator==GAL_ARITHMETIC_OP_RANDOM_FROM_HIST
+ ? ((b[rind]-halfbinwidth) + gsl_rng_uniform(rngu)*binwidth)
+ : b[rind] );
+ }
+ while(++d<df);
+
+ /* Clean up and return. */
+ gsl_rng_free(rng);
+ if(rngu) gsl_rng_free(rngu);
+ gsl_ran_discrete_free(disc);
+ if(flags & GAL_ARITHMETIC_FLAG_FREE)
+ {
+ gal_data_free(in);
+ gal_data_free(hist);
+ gal_data_free(binc);
+ }
+ return out;
+}
+
+
+
+
@@ -2419,7 +2580,24 @@ arithmetic_box_around_ellipse(gal_data_t *d1, gal_data_t
*d2,
-/* Make a new dataset. */
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/**********************************************************************/
+/**************** New datasets *****************/
+/**********************************************************************/
gal_data_t *
arithmetic_makenew(gal_data_t *sizes)
{
@@ -2474,6 +2652,113 @@ arithmetic_makenew(gal_data_t *sizes)
+gal_data_t *
+gal_arithmetic_load_col(char *str, int searchin, int ignorecase,
+ size_t minmapsize, int quietmmap)
+{
+ gal_data_t *out=NULL;
+ gal_list_str_t *colid=NULL;
+ char *c, *cf, *copy, *hdu=NULL, *filename=NULL;
+ size_t numthreads=1; /* We only want to read a single column! */
+
+ /* This is the shortest possible string (with each component being given
+ a one character value). Recall that in C, simply putting literal
+ strings after each other, will merge them together into one literal
+ string.*/
+ char *checkstr = ( GAL_ARITHMETIC_OPSTR_LOADCOL_PREFIX "a"
+ GAL_ARITHMETIC_OPSTR_LOADCOL_FILE "a" );
+
+ /* This function is called on every call of Arithmetic, so before going
+ into any further tests, first make sure the string is long enough, and
+ that it starts with the fixed format string 'FMTCOL'. */
+ if( strlen(str)<strlen(checkstr)
+ || strncmp(str, GAL_ARITHMETIC_OPSTR_LOADCOL_PREFIX,
+ GAL_ARITHMETIC_OPSTR_LOADCOL_PREFIX_LEN) ) return NULL;
+
+ /* To separate the components, we need to put '\0's within the input
+ string. But we don't want to ruin the input string (in case it isn't
+ for this purpose), so for the rest of the steps we'll make a copy of
+ the input string and work on that. */
+ gal_checkset_allocate_copy(str, ©);
+
+ /* Parse the string and separate the components. */
+ gal_list_str_add(&colid,
+ ©[GAL_ARITHMETIC_OPSTR_LOADCOL_PREFIX_LEN], 0);
+ cf=copy+strlen(str);
+ c=colid->v+1; /* colID has at least one character long, so we'll */
+ while(c<cf) /* start the parsing from the next character */
+ {
+ /* If we are on the file-name component, then we can set the end of
+ the column ID component, and set the start of the HDU. But this is
+ only valid if 'hdu' hasn't already been set. */
+ if( !strncmp(c, GAL_ARITHMETIC_OPSTR_LOADCOL_FILE,
+ GAL_ARITHMETIC_OPSTR_LOADCOL_FILE_LEN ) )
+ {
+ /* If 'filename' or 'hdu' have already been set, then this
+ doesn't conform to the format, and we should leave this
+ function. */
+ if(filename || hdu) { free(copy); return NULL; };
+
+ /* Set the current position to '\0' (to end the column name) */
+ *c='\0';
+
+ /* Set the HDU's starting pointer. */
+ filename=c+GAL_ARITHMETIC_OPSTR_LOADCOL_FILE_LEN;
+ c=filename+1; /* Similar to 'colid->v' above. */
+ }
+
+ /* If we are on a HDU component, steps are very similar to the
+ filename steps above. */
+ else if( !strncmp(c, GAL_ARITHMETIC_OPSTR_LOADCOL_HDU,
+ GAL_ARITHMETIC_OPSTR_LOADCOL_HDU_LEN) )
+ {
+ if(hdu) { free(copy); return NULL; }
+ *c='\0'; hdu=c+GAL_ARITHMETIC_OPSTR_LOADCOL_HDU_LEN;
+ c=hdu+1;
+ }
+
+ /* If there was no match with HDU or file strings, then simply
+ increment the pointer. */
+ else ++c;
+ }
+
+ /* If a file-name couldn't be identified, then return NULL. */
+ if(filename==NULL) { free(copy); return NULL; }
+
+ /* If a HDU wasn't given and the file is a FITS file, print a warning and
+ use the default "1". */
+ if(hdu==NULL && gal_fits_name_is_fits(filename))
+ error(EXIT_FAILURE, 0, "WARNING: '%s' is a FITS file, but no HDU "
+ "has been given (recall that a FITS file can contain "
+ "multiple HDUs). Please add a '-hdu-AAA' suffix to your "
+ "'"GAL_ARITHMETIC_OPSTR_LOADCOL_PREFIX"' operator to specify "
+ "the HDU (where 'AAA' is your HDU name or counter, counting "
+ "from zero); like this: '"GAL_ARITHMETIC_OPSTR_LOADCOL_PREFIX
+ "%s"GAL_ARITHMETIC_OPSTR_LOADCOL_FILE"%s"
+ GAL_ARITHMETIC_OPSTR_LOADCOL_HDU"AAA'",
+ filename, colid->v, filename);
+
+ /* Read the column from the table. */
+ out=gal_table_read(filename, hdu, NULL, colid, searchin, ignorecase,
+ numthreads, minmapsize, quietmmap, NULL);
+
+ /* Make sure that only a single column matched. */
+ if(out->next)
+ error(EXIT_FAILURE, 0, "%s: '%s' matches more than one column! "
+ "To load columns during arithmetic, it is important that "
+ "'load-col' returns only a single column", colid->v,
+ gal_fits_name_save_as_string(filename, hdu));
+
+ /* Clean up and return. */
+ gal_list_str_free(colid, 0);
+ free(copy);
+ return out;
+}
+
+
+
+
+
@@ -2636,6 +2921,10 @@ gal_arithmetic_set_operator(char *string, size_t
*num_operands)
{ op=GAL_ARITHMETIC_OP_MKNOISE_POISSON; *num_operands=2; }
else if (!strcmp(string, "mknoise-uniform"))
{ op=GAL_ARITHMETIC_OP_MKNOISE_UNIFORM; *num_operands=2; }
+ else if (!strcmp(string, "random-from-hist"))
+ { op=GAL_ARITHMETIC_OP_RANDOM_FROM_HIST; *num_operands=3; }
+ else if (!strcmp(string, "random-from-hist-raw"))
+ { op=GAL_ARITHMETIC_OP_RANDOM_FROM_HIST_RAW; *num_operands=3; }
/* The size operator */
else if (!strcmp(string, "size"))
@@ -2817,6 +3106,8 @@ gal_arithmetic_operator_string(int operator)
case GAL_ARITHMETIC_OP_MKNOISE_SIGMA: return "mknoise-sigma";
case GAL_ARITHMETIC_OP_MKNOISE_POISSON: return "mknoise-poisson";
case GAL_ARITHMETIC_OP_MKNOISE_UNIFORM: return "mknoise-uniform";
+ case GAL_ARITHMETIC_OP_RANDOM_FROM_HIST:return "random-from-hist";
+ case GAL_ARITHMETIC_OP_RANDOM_FROM_HIST_RAW:return "random-from-hist-raw";
case GAL_ARITHMETIC_OP_SIZE: return "size";
case GAL_ARITHMETIC_OP_STITCH: return "stitch";
@@ -2997,13 +3288,22 @@ gal_arithmetic(int operator, size_t numthreads, int
flags, ...)
out=arithmetic_bitwise_not(flags, d1);
break;
- /* Adding noise. */
+ /* Random steps. */
case GAL_ARITHMETIC_OP_MKNOISE_SIGMA:
case GAL_ARITHMETIC_OP_MKNOISE_POISSON:
case GAL_ARITHMETIC_OP_MKNOISE_UNIFORM:
+ case GAL_ARITHMETIC_OP_RANDOM_FROM_HIST:
+ case GAL_ARITHMETIC_OP_RANDOM_FROM_HIST_RAW:
d1 = va_arg(va, gal_data_t *);
d2 = va_arg(va, gal_data_t *);
- out=arithmetic_mknoise(operator, flags, d1, d2);
+ if(operator==GAL_ARITHMETIC_OP_RANDOM_FROM_HIST
+ || operator==GAL_ARITHMETIC_OP_RANDOM_FROM_HIST_RAW)
+ {
+ d3 = va_arg(va, gal_data_t *);
+ out=arithmetic_random_from_hist(operator, flags, d1, d2, d3);
+ }
+ else
+ out=arithmetic_mknoise(operator, flags, d1, d2);
break;
/* Size operator */
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index ac64e2f8..a82dd750 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -71,6 +71,15 @@ __BEGIN_C_DECLS /* From C++ preparations */
| GAL_ARITHMETIC_FLAG_FREE \
| GAL_ARITHMETIC_FLAG_NUMOK )
+/* Operator fixed strings. */
+#define GAL_ARITHMETIC_OPSTR_LOADCOL_PREFIX "load-col-"
+#define GAL_ARITHMETIC_OPSTR_LOADCOL_HDU "-hdu-"
+#define GAL_ARITHMETIC_OPSTR_LOADCOL_FILE "-from-"
+#define GAL_ARITHMETIC_OPSTR_LOADCOL_PREFIX_LEN
strlen(GAL_ARITHMETIC_OPSTR_LOADCOL_PREFIX)
+#define GAL_ARITHMETIC_OPSTR_LOADCOL_HDU_LEN
strlen(GAL_ARITHMETIC_OPSTR_LOADCOL_HDU)
+#define GAL_ARITHMETIC_OPSTR_LOADCOL_FILE_LEN
strlen(GAL_ARITHMETIC_OPSTR_LOADCOL_FILE)
+
+
/* Identifiers for each operator. */
enum gal_arithmetic_operators
@@ -165,6 +174,8 @@ enum gal_arithmetic_operators
GAL_ARITHMETIC_OP_MKNOISE_SIGMA,/* Fixed-sigma noise to every element. */
GAL_ARITHMETIC_OP_MKNOISE_POISSON,/* Poission noise on every element. */
GAL_ARITHMETIC_OP_MKNOISE_UNIFORM,/* Uniform noise on every element. */
+ GAL_ARITHMETIC_OP_RANDOM_FROM_HIST,/* Randoms from a histogram (uniform).*/
+ GAL_ARITHMETIC_OP_RANDOM_FROM_HIST_RAW,/* Randoms from a histogram (raw).*/
GAL_ARITHMETIC_OP_SIZE, /* Size of the dataset along an axis */
@@ -194,6 +205,10 @@ gal_arithmetic_operator_string(int operator);
int
gal_arithmetic_set_operator(char *string, size_t *num_operands);
+gal_data_t *
+gal_arithmetic_load_col(char *str, int searchin, int ignorecase,
+ size_t minmapsize, int quietmmap);
+
gal_data_t *
gal_arithmetic(int operator, size_t numthreads, int flags, ...);
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 7bb809c2: Library (arithmetic.h): loading external columns and random from hist,
Mohammad Akhlaghi <=