[gnuastro-commits] master a473254: Library (arithmetic): randomly select from uniform distribution

[Mohammad Akhlaghi]

branch: master
commit a473254cab5c411c171336235624f01c7aabe183
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Library (arithmetic): randomly select from uniform distribution
    
    Until now, the arithmetic library (and thus the arithmetic program, or
    column arithmetic in Table) only had two ways to select a random value: a
    Gaussian or Poisson distribution. But these distributions have a peak
    value, so the selected values won't be uniformly distributed (which is
    necessary in some scenarios).
    
    Witht his commit, a new 'mknoise-uniform' operator has been added to allow
    easy creation of random values that are uniformly distributed within the
    requested range.
---
 NEWS                      |  1 +
 doc/gnuastro.texi         | 19 +++++++++++++++++--
 lib/arithmetic.c          | 14 ++++++++++++--
 lib/gnuastro/arithmetic.h |  1 +
 4 files changed, 31 insertions(+), 4 deletions(-)

diff --git a/NEWS b/NEWS
index 4adb15e..648cb3f 100644
--- a/NEWS
+++ b/NEWS
@@ -65,6 +65,7 @@ See the end of the file for license conditions.
      - atanh: Inverse of hyperbolic tangent.
      - mknoise-sigma: Add Gaussian noise with the fixed sigma.
      - mknoise-poisson: Add Poisson noise with the given background.
+     - mknoise-uniform: Add uniform noise around existing value.
      - counts-to-mag: Convert counts to magnitudes with given zero point.
      - counts-to-jy: Convert counts to Janskys through a zero point based
           on AB magnitudes.
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 2a74662..2c7d6d0 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -12960,6 +12960,20 @@ This operator takes two arguments: the top/first 
popped operand is the backgroun
 Except for the noise-model, this operator is very similar to 
@code{mknoise-sigma} and the examples there apply here too.
 The main difference with @code{mknoise-sigma} is that in a Poisson 
distribution the scatter/sigma will depend on each element's value.
 
+@item mknoise-uniform
+Add uniform noise to each element of the input dataset.
+This operator takes two arguments: the top/first popped operand is the width 
of the interval, the second popped operand is the dataset that the noise should 
be added to (each element will be the center of the interval).
+The returned random values may happen to be the minimum interval value, but 
will never be the maximum.
+Except for the noise-model, this operator behaves very similar to 
@code{mknoise-sigma}, see the explanation there for more.
+
+For example with the command below, a random value will be selected between 10 
to 14 (centered on 12, which is the only input data element, with a total width 
of 4).
+
+@example
+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 size
 Size of the dataset along a given FITS/Fortran dimension (counting from 1).
 The desired dimension should be the first popped operand and the dataset must 
be the second popped operand.
@@ -26588,9 +26602,10 @@ for the rest it will be @code{GAL_TYPE_FLOAT32}.
 
 @deffn  Macro GAL_ARITHMETIC_OP_MKNOISE_SIGMA
 @deffnx Macro GAL_ARITHMETIC_OP_MKNOISE_POISSON
+@deffnx Macro GAL_ARITHMETIC_OP_MKNOISE_UNIFORM
 Add noise to the input dataset.
-Both operators take two arguments: the first is the input data set (can have 
any dimensionality or number of elements.
-The second argument is the noise specifier (a single element, of any type): 
for a fixed-sigma noise, it is the Gaussian standard deviation, for the Poisson 
noise, it is the background (see @ref{Photon counting noise}).
+These operators take two arguments: the first is the input data set (can have 
any dimensionality or number of elements.
+The second argument is the noise specifier (a single element, of any type): 
for a fixed-sigma noise, it is the Gaussian standard deviation, for the Poisson 
noise, it is the background (see @ref{Photon counting noise}) and for the 
uniform distribution it is the width of the interval around each element of the 
input dataset.
 
 By default, a separate random number generator seed will be used on each 
separate run of these operators.
 Therefore two identical runs on the same input will produce different results.
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 6ee39e6..3e49430 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -698,8 +698,9 @@ arithmetic_mknoise(int operator, int flags, gal_data_t *in, 
gal_data_t *arg)
   /* Make sure the noise identifier is positive. */
   if(arg_v<0)
     error(EXIT_FAILURE, 0, "the noise identifier (sigma for "
-          "'mknoise-sigma' or background for 'mknoise-poisson') must "
-          "be positive (it is %g)", arg_v);
+          "'mknoise-sigma', background for 'mknoise-poisson', or "
+          "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 &
@@ -733,6 +734,11 @@ arithmetic_mknoise(int operator, int flags, gal_data_t 
*in, gal_data_t *arg)
         *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 "
@@ -2125,6 +2131,8 @@ gal_arithmetic_set_operator(char *string, size_t 
*num_operands)
     { op=GAL_ARITHMETIC_OP_MKNOISE_SIGMA;     *num_operands=2; }
   else if (!strcmp(string, "mknoise-poisson"))
     { op=GAL_ARITHMETIC_OP_MKNOISE_POISSON;   *num_operands=2; }
+  else if (!strcmp(string, "mknoise-uniform"))
+    { op=GAL_ARITHMETIC_OP_MKNOISE_UNIFORM;   *num_operands=2; }
 
   /* The size operator */
   else if (!strcmp(string, "size"))
@@ -2284,6 +2292,7 @@ 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_SIZE:            return "size";
 
@@ -2448,6 +2457,7 @@ gal_arithmetic(int operator, size_t numthreads, int 
flags, ...)
     /* Adding noise. */
     case GAL_ARITHMETIC_OP_MKNOISE_SIGMA:
     case GAL_ARITHMETIC_OP_MKNOISE_POISSON:
+    case GAL_ARITHMETIC_OP_MKNOISE_UNIFORM:
       d1 = va_arg(va, gal_data_t *);
       d2 = va_arg(va, gal_data_t *);
       out=arithmetic_mknoise(operator, flags, d1, d2);
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index d8922f9..f6c2712 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -152,6 +152,7 @@ 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_SIZE,         /* Size of the dataset along an axis     */