1 /*
2  *  @file    dx-math.c
3  *  @author  Nathan Baker
4  *  @brief   Arithmetic with DX files
5  *  @version $Id$
6  */
7 
8 /*
9  *  Last update: 08/29/2016 by Leighton Wilson
10  *  Description: Added ability to read in binary DX files as input
11  */
12 
13 #include "apbs.h"
14 
15 #define IJK(i,j,k)  (((k)*(nx)*(ny))+((j)*(nx))+(i))
16 #define ERRRC 13
17 #define DXM_MAXOP 20
18 #define DXM_ISGRID 0
19 #define DXM_ISSCALAR 1
20 #define DXM_ISGRIDBIN 2
21 
22 VEMBED(rcsid="$Id$")
23 
24 typedef enum Dxmath_Opcode {
25     DXM_ADD,  /**< + */
26     DXM_SUB,  /**< - */
27     DXM_DIV,  /**< / */
28     DXM_MUL,  /**< * */
29     DXM_EQU,  /**< = */
30     DXM_EXP   /**< ^ */
31 } Dxmath_Opcode;
32 
main(int argc,char ** argv)33 int main(int argc, char **argv) {
34 
35     /* *************** VARIABLES ******************* */
36     char *input_path;
37     char *dot;
38     char *MCwhiteChars = " \t\n";
39     char *MCcommChars  = "#%";
40     char tok[VMAX_BUFSIZE];
41     char format[VMAX_BUFSIZE];
42     char gridPath[DXM_MAXOP+1][VMAX_BUFSIZE];
43     double scalar[DXM_MAXOP+1];
44     int obType[DXM_MAXOP+1];
45     int iop, numop;
46     int i, nx, ny, nz, len;
47     Dxmath_Opcode op[DXM_MAXOP];
48     Vio *sock = VNULL;
49     Vgrid *grid1 = VNULL;
50     Vgrid *grid2 = VNULL;
51 
52     char *header = "\n\n\
53     ----------------------------------------------------------------------\n\
54     ----------------------------------------------------------------------\n\
55     \n\n";
56 
57     char *usage = "\n\n\
58     ----------------------------------------------------------------------\n\
59     This driver program (like its UHBD counterpart) does simple arithmetic\n\
60     with Cartesian grid data.  It is invoked as:\n\n\
61       dx-math <path>\n\n\
62     where <path> is the path to a dx file with operations specified\n\
63     in a stack-based (RPN) manner.  For example, a command file which adds\n\
64     grid1 and grid2, multiplies the result by 5.3, adds grid4, subtracts\n\
65     99.3 from the whole thing, and writes the result on grid5 would have the\n\
66     form:\n\n\
67       grid1\n\
68       grid2 +\n\
69       5.3 *\n\
70       grid4 +\n\
71       99.3 -\n\
72       grid5 =\n\n\
73     where the file names, scalar values, and operations must be separated by\n\
74     tabs, line breaks, or white space.  Comments can be included between the\n\
75     character # and a new line (in the usual shell script fashion).\n\n\
76     If the file extension is .dxbin, the file is assumed to be in binary dx\n\
77     format. Otherwise, the file is assumed to be in standard dx format.\n\
78     Similarly, if the output file extension is .dxbin, the output grid file\n\
79     will be in binary dx format, and standard dx otherwise.\n\
80     ----------------------------------------------------------------------\n\n";
81 
82 
83     /* *************** CHECK INVOCATION ******************* */
84     Vio_start();
85     Vnm_print(1, "%s", header);
86     if (argc != 2) {
87         Vnm_print(2,"\n*** Syntax error: got %d arguments, expected 2.\n\n",
88           argc);
89         Vnm_print(2,"%s\n", usage);
90         return ERRRC;
91     }
92     input_path = argv[1];
93 
94     /* *************** OPEN INPUT FILE ******************* */
95     sock = Vio_ctor("FILE", "ASC", VNULL, input_path, "r");
96     if (sock == VNULL) {
97         Vnm_print(2, "main:  Null socket!\n");
98         return ERRRC;
99     }
100     if (Vio_accept(sock, 0) < 0) {
101         Vnm_print(2, "main:  Problem reading from socket!\n");
102         return 0;
103     }
104     Vio_setWhiteChars(sock, MCwhiteChars);
105     Vio_setCommChars(sock, MCcommChars);
106 
107     /* *************** PARSE INPUT FILE ******************* */
108     /* After reading in the first arg, we should alternate between objects and
109      * operations, starting with the objects.  For each opject, we assign a
110      * type (0 = grid path, 1 = scalar) */
111     if (Vio_scanf(sock, "%s", tok) != 1) {
112         Vnm_print(2, "main:  Ran out of tokens when parsing initial input!\n");
113         return ERRRC;
114     }
115 
116     strncpy(gridPath[0], tok, VMAX_BUFSIZE);
117     dot = strrchr(tok, '.');
118     if (dot && !Vstring_strcasecmp(dot, ".dxbin")) {
119         Vnm_print(1, "main:  Assuming %s is OpenDX binary format...\n", tok);
120         obType[0] = DXM_ISGRIDBIN;
121     } else {
122         Vnm_print(1, "main:  Assuming %s is standard OpenDX format...\n", tok);
123         obType[0] = DXM_ISGRID;
124     }
125 
126     numop = 0;
127     while (Vio_scanf(sock, "%s", tok) == 1) {
128         if (sscanf(tok, "%lg", &scalar[numop+1]) == 1) {
129             obType[numop+1] = DXM_ISSCALAR;
130         } else {
131             strncpy(gridPath[numop+1], tok, VMAX_BUFSIZE);
132 
133             dot = strrchr(tok, '.');
134             if (dot && !Vstring_strcasecmp(dot, ".dxbin")) {
135                 Vnm_print(1, "main:  Assuming %s is OpenDX binary format...\n", tok);
136                 obType[numop+1] = DXM_ISGRIDBIN;
137             } else {
138                 Vnm_print(1, "main:  Assuming %s is standard OpenDX format...\n", tok);
139                 obType[numop+1] = DXM_ISGRID;
140             }
141         }
142 
143         if (Vio_scanf(sock, "%s", tok) != 1) {
144             Vnm_print(2, "main:  Ran out of tokens when parsing input!\n");
145             Vnm_print(2, "main:  Last token = %s.\n", tok);
146             return ERRRC;
147         }
148 
149         if (strcmp(tok, "*") == 0) op[numop] = DXM_MUL;
150         else if (strcmp(tok, "+") == 0) op[numop] = DXM_ADD;
151         else if (strcmp(tok, "-") == 0) op[numop] = DXM_SUB;
152         else if (strcmp(tok, "/") == 0) op[numop] = DXM_DIV;
153         else if (strcmp(tok, "^") == 0) op[numop] = DXM_EXP;
154         else if (strcmp(tok, "=") == 0) {
155             op[numop] = DXM_EQU;
156             numop++;
157             break;
158         } else {
159             Vnm_print(2, "main:  Undefined operation '%s'!\n", tok);
160             return ERRRC;
161         }
162         numop++;
163         if (numop == DXM_MAXOP) {
164             Vnm_print(2, "main:  Exceed maximum number of operations (%d)!\n",
165               DXM_MAXOP);
166             return ERRRC;
167         }
168     }
169     Vio_acceptFree(sock);
170 
171     /* Spit out what we parsed: */
172     Vnm_print(1, "main:  Performing following operations:\n");
173     Vnm_print(1, "main:    %s \n", gridPath[0]);
174     for (iop=0; iop<numop; iop++) {
175         if (obType[iop+1] == DXM_ISGRID)
176           Vnm_print(1, "main:    %s (dx grid) ", gridPath[iop+1]);
177         else if (obType[iop+1] == DXM_ISGRIDBIN)
178           Vnm_print(1, "main:    %s (dxbin grid) ", gridPath[iop+1]);
179         else
180           Vnm_print(1, "main:    %g (scalar) ", scalar[iop+1]);
181         switch (op[iop]) {
182             case DXM_MUL:
183                 Vnm_print(1, "*\n");
184                 break;
185             case DXM_ADD:
186                 Vnm_print(1, "+\n");
187                 break;
188             case DXM_SUB:
189                 Vnm_print(1, "-\n");
190                 break;
191             case DXM_DIV:
192                 Vnm_print(1, "/\n");
193                 break;
194             case DXM_EQU:
195                 Vnm_print(1, "=\n");
196                 break;
197             case DXM_EXP:
198                 Vnm_print(1, "^\n");
199                 break;
200             default:
201                 Vnm_print(2, "\nmain:  Unknown operation (%d)!", op[iop]);
202                 return ERRRC;
203         }
204     }
205 
206     /* *************** PARSE INPUT FILE ******************* */
207     Vnm_print(1, "main:  Reading grid from %s...\n", gridPath[0]);
208     grid1 = Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
209     if (obType[0] == DXM_ISGRID) {
210         if (!Vgrid_readDX(grid1, "FILE", "ASC", VNULL, gridPath[0])) {
211             Vnm_print(2, "main:  Problem reading OpenDX-format grid from %s\n",
212                       gridPath[0]);
213             return ERRRC;
214         }
215     } else if (obType[0] == DXM_ISGRIDBIN) {
216         if (!Vgrid_readDXBIN(grid1, "FILE", "ASC", VNULL, gridPath[0])) {
217             Vnm_print(2, "main:  Problem reading OpenDX binary-format grid from %s\n",
218                       gridPath[0]);
219             return ERRRC;
220         }
221     } else {
222             Vnm_print(2, "main:  First argument must be a grid\n");
223             return ERRRC;
224     }
225 
226     nx = grid1->nx;
227     ny = grid1->ny;
228     nz = grid1->nz;
229     len = nx * ny * nz;
230 
231     for (iop=0; iop<numop-1; iop++) {
232         if (obType[iop+1] == DXM_ISGRID || obType[iop+1] == DXM_ISGRIDBIN) {
233             Vnm_print(1, "main:  Reading grid from %s...\n", gridPath[iop+1]);
234             grid2 = Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
235 
236             if (obType[iop+1] == DXM_ISGRID) {
237                 if (!Vgrid_readDX(grid2, "FILE", "ASC", VNULL, gridPath[iop+1])) {
238                     Vnm_print(2, "main:  Problem reading OpenDX-format grid from %s\n",
239                           gridPath[iop+1]);
240                     return ERRRC;
241                 }
242             } else if (obType[iop+1] == DXM_ISGRIDBIN) {
243                 if (!Vgrid_readDXBIN(grid2, "FILE", "ASC", VNULL, gridPath[iop+1])) {
244                     Vnm_print(2, "main:  Problem reading OpenDX binary-format grid \
245                           from %s\n",
246                           gridPath[iop+1]);
247                     return ERRRC;
248                 }
249             }
250 
251             if ((grid2->nx != nx) || (grid2->ny != ny) || (grid2->nz != nz)) {
252                 Vnm_print(2, "main:  Grid dimension mis-match!\n");
253                 Vnm_print(2, "main:  Grid 1 is %d x %d x %d\n", nx, ny, nz);
254                 Vnm_print(2, "main:  Grid 2 is %d x %d x %d\n",
255                   grid2->nx, grid2->ny, grid2->nz);
256                 return ERRRC;
257             }
258 
259             switch (op[iop]) {
260                 case DXM_ADD:
261                     Vnm_print(1, "main:  Adding...\n");
262                     for (i=0; i<len; i++)
263                       grid1->data[i] = grid1->data[i] + grid2->data[i];
264                     break;
265                 case DXM_MUL:
266                     Vnm_print(1, "main:  Multiplying...\n");
267                     for (i=0; i<len; i++)
268                       grid1->data[i] = grid1->data[i] * grid2->data[i];
269                     break;
270                 case DXM_SUB:
271                     Vnm_print(1, "main:  Subtracting...\n");
272                     for (i=0; i<len; i++)
273                       grid1->data[i] = grid1->data[i] - grid2->data[i];
274                     break;
275                 case DXM_DIV:
276                     Vnm_print(1, "main:  Dividing...\n");
277                     for (i=0; i<len; i++)
278                       grid1->data[i] = grid1->data[i] / grid2->data[i];
279                     break;
280                 case DXM_EXP:
281                     Vnm_print(1, "main:  Applying exponents...\n");
282                     for(i = 0; i < grid2->data[i]; i++) {
283                         grid1->data[i] *= grid1->data[i];
284                     }
285                     break;
286                 default:
287                     Vnm_print(2, "main:  Unexpected operation (%d)!\n",
288                       op[iop]);
289                     break;
290             }
291             Vgrid_dtor(&grid2);
292         } else { /* if (obType[iop+1] == DXM_ISGRID) */
293             Vnm_print(1, "main:  Loading scalar %g...\n", scalar[iop+1]);
294             switch (op[iop]) {
295                 case DXM_ADD:
296                     Vnm_print(1, "main:  Adding...\n");
297                     for (i=0; i<len; i++)
298                       grid1->data[i] = grid1->data[i] + scalar[iop+1];
299                     break;
300                 case DXM_MUL:
301                     Vnm_print(1, "main:  Multiplying...\n");
302                     for (i=0; i<len; i++)
303                       grid1->data[i] = grid1->data[i] * scalar[iop+1];
304                     break;
305                 case DXM_SUB:
306                     Vnm_print(1, "main:  Subtracting...\n");
307                     for (i=0; i<len; i++)
308                       grid1->data[i] = grid1->data[i] - scalar[iop+1];
309                     break;
310                 case DXM_DIV:
311                     Vnm_print(1, "main:  Dividing...\n");
312                     for (i=0; i<len; i++)
313                       grid1->data[i] = grid1->data[i] / scalar[iop+1];
314                     break;
315                 case DXM_EXP:
316                     Vnm_print(1, "main:  Applying exponents...\n");
317                     for(i = 0; i < scalar[iop+1]; i++) {
318                         grid1->data[i] *= grid1->data[i];
319                     }
320                     break;
321                 default:
322                     Vnm_print(2, "main:  Unexpected operation (%d)!\n",
323                       op[iop]);
324                     break;
325             }
326         } /* if (obType[iop+1] == DXM_ISGRID) */
327     } /* for (iop=0; iop<numop-1; iop++) */
328 
329     /* The last operation is the = sign, implying that we write out the grid */
330     Vnm_print(1, "main:  Writing results to %s...\n", gridPath[numop]);
331 
332     if (obType[numop] == DXM_ISGRID) {
333         Vgrid_writeDX(grid1, "FILE", "ASC", VNULL, gridPath[numop],
334                       "DXMATH RESULTS", VNULL);
335     } else if (obType[numop] == DXM_ISGRIDBIN) {
336         Vgrid_writeDXBIN(grid1, "FILE", "ASC", VNULL, gridPath[numop],
337                          "DXMATH RESULTS", VNULL);
338     } else {
339         Vnm_print(2, "main:  Last object must be output grid\n");
340         return ERRRC;
341     }
342 
343     Vnm_print(1, "main:  All done -- exiting.\n");
344     return 0;
345 }
346