[gnuastro-commits] master d0aa78e 005/125: Arithmetic operation on data structures in library

[Mohammad Akhlaghi]

branch: master
commit d0aa78e7d24c0eabf36f4ce3789432184fca45e9
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Arithmetic operation on data structures in library
    
    The operator functions in Arithmetic are mostly generally useful in any
    program for working on functions of any type. So I am now working on moving
    the main operator functions in Arithmetic into the library, while also
    designing a general purpose function for the basic jobs.
    
     - Arithmetic no longer (for the time being) needs an operands.[ch].
    
     - A new function (`gal_data_string_to_number') was added in `data.h' to
       read a number in string format into the proper type within a data
       structure.
    
     - A new function (`gal_data_alloc_number') will allocate a single element
       array and write the given number into it, this is now used by
       `gal_data_alloc_blank' and the new `gal_data_string_to_number'.
    
     - `minmapsize' is a new common argument/parameter, to keep the minimum
       number of bytes, to use mmap (in the hdd/ssd) instead of malloc (in the
       RAM). It is now incorporated into the functionst that read a FITS image.
    
     - Arithmetic's add_operand or pop_operand functions don't need a separate
       number or array argument, there is only a single data-structure.
    
     - `gal_data_alloc', now also accepts an input array. When present no array
       is allocated, and the given array will be used in the data structure.
    
     - The new `gal_data_copy' will copy one data structure into an identical
       data structure.
    
     - The new `gal_data_out_type' will return the type that should be used
       from two possibly different types.
    
     - The new `gal_data_to_same_type' will convert the two given data types to
       a given type and possibly free the input arrays.
    
     - The new `gal_data_arithmetic' function is now in charge of doing basic
       arithmetic operations on a variable number of input datatypes.
---
 bin/arithmetic/Makefile.am                         |   2 +-
 bin/arithmetic/arithmetic.c                        |  95 ++--
 bin/arithmetic/astarithmetic.conf                  |   1 +
 bin/arithmetic/main.h                              |  41 +-
 bin/arithmetic/operands.c                          |  34 +-
 bin/arithmetic/operands.h                          |   8 +-
 bin/arithmetic/ui.c                                |   4 +
 lib/Makefile.am                                    |  14 +-
 lib/commonargs.h                                   |  15 +-
 lib/commonparams.h                                 |   2 +
 lib/configfiles.h                                  |   8 +
 bin/arithmetic/operands.h => lib/data-arithmetic.c |  24 +-
 lib/data-arithmetic.h                              | 152 +++++++
 lib/{changetype.c => data-changetype.c}            |   2 +-
 lib/{changetype.h => data-changetype.h}            |   0
 lib/data.c                                         | 486 +++++++++++++++++----
 lib/fits.c                                         |  55 +--
 lib/gnuastro/data.h                                |  41 +-
 lib/gnuastro/fits.h                                |   7 +-
 19 files changed, 766 insertions(+), 225 deletions(-)

diff --git a/bin/arithmetic/Makefile.am b/bin/arithmetic/Makefile.am
index cf714b6..657a834 100644
--- a/bin/arithmetic/Makefile.am
+++ b/bin/arithmetic/Makefile.am
@@ -29,7 +29,7 @@ AM_CPPFLAGS = -I\$(top_srcdir)/bootstrapped/lib
 bin_PROGRAMS = astarithmetic
 
 astarithmetic_SOURCES = main.c main.h cite.h ui.c ui.h args.h  \
-arithmetic.c arithmetic.h operands.c operands.h operators.c operators.h
+arithmetic.c arithmetic.h operands.c operands.h
 
 astarithmetic_LDADD = $(top_builddir)/bootstrapped/lib/libgnu.la       \
 -lgnuastro
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index efb10aa..dfcef4a 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -29,6 +29,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <string.h>
 #include <stdlib.h>
 
+#include <gnuastro/data.h>
 #include <gnuastro/fits.h>
 #include <gnuastro/array.h>
 #include <gnuastro/statistics.h>
@@ -67,30 +68,48 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 void
 reversepolish(struct imgarithparams *p)
 {
-  double number;
-  gal_data_t *ddata;
   struct gal_linkedlist_stll *token;
+  gal_data_t *d1=NULL, *d2=NULL, *d3=NULL;
+  unsigned char flags = ( GAL_DATA_ARITH_INPLACE | GAL_DATA_ARITH_FREE
+                          | GAL_DATA_ARITH_NUMOK );
+
 
   /* Prepare the processing: */
   p->operands=NULL;
   p->addcounter=p->popcounter=0;
 
+
   /* Go over each input token and do the work. */
   for(token=p->tokens;token!=NULL;token=token->next)
     {
-      /* If we have a name or number, then add it to the operands
-         linked list. Otherwise, pull out two members and do the
-         specified operation on them. */
+      /* If we have a name or number, then add it to the operands linked
+         list. Otherwise, pull out two members and do the specified
+         operation on them. */
       if(gal_fits_name_is_fits(token->v))
-        add_operand(p, token->v, NOOPTNUMBER, NOOPTDATA);
-      else if(gal_checkset_str_is_double(token->v, &number))
-        add_operand(p, NOOPTFILENAME, number, NOOPTDATA);
+        add_operand(p, token->v, NULL);
+      else if( (d1=gal_data_string_to_number(token->v)) )
+        add_operand(p, NULL, d1);
       else
         {
-          if     (!strcmp(token->v, "+"))         sum(p);
-          else if(!strcmp(token->v, "-"))         subtract(p);
-          else if(!strcmp(token->v, "*"))         multiply(p);
-          else if(!strcmp(token->v, "/"))         divide(p);
+          /* Note that the first popped operand is the right/second given
+             operand on the command-line. */
+          if( !strcmp(token->v, "+") || !strcmp(token->v, "-")
+              || !strcmp(token->v, "*") || !strcmp(token->v, "/") )
+            {
+              d2=pop_operand(p, token->v);
+              d1=pop_operand(p, token->v);
+            }
+          else
+            error(EXIT_FAILURE, 0, "the argument \"%s\" could not be "
+                  "interpretted as a FITS file, number, or operator",
+                  token->v);
+          add_operand(p, NULL,
+                      gal_data_arithmetic(token->v, flags, d1, d2, d3));
+        }
+    }
+
+#if 0
+
           else if(!strcmp(token->v, "abs"))       takeabs(p);
           else if(!strcmp(token->v, "pow"))       topower(p, NULL);
           else if(!strcmp(token->v, "sqrt"))      takesqrt(p);
@@ -113,46 +132,36 @@ reversepolish(struct imgarithparams *p)
           else if(!strcmp(token->v, "not"))       notfunc(p);
           else if(!strcmp(token->v, "isblank"))   opisblank(p);
           else if(!strcmp(token->v, "where"))     where(p);
-          else
-            error(EXIT_FAILURE, 0, "the argument \"%s\" could not be "
-                  "interpretted as an operator", token->v);
-        }
-    }
+#endif
+
 
   /* If there is more than one node in the operands stack then the user has
      given too many operands which is an error. */
   if(p->operands->next!=NULL)
     error(EXIT_FAILURE, 0, "too many operands");
 
-  /* If the remaining operand is an array then save the array as a
-     FITS image, if not, simply print the floating point number. */
-  if(p->operands->data)
-    {
-      /* Internally, all arrays were double type. But the user could set
-         the the output type using the type option. So if the user has
-         asked for anything other than a double, we will have to convert
-         the arrays.*/
-      if(p->outtype==GAL_DATA_TYPE_DOUBLE)
-        ddata=p->operands->data;
-      else
-        {
-          ddata=gal_data_copy_to_new_type(p->operands->data, p->outtype);
-          gal_data_free(p->operands->data);
-        }
 
-      /* Copy the WCS structure into the final dataset and write it into
-         the output file. */
-      ddata->wcs=p->refdata.wcs;
-      gal_fits_write_img(ddata, p->cp.output, "Arithmetic", NULL,
-                         SPACK_STRING);
+  /* Copy the WCS structure into the final dataset. */
+  p->operands->data->wcs=p->refdata.wcs;
 
-      /* Clean up: */
-      gal_data_free(ddata);
-      free(p->refdata.dsize);
-      wcsfree(p->refdata.wcs);
-    }
+
+  /* Set the output type. */
+  if(p->outtype==p->operands->data->type)
+    d1=p->operands->data;
   else
-    printf("%g\n", p->operands->number);
+    {
+      d1=gal_data_copy_to_new_type(p->operands->data, p->outtype);
+      gal_data_free(p->operands->data);
+    }
+  gal_fits_write_img(p->operands->data, p->cp.output, "Arithmetic", NULL,
+                     SPACK_STRING);
+
+
+  /* Clean up: */
+  gal_data_free(d1);
+  free(p->refdata.dsize);
+  wcsfree(p->refdata.wcs);
+
 
   /* Clean up. Note that the tokens were taken from the command-line
      arguments, so the string within each token linked list must not be
diff --git a/bin/arithmetic/astarithmetic.conf 
b/bin/arithmetic/astarithmetic.conf
index bdd347f..b877df2 100644
--- a/bin/arithmetic/astarithmetic.conf
+++ b/bin/arithmetic/astarithmetic.conf
@@ -28,3 +28,4 @@
 
 # Output:
  type       float
+ minmapsize    1000000000
\ No newline at end of file
diff --git a/bin/arithmetic/main.h b/bin/arithmetic/main.h
index 3aa5000..2e4409a 100644
--- a/bin/arithmetic/main.h
+++ b/bin/arithmetic/main.h
@@ -39,37 +39,35 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #define NOTCOMMONHDU   1
 
 
+
+
+
 /* Constants: */
 #define NEGDASHREPLACE 11  /* A vertical tab (ASCII=11) for negative dash */
-#define NOOPTDATA      NULL
-#define NOOPTNUMBER    NAN
-#define NOOPTFILENAME  ""
 
 
-/* The operands can be in three ways:
 
-      1. A file name which is not read yet. When inactive, NOOPTFILENAME.
-      2. A number which is not used yet. When inactive, NOOPTNUMBER.
-      3. A dataset (operator output). When inactive, NOOPTDATA.
 
-   In every node of the operand linked list, only one of these should be
-   active. The other two should be inactive. Otherwise it will be a bug and
-   will cause problems. All the operands operate on this premise.
-*/
+
+/* In every node of the operand linked list, only one of the `filename' or
+   `data' should be non-NULL. Otherwise it will be a bug and will cause
+   problems. All the operands operate on this premise. */
 struct operand
 {
-  char       *filename;         /* !=NULL if the operand is a filename. */
-  char            *hdu;         /* !=NULL if the operand is a filename. */
-  double        number;         /* Value used if not a dataset.         */
-  gal_data_t     *data;         /* !=NULL if the operand is a dataset.  */
-  struct operand *next;         /* Pointer to next operand.             */
+  char       *filename;    /* !=NULL if the operand is a filename. */
+  char            *hdu;    /* !=NULL if the operand is a filename. */
+  gal_data_t     *data;    /* !=NULL if the operand is a dataset.  */
+  struct operand *next;    /* Pointer to next operand.             */
 };
 
 
+
+
+
 struct uiparams
 {
-  char        *maskname;  /* Name of mask image.                        */
-  char            *mhdu;  /* Mask image HDU.                            */
+  char        *maskname;   /* Name of mask image.                   */
+  char            *mhdu;   /* Mask image HDU.                       */
 
   int           typeset;
   int       masknameset;
@@ -78,6 +76,9 @@ struct uiparams
 };
 
 
+
+
+
 struct imgarithparams
 {
   /* Other structures: */
@@ -102,4 +103,8 @@ struct imgarithparams
   time_t           rawtime;  /* Starting time of the program.           */
 };
 
+
+
+
+
 #endif
diff --git a/bin/arithmetic/operands.c b/bin/arithmetic/operands.c
index 369917c..1153941 100644
--- a/bin/arithmetic/operands.c
+++ b/bin/arithmetic/operands.c
@@ -56,8 +56,7 @@ num_operands(struct imgarithparams *p)
 
 
 void
-add_operand(struct imgarithparams *p, char *filename, double number,
-            gal_data_t *data)
+add_operand(struct imgarithparams *p, char *filename, gal_data_t *data)
 {
   struct operand *newnode;
 
@@ -69,7 +68,6 @@ add_operand(struct imgarithparams *p, char *filename, double 
number,
 
   /* Fill in the values. */
   newnode->data=data;
-  newnode->number=number;
   newnode->filename=filename;
 
   if(filename != NULL && gal_fits_name_is_fits(filename))
@@ -90,11 +88,11 @@ add_operand(struct imgarithparams *p, char *filename, 
double number,
 
 
 
-void
-pop_operand(struct imgarithparams *p, double *number, gal_data_t **data,
-            char *operator)
+gal_data_t *
+pop_operand(struct imgarithparams *p, char *operator)
 {
   size_t i;
+  gal_data_t *data;
   struct uiparams *up=&p->up;
   struct operand *operands=p->operands;
   char *maskname, *mhdu, *filename, *hdu;
@@ -102,13 +100,13 @@ pop_operand(struct imgarithparams *p, double *number, 
gal_data_t **data,
   /* If the operand linked list has finished, then give an error and
      exit. */
   if(operands==NULL)
-    error(EXIT_FAILURE, 0, "not enough operands for the \"%s\" operator",
-          operator);
+    error(EXIT_FAILURE, 0, "not enough operands for the \"%s\" "
+          "operator", operator);
 
 
   /* Set the dataset. If filename is present then read the file
      and fill in the array, if not then just set the array. */
-  if(strlen(operands->filename))
+  if(operands->filename)
     {
       hdu=operands->hdu;
       filename=operands->filename;
@@ -124,8 +122,8 @@ pop_operand(struct imgarithparams *p, double *number, 
gal_data_t **data,
                           &p->refdata.wcs);
 
       /* Read the input image as a double type array. */
-      *data=gal_fits_read_to_type(filename, maskname, hdu, mhdu,
-                                  GAL_DATA_TYPE_DOUBLE);
+      data=gal_fits_read_img_hdu(filename, hdu, maskname, mhdu,
+                                 p->cp.minmapsize);
 
       /* When the reference data structure's dimensionality is non-zero, it
          means that this is not the first image read. So, write its basic
@@ -133,7 +131,7 @@ pop_operand(struct imgarithparams *p, double *number, 
gal_data_t **data,
          checks. */
       if(p->refdata.ndim)
         {
-          if(gal_data_dsize_is_different(&p->refdata, *data))
+          if(gal_data_dsize_is_different(&p->refdata, data))
             error(EXIT_FAILURE, 0, "%s (hdu=%s): has a different size "
                   "compared to previous images. All the images must be "
                   "the same size in order for Arithmetic to work",
@@ -142,7 +140,7 @@ pop_operand(struct imgarithparams *p, double *number, 
gal_data_t **data,
       else
         {
           /* Set the dimensionality. */
-          p->refdata.ndim=(*data)->ndim;
+          p->refdata.ndim=(data)->ndim;
 
           /* Allocate the dsize array. */
           errno=0;
@@ -153,7 +151,7 @@ pop_operand(struct imgarithparams *p, double *number, 
gal_data_t **data,
 
           /* Write the values into it. */
           for(i=0;i<p->refdata.ndim;++i)
-            p->refdata.dsize[i]=(*data)->dsize[i];
+            p->refdata.dsize[i]=data->dsize[i];
         }
 
       /* Free the HDU string: */
@@ -166,12 +164,10 @@ pop_operand(struct imgarithparams *p, double *number, 
gal_data_t **data,
       if(p->cp.verb) printf("%s is read.\n", filename);
     }
   else
-    *data=operands->data;
-
-  /* Set the number: */
-  *number=operands->number;
+    data=operands->data;
 
-  /* Remove this node from the queue. */
+  /* Remove this node from the queue, return the data structure. */
   p->operands=operands->next;
   free(operands);
+  return data;
 }
diff --git a/bin/arithmetic/operands.h b/bin/arithmetic/operands.h
index 4f3798d..6b6f6e0 100644
--- a/bin/arithmetic/operands.h
+++ b/bin/arithmetic/operands.h
@@ -27,12 +27,10 @@ size_t
 num_operands(struct imgarithparams *p);
 
 void
-add_operand(struct imgarithparams *p, char *filename, double number,
-            gal_data_t *data);
+add_operand(struct imgarithparams *p, char *filename, gal_data_t *data);
 
-void
-pop_operand(struct imgarithparams *p, double *number, gal_data_t **data,
-            char *operator);
+gal_data_t *
+pop_operand(struct imgarithparams *p, char *operator);
 
 
 #endif
diff --git a/bin/arithmetic/ui.c b/bin/arithmetic/ui.c
index 556b910..4a18be8 100644
--- a/bin/arithmetic/ui.c
+++ b/bin/arithmetic/ui.c
@@ -199,6 +199,10 @@ checkifset(struct imgarithparams *p)
   char comment[100];
   size_t numhdus=gal_linkedlist_num_in_stll(p->hdus);
 
+  /* Operating modes: */
+  if(p->cp.minmapsizeset==0)
+    GAL_CONFIGFILES_REPORT_NOTSET("minmapsize");
+
   /* Output parameters: */
   if(p->up.typeset==0)
     GAL_CONFIGFILES_REPORT_NOTSET("type");
diff --git a/lib/Makefile.am b/lib/Makefile.am
index e1704f5..49d36a6 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -31,10 +31,10 @@ libgnuastro_la_LDFLAGS = -version-info $(GAL_LT_VERSION)
 
 
 # Specify the library .c files
-libgnuastro_la_SOURCES = array.c box.c changetype.c checkset.c           \
-  configfiles.c data.c fits.c git.c linkedlist.c mesh.c mode.c polygon.c \
-  qsort.c spatialconvolve.c statistics.c threads.c timing.c txtarray.c   \
-  wcs.c
+libgnuastro_la_SOURCES = array.c box.c checkset.c configfiles.c data.c \
+  data-arithmetic.c data-changetype.c fits.c git.c linkedlist.c mesh.c \
+  mode.c polygon.c qsort.c spatialconvolve.c statistics.c threads.c    \
+  timing.c txtarray.c wcs.c
 
 
 # Installed headers, note that we are not blindly including all `.h' files
@@ -54,9 +54,9 @@ pkginclude_HEADERS = gnuastro/config.h $(headersdir)/array.h  
          \
 # and if they are not explicitly mentioned somewhere in the Makefile, they
 # will not distributed, so we need to explicitly tell Automake to
 # distribute them here.
-EXTRA_DIST = gnuastro.pc.in changetype.h config.h.in checkset.h         \
-  commonargs.h commonparams.h configfiles.h fixedstringmacros.h mode.h  \
-  neighbors.h timing.h $(headersdir)/README
+EXTRA_DIST = gnuastro.pc.in data-arithmetic.h data-changetype.h         \
+  config.h.in checkset.h commonargs.h commonparams.h configfiles.h      \
+  fixedstringmacros.h mode.h neighbors.h timing.h $(headersdir)/README
 
 
 # Definitions for Gnuastro's the pkg-config file (inspired from GSL's
diff --git a/lib/commonargs.h b/lib/commonargs.h
index 99425e2..83caf0d 100644
--- a/lib/commonargs.h
+++ b/lib/commonargs.h
@@ -51,7 +51,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
    a b c d e f g i j k l m n p r s t u v w x y z
    A B C E F G H I J L M O Q R T U W X Y Z
 
-   Used numbers <= 1003
+   Used numbers <= 1004
 
    You can use this above list to set short options for the different
    utilities.
@@ -131,6 +131,14 @@ static struct argp_option gal_commonargs_options[] =
       "No log file for programs which make one.",
       -1
     },
+    {
+      "minmapsize",
+      1004,
+      0,
+      0,
+      "Minimum no. bytes to map arrays to hdd/ssd.",
+      -1
+    },
 
 
 
@@ -248,6 +256,11 @@ gal_checkset_commonargs_cparse_opt(int key, char *arg, 
struct argp_state *state)
       cp->nolog=1;
       cp->nologset=1;
       break;
+    case 1004:
+      gal_checkset_sizet_l_zero(arg, &cp->minmapsize, "minmapsize", key,
+                                cp->spack, NULL, 0);
+      cp->minmapsizeset=1;
+      break;
 
     /* Input/output: */
     case 'h':
diff --git a/lib/commonparams.h b/lib/commonparams.h
index d426e00..70de2f5 100644
--- a/lib/commonparams.h
+++ b/lib/commonparams.h
@@ -47,6 +47,7 @@ struct gal_commonparams
   size_t  numthreads;  /* Number of threads to use.                 */
   int    onlydirconf;  /* Only check current directory conf. file.  */
   char  *onlyversion;  /* The string of the requested version.      */
+  size_t  minmapsize;  /* The minimum bytes necessary to use mmap.  */
   int          nolog;  /* ==1: do not make a log file.              */
 
   /* Check: */
@@ -57,6 +58,7 @@ struct gal_commonparams
   int         hduset;  /* If the input image extension is set.      */
   int      outputset;  /* If the output is set.                     */
   int       nologset;  /* If nolog is set.                          */
+  int  minmapsizeset;  /* If minmapsize is set.                     */
   int  dontdeleteset;  /* If the --dontdelete option was called.    */
   int removedirinfoset;  /* If --keepinputdir was called.           */
 };
diff --git a/lib/configfiles.h b/lib/configfiles.h
index 51ead50..7055d92 100644
--- a/lib/configfiles.h
+++ b/lib/configfiles.h
@@ -226,6 +226,13 @@ __BEGIN_C_DECLS  /* From C++ preparations */
                                      SPACK, filename, lineno);          \
         cp->nologset=1;                                                 \
       }                                                                 \
+    else if(strcmp(name, "minmapsize")==0)                              \
+      {                                                                 \
+        if(cp->minmapsizeset) continue;                                 \
+        gal_checkset_sizet_l_zero(value, &cp->minmapsize, name, key,    \
+                                  SPACK, filename, lineno);             \
+        cp->minmapsizeset=1;                                            \
+      }                                                                 \
                                                                         \
                                                                         \
     else if(strcmp(name, "dontdelete")==0)                              \
@@ -250,6 +257,7 @@ __BEGIN_C_DECLS  /* From C++ preparations */
 
 
 
+
 /* Write common options: */
 #define GAL_CONFIGFILES_PRINT_COMMONOPTIONS {                           \
     if(cp->quietset)                                                    \
diff --git a/bin/arithmetic/operands.h b/lib/data-arithmetic.c
similarity index 63%
copy from bin/arithmetic/operands.h
copy to lib/data-arithmetic.c
index 4f3798d..8c08a71 100644
--- a/bin/arithmetic/operands.h
+++ b/lib/data-arithmetic.c
@@ -1,11 +1,11 @@
 /*********************************************************************
-Arithmetic - Do arithmetic operations on images.
-Arithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
+Arithmetic operations on data structures.
+This is part of GNU Astronomy Utilities (Gnuastro) package.
 
 Original author:
      Mohammad Akhlaghi <address@hidden>
 Contributing author(s):
-Copyright (C) 2015, Free Software Foundation, Inc.
+Copyright (C) 2016, 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
@@ -20,19 +20,9 @@ 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 OPERANDS_H
-#define OPERANDS_H
+#include <config.h>
 
-size_t
-num_operands(struct imgarithparams *p);
+#include <errno.h>
+#include <error.h>
 
-void
-add_operand(struct imgarithparams *p, char *filename, double number,
-            gal_data_t *data);
-
-void
-pop_operand(struct imgarithparams *p, double *number, gal_data_t **data,
-            char *operator);
-
-
-#endif
+#include <gnuastro/data.h>
diff --git a/lib/data-arithmetic.h b/lib/data-arithmetic.h
new file mode 100644
index 0000000..c3c78c4
--- /dev/null
+++ b/lib/data-arithmetic.h
@@ -0,0 +1,152 @@
+/*********************************************************************
+Arithmetic operations on data structures.
+This is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+     Mohammad Akhlaghi <address@hidden>
+Contributing author(s):
+Copyright (C) 2016, 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_ARITHMETIC_H__
+#define __GAL_ARITHMETIC_H__
+
+
+
+
+
+#define BINARY_DEFINITIONS                                                \
+  unsigned char  *luc = l->array, *ruc = r->array, *ouc = o->array, *ucf; \
+  char           *lc  = l->array, *rc  = r->array, *oc  = o->array,  *cf; \
+  unsigned short *lus = l->array, *rus = r->array, *ous = o->array, *usf; \
+  short          *ls  = l->array, *rs  = r->array, *os  = o->array,  *sf; \
+  unsigned int   *lui = l->array, *rui = r->array, *oui = o->array, *uif; \
+  int            *li  = l->array, *ri  = r->array, *oi  = o->array, *iif; \
+  unsigned long  *lul = l->array, *rul = r->array, *oul = o->array, *ulf; \
+  long           *ll  = l->array, *rl  = r->array, *ol  = o->array,  *lf; \
+  LONGLONG       *lL  = l->array, *rL  = r->array, *oL  = o->array,  *Lf; \
+  float          *lff = l->array, *rf  = r->array, *of  = o->array,  *ff; \
+  double         *ld  = l->array, *rd  = r->array, *od  = o->array,  *df;
+
+
+
+
+#define OP_FUNC(OP,L, R)
+#define BINARY_OPERATOR_FOR_TYPE(Tf, oT, lT, rT, OP)                    \
+  Tf = oT + o->size;                                                    \
+  if(l->size==r->size) do *oT = *lT++ OP *rT++; while(++oT<Tf);         \
+  else if(l->size==1)  do *oT = *lT   OP *rT++; while(++oT<Tf);         \
+  else                 do *oT = *lT++ OP *rT;   while(++oT<Tf);         \
+  break;
+
+
+
+
+
+/* Preparations for binary operators.
+
+   Note that after `gal_data_to_same_type', we are only concerned with the
+   `l' and `s' pointers. */
+#define BINARY_INTERNAL(OP, OUT_TYPE)                                   \
+                                                                        \
+  gal_data_t *l, *r, *left, *right;                                     \
+                                                                        \
+                                                                        \
+  /* Read the variable arguments. */                                    \
+  left = va_arg(va, gal_data_t *);                                      \
+  right = va_arg(va, gal_data_t *);                                     \
+                                                                        \
+                                                                        \
+  /* Simple sanity check, then choose the common type. */               \
+  if( (flags & GAL_DATA_ARITH_NUMOK) && (left->size==1 || right->size==1 ) ) \
+    type=0;      /* Everything is ok, just a place-holder. */           \
+  else if (gal_data_dsize_is_different(left, right))                    \
+    error(EXIT_FAILURE, 0, "The datasets don't have the same "          \
+          "dimension/size");                                            \
+                                                                        \
+                                                                        \
+  /* Set the two datasets to the same type if they were different. */   \
+  type=gal_data_out_type(left, right);                                  \
+  gal_data_to_same_type(left, right, &l, &r, type,                      \
+                        flags & GAL_DATA_ARITH_FREE );                  \
+                                                                        \
+                                                                        \
+  /* Output can point to any one of the two arrays. */                  \
+  if(flags & GAL_DATA_ARITH_INPLACE)                                    \
+    o = l->size>1 ? l : r;                                              \
+  else                                                                  \
+    o = gal_data_alloc(NULL, type,                                      \
+                       l->size>1 ? l->ndim  : r->ndim,                  \
+                       l->size>1 ? l->dsize : r->dsize,                 \
+                       0, l->mmapped || r->mmapped);                    \
+                                                                        \
+                                                                        \
+  /* In block to allow definitions. */                                  \
+  {                                                                     \
+    BINARY_DEFINITIONS;                                                 \
+    switch(l->type)                                                     \
+      {                                                                 \
+      case GAL_DATA_TYPE_UCHAR:                                         \
+        BINARY_OPERATOR_FOR_TYPE(ucf, ouc, luc, ruc, OP);               \
+                                                                        \
+      case GAL_DATA_TYPE_CHAR:                                          \
+        BINARY_OPERATOR_FOR_TYPE(cf,  oc,  lc,  rc,  OP);               \
+                                                                        \
+      case GAL_DATA_TYPE_USHORT:                                        \
+        BINARY_OPERATOR_FOR_TYPE(usf, ous, lus, rus, OP);               \
+                                                                        \
+      case GAL_DATA_TYPE_SHORT:                                         \
+        BINARY_OPERATOR_FOR_TYPE(sf,  os,  ls,  rs,  OP);               \
+                                                                        \
+      case GAL_DATA_TYPE_UINT:                                          \
+        BINARY_OPERATOR_FOR_TYPE(uif, oui, lui, rui, OP);               \
+                                                                        \
+      case GAL_DATA_TYPE_INT:                                           \
+        BINARY_OPERATOR_FOR_TYPE(iif, oi,  li,  ri,  OP);               \
+                                                                        \
+      case GAL_DATA_TYPE_ULONG:                                         \
+        BINARY_OPERATOR_FOR_TYPE(ulf, oul, lul, rul, OP);               \
+                                                                        \
+      case GAL_DATA_TYPE_LONG:                                          \
+        BINARY_OPERATOR_FOR_TYPE(lf,  ol,  ll,  rl,  OP);               \
+                                                                        \
+      case GAL_DATA_TYPE_LONGLONG:                                      \
+        BINARY_OPERATOR_FOR_TYPE(Lf,  oL,  lL,  rL,  OP);               \
+                                                                        \
+      case GAL_DATA_TYPE_FLOAT:                                         \
+        BINARY_OPERATOR_FOR_TYPE(ff, of,  lff,  rf,  OP);               \
+                                                                        \
+      case GAL_DATA_TYPE_DOUBLE:                                        \
+        BINARY_OPERATOR_FOR_TYPE(df,  od,  ld,  rd,  OP);               \
+                                                                        \
+      default:                                                          \
+        error(EXIT_FAILURE, 0, "type %d not recognized in "             \
+              "data-arithmetic", type);                                 \
+      }                                                                 \
+  }                                                                     \
+                                                                        \
+    /* Clean up. */                                                     \
+    if(flags & GAL_DATA_ARITH_FREE)                                     \
+    {                                                                   \
+      if(o==l) gal_data_free(r);                                        \
+      else if(o==r) gal_data_free(l);                                   \
+      else {gal_data_free(l); gal_data_free(r);}                        \
+    }
+
+
+
+
+
+#endif
diff --git a/lib/changetype.c b/lib/data-changetype.c
similarity index 99%
rename from lib/changetype.c
rename to lib/data-changetype.c
index 21f3220..81ef3de 100644
--- a/lib/changetype.c
+++ b/lib/data-changetype.c
@@ -27,7 +27,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gnuastro/data.h>
 
-#include <changetype.h>
+#include <data-changetype.h>
 
 
 
diff --git a/lib/changetype.h b/lib/data-changetype.h
similarity index 100%
rename from lib/changetype.h
rename to lib/data-changetype.h
diff --git a/lib/data.c b/lib/data.c
index 68a1230..61b130c 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -25,6 +25,9 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <errno.h>
 #include <error.h>
 #include <fcntl.h>
+#include <float.h>
+#include <ctype.h>
+#include <stdarg.h>
 #include <unistd.h>
 #include <stdlib.h>
 #include <string.h>
@@ -34,7 +37,8 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/data.h>
 
 #include <checkset.h>
-#include <changetype.h>
+#include <data-changetype.h>
+#include <data-arithmetic.h>
 
 
 
@@ -187,6 +191,103 @@ gal_data_calloc_array(int type, size_t size)
 
 
 
+/* Allocate space for one blank value of the given type and put the value
+   in it. */
+void *
+gal_data_alloc_number(int type, void *number)
+{
+  /* Define the pointers. */
+  void *allocated;
+
+  /* Allocate the space for the blank value: */
+  allocated=gal_data_malloc_array(type, 1);
+
+  /* Put the blank value into it. */
+  errno=0;
+  switch(type)
+    {
+    case GAL_DATA_TYPE_BIT:
+      error(EXIT_FAILURE, 0, "Currently Gnuastro doesn't support blank "
+            "values for `GAL_DATA_TYPE_BIT', please get in touch with "
+            "us to see how we can implement it.");
+
+    case GAL_DATA_TYPE_UCHAR:
+      *(unsigned char *)allocated=*(unsigned char *)number;
+      break;
+
+      /* CFITSIO says "int for keywords, char for table columns". Here we
+         are only assuming table columns. So in practice this also applies
+         to TSBYTE.*/
+    case GAL_DATA_TYPE_CHAR: case GAL_DATA_TYPE_LOGICAL:
+      *(char *)allocated=*(char *)number;
+      break;
+
+    case GAL_DATA_TYPE_STRING:
+      *(unsigned char **)allocated=*(unsigned char **)number;
+      break;
+
+    case GAL_DATA_TYPE_USHORT:
+      *(unsigned short *)allocated=*(unsigned short *)number;
+      break;
+
+    case GAL_DATA_TYPE_SHORT:
+      *(short *)allocated=*(short *)number;
+      break;
+
+    case GAL_DATA_TYPE_UINT:
+      *(unsigned int *)allocated=*(unsigned int *)number;
+      break;
+
+    case GAL_DATA_TYPE_INT:
+      *(int *)allocated=*(int *)number;
+      break;
+
+    case GAL_DATA_TYPE_ULONG:
+      *(unsigned long *)allocated=*(unsigned long *)number;
+      break;
+
+    case GAL_DATA_TYPE_LONG:
+      *(long *)allocated=*(long *)number;
+      break;
+
+    case GAL_DATA_TYPE_LONGLONG:
+      *(LONGLONG *)allocated=*(LONGLONG *)number;
+      break;
+
+    case GAL_DATA_TYPE_FLOAT:
+      *(float *)allocated=*(float *)number;
+      break;
+
+    case GAL_DATA_TYPE_DOUBLE:
+      *(double *)allocated=*(double *)number;
+      break;
+
+    case GAL_DATA_TYPE_COMPLEX:
+      GSL_COMPLEX_P_REAL(((gsl_complex_float *)allocated)) =
+        GSL_COMPLEX_P_REAL(((gsl_complex_float *)number));
+      GSL_COMPLEX_P_IMAG(((gsl_complex_float *)allocated)) =
+        GSL_COMPLEX_P_IMAG(((gsl_complex_float *)number));
+      break;
+
+    case GAL_DATA_TYPE_DCOMPLEX:
+      GSL_COMPLEX_P_REAL(((gsl_complex *)allocated)) =
+        GSL_COMPLEX_P_REAL(((gsl_complex *)number));
+      GSL_COMPLEX_P_IMAG(((gsl_complex *)allocated)) =
+        GSL_COMPLEX_P_IMAG(((gsl_complex *)number));
+      break;
+
+    default:
+      error(EXIT_FAILURE, 0, "type value of %d not recognized in "
+            "`gal_data_alloc_number'", type);
+    }
+
+  return allocated;
+}
+
+
+
+
+
 void
 gal_data_mmap(gal_data_t *data)
 {
@@ -240,12 +341,17 @@ gal_data_mmap(gal_data_t *data)
 
 
 
+/* Allocate a data structure based on the given parameters. If you want to
+   force the array into the hdd/ssd (mmap it), then set minmapsize=-1
+   (largest possible size_t value), in this way, no file will be larger. */
 gal_data_t *
-gal_data_alloc(int type, size_t ndim, long *dsize, int clear, int map)
+gal_data_alloc(void *array, int type, size_t ndim, long *dsize,
+               int clear, size_t minmapsize)
 {
   size_t i;
   gal_data_t *out;
 
+
   /* Allocate the space for the actual structure. */
   errno=0;
   out=malloc(sizeof *out);
@@ -258,6 +364,7 @@ gal_data_alloc(int type, size_t ndim, long *dsize, int 
clear, int map)
   out->ndim=ndim;
   out->type=type;
 
+
   /* Initialize the other values */
   out->anyblank=0;
   out->wcs=NULL;
@@ -271,6 +378,7 @@ gal_data_alloc(int type, size_t ndim, long *dsize, int 
clear, int map)
           ndim*sizeof *out->dsize);
 
 
+
   /* Fill in the dsize values: */
   out->size=1;
   for(i=0;i<ndim;++i)
@@ -285,22 +393,29 @@ gal_data_alloc(int type, size_t ndim, long *dsize, int 
clear, int map)
       out->size *= ( out->dsize[i] = dsize[i] );
     }
 
+
   /* Allocate space for the array, clear it if necessary: */
-  if(map)
-    gal_data_mmap(out);
+  if(array)
+    out->array=array;
   else
     {
-      /* Allocate the space for the array. */
-      if(clear)
-        out->array = gal_data_calloc_array(out->type, out->size);
+      if( gal_data_sizeof(type)*out->size  > minmapsize )
+        gal_data_mmap(out);
       else
-        out->array = gal_data_malloc_array(out->type, out->size);
-
-      /* Set the values. */
-      out->mmapped=0;
-      out->mmapname=NULL;
+        {
+          /* Allocate the space for the array. */
+          if(clear)
+            out->array = gal_data_calloc_array(out->type, out->size);
+          else
+            out->array = gal_data_malloc_array(out->type, out->size);
+
+          /* Set the values. */
+          out->mmapped=0;
+          out->mmapname=NULL;
+        }
     }
 
+
   /* Return the final structure. */
   return out;
 }
@@ -356,34 +471,26 @@ gal_data_free(gal_data_t *data)
 /*************************************************************
  **************          Blank data            ***************
  *************************************************************/
-
-/* Allocate space for one blank value of the given type and put the value
-   in it. */
 void *
 gal_data_alloc_blank(int type)
 {
   /* Define the pointers. */
-  void *allocated;
-  unsigned char *b;
-  char *c;
-  char **str;
-  unsigned short *us;
-  short *s;
-  unsigned int *ui;
-  int *i;
-  unsigned long *ul;
-  long *l;
-  LONGLONG *L;
-  float *f;
-  double *d;
-  gsl_complex_float *cx;
-  gsl_complex *dcx;
-
-  /* Allocate the space for the blank value: */
-  allocated=gal_data_malloc_array(type, 1);
+  unsigned char     uc = GAL_DATA_BLANK_UCHAR;
+  char               c = GAL_DATA_BLANK_CHAR;
+  char            *str = GAL_DATA_BLANK_STRING;
+  unsigned short    us = GAL_DATA_BLANK_USHORT;
+  short              s = GAL_DATA_BLANK_SHORT;
+  unsigned int      ui = GAL_DATA_BLANK_UINT;
+  int                i = GAL_DATA_BLANK_INT;
+  unsigned long     ul = GAL_DATA_BLANK_ULONG;
+  long               l = GAL_DATA_BLANK_LONG;
+  LONGLONG           L = GAL_DATA_BLANK_LONGLONG;
+  float              f = GAL_DATA_BLANK_FLOAT;
+  double             d = GAL_DATA_BLANK_DOUBLE;
+  gsl_complex_float cx;
+  gsl_complex      dcx;
 
   /* Put the blank value into it. */
-  errno=0;
   switch(type)
     {
     case GAL_DATA_TYPE_BIT:
@@ -392,77 +499,51 @@ gal_data_alloc_blank(int type)
             "us to see how we can implement it.");
 
     case GAL_DATA_TYPE_UCHAR:
-      b=allocated;
-      *b=GAL_DATA_BLANK_UCHAR;
-      return b;
+      return gal_data_alloc_number(type, &uc);
 
       /* CFITSIO says "int for keywords, char for table columns". Here we
          are only assuming table columns. So in practice this also applies
          to TSBYTE.*/
     case GAL_DATA_TYPE_CHAR: case GAL_DATA_TYPE_LOGICAL:
-      c=allocated;
-      *c=GAL_DATA_BLANK_CHAR;
-      return c;
+      return gal_data_alloc_number(type, &c);
 
     case GAL_DATA_TYPE_STRING:
-      str=allocated;
-      *str=GAL_DATA_BLANK_STRING;
-      return str;
+      return gal_data_alloc_number(type, &str);
 
     case GAL_DATA_TYPE_USHORT:
-      us=allocated;
-      *us=GAL_DATA_BLANK_USHORT;
-      return us;
+      return gal_data_alloc_number(type, &us);
 
     case GAL_DATA_TYPE_SHORT:
-      s=allocated;
-      *s=GAL_DATA_BLANK_SHORT;
-      return s;
+      return gal_data_alloc_number(type, &s);
 
     case GAL_DATA_TYPE_UINT:
-      ui=allocated;
-      *ui=GAL_DATA_BLANK_UINT;
-      return ui;
+      return gal_data_alloc_number(type, &ui);
 
     case GAL_DATA_TYPE_INT:
-      i=allocated;
-      *i=GAL_DATA_BLANK_INT;
-      return i;
+      return gal_data_alloc_number(type, &i);
 
     case GAL_DATA_TYPE_ULONG:
-      ul=allocated;
-      *ul=GAL_DATA_BLANK_ULONG;
-      return ul;
+      return gal_data_alloc_number(type, &ul);
 
     case GAL_DATA_TYPE_LONG:
-      l=allocated;
-      *l=GAL_DATA_BLANK_LONG;
-      return l;
+      return gal_data_alloc_number(type, &l);
 
     case GAL_DATA_TYPE_LONGLONG:
-      L=allocated;
-      *L=GAL_DATA_BLANK_LONGLONG;
-      return L;
+      return gal_data_alloc_number(type, &L);
 
     case GAL_DATA_TYPE_FLOAT:
-      f=allocated;
-      *f=GAL_DATA_BLANK_FLOAT;
-      return f;
+      return gal_data_alloc_number(type, &f);
 
     case GAL_DATA_TYPE_DOUBLE:
-      d=allocated;
-      *d=GAL_DATA_BLANK_DOUBLE;
-      return d;
+      return gal_data_alloc_number(type, &d);
 
     case GAL_DATA_TYPE_COMPLEX:
-      cx=allocated;
-      GSL_SET_COMPLEX(cx, GAL_DATA_BLANK_FLOAT, GAL_DATA_BLANK_FLOAT);
-      return cx;
+      GSL_SET_COMPLEX(&cx, GAL_DATA_BLANK_FLOAT, GAL_DATA_BLANK_FLOAT);
+      return gal_data_alloc_number(type, &cx);
 
     case GAL_DATA_TYPE_DCOMPLEX:
-      dcx=allocated;
-      GSL_SET_COMPLEX(dcx, GAL_DATA_BLANK_DOUBLE, GAL_DATA_BLANK_DOUBLE);
-      return dcx;
+      GSL_SET_COMPLEX(&dcx, GAL_DATA_BLANK_DOUBLE, GAL_DATA_BLANK_DOUBLE);
+      return gal_data_alloc_number(type, &dcx);
 
     default:
       error(EXIT_FAILURE, 0, "type value of %d not recognized in "
@@ -790,12 +871,12 @@ gal_data_blank_to_value(gal_data_t *data, void *value)
 
 
 /*************************************************************
- **************          Change type           ***************
+ **************             Copy           ***************
  *************************************************************/
-int
-gal_data_out_type(gal_data_t *first, gal_data_t *second)
+gal_data_t *
+gal_data_copy(gal_data_t *in)
 {
-  return first->type > second->type ? first->type : second->type;
+  return gal_data_copy_to_new_type(in, in->type);
 }
 
 
@@ -808,7 +889,19 @@ gal_data_copy_to_new_type(gal_data_t *in, int newtype)
   gal_data_t *out;
 
   /* Allocate space for the output type */
-  out=gal_data_alloc(newtype, in->ndim, in->dsize, 0, in->mmapped);
+  out=gal_data_alloc(NULL, newtype, in->ndim, in->dsize, 0, in->mmapped);
+
+  /* Copy the WCS structure, we need to have a blank WCS structure
+     allocated outside of WCSLIB, then copy the contents. */
+  if(in->wcs)
+    {
+      errno=0;
+      out->wcs=malloc(sizeof *out->wcs);
+      if(out->wcs==NULL)
+        error(EXIT_FAILURE, errno, "%zu bytes for out->wcs in "
+              "gal_data_copy_to_new_type", sizeof *out->wcs);
+      wcscopy(1, in->wcs, out->wcs);
+    }
 
   /* Fill in the output data array while doing the conversion */
   switch(newtype)
@@ -865,3 +958,234 @@ gal_data_copy_to_new_type(gal_data_t *in, int newtype)
   /* Return the created array */
   return out;
 }
+
+
+
+
+
+int
+gal_data_out_type(gal_data_t *first, gal_data_t *second)
+{
+  return first->type > second->type ? first->type : second->type;
+}
+
+
+
+
+
+/* The two input `f' and `s' datasets can be any type. But `of' and `os'
+   will have type `type', if freeinputs is non-zero, then the input arrays
+   will be freed if they needed to be changed to a new type. */
+void
+gal_data_to_same_type(gal_data_t *f,   gal_data_t *s,
+                      gal_data_t **of, gal_data_t **os,
+                      int type, int freeinputs)
+{
+  /* Change first dataset into the new type if necessary. */
+  if( f->type != type )
+    {
+      *of=gal_data_copy_to_new_type(f, type);
+      if(freeinputs)
+        gal_data_free(f);
+    }
+  else
+    *of=f;
+
+  /* Change second dataset into the new type if necessary. */
+  if( s->type != type )
+    {
+      *os=gal_data_copy_to_new_type(s, type);
+      if(freeinputs)
+        gal_data_free(s);
+    }
+  else
+    *os=s;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*************************************************************
+ **************              Read              ***************
+ *************************************************************/
+/* If the data structure was correctly created (the string was a number),
+   then return its pointer. Otherwise, return NULL. */
+gal_data_t *
+gal_data_string_to_number(char *string)
+{
+  int type;
+  long dsize[1]={1};
+  int fnz=-1, lnz=0;     /* `F'irst (or `L'ast) `N'on-`Z'ero. */
+  void *ptr, *numarr;
+  char *tailptr, *cp;
+
+  /* Define the pointers. */
+  unsigned char     uc;
+  char               c;
+  unsigned short    us;
+  short              s;
+  unsigned int      ui;
+  int                i;
+  unsigned long     ul;
+  long               l;
+  LONGLONG           L;
+  float              f;
+  double             d;
+
+  /* First see if the number is a double. */
+  errno=0;
+  d=strtod(string, &tailptr);
+  if(*tailptr!='\0')
+    return NULL;
+
+  /* See if the number is actually an integer: */
+  if( ceil(d) == d )
+    {
+      /* If the number is negative, put it in the signed types (based on
+         its value). If its zero or positive, then put it in the unsigned
+         types. */
+      if( d < 0 )
+        {
+          if     (d>CHAR_MIN)   {c=d; ptr=&c; type=GAL_DATA_TYPE_CHAR;}
+          else if(d>SHRT_MIN)   {s=d; ptr=&s; type=GAL_DATA_TYPE_SHORT;}
+          else if(d>INT_MIN)    {i=d; ptr=&i; type=GAL_DATA_TYPE_INT;}
+          else if(d>LONG_MIN)   {l=d; ptr=&l; type=GAL_DATA_TYPE_LONG;}
+          else                  {L=d; ptr=&L; type=GAL_DATA_TYPE_LONGLONG;}
+        }
+      else
+        {
+          if     (d<=UCHAR_MAX) {uc=d; ptr=&uc; type=GAL_DATA_TYPE_UCHAR;}
+          else if(d<=USHRT_MAX) {us=d; ptr=&us; type=GAL_DATA_TYPE_USHORT;}
+          else if(d<=UINT_MAX)  {ui=d; ptr=&ui; type=GAL_DATA_TYPE_UINT;}
+          else if(d<=ULONG_MAX) {ul=d; ptr=&ul; type=GAL_DATA_TYPE_ULONG;}
+          else                  {L=d;  ptr=&L;  type=GAL_DATA_TYPE_LONGLONG;}
+        }
+    }
+  else
+    {
+      /* The maximum number of decimal digits to store in float or double
+         precision floating point are:
+
+         float:  23 mantissa bits + 1 hidden bit: log(224)÷log(10) = 7.22
+         double: 52 mantissa bits + 1 hidden bit: log(253)÷log(10) = 15.95
+
+         FLT_DIG (at least 6 in ISO C) keeps the number of digits (not zero
+         before or after) that can be represented by a single precision
+         floating point number. If there are more digits, then we should
+         store the value as a double precision.
+
+         Note that the number can have non-digit characters that we don't
+         want, like: `.', `e', `E', `,'. */
+      for(cp=string;*cp!='\0';++cp)
+        if(isdigit(*cp) && *cp!='0' && fnz==-1)
+          fnz=cp-string;
+
+      /* In the previous loop, we went to the end of the string, so `cp'
+         now points to its `\0'. We just have to iterate backwards! */
+      for(;cp!=string;--cp)
+        if(isdigit(*cp) && *cp!='0')
+          {
+            lnz=cp-string;
+            break;
+          }
+
+      /* Calculate the number of decimal digits and decide if it the number
+         should be a float or a double. */
+      if( lnz-fnz < FLT_DIG || ( d<FLT_MAX && d>FLT_MIN ) )
+        { f=d; ptr=&f; type=GAL_DATA_TYPE_FLOAT; }
+      else
+        {      ptr=&d; type=GAL_DATA_TYPE_DOUBLE; }
+    }
+
+  /* Return the pointer to the data structure. */
+  numarr=gal_data_alloc_number(type, ptr);
+  return gal_data_alloc(numarr, type, 1, dsize, 0, 0);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*************************************************************
+ **************           Arithmetic           ***************
+ *************************************************************/
+gal_data_t *
+gal_data_arithmetic(char *operator, unsigned char flags, ...)
+{
+  int type;
+  va_list va;
+  gal_data_t *o=NULL;
+
+  /* Prepare the variable arguments (starting after the flags argument). */
+  va_start(va, flags);
+
+  /* Depending on the operator do the job: */
+  if      (!strcmp(operator, "+"))   { BINARY_INTERNAL(+, 0); }
+  else if (!strcmp(operator, "-"))   { BINARY_INTERNAL(-, 0); }
+  else if (!strcmp(operator, "*"))   { BINARY_INTERNAL(*, 0); }
+  else if (!strcmp(operator, "/"))   { BINARY_INTERNAL(/, 0); }
+
+#if 0
+  else if(!strcmp(operator, "abs"))       takeabs(p);
+  else if(!strcmp(operator, "pow"))       topower(p, NULL);
+  else if(!strcmp(operator, "sqrt"))      takesqrt(p);
+  else if(!strcmp(operator, "log"))       takelog(p);
+  else if(!strcmp(operator, "log10"))     takelog10(p);
+  else if(!strcmp(operator, "minvalue"))  findmin(p);
+  else if(!strcmp(operator, "maxvalue"))  findmax(p);
+  else if(!strcmp(operator, "min")
+          || !strcmp(operator, "max")
+          || !strcmp(operator, "average")
+          || !strcmp(operator, "median")) alloppixs(p, operator);
+  else if(!strcmp(operator, "lt")
+          || !strcmp(operator, "le")
+          || !strcmp(operator, "gt")
+          || !strcmp(operator, "ge")
+          || !strcmp(operator, "eq")
+          || !strcmp(operator, "neq"))    conditionals(p, operator);
+  else if(!strcmp(operator, "and")
+          || !strcmp(operator, "or"))     andor(p, operator);
+  else if(!strcmp(operator, "not"))       notfunc(p);
+  else if(!strcmp(operator, "isblank"))   opisblank(p);
+  else if(!strcmp(operator, "where"))     where(p);
+#endif
+
+  else
+    error(EXIT_FAILURE, 0, "the argument \"%s\" could not be "
+          "interpretted as an operator", operator);
+
+  /* End the variable argument structure and return. */
+  va_end(va);
+  return o;
+}
diff --git a/lib/fits.c b/lib/fits.c
index 852876b..dfd8836 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -1198,14 +1198,15 @@ gal_fits_read_wcs(char *filename, char *hdu, size_t 
hstartwcs,
    function. The value that is placed for those pixels is defined by
    the macros in fitsarrayvv.h and depends on the type of the data.*/
 gal_data_t *
-gal_fits_read_img_hdu(char *filename, char *hdu)
+gal_fits_read_img_hdu(char *filename, char *hdu, char *maskname,
+                      char *mhdu, size_t minmapsize)
 {
   void *blank;
   size_t i, ndim;
   fitsfile *fptr;
-  gal_data_t *data;
   int status=0, type;
   long *fpixel, *dsize;
+  gal_data_t *img, *mask;
 
 
   /* Check HDU for realistic conditions: */
@@ -1236,15 +1237,14 @@ gal_fits_read_img_hdu(char *filename, char *hdu)
 
 
   /* Allocate the space for the array and for the blank values. */
-  printf("alloc for %s\n", filename);
-  data=gal_data_alloc(type, (long)ndim, dsize, 0, 1);
+  img=gal_data_alloc(NULL, type, (long)ndim, dsize, 0, minmapsize);
   blank=gal_data_alloc_blank(type);
   free(dsize);
 
 
   /* Read the image into the allocated array: */
   fits_read_pix(fptr, gal_fits_type_to_datatype(type), fpixel,
-                data->size, blank, data->array, &data->anyblank, &status);
+                img->size, blank, img->array, &img->anyblank, &status);
   if(status) gal_fits_io_error(status, NULL);
   free(fpixel);
   free(blank);
@@ -1253,7 +1253,23 @@ gal_fits_read_img_hdu(char *filename, char *hdu)
   /* Close the FITS file, and return the data pointer */
   fits_close_file(fptr, &status);
   gal_fits_io_error(status, NULL);
-  return data;
+
+
+  /* If a mask was specified, read it into the mask data structure, then
+     set all the corresponding pixels of the input image to NaN. */
+  if(maskname)
+    {
+      /* Read the mask HDU. */
+      mask=gal_fits_read_img_hdu(maskname, mhdu, NULL, NULL, minmapsize);
+
+      /* Apply the mask on the input. */
+      gal_data_apply_mask(img, mask);
+
+      /* Free the mask space. */
+      gal_data_free(mask);
+    }
+
+  return img;
 }
 
 
@@ -1268,13 +1284,13 @@ gal_fits_read_img_hdu(char *filename, char *hdu)
    original input bitpix will be stored so if you want to, you can
    return it back to the input type if you please. */
 gal_data_t *
-gal_fits_read_to_type(char *inputname, char *maskname, char *inhdu,
-                      char *mhdu, int type)
+gal_fits_read_to_type(char *inputname, char *inhdu, char *maskname,
+                      char *mhdu, int type, size_t minmapsize)
 {
-  gal_data_t *in, *mask, *converted;
+  gal_data_t *in, *converted;
 
   /* Read the specified input image HDU. */
-  in=gal_fits_read_img_hdu(inputname, inhdu);
+  in=gal_fits_read_img_hdu(inputname, inhdu, maskname, mhdu, minmapsize);
 
   /* If the input had another type, convert it to float. */
   if(in->type!=type)
@@ -1284,20 +1300,7 @@ gal_fits_read_to_type(char *inputname, char *maskname, 
char *inhdu,
       in=converted;
     }
 
-  /* If a mask was specified, read it as a float image, then set all
-     the corresponding pixels of the input image to NaN. */
-  if(maskname)
-    {
-      /* Read the mask HDU. */
-      mask=gal_fits_read_img_hdu(maskname, mhdu);
-
-      /* Apply the mask on the input. */
-      gal_data_apply_mask(in, mask);
-
-      /* Free the mask space */
-      gal_data_free(mask);
-    }
-
+  /* Return the final structure. */
   return in;
 }
 
@@ -1316,8 +1319,8 @@ gal_fits_read_float_kernel(char *inputname, char *inhdu, 
float **outkernel,
   float *f, *fp, tmp;
 
   /* Read the image as a float */
-  kernel=gal_fits_read_to_type(inputname, NULL, inhdu, NULL,
-                                GAL_DATA_TYPE_FLOAT);
+  kernel=gal_fits_read_to_type(inputname, inhdu, NULL, NULL,
+                               GAL_DATA_TYPE_FLOAT, -1);
 
   /* Check if the size along each dimension of the kernel is an odd
      number. If they are all an odd number, then the for each dimension,
diff --git a/lib/gnuastro/data.h b/lib/gnuastro/data.h
index ea1efa8..f37e10b 100644
--- a/lib/gnuastro/data.h
+++ b/lib/gnuastro/data.h
@@ -59,6 +59,11 @@ __BEGIN_C_DECLS  /* From C++ preparations */
 /* The maximum dimensionality of datasets. */
 #define GAL_DATA_MAXDIM    999
 
+/* Arithmetic macros (powers of 2). */
+#define GAL_DATA_ARITH_INPLACE  1
+#define GAL_DATA_ARITH_FREE     2
+#define GAL_DATA_ARITH_NUMOK    4
+
 /* Blank values: Note that for the unsigned types or small types (like
    char), the maximum value is considered as a blank value, since the
    minimum value of an unsigned type is zero and zero is often meaningful
@@ -147,7 +152,8 @@ void *
 gal_data_calloc_array(int type, size_t size);
 
 gal_data_t *
-gal_data_alloc(int type, size_t ndim, long *dsize, int clear, int map);
+gal_data_alloc(void *array, int type, size_t ndim, long *dsize,
+               int clear, size_t minmapsize);
 
 void
 gal_data_free(gal_data_t *data);
@@ -174,13 +180,42 @@ gal_data_blank_to_value(gal_data_t *data, void *value);
 
 
 /*************************************************************
- **************         Convert types          ***************
+ **************             Copy               ***************
  *************************************************************/
+gal_data_t *
+gal_data_copy(gal_data_t *in);
+
+gal_data_t *
+gal_data_copy_to_new_type(gal_data_t *in, int newtype);
+
 int
 gal_data_out_type(gal_data_t *first, gal_data_t *second);
 
+void
+gal_data_to_same_type(gal_data_t *f, gal_data_t *s,
+                      gal_data_t **of, gal_data_t **os,
+                      int type, int freeinputs);
+
+
+
+
+
+/*************************************************************
+ **************              Read              ***************
+ *************************************************************/
 gal_data_t *
-gal_data_copy_to_new_type(gal_data_t *in, int newtype);
+gal_data_string_to_number(char *string);
+
+
+
+
+
+/*************************************************************
+ **************           Arithmetic           ***************
+ *************************************************************/
+gal_data_t *
+gal_data_arithmetic(char *operator, unsigned char flags, ...);
+
 
 
 
diff --git a/lib/gnuastro/fits.h b/lib/gnuastro/fits.h
index 6d99114..de2c7c0 100644
--- a/lib/gnuastro/fits.h
+++ b/lib/gnuastro/fits.h
@@ -228,11 +228,12 @@ gal_fits_read_wcs(char *filename, char *hdu, size_t 
hstartwcs,
  ******************     Array functions      *****************
  *************************************************************/
 gal_data_t *
-gal_fits_read_img_hdu(char *filename, char *hdu);
+gal_fits_read_img_hdu(char *filename, char *hdu, char *maskname,
+                      char *mhdu, size_t minmapsize);
 
 gal_data_t *
-gal_fits_read_to_type(char *inputname, char *maskname, char *inhdu,
-                      char *mhdu, int type);
+gal_fits_read_to_type(char *inputname, char *inhdu, char *maskname,
+                      char *mhdu, int type, size_t minmapsize);
 
 gal_data_t *
 gal_fits_read_float_kernel(char *inputname, char *inhdu, float **outkernel,