1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
11 *
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
16 *
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
21 *
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 *
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
34 *
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
37 */
38 #include "gmxpre.h"
39
40 #include "xvgr.h"
41
42 #include <cassert>
43 #include <cctype>
44 #include <cstring>
45
46 #include <string>
47
48 #include "gromacs/fileio/gmxfio.h"
49 #include "gromacs/fileio/oenv.h"
50 #include "gromacs/math/vec.h"
51 #include "gromacs/utility/basedefinitions.h"
52 #include "gromacs/utility/binaryinformation.h"
53 #include "gromacs/utility/coolstuff.h"
54 #include "gromacs/utility/cstringutil.h"
55 #include "gromacs/utility/exceptions.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/futil.h"
58 #include "gromacs/utility/smalloc.h"
59 #include "gromacs/utility/sysinfo.h"
60
output_env_get_print_xvgr_codes(const gmx_output_env_t * oenv)61 gmx_bool output_env_get_print_xvgr_codes(const gmx_output_env_t* oenv)
62 {
63 XvgFormat xvgFormat = output_env_get_xvg_format(oenv);
64
65 return (xvgFormat == XvgFormat::Xmgrace || xvgFormat == XvgFormat::Xmgr);
66 }
67
xvgrstr(const std::string & gmx,const gmx_output_env_t * oenv,char * buf,int buflen)68 static char* xvgrstr(const std::string& gmx, const gmx_output_env_t* oenv, char* buf, int buflen)
69 {
70 /* Supported greek letter names and corresponding xmgrace/xmgr symbols */
71 const char* sym[] = { "beta", "chi", "delta", "eta", "lambda", "mu",
72 "omega", "phi", "psi", "rho", "theta", nullptr };
73 const char symc[] = { 'b', 'c', 'd', 'h', 'l', 'm', 'w', 'f', 'y', 'r', 'q', '\0' };
74 gmx_bool bXVGR;
75 int g, b, i;
76 char c;
77
78 XvgFormat xvgf = output_env_get_xvg_format(oenv);
79 bXVGR = (xvgf == XvgFormat::Xmgrace || xvgf == XvgFormat::Xmgr);
80
81 g = 0;
82 b = 0;
83 while (gmx[g] != '\0')
84 {
85 /* Check with the largest string we have ("lambda"), add one for \0 */
86 if (b + 6 + 1 >= buflen)
87 {
88 gmx_fatal(FARGS,
89 "Output buffer length in xvgstr (%d) too small to process xvg input string "
90 "'%s'",
91 buflen, gmx.c_str());
92 }
93 if (gmx[g] == '\\')
94 {
95 g++;
96 if (gmx[g] == 's')
97 {
98 /* Subscript */
99 if (bXVGR)
100 {
101 buf[b++] = '\\';
102 buf[b++] = 's';
103 }
104 else
105 {
106 buf[b++] = '_';
107 }
108 g++;
109 }
110 else if (gmx[g] == 'S')
111 {
112 /* Superscript */
113 if (bXVGR)
114 {
115 buf[b++] = '\\';
116 buf[b++] = 'S';
117 }
118 else
119 {
120 buf[b++] = '^';
121 }
122 g++;
123 }
124 else if (gmx[g] == 'N')
125 {
126 /* End sub/superscript */
127 if (bXVGR)
128 {
129 buf[b++] = '\\';
130 buf[b++] = 'N';
131 }
132 else
133 {
134 if (gmx[g + 1] != ' ')
135 {
136 buf[b++] = ' ';
137 }
138 }
139 g++;
140 }
141 else if (gmx[g] == '4')
142 {
143 /* Backward compatibility for xmgr normal font "\4" */
144 switch (xvgf)
145 {
146 case XvgFormat::Xmgrace: sprintf(buf + b, "%s", "\\f{}"); break;
147 case XvgFormat::Xmgr: sprintf(buf + b, "%s", "\\4"); break;
148 default: buf[b] = '\0'; break;
149 }
150 g++;
151 b = strlen(buf);
152 }
153 else if (gmx[g] == '8')
154 {
155 /* Backward compatibility for xmgr symbol font "\8" */
156 switch (xvgf)
157 {
158 case XvgFormat::Xmgrace: sprintf(buf + b, "%s", "\\x"); break;
159 case XvgFormat::Xmgr: sprintf(buf + b, "%s", "\\8"); break;
160 default: buf[b] = '\0'; break;
161 }
162 g++;
163 b = std::strlen(buf);
164 }
165 else
166 {
167 /* Check for special symbol */
168 i = 0;
169 while (sym[i] != nullptr && gmx_strncasecmp(sym[i], &gmx[g], std::strlen(sym[i])) != 0)
170 {
171 i++;
172 }
173 if (sym[i] != nullptr)
174 {
175 c = symc[i];
176 if (std::isupper(gmx[g]))
177 {
178 c = std::toupper(c);
179 }
180 switch (xvgf)
181 {
182 case XvgFormat::Xmgrace:
183 sprintf(buf + b, "%s%c%s", "\\x", c, "\\f{}");
184 break;
185 case XvgFormat::Xmgr: sprintf(buf + b, "%s%c%s", "\\8", c, "\\4"); break;
186 default:
187 std::strncat(buf + b, &gmx[g], std::strlen(sym[i]));
188 b += std::strlen(sym[i]);
189 if (gmx[g + std::strlen(sym[i])] != ' ')
190 {
191 buf[b++] = ' ';
192 }
193 buf[b] = '\0';
194 break;
195 }
196 g += std::strlen(sym[i]);
197 b = std::strlen(buf);
198 }
199 else
200 {
201 /* Unknown escape sequence, this should not happen.
202 * We do not generate a fatal error, since that might
203 * stop otherwise functioning code from working.
204 * Copy the backslash to the output and continue processing.
205 */
206 buf[b++] = '\\';
207 }
208 }
209 }
210 else
211 {
212 buf[b++] = gmx[g++];
213 }
214 }
215
216 buf[b++] = '\0';
217
218 return buf;
219 }
220
xvgr_header(FILE * fp,const char * title,const std::string & xaxis,const std::string & yaxis,int exvg_graph_type,const gmx_output_env_t * oenv)221 void xvgr_header(FILE* fp,
222 const char* title,
223 const std::string& xaxis,
224 const std::string& yaxis,
225 int exvg_graph_type,
226 const gmx_output_env_t* oenv)
227 {
228 char buf[STRLEN];
229
230 if (output_env_get_print_xvgr_codes(oenv))
231 {
232 fprintf(fp, "# This file was created %s", gmx_format_current_time().c_str());
233 try
234 {
235 gmx::BinaryInformationSettings settings;
236 settings.generatedByHeader(true);
237 settings.linePrefix("# ");
238 gmx::printBinaryInformation(fp, output_env_get_program_context(oenv), settings);
239 }
240 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
241 fprintf(fp, "# %s is part of G R O M A C S:\n#\n", output_env_get_program_display_name(oenv));
242 fprintf(fp, "# %s\n#\n", gmx::bromacs().c_str());
243 fprintf(fp, "@ title \"%s\"\n", xvgrstr(title, oenv, buf, STRLEN));
244 fprintf(fp, "@ xaxis label \"%s\"\n", xvgrstr(xaxis, oenv, buf, STRLEN));
245 fprintf(fp, "@ yaxis label \"%s\"\n", xvgrstr(yaxis, oenv, buf, STRLEN));
246 switch (exvg_graph_type)
247 {
248 case exvggtXNY:
249 if (output_env_get_xvg_format(oenv) == XvgFormat::Xmgr)
250 {
251 fprintf(fp, "@TYPE nxy\n");
252 }
253 else
254 {
255 fprintf(fp, "@TYPE xy\n");
256 }
257 break;
258 case exvggtXYDY: fprintf(fp, "@TYPE xydy\n"); break;
259 case exvggtXYDYDY: fprintf(fp, "@TYPE xydydy\n"); break;
260 }
261 }
262 }
263
xvgropen_type(const char * fn,const char * title,const std::string & xaxis,const std::string & yaxis,int exvg_graph_type,const gmx_output_env_t * oenv)264 FILE* xvgropen_type(const char* fn,
265 const char* title,
266 const std::string& xaxis,
267 const std::string& yaxis,
268 int exvg_graph_type,
269 const gmx_output_env_t* oenv)
270 {
271 FILE* fp;
272
273 fp = gmx_fio_fopen(fn, "w");
274
275 xvgr_header(fp, title, xaxis, yaxis, exvg_graph_type, oenv);
276
277 return fp;
278 }
279
xvgropen(const char * fn,const char * title,const std::string & xaxis,const std::string & yaxis,const gmx_output_env_t * oenv)280 FILE* xvgropen(const char* fn,
281 const char* title,
282 const std::string& xaxis,
283 const std::string& yaxis,
284 const gmx_output_env_t* oenv)
285 {
286 return xvgropen_type(fn, title, xaxis, yaxis, exvggtXNY, oenv);
287 }
288
xvgrclose(FILE * fp)289 void xvgrclose(FILE* fp)
290 {
291 gmx_fio_fclose(fp);
292 }
293
xvgr_subtitle(FILE * out,const char * subtitle,const gmx_output_env_t * oenv)294 void xvgr_subtitle(FILE* out, const char* subtitle, const gmx_output_env_t* oenv)
295 {
296 char buf[STRLEN];
297
298 if (output_env_get_print_xvgr_codes(oenv))
299 {
300 fprintf(out, "@ subtitle \"%s\"\n", xvgrstr(subtitle, oenv, buf, STRLEN));
301 }
302 }
303
xvgr_view(FILE * out,real xmin,real ymin,real xmax,real ymax,const gmx_output_env_t * oenv)304 void xvgr_view(FILE* out, real xmin, real ymin, real xmax, real ymax, const gmx_output_env_t* oenv)
305 {
306 if (output_env_get_print_xvgr_codes(oenv))
307 {
308 fprintf(out, "@ view %g, %g, %g, %g\n", xmin, ymin, xmax, ymax);
309 }
310 }
311
xvgr_world(FILE * out,real xmin,real ymin,real xmax,real ymax,const gmx_output_env_t * oenv)312 void xvgr_world(FILE* out, real xmin, real ymin, real xmax, real ymax, const gmx_output_env_t* oenv)
313 {
314 if (output_env_get_print_xvgr_codes(oenv))
315 {
316 fprintf(out,
317 "@ world xmin %g\n"
318 "@ world ymin %g\n"
319 "@ world xmax %g\n"
320 "@ world ymax %g\n",
321 xmin, ymin, xmax, ymax);
322 }
323 }
324
stringIsEmpty(const std::string & s)325 static bool stringIsEmpty(const std::string& s)
326 {
327 return s.empty();
328 }
329
stringIsEmpty(const char * s)330 static bool stringIsEmpty(const char* s)
331 {
332 return (s == nullptr || s[0] == '\0');
333 }
334
335 template<typename T>
xvgr_legend(FILE * out,int nsets,const T * setname,const gmx_output_env_t * oenv)336 static void xvgr_legend(FILE* out, int nsets, const T* setname, const gmx_output_env_t* oenv)
337 {
338 int i;
339 char buf[STRLEN];
340
341 if (output_env_get_print_xvgr_codes(oenv))
342 {
343 xvgr_view(out, 0.15, 0.15, 0.75, 0.85, oenv);
344 fprintf(out, "@ legend on\n");
345 fprintf(out, "@ legend box on\n");
346 fprintf(out, "@ legend loctype view\n");
347 fprintf(out, "@ legend %g, %g\n", 0.78, 0.8);
348 fprintf(out, "@ legend length %d\n", 2);
349 for (i = 0; (i < nsets); i++)
350 {
351 if (!stringIsEmpty(setname[i]))
352 {
353 if (output_env_get_xvg_format(oenv) == XvgFormat::Xmgr)
354 {
355 fprintf(out, "@ legend string %d \"%s\"\n", i, xvgrstr(setname[i], oenv, buf, STRLEN));
356 }
357 else
358 {
359 fprintf(out, "@ s%d legend \"%s\"\n", i, xvgrstr(setname[i], oenv, buf, STRLEN));
360 }
361 }
362 }
363 }
364 }
365
xvgrLegend(FILE * out,const std::vector<std::string> & setNames,const struct gmx_output_env_t * oenv)366 void xvgrLegend(FILE* out, const std::vector<std::string>& setNames, const struct gmx_output_env_t* oenv)
367 {
368 xvgr_legend(out, setNames.size(), setNames.data(), oenv);
369 }
xvgr_legend(FILE * out,int nsets,const char * const * setnames,const struct gmx_output_env_t * oenv)370 void xvgr_legend(FILE* out, int nsets, const char* const* setnames, const struct gmx_output_env_t* oenv)
371 {
372 xvgr_legend<const char*>(out, nsets, setnames, oenv);
373 }
374
xvgr_new_dataset(FILE * out,int nr_first,int nsets,const char ** setname,const gmx_output_env_t * oenv)375 void xvgr_new_dataset(FILE* out, int nr_first, int nsets, const char** setname, const gmx_output_env_t* oenv)
376 {
377 int i;
378 char buf[STRLEN];
379
380 if (output_env_get_print_xvgr_codes(oenv))
381 {
382 fprintf(out, "@\n");
383 for (i = 0; (i < nsets); i++)
384 {
385 if (setname[i])
386 {
387 if (output_env_get_xvg_format(oenv) == XvgFormat::Xmgr)
388 {
389 fprintf(out, "@ legend string %d \"%s\"\n", i + nr_first,
390 xvgrstr(setname[i], oenv, buf, STRLEN));
391 }
392 else
393 {
394 fprintf(out, "@ s%d legend \"%s\"\n", i + nr_first,
395 xvgrstr(setname[i], oenv, buf, STRLEN));
396 }
397 }
398 }
399 }
400 else
401 {
402 fprintf(out, "\n");
403 }
404 }
405
xvgr_line_props(FILE * out,int NrSet,int LineStyle,int LineColor,const gmx_output_env_t * oenv)406 void xvgr_line_props(FILE* out, int NrSet, int LineStyle, int LineColor, const gmx_output_env_t* oenv)
407 {
408 if (output_env_get_print_xvgr_codes(oenv))
409 {
410 fprintf(out, "@ with g0\n");
411 fprintf(out, "@ s%d linestyle %d\n", NrSet, LineStyle);
412 fprintf(out, "@ s%d color %d\n", NrSet, LineColor);
413 }
414 }
415
416 static const char* LocTypeStr[] = { "view", "world" };
417 static const char* BoxFillStr[] = { "none", "color", "pattern" };
418
xvgr_box(FILE * out,int LocType,real xmin,real ymin,real xmax,real ymax,int LineStyle,int LineWidth,int LineColor,int BoxFill,int BoxColor,int BoxPattern,const gmx_output_env_t * oenv)419 void xvgr_box(FILE* out,
420 int LocType,
421 real xmin,
422 real ymin,
423 real xmax,
424 real ymax,
425 int LineStyle,
426 int LineWidth,
427 int LineColor,
428 int BoxFill,
429 int BoxColor,
430 int BoxPattern,
431 const gmx_output_env_t* oenv)
432 {
433 if (output_env_get_print_xvgr_codes(oenv))
434 {
435 fprintf(out, "@with box\n");
436 fprintf(out, "@ box on\n");
437 fprintf(out, "@ box loctype %s\n", LocTypeStr[LocType]);
438 fprintf(out, "@ box %g, %g, %g, %g\n", xmin, ymin, xmax, ymax);
439 fprintf(out, "@ box linestyle %d\n", LineStyle);
440 fprintf(out, "@ box linewidth %d\n", LineWidth);
441 fprintf(out, "@ box color %d\n", LineColor);
442 fprintf(out, "@ box fill %s\n", BoxFillStr[BoxFill]);
443 fprintf(out, "@ box fill color %d\n", BoxColor);
444 fprintf(out, "@ box fill pattern %d\n", BoxPattern);
445 fprintf(out, "@box def\n");
446 }
447 }
448
449 /* reads a line into ptr, adjusting len and renewing ptr if neccesary */
fgets3(FILE * fp,char ** ptr,int * len,int maxlen)450 static char* fgets3(FILE* fp, char** ptr, int* len, int maxlen)
451 {
452 int len_remaining = *len; /* remaining amount of allocated bytes in buf */
453 int curp = 0; /* current position in buf to read into */
454
455 do
456 {
457 if (len_remaining < 2)
458 {
459 if (*len + STRLEN < maxlen)
460 {
461 /* This line is longer than len characters, let's increase len! */
462 *len += STRLEN;
463 len_remaining += STRLEN;
464 srenew(*ptr, *len);
465 }
466 else
467 {
468 /*something is wrong, we'll just keep reading and return NULL*/
469 len_remaining = STRLEN;
470 curp = 0;
471 }
472 }
473 if (fgets(*ptr + curp, len_remaining, fp) == nullptr)
474 {
475 /* if last line, skip */
476 return nullptr;
477 }
478 curp += len_remaining - 1; /* overwrite the nul char in next iteration */
479 len_remaining = 1;
480 } while ((std::strchr(*ptr, '\n') == nullptr) && (feof(fp) == 0));
481
482 if (*len + STRLEN >= maxlen)
483 {
484 return nullptr; /* this line was too long */
485 }
486
487 if (feof(fp))
488 {
489 /* We reached EOF before '\n', skip this last line. */
490 return nullptr;
491 }
492 {
493 /* now remove newline */
494 int slen = std::strlen(*ptr);
495 if ((*ptr)[slen - 1] == '\n')
496 {
497 (*ptr)[slen - 1] = '\0';
498 }
499 }
500
501 return *ptr;
502 }
503
wordcount(char * ptr)504 static int wordcount(char* ptr)
505 {
506 int i, n = 0, is[2];
507 int cur = 0;
508 #define prev (1 - cur)
509
510 if (nullptr != ptr)
511 {
512 for (i = 0; (ptr[i] != '\0'); i++)
513 {
514 is[cur] = std::isspace(ptr[i]);
515 if (((0 == i) && !is[cur]) || ((i > 0) && (!is[cur] && is[prev])))
516 {
517 n++;
518 }
519 cur = prev;
520 }
521 }
522 return n;
523 }
524
read_xvgr_string(const char * line)525 static char* read_xvgr_string(const char* line)
526 {
527 const char *ptr0, *ptr1;
528 char* str;
529
530 ptr0 = std::strchr(line, '"');
531 if (ptr0 != nullptr)
532 {
533 ptr0++;
534 ptr1 = std::strchr(ptr0, '"');
535 if (ptr1 != nullptr)
536 {
537 str = gmx_strdup(ptr0);
538 str[ptr1 - ptr0] = '\0';
539 }
540 else
541 {
542 str = gmx_strdup("");
543 }
544 }
545 else
546 {
547 str = gmx_strdup("");
548 }
549
550 return str;
551 }
552
read_xvg_legend(const char * fn,double *** y,int * ny,char ** subtitle,char *** legend)553 int read_xvg_legend(const char* fn, double*** y, int* ny, char** subtitle, char*** legend)
554 {
555 FILE* fp;
556 char* ptr;
557 char* base = nullptr;
558 char* fmt = nullptr;
559 int k, line = 0, nny, nx, maxx, rval, legend_nalloc, set, nchar;
560 double lf;
561 double** yy = nullptr;
562 char* tmpbuf;
563 int len = STRLEN;
564 *ny = 0;
565 nny = 0;
566 nx = 0;
567 maxx = 0;
568 fp = gmx_fio_fopen(fn, "r");
569
570 snew(tmpbuf, len);
571 if (subtitle != nullptr)
572 {
573 *subtitle = nullptr;
574 }
575 legend_nalloc = 0;
576 if (legend != nullptr)
577 {
578 *legend = nullptr;
579 }
580
581 while ((ptr = fgets3(fp, &tmpbuf, &len, 10 * STRLEN)) != nullptr && ptr[0] != '&')
582 {
583 line++;
584 trim(ptr);
585 if (ptr[0] == '@')
586 {
587 if (legend != nullptr)
588 {
589 ptr++;
590 trim(ptr);
591 set = -1;
592 if (std::strncmp(ptr, "subtitle", 8) == 0)
593 {
594 ptr += 8;
595 if (subtitle != nullptr)
596 {
597 *subtitle = read_xvgr_string(ptr);
598 }
599 }
600 else if (std::strncmp(ptr, "legend string", 13) == 0)
601 {
602 ptr += 13;
603 sscanf(ptr, "%d%n", &set, &nchar);
604 ptr += nchar;
605 }
606 else if (ptr[0] == 's')
607 {
608 ptr++;
609 sscanf(ptr, "%d%n", &set, &nchar);
610 ptr += nchar;
611 trim(ptr);
612 if (std::strncmp(ptr, "legend", 6) == 0)
613 {
614 ptr += 6;
615 }
616 else
617 {
618 set = -1;
619 }
620 }
621 if (set >= 0)
622 {
623 if (set >= legend_nalloc)
624 {
625 legend_nalloc = set + 1;
626 srenew(*legend, legend_nalloc);
627 (*legend)[set] = read_xvgr_string(ptr);
628 }
629 }
630 }
631 }
632 else if (ptr[0] != '#')
633 {
634 if (nny == 0)
635 {
636 (*ny) = nny = wordcount(ptr);
637 /* fprintf(stderr,"There are %d columns in your file\n",nny);*/
638 if (nny == 0)
639 {
640 return 0;
641 }
642 snew(yy, nny);
643 snew(fmt, 3 * nny + 1);
644 snew(base, 3 * nny + 1);
645 }
646 /* Allocate column space */
647 if (nx >= maxx)
648 {
649 maxx += 1024;
650 for (k = 0; (k < nny); k++)
651 {
652 srenew(yy[k], maxx);
653 }
654 }
655 /* Initiate format string */
656 fmt[0] = '\0';
657 base[0] = '\0';
658
659 /* fprintf(stderr,"ptr='%s'\n",ptr);*/
660 for (k = 0; (k < nny); k++)
661 {
662 std::strcpy(fmt, base);
663 std::strcat(fmt, "%lf");
664 rval = sscanf(ptr, fmt, &lf);
665 /* fprintf(stderr,"rval = %d\n",rval);*/
666 if ((rval == EOF) || (rval == 0))
667 {
668 break;
669 }
670 yy[k][nx] = lf;
671 srenew(fmt, 3 * (nny + 1) + 1);
672 srenew(base, 3 * nny + 1);
673 std::strcat(base, "%*s");
674 }
675 if (k != nny)
676 {
677 fprintf(stderr, "Only %d columns on line %d in file %s\n", k, line, fn);
678 for (; (k < nny); k++)
679 {
680 yy[k][nx] = 0.0;
681 }
682 }
683 nx++;
684 }
685 }
686 gmx_fio_fclose(fp);
687
688 *y = yy;
689 sfree(tmpbuf);
690 sfree(base);
691 sfree(fmt);
692
693 if (legend_nalloc > 0)
694 {
695 if (*ny - 1 > legend_nalloc)
696 {
697 assert(legend);
698 srenew(*legend, *ny - 1);
699 for (set = legend_nalloc; set < *ny - 1; set++)
700 {
701 (*legend)[set] = nullptr;
702 }
703 }
704 }
705
706 return nx;
707 }
708
read_xvg(const char * fn,double *** y,int * ny)709 int read_xvg(const char* fn, double*** y, int* ny)
710 {
711 gmx::MultiDimArray<std::vector<double>, gmx::dynamicExtents2D> xvgData =
712 readXvgData(std::string(fn));
713
714 int numColumns = xvgData.extent(0);
715 int numRows = xvgData.extent(1);
716
717 double** yy = nullptr;
718 snew(yy, numColumns);
719 for (int column = 0; column < numColumns; column++)
720 {
721 snew(yy[column], numRows);
722 for (int row = 0; row < numRows; row++)
723 {
724 yy[column][row] = xvgData.asConstView()[column][row];
725 }
726 }
727
728 *y = yy;
729 *ny = numColumns;
730 int nx = numRows;
731 return nx;
732 }
733
readXvgData(const std::string & fn)734 gmx::MultiDimArray<std::vector<double>, gmx::dynamicExtents2D> readXvgData(const std::string& fn)
735 {
736 FILE* fp = gmx_fio_fopen(fn.c_str(), "r");
737 char* ptr;
738 char* base = nullptr;
739 char* fmt = nullptr;
740 char* tmpbuf;
741 int len = STRLEN;
742
743 //! This only gets properly set after the first line of data is read
744 int numColumns = 0;
745 int numRows = 0;
746 snew(tmpbuf, len);
747 std::vector<double> xvgData;
748
749 for (int line = 0; (ptr = fgets3(fp, &tmpbuf, &len, 10 * STRLEN)) != nullptr && ptr[0] != '&'; ++line)
750 {
751 trim(ptr);
752 if (ptr[0] == '@' || ptr[0] == '#')
753 {
754 continue;
755 }
756 ++numRows;
757 if (numColumns == 0)
758 {
759 numColumns = wordcount(ptr);
760 if (numColumns == 0)
761 {
762 return {}; // There are no columns and hence no data to process
763 }
764 snew(fmt, 3 * numColumns + 1);
765 snew(base, 3 * numColumns + 1);
766 }
767 /* Initiate format string */
768 fmt[0] = '\0';
769 base[0] = '\0';
770 int columnCount = 0;
771 for (columnCount = 0; (columnCount < numColumns); columnCount++)
772 {
773 double lf;
774 std::strcpy(fmt, base);
775 std::strcat(fmt, "%lf");
776 int rval = sscanf(ptr, fmt, &lf);
777 if ((rval == EOF) || (rval == 0))
778 {
779 break;
780 }
781 xvgData.push_back(lf);
782 srenew(fmt, 3 * (numColumns + 1) + 1);
783 srenew(base, 3 * numColumns + 1);
784 std::strcat(base, "%*s");
785 }
786
787 if (columnCount != numColumns)
788 {
789 fprintf(stderr, "Only %d columns on line %d in file %s\n", columnCount, line, fn.c_str());
790 for (; (columnCount < numColumns); columnCount++)
791 {
792 xvgData.push_back(0.0);
793 }
794 }
795 }
796 gmx_fio_fclose(fp);
797
798 sfree(tmpbuf);
799 sfree(base);
800 sfree(fmt);
801
802 gmx::MultiDimArray<std::vector<double>, gmx::dynamicExtents2D> xvgDataAsArray(numRows, numColumns);
803 std::copy(std::begin(xvgData), std::end(xvgData), begin(xvgDataAsArray.asView()));
804
805 gmx::MultiDimArray<std::vector<double>, gmx::dynamicExtents2D> xvgDataAsArrayTransposed(
806 numColumns, numRows);
807 for (std::ptrdiff_t row = 0; row < numRows; ++row)
808 {
809 for (std::ptrdiff_t column = 0; column < numColumns; ++column)
810 {
811 xvgDataAsArrayTransposed(column, row) = xvgDataAsArray(row, column);
812 }
813 }
814
815 return xvgDataAsArrayTransposed;
816 }
817
write_xvg(const char * fn,const char * title,int nx,int ny,real ** y,const char ** leg,const gmx_output_env_t * oenv)818 void write_xvg(const char* fn, const char* title, int nx, int ny, real** y, const char** leg, const gmx_output_env_t* oenv)
819 {
820 FILE* fp;
821 int i, j;
822
823 fp = xvgropen(fn, title, "X", "Y", oenv);
824 if (leg)
825 {
826 xvgr_legend(fp, ny - 1, leg, oenv);
827 }
828 for (i = 0; (i < nx); i++)
829 {
830 for (j = 0; (j < ny); j++)
831 {
832 fprintf(fp, " %12.5e", y[j][i]);
833 }
834 fprintf(fp, "\n");
835 }
836 xvgrclose(fp);
837 }
838
read_xvg_time(const char * fn,gmx_bool bHaveT,gmx_bool bTB,real tb,gmx_bool bTE,real te,int nsets_in,int * nset,int * nval,real * dt,real ** t)839 real** read_xvg_time(const char* fn,
840 gmx_bool bHaveT,
841 gmx_bool bTB,
842 real tb,
843 gmx_bool bTE,
844 real te,
845 int nsets_in,
846 int* nset,
847 int* nval,
848 real* dt,
849 real** t)
850 {
851 FILE* fp;
852 #define MAXLINELEN 16384
853 char line0[MAXLINELEN];
854 char* line;
855 int t_nalloc, *val_nalloc, a, narg, n, sin, set, nchar;
856 double dbl;
857 gmx_bool bEndOfSet, bTimeInRange, bFirstLine = TRUE;
858 real** val;
859
860 t_nalloc = 0;
861 *t = nullptr;
862 val = nullptr;
863 val_nalloc = nullptr;
864 *nset = 0;
865 *dt = 0;
866 fp = gmx_fio_fopen(fn, "r");
867 for (sin = 0; sin < nsets_in; sin++)
868 {
869 if (nsets_in == 1)
870 {
871 narg = 0;
872 }
873 else
874 {
875 narg = bHaveT ? 2 : 1;
876 }
877 n = 0;
878 bEndOfSet = FALSE;
879 while (!bEndOfSet && fgets(line0, MAXLINELEN, fp))
880 {
881 line = line0;
882 /* Remove whitespace */
883 while (line[0] == ' ' || line[0] == '\t')
884 {
885 line++;
886 }
887 bEndOfSet = (line[0] == '&');
888 if (line[0] != '#' && line[0] != '@' && line[0] != '\n' && !bEndOfSet)
889 {
890 if (bFirstLine && bHaveT)
891 {
892 /* Check the first line that should contain data */
893 a = sscanf(line, "%lf%lf", &dbl, &dbl);
894 if (a == 0)
895 {
896 gmx_fatal(FARGS, "Expected a number in %s on line:\n%s", fn, line0);
897 }
898 else if (a == 1)
899 {
900 fprintf(stderr,
901 "Found only 1 number on line, "
902 "assuming no time is present.\n");
903 bHaveT = FALSE;
904 if (nsets_in > 1)
905 {
906 narg = 1;
907 }
908 }
909 }
910
911 a = 0;
912 bTimeInRange = TRUE;
913 while ((a < narg || (nsets_in == 1 && n == 0))
914 && sscanf(line, "%lf%n", &dbl, &nchar) == 1 && bTimeInRange)
915 {
916 /* Use set=-1 as the time "set" */
917 if (sin)
918 {
919 if (!bHaveT || (a > 0))
920 {
921 set = sin;
922 }
923 else
924 {
925 set = -1;
926 }
927 }
928 else
929 {
930 if (!bHaveT)
931 {
932 set = a;
933 }
934 else
935 {
936 set = a - 1;
937 }
938 }
939 if (set == -1 && ((bTB && dbl < tb) || (bTE && dbl > te)))
940 {
941 bTimeInRange = FALSE;
942 }
943
944 if (bTimeInRange)
945 {
946 if (n == 0)
947 {
948 if (nsets_in == 1)
949 {
950 narg++;
951 }
952 if (set >= 0)
953 {
954 *nset = set + 1;
955 srenew(val, *nset);
956 srenew(val_nalloc, *nset);
957 val_nalloc[set] = 0;
958 val[set] = nullptr;
959 }
960 }
961 if (set == -1)
962 {
963 if (sin == 0)
964 {
965 if (n >= t_nalloc)
966 {
967 t_nalloc = over_alloc_small(n);
968 srenew(*t, t_nalloc);
969 }
970 (*t)[n] = dbl;
971 }
972 /* else we should check the time of the next sets with set 0 */
973 }
974 else
975 {
976 if (n >= val_nalloc[set])
977 {
978 val_nalloc[set] = over_alloc_small(n);
979 srenew(val[set], val_nalloc[set]);
980 }
981 val[set][n] = static_cast<real>(dbl);
982 }
983 }
984 a++;
985 line += nchar;
986 }
987 if (line0[strlen(line0) - 1] != '\n')
988 {
989 fprintf(stderr, "File %s does not end with a newline, ignoring the last line\n", fn);
990 }
991 else if (bTimeInRange)
992 {
993 if (a == 0)
994 {
995 fprintf(stderr, "Ignoring invalid line in %s:\n%s", fn, line0);
996 }
997 else
998 {
999 if (a != narg)
1000 {
1001 fprintf(stderr,
1002 "Invalid line in %s:\n%s"
1003 "Using zeros for the last %d sets\n",
1004 fn, line0, narg - a);
1005 }
1006 n++;
1007 }
1008 }
1009 if (a > 0)
1010 {
1011 bFirstLine = FALSE;
1012 }
1013 }
1014 }
1015 if (sin == 0)
1016 {
1017 *nval = n;
1018 if (!bHaveT)
1019 {
1020 snew(*t, n);
1021 for (a = 0; a < n; a++)
1022 {
1023 (*t)[a] = a;
1024 }
1025 }
1026 if (n > 1)
1027 {
1028 *dt = static_cast<real>((*t)[n - 1] - (*t)[0]) / (n - 1.0);
1029 }
1030 else
1031 {
1032 *dt = 1;
1033 }
1034 }
1035 else
1036 {
1037 if (n < *nval)
1038 {
1039 fprintf(stderr, "Set %d is shorter (%d) than the previous set (%d)\n", sin + 1, n, *nval);
1040 *nval = n;
1041 fprintf(stderr, "Will use only the first %d points of every set\n", *nval);
1042 }
1043 }
1044 }
1045 gmx_fio_fclose(fp);
1046
1047 sfree(val_nalloc);
1048
1049 return val;
1050 }
1051