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