Logo Search packages:      
Sourcecode: h5utils version File versions  Download package

h5math.c

/* Copyright (c) 1999-2009 Massachusetts Institute of Technology
 *
 * Permission is hereby granted, free of charge, to any person obtaining
 * a copy of this software and associated documentation files (the
 * "Software"), to deal in the Software without restriction, including
 * without limitation the rights to use, copy, modify, merge, publish,
 * distribute, sublicense, and/or sell copies of the Software, and to
 * permit persons to whom the Software is furnished to do so, subject to
 * the following conditions:
 * 
 * The above copyright notice and this permission notice shall be
 * included in all copies or substantial portions of the Software.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
 * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
 * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
 * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>

#include <unistd.h>

#include "config.h"
#include "arrayh5.h"
#include "copyright.h"
#include "h5utils.h"

#include <matheval.h>

#define CHECK(cond, msg) { if (!(cond)) { fprintf(stderr, "h5math error: %s\n", msg); exit(EXIT_FAILURE); } }

const char default_data_name[] = "h5math";

void usage(FILE *f)
{
     fprintf(f, "Usage: h5math [options] <output-filename> [<filenames>]\n"
           "Options:\n"
           "         -h : this help message\n"
             "         -V : print version number and copyright\n"
           "         -v : verbose output\n"
           "         -a : append to existing hdf5 file\n"
           "  -n <size> : output array dimensions [ default: from input ]\n"
           "  -f <file> : read expression to evaluate from file [ default: stdin ]\n"
           "  -e <expr> : evaluate <expr> to output\n"
           "    -x <ix> : take x=<ix> slice of data\n"
           "    -y <iy> : take y=<iy> slice of data\n"
           "    -z <iz> : take z=<iz> slice of data\n"
           "    -t <it> : take t=<it> slice of data's last dimension\n"
           "         -0 : use dataset center as origin for -x/-y/-z\n"
           "     -r <r> : use resolution <r> for xyz coordinate units in expression\n"
           "  -d <name> : use dataset <name> in the input/output files\n"
           "              [ default: first dataset/%s ]\n"
           "              -- you can also specify a dataset via <filename>:<name>\n",
           default_data_name
        );
}

#define MAX_RANK 10

int main(int argc, char **argv)
{
     arrayh5 *a, ao;
     int i, n;
     int rank = -1, dims[MAX_RANK];
     extern char *optarg;
     extern int optind;
     int c;
     int slicedim[4] = {NO_SLICE_DIM,NO_SLICE_DIM,NO_SLICE_DIM,NO_SLICE_DIM};
     int islice[4], center_slice[4] = {0,0,0,0};
     int verbose = 0;
     int append = 0;
     char *expr_string = 0, *expr_filename = 0;
     char *data_name = 0;
     char *out_fname, *out_dname;
     char **eval_vars;
     int eval_nvars;
     char **vars;
     double *vals;
     void *evaluator;
     double res = 1.0;
     int nx, ny, nz, nt, nr, ix, iy, iz, it, ir;
     double cx, cy, cz;

     while ((c = getopt(argc, argv, "hVvan:f:e:x:y:z:t:0d:r:")) != -1)
        switch (c) {
            case 'h':
               usage(stdout);
               return EXIT_SUCCESS;
            case 'V':
               printf("h5totxt " PACKAGE_VERSION " by Steven G. Johnson\n" 
                    COPYRIGHT);
               return EXIT_SUCCESS;
            case 'v':
               verbose = 1;
               break;
            case 'a':
               append = 1;
               break;
            case 'n':
            {
               int pos = 0;
               rank = 0;
               while (isdigit(optarg[pos])) {
                  CHECK(rank < MAX_RANK,
                        "Rank too big in -n argument!\n");
                  dims[rank] = 0;
                  while (isdigit(optarg[pos])) {
                       dims[rank] = dims[rank]*10 + optarg[pos]-'0';
                       ++pos;
                  }
                  ++rank;
                  if (optarg[pos] == 'x' || optarg[pos] == 'X'
                      || optarg[pos] == '*')
                       ++pos;
               }
               CHECK(rank > 0 && !optarg[pos], "Invalid -n argument; should be e.g. 23x34 or 10x10x10\n");
               break;
            }
            case 'f':
               free(expr_filename);
               expr_filename = my_strdup(optarg);
               break;
            case 'e':
               free(expr_string);
               expr_string = my_strdup(optarg);
               break;
            case 'x':
               islice[0] = atoi(optarg);
               slicedim[0] = 0;
               break;
            case 'y':
               islice[1] = atoi(optarg);
               slicedim[1] = 1;
               break;
            case 'z':
               islice[2] = atoi(optarg);
               slicedim[2] = 2;
               break;
            case 't':
               islice[3] = atoi(optarg);
               slicedim[3] = LAST_SLICE_DIM;
               break;
              case '0':
                   center_slice[0] = center_slice[1] = center_slice[2] = 1;
                   break;
            case 'r':
               res = atof(optarg);
               break;
            case 'd':
               free(data_name);
               data_name = my_strdup(optarg);
               break;            
            default:
               fprintf(stderr, "Invalid argument -%c\n", c);
               usage(stderr);
               return EXIT_FAILURE;
        }
     if (optind == argc) {  /* no parameters left */
        usage(stderr);
        return EXIT_FAILURE;
     }

     out_fname = split_fname(argv[optind], &out_dname);
     if (!out_dname[0]) {
        if (data_name)
             out_dname = data_name;
        else
             out_dname = (char *) default_data_name;
     }
     optind++;

     n = argc - optind;
     a = (arrayh5 *) malloc(sizeof(arrayh5) * n);
     CHECK(a, "out of memory");

     for (i = 0; i < n; ++i) {
        int err;
        char *fname, *dname;

          fname = split_fname(argv[i + optind], &dname);
          if (!dname[0])
               dname = data_name;

        err = arrayh5_read(&a[i], fname, dname, NULL,
                             4, slicedim, islice, center_slice);
          CHECK(!err, arrayh5_read_strerror[err]);

        CHECK(!i || arrayh5_conformant(a[i], a[i-1]),
            "all input arrays must have the same dimensions");

        if (verbose)
             printf("read variable d%d: dataset \"%s\" in file \"%s\"\n",
                  i + 1, dname ? dname : "<first>", fname);
        
        free(fname);
     }

     if (rank >= 0) {
        ao = arrayh5_create(rank, dims);
        CHECK(!n || arrayh5_conformant(ao, a[0]),
            "-n dimensions must be same as those of input arrays");
     }
     else if (n)
        ao = arrayh5_clone(a[0]);
     else
        CHECK(0, "output size must be specified with -n if no input arrays");

     if (verbose) {
        printf("rank-%d array dimensions: ", ao.rank);
        if (!ao.rank) printf("1\n");
        for (i = 0; i < ao.rank; ++i)
             printf("%s%d", i ? "x" : "", ao.dims[i]);
        printf("\n");
     }

     nx = ao.rank >= 1 ? ao.dims[0] : 1;
     ny = ao.rank >= 2 ? ao.dims[1] : 1;
     nz = ao.rank >= 3 ? ao.dims[2] : 1;
     nt = ao.rank >= 4 ? ao.dims[3] : 1;
     for (nr = 1, i = 4; i < ao.rank; ++i)
        nr *= ao.dims[i];
     cx = center_slice[0] ? (nx - 1) * 0.5 : 0.0;
     cy = center_slice[1] ? (ny - 1) * 0.5 : 0.0;
     cz = center_slice[2] ? (nz - 1) * 0.5 : 0.0;

     vars = (char **) malloc(sizeof(char *) * (n + 4));
     CHECK(vars, "out of memory");
     vals = (double *) malloc(sizeof(double) * (n + 4));
     CHECK(vals, "out of memory");
     for (i = 0; i < n; ++i) {
        vars[i] = my_strdup("dxxxxxxxxxxxx");
#ifdef HAVE_SNPRINTF
        snprintf(vars[i], 14, "d%d", i + 1);
#else
        sprintf(vars[i], "d%d", i + 1);
#endif
        vals[i] = 0.0;
     }
     vars[n+0] = strdup("x"); vals[n+0] = 0.0;
     vars[n+1] = strdup("y"); vals[n+1] = 0.0;
     vars[n+2] = strdup("z"); vals[n+2] = 0.0;
     vars[n+3] = strdup("t"); vals[n+3] = 0.0;
     
     if (!expr_string) {
        char buf[1024] = "";
        int len;
        FILE *f = expr_filename ? fopen(expr_filename, "r") : stdin;
        CHECK(f, "unable to open expression file");

        if (verbose && f == stdin)
             printf("Enter expression to write to %s:\n", out_fname);
        
        fgets(buf, 1024, f);
        expr_string = my_strdup(buf);
        len = strlen(buf) + 1;
        while (fgets(buf, 1024, f)) {
             len += strlen(buf);
             expr_string = (char *) realloc(expr_string, len);
             strcat(expr_string, buf);
        }
        for (ix = 0; ix < len; ++ix)
             if (expr_string[ix] == '\n')
                expr_string[ix] = ' '; /* matheval chokes on newlines */

        if (expr_filename) fclose(f);
     }

     CHECK(evaluator = evaluator_create(expr_string),
         "error parsing symbolic expression");

     evaluator_get_variables(evaluator, &eval_vars, &eval_nvars);
     for (ix = 0; ix < eval_nvars; ++ix) {
        for (iy = 0; iy < n + 4 && strcmp(eval_vars[ix], vars[iy]); ++iy)
             ;
        if (iy == n + 4) {
             fprintf(stderr, "h5math error: unrecognized variable \"%s\"\n",
                   eval_vars[ix]);
             exit(EXIT_FAILURE);
        }
     }
     
     if (verbose) {
        char *buf = evaluator_get_string(evaluator);
        printf("Evaluating expression: %s\n", buf);
     }

     for (ix = 0; ix < nx; ++ix)
     for (iy = 0; iy < ny; ++iy)
     for (iz = 0; iz < nz; ++iz)
     for (it = 0; it < nt; ++it)
     for (ir = 0; ir < nr; ++ir) {
        int idx = ir + nr * (it + nt * (iz + nz * (iy + ny * ix)));
        for (i = 0; i < n; ++i)
             vals[i] = a[i].data[idx];
        vals[n+0] = (ix - cx) / res;
        vals[n+1] = (iy - cy) / res;
        vals[n+2] = (iz - cz) / res;
        vals[n+3] = ao.rank >= 4 ? it : 
             (ao.rank >= 3 ? iz : (ao.rank >= 2 ? iy : ix));
        ao.data[idx] = evaluator_evaluate(evaluator, n+4, vars, vals);
     }

     if (verbose)
        printf("Writing data to \"%s\" in \"%s\"...\n", 
             out_dname ? out_dname : "<first>", out_fname);
     arrayh5_write(ao, out_fname, out_dname, append);

     free(vals);
     for (i = 0; i < n+4; ++i) free(vars[i]);
     free(vars);
     arrayh5_destroy(ao);
     for (i = 0; i < n; ++i) arrayh5_destroy(a[i]);
     free(a);
     free(out_fname);
     free(expr_filename);
     free(expr_string);
     free(data_name);

     return EXIT_SUCCESS;
}

Generated by  Doxygen 1.6.0   Back to index