Modernize t_matrix, t_mapping, t_rgb
authorMark Abraham <mark.j.abraham@gmail.com>
Mon, 12 Aug 2019 09:00:55 +0000 (11:00 +0200)
committerMark Abraham <mark.j.abraham@gmail.com>
Mon, 19 Aug 2019 19:04:16 +0000 (21:04 +0200)
This supports removing some manual string handling that regularly
troubles e.g. the static analyzer.

Made more use of vector, ArrayRef and ssize.

Stopped installing matio.h, because the transformation introduced
here was going to lead to many more files installed, and we anyway
plan to remove the serendipitous API that we have.

Refs #2899

Change-Id: Ia4c6a37d4f25dc2de93b06b2ffc46e2ddab96f39

14 files changed:
src/gromacs/fileio/CMakeLists.txt
src/gromacs/fileio/matio.cpp
src/gromacs/fileio/matio.h
src/gromacs/fileio/rgb.h
src/gromacs/gmxana/gmx_cluster.cpp
src/gromacs/gmxana/gmx_do_dssp.cpp
src/gromacs/gmxana/gmx_hbond.cpp
src/gromacs/gmxana/gmx_xpm2ps.cpp
src/gromacs/math/densityfit.cpp
src/gromacs/math/densityfittingforce.cpp
src/gromacs/math/gausstransform.cpp
src/gromacs/math/gausstransform.h
src/gromacs/math/matrix.h
src/gromacs/math/multidimarray.h

index 2b6194ca8e16c9fbc620ec79488b4c52a0ccf141..b330d2d41e0681f2a3d5db0f14d98d94de3bca03 100644 (file)
@@ -51,7 +51,6 @@ gmx_install_headers(
     enxio.h
     filetypes.h
     gmxfio.h
-    matio.h
     mtxio.h
     oenv.h
     pdbio.h
index 129ae67d63d91155491d8f34ddc31ceb363838ab..2eb0897281128df0155e3f6ed1d226219bdd9fdd 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2017,2018,2019, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -44,6 +44,7 @@
 #include <cstring>
 
 #include <algorithm>
+#include <regex>
 #include <string>
 
 #include "gromacs/fileio/gmxfio.h"
@@ -57,6 +58,8 @@
 #include "gromacs/utility/programcontext.h"
 #include "gromacs/utility/smalloc.h"
 
+using namespace gmx;
+
 static const char mapper[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789!@#$%^&*()-_=+{}|;:',<.>/?";
 #define NMAP static_cast<long int>(strlen(mapper))
 
@@ -100,33 +103,27 @@ void done_matrix(int nx, real ***m)
     *m = nullptr;
 }
 
-gmx_bool matelmt_cmp(t_xpmelmt e1, t_xpmelmt e2)
+static bool operator==(t_xpmelmt e1, t_xpmelmt e2)
 {
     return (e1.c1 == e2.c1) && (e1.c2 == e2.c2);
 }
 
-t_matelmt searchcmap(int n, t_mapping map[], t_xpmelmt c)
+//! Return the index of the first element that matches \c c, or -1 if not found.
+t_matelmt searchcmap(ArrayRef<const t_mapping> map, t_xpmelmt c)
 {
-    int i;
-
-    for (i = 0; (i < n); i++)
-    {
-        if (matelmt_cmp(map[i].code, c))
-        {
-            return i;
-        }
-    }
-
-    return -1;
+    auto findIt = std::find_if(map.begin(), map.end(), [&c](const auto &m)
+                               { return (m.code == c); });
+    return (findIt == map.end()) ? -1 : (findIt - map.begin());
 }
 
-int getcmap(FILE *in, const char *fn, t_mapping **map)
+//! Read the mapping table from in, return number of entries
+static std::vector<t_mapping> getcmap(FILE *in, const char *fn)
 {
-    int        i, n;
-    char       line[STRLEN];
-    char       code[STRLEN], desc[STRLEN];
-    double     r, g, b;
-    t_mapping *m;
+    int                    i, n;
+    char                   line[STRLEN];
+    char                   code[STRLEN], desc[STRLEN];
+    double                 r, g, b;
+    std::vector<t_mapping> m;
 
     if (fgets2(line, STRLEN-1, in) == nullptr)
     {
@@ -134,7 +131,7 @@ int getcmap(FILE *in, const char *fn, t_mapping **map)
                   "(just wanted to read number of entries)", fn);
     }
     sscanf(line, "%d", &n);
-    snew(m, n);
+    m.resize(n);
     for (i = 0; (i < n); i++)
     {
         if (fgets2(line, STRLEN-1, in) == nullptr)
@@ -150,15 +147,14 @@ int getcmap(FILE *in, const char *fn, t_mapping **map)
         m[i].rgb.g   = g;
         m[i].rgb.b   = b;
     }
-    *map = m;
 
-    return n;
+    return m;
 }
 
-int readcmap(const char *fn, t_mapping **map)
+std::vector<t_mapping> readcmap(const char *fn)
 {
-    gmx::FilePtr in = gmx::openLibraryFile(fn);
-    return getcmap(in.get(), fn, map);
+    FilePtr in = openLibraryFile(fn);
+    return getcmap(in.get(), fn);
 }
 
 void printcmap(FILE *out, int n, t_mapping map[])
@@ -254,50 +250,43 @@ static char *line2string(char **line)
     return *line;
 }
 
-static void parsestring(char *line, const char *label, char *string)
+//! If a label named \c label is found in \c line, return it. Otherwise return empty string.
+static std::string findLabelInLine(const std::string &line,
+                                   const std::string &label)
 {
-    if (strstr(line, label))
+    std::regex  re(".*" + label + "\"(.*)\"");
+    std::smatch match;
+    if (std::regex_search(line, match, re) && match.size() > 1)
     {
-        if (strstr(line, label) < strchr(line, '\"'))
-        {
-            line2string(&line);
-            strcpy(string, line);
-        }
+        return match.str(1);
     }
+    return std::string();
 }
 
-static void read_xpm_entry(FILE *in, t_matrix *mm)
+//! Read and return a matrix from \c in
+static t_matrix read_xpm_entry(FILE *in)
 {
-    t_mapping   *map;
     char        *line_buf = nullptr, *line = nullptr, *str, buf[256] = {0};
-    int          i, m, col_len, nch = 0, n_axis_x, n_axis_y, llmax;
+    int          i, m, col_len, nch = 0, llmax;
     int          llalloc = 0;
     unsigned int r, g, b;
     double       u;
     gmx_bool     bGetOnWithIt, bSetLine;
     t_xpmelmt    c;
 
-    mm->flags      = 0;
-    mm->title[0]   = 0;
-    mm->legend[0]  = 0;
-    mm->label_x[0] = 0;
-    mm->label_y[0] = 0;
-    mm->matrix     = nullptr;
-    mm->axis_x     = nullptr;
-    mm->axis_y     = nullptr;
-    mm->bDiscrete  = FALSE;
+    t_matrix     mm;
 
     llmax = STRLEN;
 
     while ((nullptr != fgetline(&line_buf, llmax, &llalloc, in)) &&
            (std::strncmp(line_buf, "static", 6) != 0))
     {
-        line = line_buf;
-        parsestring(line, "title", (mm->title));
-        parsestring(line, "legend", (mm->legend));
-        parsestring(line, "x-label", (mm->label_x));
-        parsestring(line, "y-label", (mm->label_y));
-        parsestring(line, "type", buf);
+        std::string lineString = line_buf;
+        mm.title   = findLabelInLine(lineString, "title");
+        mm.legend  = findLabelInLine(lineString, "legend");
+        mm.label_x = findLabelInLine(lineString, "x-label");
+        mm.label_y = findLabelInLine(lineString, "y-label");
+        findLabelInLine(lineString, "type"); // discard the returned string
     }
 
     if (!line || strncmp(line, "static", 6) != 0)
@@ -307,16 +296,17 @@ static void read_xpm_entry(FILE *in, t_matrix *mm)
 
     if (buf[0] && (gmx_strcasecmp(buf, "Discrete") == 0))
     {
-        mm->bDiscrete = TRUE;
+        mm.bDiscrete = TRUE;
     }
 
     if (debug)
     {
         fprintf(debug, "%s %s %s %s\n",
-                mm->title, mm->legend, mm->label_x, mm->label_y);
+                mm.title.c_str(), mm.legend.c_str(), mm.label_x.c_str(), mm.label_y.c_str());
     }
 
     /* Read sizes */
+    int nmap = 0;
     bGetOnWithIt = FALSE;
     while (!bGetOnWithIt && (nullptr != fgetline(&line_buf, llmax, &llalloc, in)))
     {
@@ -329,23 +319,23 @@ static void read_xpm_entry(FILE *in, t_matrix *mm)
         if  (line[0] == '\"')
         {
             line2string(&line);
-            sscanf(line, "%d %d %d %d", &(mm->nx), &(mm->ny), &(mm->nmap), &nch);
+            sscanf(line, "%d %d %d %d", &(mm.nx), &(mm.ny), &nmap, &nch);
             if (nch > 2)
             {
-                gmx_fatal(FARGS, "Sorry can only read xpm's with at most 2 caracters per pixel\n");
+                gmx_fatal(FARGS, "Sorry can only read xpm's with at most 2 characters per pixel\n");
             }
-            if (mm->nx <= 0 || mm->ny <= 0)
+            if (mm.nx <= 0 || mm.ny <= 0)
             {
                 gmx_fatal(FARGS, "Dimensions of xpm-file have to be larger than 0\n");
             }
-            llmax        = std::max(STRLEN, mm->nx+10);
+            llmax        = std::max(STRLEN, mm.nx+10);
             bGetOnWithIt = TRUE;
         }
     }
     if (debug)
     {
-        fprintf(debug, "mm->nx %d mm->ny %d mm->nmap %d nch %d\n",
-                mm->nx, mm->ny, mm->nmap, nch);
+        fprintf(debug, "mm.nx %d mm.ny %d nmap %d nch %d\n",
+                mm.nx, mm.ny, nmap, nch);
     }
     if (nch == 0)
     {
@@ -353,23 +343,23 @@ static void read_xpm_entry(FILE *in, t_matrix *mm)
     }
 
     /* Read color map */
-    snew(map, mm->nmap);
+    mm.map.resize(nmap);
     m = 0;
-    while ((m < mm->nmap) && (nullptr != fgetline(&line_buf, llmax, &llalloc, in)))
+    while ((m < nmap) && (nullptr != fgetline(&line_buf, llmax, &llalloc, in)))
     {
         line = std::strchr(line_buf, '\"');
         if  (line)
         {
             line++;
             /* Read xpm color map entry */
-            map[m].code.c1 = line[0];
+            mm.map[m].code.c1 = line[0];
             if (nch == 1)
             {
-                map[m].code.c2 = 0;
+                mm.map[m].code.c2 = 0;
             }
             else
             {
-                map[m].code.c2 = line[1];
+                mm.map[m].code.c2 = line[1];
             }
             line += nch;
             str   = std::strchr(line, '#');
@@ -384,16 +374,16 @@ static void read_xpm_entry(FILE *in, t_matrix *mm)
                 if (col_len == 6)
                 {
                     sscanf(line, "%*s #%2x%2x%2x", &r, &g, &b);
-                    map[m].rgb.r = r/255.0;
-                    map[m].rgb.g = g/255.0;
-                    map[m].rgb.b = b/255.0;
+                    mm.map[m].rgb.r = r/255.0;
+                    mm.map[m].rgb.g = g/255.0;
+                    mm.map[m].rgb.b = b/255.0;
                 }
                 else if (col_len == 12)
                 {
                     sscanf(line, "%*s #%4x%4x%4x", &r, &g, &b);
-                    map[m].rgb.r = r/65535.0;
-                    map[m].rgb.g = g/65535.0;
-                    map[m].rgb.b = b/65535.0;
+                    mm.map[m].rgb.r = r/65535.0;
+                    mm.map[m].rgb.g = g/65535.0;
+                    mm.map[m].rgb.b = b/65535.0;
                 }
                 else
                 {
@@ -412,26 +402,23 @@ static void read_xpm_entry(FILE *in, t_matrix *mm)
                     gmx_file("Unsupported or invalid colormap in X PixMap");
                 }
                 fprintf(stderr, "Using white for color \"%s", str);
-                map[m].rgb.r = 1;
-                map[m].rgb.g = 1;
-                map[m].rgb.b = 1;
+                mm.map[m].rgb.r = 1;
+                mm.map[m].rgb.g = 1;
+                mm.map[m].rgb.b = 1;
             }
             line = std::strchr(line, '\"');
             line++;
             line2string(&line);
-            map[m].desc = gmx_strdup(line);
+            mm.map[m].desc = gmx_strdup(line);
             m++;
         }
     }
-    if  (m != mm->nmap)
+    if  (m != nmap)
     {
-        gmx_fatal(FARGS, "Number of read colors map entries (%d) does not match the number in the header (%d)", m, mm->nmap);
+        gmx_fatal(FARGS, "Number of read colors map entries (%d) does not match the number in the header (%d)", m, nmap);
     }
-    mm->map = map;
 
     /* Read axes, if there are any */
-    n_axis_x     = 0;
-    n_axis_y     = 0;
     bSetLine     = FALSE;
     do
     {
@@ -444,22 +431,19 @@ static void read_xpm_entry(FILE *in, t_matrix *mm)
         {
             line = std::strstr(line, "x-axis");
             skipstr(line);
-            if (mm->axis_x == nullptr)
-            {
-                snew(mm->axis_x, mm->nx + 1);
-            }
+            mm.axis_x.resize(0);
+            mm.axis_x.reserve(mm.nx + 1);
             while (sscanf(line, "%lf", &u) == 1)
             {
-                if (n_axis_x > mm->nx)
+                if (ssize(mm.axis_x) > mm.nx)
                 {
-                    gmx_fatal(FARGS, "Too many x-axis labels in xpm (max %d)", mm->nx);
+                    gmx_fatal(FARGS, "Too many x-axis labels in xpm (max %d)", mm.nx);
                 }
-                else if (n_axis_x == mm->nx)
+                else if (ssize(mm.axis_x) == mm.nx)
                 {
-                    mm->flags |= MAT_SPATIAL_X;
+                    mm.flags |= MAT_SPATIAL_X;
                 }
-                mm->axis_x[n_axis_x] = u;
-                n_axis_x++;
+                mm.axis_x.push_back(u);
                 skipstr(line);
             }
         }
@@ -467,22 +451,19 @@ static void read_xpm_entry(FILE *in, t_matrix *mm)
         {
             line = std::strstr(line, "y-axis");
             skipstr(line);
-            if (mm->axis_y == nullptr)
-            {
-                snew(mm->axis_y, mm->ny + 1);
-            }
+            mm.axis_y.resize(0);
+            mm.axis_y.reserve(mm.ny + 1);
             while (sscanf(line, "%lf", &u) == 1)
             {
-                if (n_axis_y > mm->ny)
+                if (ssize(mm.axis_y) > mm.ny)
                 {
-                    gmx_file("Too many y-axis labels in xpm");
+                    gmx_fatal(FARGS, "Too many y-axis labels in xpm (max %d)", mm.ny);
                 }
-                else if (n_axis_y == mm->ny)
+                else if (ssize(mm.axis_y) == mm.ny)
                 {
-                    mm->flags |= MAT_SPATIAL_Y;
+                    mm.flags |= MAT_SPATIAL_Y;
                 }
-                mm->axis_y[n_axis_y] = u;
-                n_axis_y++;
+                mm.axis_y.push_back(u);
                 skipstr(line);
             }
         }
@@ -490,12 +471,8 @@ static void read_xpm_entry(FILE *in, t_matrix *mm)
     while ((line[0] != '\"') && (nullptr != fgetline(&line_buf, llmax, &llalloc, in)));
 
     /* Read matrix */
-    snew(mm->matrix, mm->nx);
-    for (i = 0; i < mm->nx; i++)
-    {
-        snew(mm->matrix[i], mm->ny);
-    }
-    m        = mm->ny-1;
+    mm.matrix.resize(mm.nx, mm.ny);
+    int rowIndex = mm.ny - 1;
     bSetLine = FALSE;
     do
     {
@@ -504,9 +481,9 @@ static void read_xpm_entry(FILE *in, t_matrix *mm)
             line = line_buf;
         }
         bSetLine = TRUE;
-        if (m%(1+mm->ny/100) == 0)
+        if (rowIndex%(1+mm.ny/100) == 0)
         {
-            fprintf(stderr, "%3d%%\b\b\b\b", (100*(mm->ny-m))/mm->ny);
+            fprintf(stderr, "%3d%%\b\b\b\b", (100*(mm.ny-rowIndex))/mm.ny);
         }
         while ((line[0] != '\"') && (line[0] != '\0'))
         {
@@ -514,12 +491,12 @@ static void read_xpm_entry(FILE *in, t_matrix *mm)
         }
         if (line[0] != '\"')
         {
-            gmx_fatal(FARGS, "Not enough caracters in row %d of the matrix\n", m+1);
+            gmx_fatal(FARGS, "Not enough characters in row %d of the matrix\n", rowIndex+1);
         }
         else
         {
             line++;
-            for (i = 0; i < mm->nx; i++)
+            for (i = 0; i < mm.nx; i++)
             {
                 c.c1 = line[nch*i];
                 if (nch == 1)
@@ -530,70 +507,62 @@ static void read_xpm_entry(FILE *in, t_matrix *mm)
                 {
                     c.c2 = line[nch*i+1];
                 }
-                mm->matrix[i][m] = searchcmap(mm->nmap, mm->map, c);
+                mm.matrix(i, rowIndex) = searchcmap(mm.map, c);
             }
-            m--;
+            rowIndex--;
         }
     }
-    while ((m >= 0) && (nullptr != fgetline(&line_buf, llmax, &llalloc, in)));
-    if (m >= 0)
+    while ((rowIndex >= 0) && (nullptr != fgetline(&line_buf, llmax, &llalloc, in)));
+    if (rowIndex >= 0)
     {
         gmx_incons("Not enough rows in the matrix");
     }
 
     sfree(line_buf);
+    return mm;
 }
 
-int read_xpm_matrix(const char *fnm, t_matrix **mat)
+std::vector<t_matrix> read_xpm_matrix(const char *fnm)
 {
     FILE *in;
-    char *line = nullptr;
-    int   nmat;
+    char *line    = nullptr;
     int   llalloc = 0;
 
     in = gmx_fio_fopen(fnm, "r");
 
-    nmat = 0;
+    std::vector<t_matrix> mat;
     while (nullptr != fgetline(&line, STRLEN, &llalloc, in))
     {
         if (std::strstr(line, "/* XPM */"))
         {
-            srenew(*mat, nmat+1);
-            read_xpm_entry(in, &(*mat)[nmat]);
-            nmat++;
+            mat.emplace_back(read_xpm_entry(in));
         }
     }
     gmx_fio_fclose(in);
 
-    if (nmat == 0)
+    if (mat.empty())
     {
         gmx_file("Invalid XPixMap");
     }
 
     sfree(line);
 
-    return nmat;
+    return mat;
 }
 
 real **matrix2real(t_matrix *in, real **out)
 {
-    t_mapping *map;
-    double     tmp;
-    real      *rmap;
-    int        i, j, nmap;
+    double            tmp;
 
-    nmap = in->nmap;
-    map  = in->map;
-    snew(rmap, nmap);
+    std::vector<real> rmap(in->map.size());
 
-    for (i = 0; i < nmap; i++)
+    for (gmx::index i = 0; i != ssize(in->map); ++i)
     {
-        if ((map[i].desc == nullptr) || (sscanf(map[i].desc, "%lf", &tmp) != 1))
+        if ((in->map[i].desc == nullptr) || (sscanf(in->map[i].desc, "%lf", &tmp) != 1))
         {
             fprintf(stderr, "Could not convert matrix to reals,\n"
-                    "color map entry %d has a non-real description: \"%s\"\n",
-                    i, map[i].desc);
-            sfree(rmap);
+                    "color map entry %zd has a non-real description: \"%s\"\n",
+                    i, in->map[i].desc);
             return nullptr;
         }
         rmap[i] = tmp;
@@ -602,38 +571,36 @@ real **matrix2real(t_matrix *in, real **out)
     if (out == nullptr)
     {
         snew(out, in->nx);
-        for (i = 0; i < in->nx; i++)
+        for (int i = 0; i < in->nx; i++)
         {
             snew(out[i], in->ny);
         }
     }
-    for (i = 0; i < in->nx; i++)
+    for (int i = 0; i < in->nx; i++)
     {
-        for (j = 0; j < in->ny; j++)
+        for (int j = 0; j < in->ny; j++)
         {
-            out[i][j] = rmap[in->matrix[i][j]];
+            out[i][j] = rmap[in->matrix(i, j)];
         }
     }
 
-    sfree(rmap);
-
-    fprintf(stderr, "Converted a %dx%d matrix with %d levels to reals\n",
-            in->nx, in->ny, nmap);
+    fprintf(stderr, "Converted a %dx%d matrix with %zu levels to reals\n",
+            in->nx, in->ny, in->map.size());
 
     return out;
 }
 
 static void write_xpm_header(FILE *out,
-                             const char *title, const char *legend,
-                             const char *label_x, const char *label_y,
+                             const std::string &title, const std::string &legend,
+                             const std::string &label_x, const std::string &label_y,
                              gmx_bool bDiscrete)
 {
     fprintf(out,  "/* XPM */\n");
     fprintf(out,  "/* This file can be converted to EPS by the GROMACS program xpm2ps */\n");
-    fprintf(out,  "/* title:   \"%s\" */\n", title);
-    fprintf(out,  "/* legend:  \"%s\" */\n", legend);
-    fprintf(out,  "/* x-label: \"%s\" */\n", label_x);
-    fprintf(out,  "/* y-label: \"%s\" */\n", label_y);
+    fprintf(out,  "/* title:   \"%s\" */\n", title.c_str());
+    fprintf(out,  "/* legend:  \"%s\" */\n", legend.c_str());
+    fprintf(out,  "/* x-label: \"%s\" */\n", label_x.c_str());
+    fprintf(out,  "/* y-label: \"%s\" */\n", label_y.c_str());
     if (bDiscrete)
     {
         fprintf(out, "/* type:    \"Discrete\" */\n");
@@ -841,27 +808,26 @@ static void write_xpm_map(FILE *out, int n_x, int n_y, int *nlevels,
     }
 }
 
-static void write_xpm_axis(FILE *out, const char *axis, gmx_bool bSpatial,
-                           int n, real *label)
+static void writeXpmAxis(FILE *out, const char *axis,
+                         ArrayRef<const real> label)
 {
-    int i;
-
-    if (label)
+    if (label.empty())
+    {
+        return;
+    }
+    for (gmx::index i = 0; i != ssize(label); ++i)
     {
-        for (i = 0; i < (bSpatial ? n+1 : n); i++)
+        if (i % 80 == 0)
         {
-            if (i % 80 == 0)
+            if (i)
             {
-                if (i)
-                {
-                    fprintf(out, "*/\n");
-                }
-                fprintf(out, "/* %s-axis:  ", axis);
+                fprintf(out, "*/\n");
             }
-            fprintf(out, "%g ", label[i]);
+            fprintf(out, "/* %s-axis:  ", axis);
         }
-        fprintf(out, "*/\n");
+        fprintf(out, "%g ", label[i]);
     }
+    fprintf(out, "*/\n");
 }
 
 static void write_xpm_data(FILE *out, int n_x, int n_y, real **mat,
@@ -880,7 +846,7 @@ static void write_xpm_data(FILE *out, int n_x, int n_y, real **mat,
         fprintf(out, "\"");
         for (i = 0; (i < n_x); i++)
         {
-            c = gmx::roundToInt((mat[i][j]-lo)*invlevel);
+            c = roundToInt((mat[i][j]-lo)*invlevel);
             if (c < 0)
             {
                 c = 0;
@@ -930,11 +896,11 @@ static void write_xpm_data3(FILE *out, int n_x, int n_y, real **mat,
         {
             if (mat[i][j] >= mid)
             {
-                c = nmid+gmx::roundToInt((mat[i][j]-mid)*invlev_hi);
+                c = nmid+roundToInt((mat[i][j]-mid)*invlev_hi);
             }
             else if (mat[i][j] >= lo)
             {
-                c = gmx::roundToInt((mat[i][j]-lo)*invlev_lo);
+                c = roundToInt((mat[i][j]-lo)*invlev_lo);
             }
             else
             {
@@ -990,7 +956,7 @@ static void write_xpm_data_split(FILE *out, int n_x, int n_y, real **mat,
         {
             if (i < j)
             {
-                c = nlevel_bot+gmx::roundToInt((mat[i][j]-lo_top)*invlev_top);
+                c = nlevel_bot+roundToInt((mat[i][j]-lo_top)*invlev_top);
                 if ((c < nlevel_bot) || (c >= nlevel_bot+nlevel_top))
                 {
                     gmx_fatal(FARGS, "Range checking i = %d, j = %d, c = %d, bot = %d, top = %d matrix[i,j] = %f", i, j, c, nlevel_bot, nlevel_top, mat[i][j]);
@@ -998,7 +964,7 @@ static void write_xpm_data_split(FILE *out, int n_x, int n_y, real **mat,
             }
             else if (i > j)
             {
-                c = gmx::roundToInt((mat[i][j]-lo_bot)*invlev_bot);
+                c = roundToInt((mat[i][j]-lo_bot)*invlev_bot);
                 if ((c < 0) || (c >= nlevel_bot+nlevel_bot))
                 {
                     gmx_fatal(FARGS, "Range checking i = %d, j = %d, c = %d, bot = %d, top = %d matrix[i,j] = %f", i, j, c, nlevel_bot, nlevel_top, mat[i][j]);
@@ -1024,9 +990,6 @@ static void write_xpm_data_split(FILE *out, int n_x, int n_y, real **mat,
 
 void write_xpm_m(FILE *out, t_matrix m)
 {
-    /* Writes a t_matrix struct to .xpm file */
-
-    int           i, j;
     gmx_bool      bOneChar;
     t_xpmelmt     c;
 
@@ -1034,19 +997,19 @@ void write_xpm_m(FILE *out, t_matrix m)
     write_xpm_header(out, m.title, m.legend, m.label_x, m.label_y,
                      m.bDiscrete);
     fprintf(out, "static char *gromacs_xpm[] = {\n");
-    fprintf(out, "\"%d %d   %d %d\",\n", m.nx, m.ny, m.nmap, bOneChar ? 1 : 2);
-    for (i = 0; (i < m.nmap); i++)
+    fprintf(out, "\"%d %d   %zu %d\",\n", m.nx, m.ny, m.map.size(), bOneChar ? 1 : 2);
+    for (const auto &map : m.map)
     {
         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%s\" */,\n",
-                m.map[i].code.c1,
-                bOneChar ? ' ' : m.map[i].code.c2,
-                static_cast<unsigned int>(round(m.map[i].rgb.r*255)),
-                static_cast<unsigned int>(round(m.map[i].rgb.g*255)),
-                static_cast<unsigned int>(round(m.map[i].rgb.b*255)), m.map[i].desc);
+                map.code.c1,
+                bOneChar ? ' ' : map.code.c2,
+                static_cast<unsigned int>(round(map.rgb.r*255)),
+                static_cast<unsigned int>(round(map.rgb.g*255)),
+                static_cast<unsigned int>(round(map.rgb.b*255)), map.desc);
     }
-    write_xpm_axis(out, "x", (m.flags & MAT_SPATIAL_X) != 0u, m.nx, m.axis_x);
-    write_xpm_axis(out, "y", (m.flags & MAT_SPATIAL_Y) != 0u, m.ny, m.axis_y);
-    for (j = m.ny-1; (j >= 0); j--)
+    writeXpmAxis(out, "x", m.axis_x);
+    writeXpmAxis(out, "y", m.axis_y);
+    for (int j = m.ny-1; (j >= 0); j--)
     {
         if (j%(1+m.ny/100) == 0)
         {
@@ -1055,16 +1018,16 @@ void write_xpm_m(FILE *out, t_matrix m)
         fprintf(out, "\"");
         if (bOneChar)
         {
-            for (i = 0; (i < m.nx); i++)
+            for (int i = 0; (i < m.nx); i++)
             {
-                fprintf(out, "%c", m.map[m.matrix[i][j]].code.c1);
+                fprintf(out, "%c", m.map[m.matrix(i, j)].code.c1);
             }
         }
         else
         {
-            for (i = 0; (i < m.nx); i++)
+            for (int i = 0; (i < m.nx); i++)
             {
-                c = m.map[m.matrix[i][j]].code;
+                c = m.map[m.matrix(i, j)].code;
                 fprintf(out, "%c%c", c.c1, c.c2);
             }
         }
@@ -1095,10 +1058,10 @@ void write_xpm3(FILE *out, unsigned int flags,
         gmx_fatal(FARGS, "hi (%g) <= lo (%g)", hi, lo);
     }
 
-    write_xpm_header(out, title.c_str(), legend.c_str(), label_x.c_str(), label_y.c_str(), FALSE);
+    write_xpm_header(out, title, legend, label_x, label_y, FALSE);
     write_xpm_map3(out, n_x, n_y, nlevels, lo, mid, hi, rlo, rmid, rhi);
-    write_xpm_axis(out, "x", (flags & MAT_SPATIAL_X) != 0u, n_x, axis_x);
-    write_xpm_axis(out, "y", (flags & MAT_SPATIAL_Y) != 0u, n_y, axis_y);
+    writeXpmAxis(out, "x", ArrayRef<real>(axis_x, axis_x + n_x + ((flags & MAT_SPATIAL_X) != 0u ? 1 : 0)));
+    writeXpmAxis(out, "y", ArrayRef<real>(axis_y, axis_y + n_y + ((flags & MAT_SPATIAL_Y) != 0u ? 1 : 0)));
     write_xpm_data3(out, n_x, n_y, mat, lo, mid, hi, *nlevels);
 }
 
@@ -1130,11 +1093,11 @@ void write_xpm_split(FILE *out, unsigned int flags,
         gmx_impl("Can not plot more than 16 discrete colors");
     }
 
-    write_xpm_header(out, title.c_str(), legend.c_str(), label_x.c_str(), label_y.c_str(), FALSE);
+    write_xpm_header(out, title, legend, label_x, label_y, FALSE);
     write_xpm_map_split(out, n_x, n_y, nlevel_top, lo_top, hi_top, rlo_top, rhi_top,
                         bDiscreteColor, nlevel_bot, lo_bot, hi_bot, rlo_bot, rhi_bot);
-    write_xpm_axis(out, "x", (flags & MAT_SPATIAL_X) != 0u, n_x, axis_x);
-    write_xpm_axis(out, "y", (flags & MAT_SPATIAL_Y) != 0u, n_y, axis_y);
+    writeXpmAxis(out, "x", ArrayRef<real>(axis_x, axis_x + n_x + ((flags & MAT_SPATIAL_X) != 0u ? 1 : 0)));
+    writeXpmAxis(out, "y", ArrayRef<real>(axis_y, axis_y + n_y + ((flags & MAT_SPATIAL_Y) != 0u ? 1 : 0)));
     write_xpm_data_split(out, n_x, n_y, mat, lo_top, hi_top, *nlevel_top,
                          lo_bot, hi_bot, *nlevel_bot);
 }
@@ -1167,9 +1130,9 @@ void write_xpm(FILE *out, unsigned int flags,
         gmx_fatal(FARGS, "hi (%f) <= lo (%f)", hi, lo);
     }
 
-    write_xpm_header(out, title.c_str(), legend.c_str(), label_x.c_str(), label_y.c_str(), FALSE);
+    write_xpm_header(out, title, legend, label_x, label_y, FALSE);
     write_xpm_map(out, n_x, n_y, nlevels, lo, hi, rlo, rhi);
-    write_xpm_axis(out, "x", (flags & MAT_SPATIAL_X) != 0u, n_x, axis_x);
-    write_xpm_axis(out, "y", (flags & MAT_SPATIAL_Y) != 0u, n_y, axis_y);
+    writeXpmAxis(out, "x", ArrayRef<real>(axis_x, axis_x + n_x + ((flags & MAT_SPATIAL_X) != 0u ? 1 : 0)));
+    writeXpmAxis(out, "y", ArrayRef<real>(axis_y, axis_y + n_y + ((flags & MAT_SPATIAL_Y) != 0u ? 1 : 0)));
     write_xpm_data(out, n_x, n_y, mat, lo, hi, *nlevels);
 }
index a06058370328145a034abae8abc789f2532a8f64..5de0335b2f100fac53bea46bb6df776799cd090c 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2017,2019, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
 #include <stdio.h>
 
 #include <string>
+#include <vector>
 
 #include "gromacs/fileio/rgb.h"
+#include "gromacs/math/multidimarray.h"
 #include "gromacs/utility/basedefinitions.h"
 #include "gromacs/utility/real.h"
 
-typedef struct {
-    char c1; /* should all be non-zero (and printable and not '"') */
-    char c2; /*
-              * should all be zero (single char color names: smaller xpm's)
-              * or should all be non-zero (double char color names: more colors)
-              */
-} t_xpmelmt;
-
+/*! \brief Models an XPM element
+ *
+ * \libinternal */
+struct t_xpmelmt
+{
+    //! Should all be non-zero (and printable and not '"')
+    char c1 = 0;
+    /*! \brief Should all be zero (single char color names: smaller xpm's)
+     * or should all be non-zero (double char color names: more colors). */
+    char c2 = 0;
+};
+
+//! Type of a matrix element
 typedef short t_matelmt;
 
-typedef struct {
-    t_xpmelmt   code; /* see comment for t_xpmelmt */
-    const char *desc;
+/*! \brief Maps an XPM element to an RGB color and a string description.
+ *
+ * \libinternal */
+struct t_mapping
+{
+    //! XPM element code
+    t_xpmelmt   code;
+    //! Description
+    const char *desc = nullptr;
+    //! RGB color
     t_rgb       rgb;
-} t_mapping;
+};
 
 #define MAT_SPATIAL_X (1<<0)
 #define MAT_SPATIAL_Y (1<<1)
@@ -69,44 +83,50 @@ typedef struct {
  * when set, there are n+1 axis ticks at the edges of the elements.
  */
 
-typedef struct {
-    unsigned int flags; /* The possible flags are defined above */
-    int          nx, ny;
-    int          y0;
-    char         title[256];
-    char         legend[256];
-    char         label_x[256];
-    char         label_y[256];
-    gmx_bool     bDiscrete;
-    real        *axis_x;
-    real        *axis_y;
-    t_matelmt  **matrix;
-    int          nmap;
-    t_mapping   *map;
-} t_matrix;
-/* title      matrix title
- * legend     label for the continuous legend
- * label_x    label for the x-axis
- * label_y    label for the y-axis
- * nx, ny     size of the matrix
- * axis_x[]   the x-ticklabels
- * axis_y[]   the y-ticklables
- * *matrix[]  element x,y is matrix[x][y]
- * nmap       number of color levels for the output(?)
- */
+//! Convenience typedef
+template <typename T>
+using DynamicMatrix2D = gmx::MultiDimArray < std::vector<T>, gmx::extents < gmx::dynamic_extent, gmx::dynamic_extent>>;
 
-gmx_bool matelmt_cmp(t_xpmelmt e1, t_xpmelmt e2);
-
-t_matelmt searchcmap(int n, t_mapping map[], t_xpmelmt c);
-/* Seach in the map for code 'c' and return entry number.
- * return -1 if not found
- */
-
-int getcmap(FILE *in, const char *fn, t_mapping **map);
-/* Read the mapping table from in, return number of entries */
-
-int readcmap(const char *fn, t_mapping **map);
-/* Read the mapping table from fn, return number of entries */
+/*! \brief A matrix of integers, plus supporting values, such as used in XPM output
+ *
+ * \libinternal
+ *
+ * \todo nx and ny should probably be replaced by operations on the
+ * extent of matrix, but currently there is only limited ability to
+ * iterate over contents of matrix. */
+struct t_matrix
+{
+    //! Defines if x and y are spatial dimensions. See comments on MAT_*.
+    unsigned int               flags = 0;
+    //! Size of the matrix in x
+    int                        nx = 0;
+    //! Size of the matrix in y
+    int                        ny = 0;
+    //! Matrix title
+    std::string                title;
+    //! Label for the continuous legend
+    std::string                legend;
+    //! Label for the x-axis
+    std::string                label_x;
+    //! Label for the y-axis
+    std::string                label_y;
+    //! Whether the quantity described is discrete or continuous.
+    bool                       bDiscrete = false;
+    //! The x-ticklabels
+    std::vector<real>          axis_x;
+    //! The y-ticklabels
+    std::vector<real>          axis_y;
+    //! Element x,y is matrix(x, y)
+    DynamicMatrix2D<t_matelmt> matrix;
+    //! Color levels for the output(?)
+    std::vector<t_mapping>     map;
+};
+
+//! Seach in the \c map for code \c c and return its entry, or -1 if not found.
+t_matelmt searchcmap(gmx::ArrayRef<const t_mapping> map, t_xpmelmt c);
+
+//! Read the mapping table from fn, return number of entries
+std::vector<t_mapping> readcmap(const char *fn);
 
 void printcmap(FILE *out, int n, t_mapping map[]);
 /* print mapping table to out */
@@ -114,8 +134,8 @@ void printcmap(FILE *out, int n, t_mapping map[]);
 void writecmap(const char *fn, int n, t_mapping map[]);
 /* print mapping table to fn */
 
-int read_xpm_matrix(const char *fnm, t_matrix **mat);
-/* Reads a number of matrices from .xpm file fnm and returns this number */
+//! Reads and returns a number of matrices from .xpm file \c fnm.
+std::vector<t_matrix> read_xpm_matrix(const char *fnm);
 
 real **matrix2real(t_matrix *in, real **out);
 /* Converts an matrix in a t_matrix struct to a matrix of reals
index 37d375a9c13c35dcd7b119549e7eb1d9baeb26d5..0a273e0c31f0e470f92a77edff94070b6e7f10f5 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2011,2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2011,2014,2015,2016,2019, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
 #ifndef GMX_FILEIO_RGB_H
 #define GMX_FILEIO_RGB_H
 
-typedef struct t_rgb {
-    double r, g, b;
-} t_rgb;
+struct t_rgb
+{
+    double r = 0;
+    double g = 0;
+    double b = 0;
+};
 
 #endif
index 9fe99f540fc317ce724bda6e52d7b81fc79316d4..619ad7c214ca79239e395755cdf89daee01bf9b5 100644 (file)
@@ -1313,12 +1313,6 @@ static void convert_mat(t_matrix *mat, t_mat *rms)
 
     rms->n1 = mat->nx;
     matrix2real(mat, rms->mat);
-    /* free input xpm matrix data */
-    for (i = 0; i < mat->nx; i++)
-    {
-        sfree(mat->matrix[i]);
-    }
-    sfree(mat->matrix);
 
     for (i = 0; i < mat->nx; i++)
     {
@@ -1428,7 +1422,6 @@ int gmx_cluster(int argc, char *argv[])
     t_topology         top;
     int                ePBC;
     t_atoms            useatoms;
-    t_matrix          *readmat = nullptr;
     real              *eigenvectors;
 
     int                isize = 0, ifsize = 0, iosize = 0;
@@ -1720,10 +1713,11 @@ int gmx_cluster(int argc, char *argv[])
         }
     }
 
+    std::vector<t_matrix> readmat;
     if (bReadMat)
     {
         fprintf(stderr, "Reading rms distance matrix ");
-        read_xpm_matrix(opt2fn("-dm", NFILE, fnm), &readmat);
+        readmat = read_xpm_matrix(opt2fn("-dm", NFILE, fnm));
         fprintf(stderr, "\n");
         if (readmat[0].nx != readmat[0].ny)
         {
@@ -1738,7 +1732,7 @@ int gmx_cluster(int argc, char *argv[])
 
         nf = readmat[0].nx;
         sfree(time);
-        time        = readmat[0].axis_x;
+        time        = readmat[0].axis_x.data();
         time_invfac = output_env_get_time_invfactor(oenv);
         for (i = 0; i < nf; i++)
         {
@@ -1748,7 +1742,7 @@ int gmx_cluster(int argc, char *argv[])
         rms = init_mat(readmat[0].nx, method == m_diagonalize);
         convert_mat(&(readmat[0]), rms);
 
-        nlevels = readmat[0].nmap;
+        nlevels = gmx::ssize(readmat[0].map);
     }
     else   /* !bReadMat */
     {
@@ -1959,7 +1953,7 @@ int gmx_cluster(int argc, char *argv[])
     if (bReadMat)
     {
         write_xpm(fp, 0, readmat[0].title, readmat[0].legend, readmat[0].label_x,
-                  readmat[0].label_y, nf, nf, readmat[0].axis_x, readmat[0].axis_y,
+                  readmat[0].label_y, nf, nf, readmat[0].axis_x.data(), readmat[0].axis_y.data(),
                   rms->mat, 0.0, rms->maxrms, rlo_top, rhi_top, &nlevels);
     }
     else
index a466b222c2989d9a63c0a4cb7227eb53c1d971dc..e0ab634dc4d7912d5ebeaa2bf29da3e490223cc7 100644 (file)
@@ -42,6 +42,7 @@
 #include <cstring>
 
 #include <algorithm>
+#include <numeric>
 
 #include "gromacs/commandline/pargs.h"
 #include "gromacs/commandline/viewit.h"
@@ -104,13 +105,11 @@ static int strip_dssp(FILE *tapeout, int nres,
 {
     static gmx_bool bFirst = TRUE;
     static char    *ssbuf;
-    static int      xsize, frame;
     char            buf[STRLEN+1];
     char            SSTP;
-    int             i, nr, iacc, nresidues;
+    int             nr, iacc, nresidues;
     int             naccf, naccb; /* Count hydrophobic and hydrophilic residues */
     real            iaccf, iaccb;
-    t_xpmelmt       c;
 
 
     /* Skip header */
@@ -178,39 +177,27 @@ static int strip_dssp(FILE *tapeout, int nres,
             fprintf(stderr, "%d residues were classified as hydrophobic and %d as hydrophilic.\n", naccb, naccf);
         }
 
-        sprintf(mat->title, "Secondary structure");
-        mat->legend[0] = 0;
-        sprintf(mat->label_x, "%s", output_env_get_time_label(oenv).c_str());
-        sprintf(mat->label_y, "Residue");
-        mat->bDiscrete = TRUE;
+        mat->title     = "Secondary structure";
+        mat->legend    = "";
+        mat->label_x   = output_env_get_time_label(oenv);
+        mat->label_y   = "Residue";
+        mat->bDiscrete = true;
         mat->ny        = nr;
-        snew(mat->axis_y, nr);
-        for (i = 0; i < nr; i++)
-        {
-            mat->axis_y[i] = i+1;
-        }
-        mat->axis_x = nullptr;
-        mat->matrix = nullptr;
-        xsize       = 0;
-        frame       = 0;
-        bFirst      = FALSE;
-    }
-    if (frame >= xsize)
-    {
-        xsize += 10;
-        srenew(mat->axis_x, xsize);
-        srenew(mat->matrix, xsize);
+        mat->axis_y.resize(nr);
+        std::iota(mat->axis_y.begin(), mat->axis_y.end(), 1);
+        mat->axis_x.resize(0);
+        mat->matrix.resize(0, 0);
+        bFirst      = false;
     }
-    mat->axis_x[frame] = t;
-    snew(mat->matrix[frame], nr);
-    c.c2 = 0;
-    for (i = 0; i < nr; i++)
+    mat->axis_x.push_back(t);
+    mat->matrix.resize(mat->matrix.extent(0), nr);
+    mat->nx = mat->matrix.extent(0);
+    auto columnIndex = mat->nx - 1;
+    for (int i = 0; i < nr; i++)
     {
-        c.c1                  = ssbuf[i];
-        mat->matrix[frame][i] = std::max(static_cast<t_matelmt>(0), searchcmap(mat->nmap, mat->map, c));
+        t_xpmelmt c = {ssbuf[i], 0};
+        mat->matrix(columnIndex, i) = std::max(static_cast<t_matelmt>(0), searchcmap(mat->map, c));
     }
-    frame++;
-    mat->nx = frame;
 
     if (fTArea)
     {
@@ -327,44 +314,35 @@ static void norm_acc(t_atoms *atoms, int nres,
 
 static void prune_ss_legend(t_matrix *mat)
 {
-    gmx_bool  *present;
-    int       *newnum;
-    int        i, r, f, newnmap;
-    t_mapping *newmap;
+    std::vector<bool>      isPresent(mat->map.size());
+    std::vector<int>       newnum(mat->map.size());
+    std::vector<t_mapping> newmap;
 
-    snew(present, mat->nmap);
-    snew(newnum, mat->nmap);
-
-    for (f = 0; f < mat->nx; f++)
+    for (int f = 0; f < mat->nx; f++)
     {
-        for (r = 0; r < mat->ny; r++)
+        for (int r = 0; r < mat->ny; r++)
         {
-            present[mat->matrix[f][r]] = TRUE;
+            isPresent[mat->matrix(f, r)] = true;
         }
     }
 
-    newnmap = 0;
-    newmap  = nullptr;
-    for (i = 0; i < mat->nmap; i++)
+    for (size_t i = 0; i < mat->map.size(); i++)
     {
         newnum[i] = -1;
-        if (present[i])
+        if (isPresent[i])
         {
-            newnum[i] = newnmap;
-            newnmap++;
-            srenew(newmap, newnmap);
-            newmap[newnmap-1] = mat->map[i];
+            newnum[i] = newmap.size();
+            newmap.emplace_back(mat->map[i]);
         }
     }
-    if (newnmap != mat->nmap)
+    if (newmap.size() != mat->map.size())
     {
-        mat->nmap = newnmap;
-        mat->map  = newmap;
-        for (f = 0; f < mat->nx; f++)
+        std::swap(mat->map, newmap);
+        for (int f = 0; f < mat->nx; f++)
         {
-            for (r = 0; r < mat->ny; r++)
+            for (int r = 0; r < mat->ny; r++)
             {
-                mat->matrix[f][r] = newnum[mat->matrix[f][r]];
+                mat->matrix(f, r) = newnum[mat->matrix(f, r)];
             }
         }
     }
@@ -392,7 +370,7 @@ static void write_sas_mat(const char *fn, real **accr, int nframe, int nres, t_m
         nlev = static_cast<int>(hi-lo+1);
         write_xpm(fp, 0, "Solvent Accessible Surface", "Surface (A^2)",
                   "Time", "Residue Index", nframe, nres,
-                  mat->axis_x, mat->axis_y, accr, lo, hi, rlo, rhi, &nlev);
+                  mat->axis_x.data(), mat->axis_y.data(), accr, lo, hi, rlo, rhi, &nlev);
         gmx_ffclose(fp);
     }
 }
@@ -400,20 +378,20 @@ static void write_sas_mat(const char *fn, real **accr, int nframe, int nres, t_m
 static void analyse_ss(const char *outfile, t_matrix *mat, const char *ss_string,
                        const gmx_output_env_t *oenv)
 {
-    FILE        *fp;
-    t_mapping   *map;
-    int          f, r, *count, *total, ss_count, total_count;
-    size_t       s;
-    const char** leg;
+    FILE                    *fp;
+    int                      ss_count, total_count;
 
-    map = mat->map;
-    snew(count, mat->nmap);
-    snew(total, mat->nmap);
-    snew(leg, mat->nmap+1);
-    leg[0] = "Structure";
-    for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
+    gmx::ArrayRef<t_mapping> map = mat->map;
+    std::vector<int>         count(map.size());
+    std::vector<int>         total(map.size(), 0);
+    // This copying would not be necessary if xvgr_legend could take a
+    // view of string views
+    std::vector<std::string> leg;
+    leg.reserve(map.size()+1);
+    leg.emplace_back("Structure");
+    for (const auto &m : map)
     {
-        leg[s+1] = gmx_strdup(map[s].desc);
+        leg.emplace_back(m.desc);
     }
 
     fp = xvgropen(outfile, "Secondary Structure",
@@ -422,41 +400,37 @@ static void analyse_ss(const char *outfile, t_matrix *mat, const char *ss_string
     {
         fprintf(fp, "@ subtitle \"Structure = ");
     }
-    for (s = 0; s < std::strlen(ss_string); s++)
+    for (size_t s = 0; s < std::strlen(ss_string); s++)
     {
         if (s > 0)
         {
             fprintf(fp, " + ");
         }
-        for (f = 0; f < mat->nmap; f++)
+        for (const auto &m : map)
         {
-            if (ss_string[s] == map[f].code.c1)
+            if (ss_string[s] == m.code.c1)
             {
-                fprintf(fp, "%s", map[f].desc);
+                fprintf(fp, "%s", m.desc);
             }
         }
     }
     fprintf(fp, "\"\n");
-    xvgr_legend(fp, mat->nmap+1, leg, oenv);
+    xvgrLegend(fp, leg, oenv);
 
     total_count = 0;
-    for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
-    {
-        total[s] = 0;
-    }
-    for (f = 0; f < mat->nx; f++)
+    for (int f = 0; f < mat->nx; f++)
     {
         ss_count = 0;
-        for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
+        for (auto &c : count)
         {
-            count[s] = 0;
+            c = 0;
         }
-        for (r = 0; r < mat->ny; r++)
+        for (int r = 0; r < mat->ny; r++)
         {
-            count[mat->matrix[f][r]]++;
-            total[mat->matrix[f][r]]++;
+            count[mat->matrix(f, r)]++;
+            total[mat->matrix(f, r)]++;
         }
-        for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
+        for (gmx::index s = 0; s != gmx::ssize(map); ++s)
         {
             if (std::strchr(ss_string, map[s].code.c1))
             {
@@ -465,31 +439,29 @@ static void analyse_ss(const char *outfile, t_matrix *mat, const char *ss_string
             }
         }
         fprintf(fp, "%8g %5d", mat->axis_x[f], ss_count);
-        for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
+        for (const auto &c : count)
         {
-            fprintf(fp, " %5d", count[s]);
+            fprintf(fp, " %5d", c);
         }
         fprintf(fp, "\n");
     }
     /* now print column totals */
     fprintf(fp, "%-8s %5d", "# Totals", total_count);
-    for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
+    for (const auto &t : total)
     {
-        fprintf(fp, " %5d", total[s]);
+        fprintf(fp, " %5d", t);
     }
     fprintf(fp, "\n");
 
     /* now print probabilities */
     fprintf(fp, "%-8s %5.2f", "# SS pr.", total_count / static_cast<real>(mat->nx * mat->ny));
-    for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
+    for (const auto &t : total)
     {
-        fprintf(fp, " %5.2f", total[s] / static_cast<real>(mat->nx * mat->ny));
+        fprintf(fp, " %5.2f", t / static_cast<real>(mat->nx * mat->ny));
     }
     fprintf(fp, "\n");
 
     xvgrclose(fp);
-    sfree(leg);
-    sfree(count);
 }
 
 int gmx_do_dssp(int argc, char *argv[])
@@ -554,10 +526,10 @@ int gmx_do_dssp(int argc, char *argv[])
     int                nres, nr0, naccr, nres_plus_separators;
     gmx_bool          *bPhbres, bDoAccSurf;
     real               t;
-    int                i, j, natoms, nframe = 0;
+    int                natoms, nframe = 0;
     matrix             box = {{0}};
     int                gnx;
-    char              *grpnm, *ss_str;
+    char              *grpnm;
     int               *index;
     rvec              *xp, *x;
     int               *average_area;
@@ -602,7 +574,7 @@ int gmx_do_dssp(int argc, char *argv[])
     get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &gnx, &index, &grpnm);
     nres = 0;
     nr0  = -1;
-    for (i = 0; (i < gnx); i++)
+    for (int i = 0; (i < gnx); i++)
     {
         if (atoms->atom[index[i]].resind != nr0)
         {
@@ -687,8 +659,7 @@ int gmx_do_dssp(int argc, char *argv[])
         fTArea = nullptr;
     }
 
-    mat.map  = nullptr;
-    mat.nmap = readcmap(opt2fn("-map", NFILE, fnm), &(mat.map));
+    mat.map = readcmap(opt2fn("-map", NFILE, fnm));
 
     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
     if (natoms > atoms->nr)
@@ -714,7 +685,7 @@ int gmx_do_dssp(int argc, char *argv[])
         {
             naccr += 10;
             srenew(accr, naccr);
-            for (i = naccr-10; i < naccr; i++)
+            for (int i = naccr-10; i < naccr; i++)
             {
                 snew(accr[i], 2*atoms->nres-1);
             }
@@ -773,19 +744,17 @@ int gmx_do_dssp(int argc, char *argv[])
     if (opt2bSet("-ssdump", NFILE, fnm))
     {
         ss = opt2FILE("-ssdump", NFILE, fnm, "w");
-        snew(ss_str, nres+1);
         fprintf(ss, "%d\n", nres);
-        for (j = 0; j < mat.nx; j++)
+        for (gmx::index j = 0; j != mat.matrix.extent(0); ++j)
         {
-            for (i = 0; (i < mat.ny); i++)
+            auto row = mat.matrix.asView()[j];
+            for (gmx::index i = 0; i != row.extent(0); ++i)
             {
-                ss_str[i] = mat.map[mat.matrix[j][i]].code.c1;
+                fputc(mat.map[row[i]].code.c1, ss);
             }
-            ss_str[i] = '\0';
-            fprintf(ss, "%s\n", ss_str);
+            fputc('\n', ss);
         }
         gmx_ffclose(ss);
-        sfree(ss_str);
     }
     analyse_ss(fnSCount, &mat, ss_string, oenv);
 
@@ -793,7 +762,7 @@ int gmx_do_dssp(int argc, char *argv[])
     {
         write_sas_mat(fnArea, accr, nframe, nres_plus_separators, &mat);
 
-        for (i = 0; i < atoms->nres; i++)
+        for (int i = 0; i < atoms->nres; i++)
         {
             av_area[i] = (average_area[i] / static_cast<real>(nframe));
         }
@@ -804,7 +773,7 @@ int gmx_do_dssp(int argc, char *argv[])
         {
             acc = xvgropen(fnAArea, "Average Accessible Area",
                            "Residue", "A\\S2", oenv);
-            for (i = 0; (i < nres); i++)
+            for (int i = 0; (i < nres); i++)
             {
                 fprintf(acc, "%5d  %10g %10g\n", i+1, av_area[i], norm_av_area[i]);
             }
index 88803c7cf937795752815451a1df04cc24cbcc48..f4dc3d3701e4c978b067362b1ec5ee5389ca5f0e 100644 (file)
@@ -43,6 +43,7 @@
 #include <cstring>
 
 #include <algorithm>
+#include <numeric>
 
 #include "gromacs/commandline/pargs.h"
 #include "gromacs/commandline/viewit.h"
@@ -3326,17 +3327,17 @@ int gmx_hbond(int argc, char *argv[])
         {
             t_matrix mat;
             int      id, ia, hh, x, y;
-            mat.flags = mat.y0 = 0;
+            mat.flags = 0;
 
             if ((nframes > 0) && (hb->nrhb > 0))
             {
                 mat.nx = nframes;
                 mat.ny = hb->nrhb;
 
-                snew(mat.matrix, mat.nx);
-                for (x = 0; (x < mat.nx); x++)
+                mat.matrix.resize(mat.nx, mat.ny);
+                for (auto &value : mat.matrix.toArrayRef())
                 {
-                    snew(mat.matrix[x], mat.ny);
+                    value = 0;
                 }
                 y = 0;
                 for (id = 0; (id < hb->d.nrd); id++)
@@ -3353,7 +3354,7 @@ int gmx_hbond(int argc, char *argv[])
                                     {
                                         int nn0 = hb->hbmap[id][ia]->n0;
                                         range_check(y, 0, mat.ny);
-                                        mat.matrix[x+nn0][y] = static_cast<t_matelmt>(is_hb(hb->hbmap[id][ia]->h[hh], x));
+                                        mat.matrix(x+nn0, y) = static_cast<t_matelmt>(is_hb(hb->hbmap[id][ia]->h[hh], x));
                                     }
                                     y++;
                                 }
@@ -3361,36 +3362,25 @@ int gmx_hbond(int argc, char *argv[])
                         }
                     }
                 }
-                mat.axis_x = hb->time;
-                snew(mat.axis_y, mat.ny);
-                for (j = 0; j < mat.ny; j++)
+                std::copy(hb->time, hb->time + mat.nx, mat.axis_x.begin());
+                mat.axis_y.resize(mat.ny);
+                std::iota(mat.axis_y.begin(), mat.axis_y.end(), 0);
+                mat.title = (bContact ? "Contact Existence Map" :
+                             "Hydrogen Bond Existence Map");
+                mat.legend    = bContact ? "Contacts" : "Hydrogen Bonds";
+                mat.label_x   = output_env_get_xvgr_tlabel(oenv);
+                mat.label_y   = bContact ? "Contact Index" : "Hydrogen Bond Index";
+                mat.bDiscrete = true;
+                mat.map.resize(2);
+                for (auto &m : mat.map)
                 {
-                    mat.axis_y[j] = j;
-                }
-                sprintf(mat.title, bContact ? "Contact Existence Map" :
-                        "Hydrogen Bond Existence Map");
-                sprintf(mat.legend, bContact ? "Contacts" : "Hydrogen Bonds");
-                sprintf(mat.label_x, "%s", output_env_get_xvgr_tlabel(oenv).c_str());
-                sprintf(mat.label_y, bContact ? "Contact Index" : "Hydrogen Bond Index");
-                mat.bDiscrete = TRUE;
-                mat.nmap      = 2;
-                snew(mat.map, mat.nmap);
-                for (i = 0; i < mat.nmap; i++)
-                {
-                    mat.map[i].code.c1 = hbmap[i];
-                    mat.map[i].desc    = hbdesc[i];
-                    mat.map[i].rgb     = hbrgb[i];
+                    m.code.c1 = hbmap[i];
+                    m.desc    = hbdesc[i];
+                    m.rgb     = hbrgb[i];
                 }
                 fp = opt2FILE("-hbm", NFILE, fnm, "w");
                 write_xpm_m(fp, mat);
                 gmx_ffclose(fp);
-                for (x = 0; x < mat.nx; x++)
-                {
-                    sfree(mat.matrix[x]);
-                }
-                sfree(mat.axis_y);
-                sfree(mat.matrix);
-                sfree(mat.map);
             }
             else
             {
index 1890cfb0455f0f17cea3f4b86c2d877be697244b..7aca5d4bf47766d36712881eaaf79acaf1ca5e1b 100644 (file)
@@ -43,6 +43,7 @@
 #include <cstring>
 
 #include <algorithm>
+#include <numeric>
 #include <string>
 
 #include "gromacs/commandline/pargs.h"
@@ -53,6 +54,7 @@
 #include "gromacs/fileio/warninp.h"
 #include "gromacs/fileio/writeps.h"
 #include "gromacs/gmxana/gmx_ana.h"
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/arraysize.h"
 #include "gromacs/utility/cstringutil.h"
 #include "gromacs/utility/exceptions.h"
@@ -195,48 +197,10 @@ static t_rgb white = { 1, 1, 1 };
 /* this must correspond to *colors[] in get_params */
 static t_rgb *linecolors[] = { nullptr, &black, &white, nullptr };
 
-static gmx_bool diff_maps(int nmap1, t_mapping *map1, int nmap2, t_mapping *map2)
+static void leg_discrete(t_psdata *ps, real x0, real y0, const std::string &label,
+                         real fontsize, char *font,
+                         gmx::ArrayRef<const t_mapping> map)
 {
-    int      i;
-    gmx_bool bDiff, bColDiff = FALSE;
-
-    if (nmap1 != nmap2)
-    {
-        bDiff = TRUE;
-    }
-    else
-    {
-        bDiff = FALSE;
-        for (i = 0; i < nmap1; i++)
-        {
-            if (!matelmt_cmp(map1[i].code, map2[i].code))
-            {
-                bDiff = TRUE;
-            }
-            if (std::strcmp(map1[i].desc, map2[i].desc) != 0)
-            {
-                bDiff = TRUE;
-            }
-            if ((map1[i].rgb.r != map2[i].rgb.r) ||
-                (map1[i].rgb.g != map2[i].rgb.g) ||
-                (map1[i].rgb.b != map2[i].rgb.b))
-            {
-                bColDiff = TRUE;
-            }
-        }
-        if (!bDiff && bColDiff)
-        {
-            fprintf(stderr, "Warning: two colormaps differ only in RGB value, using one colormap.\n");
-        }
-    }
-
-    return bDiff;
-}
-
-static void leg_discrete(t_psdata *ps, real x0, real y0, char *label,
-                         real fontsize, char *font, int nmap, t_mapping map[])
-{
-    int   i;
     real  yhh;
     real  boxhh;
 
@@ -245,87 +209,88 @@ static void leg_discrete(t_psdata *ps, real x0, real y0, char *label,
     ps_rgb(ps, BLACK);
     ps_strfont(ps, font, fontsize);
     yhh = y0+fontsize+3*DDD;
-    if (std::strlen(label) > 0)
+    if (!label.empty())
     {
         ps_ctext(ps, x0, yhh, label, eXLeft);
     }
     ps_moveto(ps, x0, y0);
-    for (i = 0; (i < nmap); i++)
+    for (const auto &m : map)
     {
         ps_setorigin(ps);
-        ps_rgb(ps, &(map[i].rgb));
+        ps_rgb(ps, &m.rgb);
         ps_fillbox(ps, DDD, DDD, DDD+fontsize, boxhh-DDD);
         ps_rgb(ps, BLACK);
         ps_box(ps, DDD, DDD, DDD+fontsize, boxhh-DDD);
-        ps_ctext(ps, boxhh+2*DDD, fontsize/3, map[i].desc, eXLeft);
+        ps_ctext(ps, boxhh+2*DDD, fontsize/3, m.desc, eXLeft);
         ps_unsetorigin(ps);
         ps_moverel(ps, DDD, -fontsize/3);
     }
 }
 
-static void leg_continuous(t_psdata *ps, real x0, real x, real y0, char *label,
+static void leg_continuous(t_psdata *ps, real x0, real x, real y0, const std::string &label,
                            real fontsize, char *font,
-                           int nmap, t_mapping map[],
+                           gmx::ArrayRef<const t_mapping> map,
                            int mapoffset)
 {
-    int   i;
-    real  xx0;
-    real  yhh, boxxh, boxyh;
+    real       xx0;
+    real       yhh, boxxh, boxyh;
+    gmx::index mapIndex = gmx::ssize(map) - mapoffset;
 
     boxyh = fontsize;
     if (x < 8*fontsize)
     {
         x = 8*fontsize;
     }
-    boxxh = x/(nmap-mapoffset);
+    boxxh = x/mapIndex;
     if (boxxh > fontsize)
     {
         boxxh = fontsize;
     }
 
-    GMX_RELEASE_ASSERT(map != nullptr, "NULL map array provided to leg_continuous()");
+    GMX_RELEASE_ASSERT(!map.empty(), "NULL map array provided to leg_continuous()");
 
     /* LANDSCAPE */
-    xx0 = x0-((nmap-mapoffset)*boxxh)/2.0;
+    xx0 = x0-(mapIndex*boxxh)/2.0;
 
-    for (i = 0; (i < nmap-mapoffset); i++)
+    for (gmx::index i = 0; (i < mapIndex); i++)
     {
         ps_rgb(ps, &(map[i+mapoffset].rgb));
         ps_fillbox(ps, xx0+i*boxxh, y0, xx0+(i+1)*boxxh, y0+boxyh);
     }
     ps_strfont(ps, font, fontsize);
     ps_rgb(ps, BLACK);
-    ps_box(ps, xx0, y0, xx0+(nmap-mapoffset)*boxxh, y0+boxyh);
+    ps_box(ps, xx0, y0, xx0+mapIndex*boxxh, y0+boxyh);
 
     yhh = y0+boxyh+3*DDD;
     ps_ctext(ps, xx0+boxxh/2, yhh, map[0].desc, eXCenter);
-    if (std::strlen(label) > 0)
+    if (!label.empty())
     {
         ps_ctext(ps, x0, yhh, label, eXCenter);
     }
-    ps_ctext(ps, xx0+((nmap-mapoffset)*boxxh)
-             - boxxh/2, yhh, map[nmap-1].desc, eXCenter);
+    ps_ctext(ps, xx0+(mapIndex*boxxh)
+             - boxxh/2, yhh, map[map.size()-1].desc, eXCenter);
 }
 
-static void leg_bicontinuous(t_psdata *ps, real x0, real x, real y0, char *label1,
-                             char *label2, real fontsize, char *font,
-                             int nmap1, t_mapping map1[], int nmap2, t_mapping map2[])
+static void leg_bicontinuous(t_psdata *ps, real x0, real x, real y0, const std::string &label1,
+                             const std::string &label2, real fontsize, char *font,
+                             gmx::ArrayRef<const t_mapping> map1,
+                             gmx::ArrayRef<const t_mapping> map2)
 {
     real xx1, xx2, x1, x2;
 
-    x1  = x/(nmap1+nmap2)*nmap1; /* width of legend 1 */
-    x2  = x/(nmap1+nmap2)*nmap2; /* width of legend 2 */
-    xx1 = x0-(x2/2.0)-fontsize;  /* center of legend 1 */
-    xx2 = x0+(x1/2.0)+fontsize;  /* center of legend 2 */
-    x1 -= fontsize/2;            /* adjust width */
-    x2 -= fontsize/2;            /* adjust width */
-    leg_continuous(ps, xx1, x1, y0, label1, fontsize, font, nmap1, map1, 0);
-    leg_continuous(ps, xx2, x2, y0, label2, fontsize, font, nmap2, map2, 0);
+    x1  = x/(map1.size()+map2.size())*map1.size(); /* width of legend 1 */
+    x2  = x/(map1.size()+map2.size())*map2.size(); /* width of legend 2 */
+    xx1 = x0-(x2/2.0)-fontsize;                    /* center of legend 1 */
+    xx2 = x0+(x1/2.0)+fontsize;                    /* center of legend 2 */
+    x1 -= fontsize/2;                              /* adjust width */
+    x2 -= fontsize/2;                              /* adjust width */
+    leg_continuous(ps, xx1, x1, y0, label1, fontsize, font, map1, 0);
+    leg_continuous(ps, xx2, x2, y0, label2, fontsize, font, map2, 0);
 }
 
-static real box_height(t_matrix *mat, t_psrec *psr)
+static real box_height(const t_matrix &mat, t_psrec *psr)
 {
-    return mat->ny*psr->yboxsize;
+    return mat.ny*psr->yboxsize;
 }
 
 static real box_dh(t_psrec *psr)
@@ -333,7 +298,6 @@ static real box_dh(t_psrec *psr)
     return psr->boxspacing;
 }
 
-#define IS_ONCE (i == nmat-1)
 static real box_dh_top(gmx_bool bOnce, t_psrec *psr)
 {
     real dh;
@@ -361,14 +325,13 @@ static gmx_bool box_do_all_x_min_ticks(t_psrec *psr)
 }
 
 static void draw_boxes(t_psdata *ps, real x0, real y0, real w,
-                       int nmat, t_matrix mat[], t_psrec *psr)
+                       gmx::ArrayRef<t_matrix> mat, t_psrec *psr)
 {
     char     buf[128];
-    char    *mylab;
     real     xxx;
     char   **xtick, **ytick;
     real     xx, yy, dy, xx00, yy00, offset_x, offset_y;
-    int      i, j, x, y, ntx, nty;
+    int      x, ntx, nty;
     size_t   strlength;
 
     /* Only necessary when there will be no y-labels */
@@ -378,70 +341,71 @@ static void draw_boxes(t_psdata *ps, real x0, real y0, real w,
     ps_rgb(ps, BLACK);
     ps_linewidth(ps, static_cast<int>(psr->boxlinewidth));
     yy00 = y0;
-    for (i = 0; (i < nmat); i++)
+    for (auto m = mat.begin(); m != mat.end(); ++m)
     {
-        dy = box_height(&(mat[i]), psr);
+        dy = box_height(*m, psr);
         ps_box(ps, x0-1, yy00-1, x0+w+1, yy00+dy+1);
-        yy00 += dy+box_dh(psr)+box_dh_top(IS_ONCE, psr);
+        yy00 += dy+box_dh(psr)+box_dh_top(m+1 == mat.end(), psr);
     }
 
     /* Draw the ticks on the axes */
     ps_linewidth(ps, static_cast<int>(psr->ticklinewidth));
     xx00 = x0-1;
     yy00 = y0-1;
-    for (i = 0; (i < nmat); i++)
+    auto halfway = mat.begin() + (mat.size()/2);
+    for (auto m = mat.begin(); m != mat.end(); ++m)
     {
-        if (mat[i].flags & MAT_SPATIAL_X)
+        if (m->flags & MAT_SPATIAL_X)
         {
-            ntx      = mat[i].nx + 1;
+            ntx      = m->nx + 1;
             offset_x = 0.1;
         }
         else
         {
-            ntx      = mat[i].nx;
+            ntx      = m->nx;
             offset_x = 0.6;
         }
-        if (mat[i].flags & MAT_SPATIAL_Y)
+        if (m->flags & MAT_SPATIAL_Y)
         {
-            nty      = mat[i].ny + 1;
+            nty      = m->ny + 1;
             offset_y = 0.1;
         }
         else
         {
-            nty      = mat[i].ny;
+            nty      = m->ny;
             offset_y = 0.6;
         }
         snew(xtick, ntx);
-        for (j = 0; (j < ntx); j++)
+        for (int j = 0; (j < ntx); j++)
         {
-            sprintf(buf, "%g", mat[i].axis_x[j]);
+            sprintf(buf, "%g", m->axis_x[j]);
             xtick[j] = gmx_strdup(buf);
         }
         ps_strfont(ps, psr->X.tickfont, psr->X.tickfontsize);
         for (x = 0; (x < ntx); x++)
         {
             xx = xx00 + (x + offset_x)*psr->xboxsize;
-            if ( ( bRmod(mat[i].axis_x[x], psr->X.offset, psr->X.major) ||
+            if ( ( bRmod(m->axis_x[x], psr->X.offset, psr->X.major) ||
                    (psr->X.first && (x == 0))) &&
-                 ( (i == 0) || box_do_all_x_maj_ticks(psr) ) )
+                 ( m == mat.begin() || box_do_all_x_maj_ticks(psr) ) )
             {
                 /* Longer tick marks */
                 ps_line (ps, xx, yy00, xx, yy00-psr->X.majorticklen);
                 /* Plot label on lowest graph only */
-                if (i == 0)
+                if (m == mat.begin())
                 {
                     ps_ctext(ps, xx,
                              yy00-DDD-psr->X.majorticklen-psr->X.tickfontsize*0.8,
                              xtick[x], eXCenter);
                 }
             }
-            else if (bRmod(mat[i].axis_x[x], psr->X.offset, psr->X.minor) &&
-                     ( (i == 0) || box_do_all_x_min_ticks(psr) ) )
+            else if (bRmod(m->axis_x[x], psr->X.offset, psr->X.minor) &&
+                     ( (m == mat.begin()) || box_do_all_x_min_ticks(psr) ) )
             {
                 /* Shorter tick marks */
                 ps_line(ps, xx, yy00, xx, yy00-psr->X.minorticklen);
             }
-            else if (bRmod(mat[i].axis_x[x], psr->X.offset, psr->X.major) )
+            else if (bRmod(m->axis_x[x], psr->X.offset, psr->X.major) )
             {
                 /* Even shorter marks, only each X.major */
                 ps_line(ps, xx, yy00, xx, yy00-(psr->boxspacing/2));
@@ -449,16 +413,16 @@ static void draw_boxes(t_psdata *ps, real x0, real y0, real w,
         }
         ps_strfont(ps, psr->Y.tickfont, psr->Y.tickfontsize);
         snew(ytick, nty);
-        for (j = 0; (j < nty); j++)
+        for (int j = 0; (j < nty); j++)
         {
-            sprintf(buf, "%g", mat[i].axis_y[j]);
+            sprintf(buf, "%g", m->axis_y[j]);
             ytick[j] = gmx_strdup(buf);
         }
 
-        for (y = 0; (y < nty); y++)
+        for (int y = 0; (y < nty); y++)
         {
             yy = yy00 + (y + offset_y)*psr->yboxsize;
-            if (bRmod(mat[i].axis_y[y], psr->Y.offset, psr->Y.major) ||
+            if (bRmod(m->axis_y[y], psr->Y.offset, psr->Y.major) ||
                 (psr->Y.first && (y == 0)))
             {
                 /* Major ticks */
@@ -467,7 +431,7 @@ static void draw_boxes(t_psdata *ps, real x0, real y0, real w,
                 ps_ctext(ps, xx00-psr->Y.majorticklen-DDD,
                          yy-psr->Y.tickfontsize/3.0, ytick[y], eXRight);
             }
-            else if (bRmod(mat[i].axis_y[y], psr->Y.offset, psr->Y.minor) )
+            else if (bRmod(m->axis_y[y], psr->Y.offset, psr->Y.minor) )
             {
                 /* Minor ticks */
                 ps_line(ps, xx00, yy, xx00-psr->Y.minorticklen, yy);
@@ -477,30 +441,32 @@ static void draw_boxes(t_psdata *ps, real x0, real y0, real w,
         sfree(ytick);
 
         /* Label on Y-axis */
-        if (!psr->bYonce || i == nmat/2)
+        if (!psr->bYonce || m == halfway)
         {
+            std::string mylab;
             if (strlen(psr->Y.label) > 0)
             {
                 mylab = psr->Y.label;
             }
             else
             {
-                mylab = mat[i].label_y;
+                mylab = m->label_y;
             }
-            if (strlen(mylab) > 0)
+            if (!mylab.empty())
             {
                 ps_strfont(ps, psr->Y.font, psr->Y.fontsize);
                 ps_flip(ps, TRUE);
                 xxx = x0-psr->X.majorticklen-psr->X.tickfontsize*strlength-DDD;
-                ps_ctext(ps, yy00+box_height(&mat[i], psr)/2.0, 612.5-xxx,
+                ps_ctext(ps, yy00+box_height(*m, psr)/2.0, 612.5-xxx,
                          mylab, eXCenter);
                 ps_flip(ps, FALSE);
             }
         }
 
-        yy00 += box_height(&(mat[i]), psr)+box_dh(psr)+box_dh_top(IS_ONCE, psr);
+        yy00 += box_height(*m, psr)+box_dh(psr)+box_dh_top(m+1 == mat.end(), psr);
     }
     /* Label on X-axis */
+    std::string mylab;
     if (strlen(psr->X.label) > 0)
     {
         mylab = psr->X.label;
@@ -509,7 +475,7 @@ static void draw_boxes(t_psdata *ps, real x0, real y0, real w,
     {
         mylab = mat[0].label_x;
     }
-    if (strlen(mylab) > 0)
+    if (!mylab.empty())
     {
         ps_strfont(ps, psr->X.font, psr->X.fontsize);
         ps_ctext(ps, x0+w/2, y0-DDD-psr->X.majorticklen-psr->X.tickfontsize*FUDGE-
@@ -518,28 +484,28 @@ static void draw_boxes(t_psdata *ps, real x0, real y0, real w,
 }
 
 static void draw_zerolines(t_psdata *out, real x0, real y0, real w,
-                           int nmat, t_matrix mat[], t_psrec *psr)
+                           gmx::ArrayRef<t_matrix> mat, t_psrec *psr)
 {
     real   xx, yy, dy, xx00, yy00;
-    int    i, x, y;
+    int    x, y;
 
     xx00 = x0-1.5;
     yy00 = y0-1.5;
     ps_linewidth(out, static_cast<int>(psr->zerolinewidth));
-    for (i = 0; (i < nmat); i++)
+    for (auto m = mat.begin(); m != mat.end(); ++m)
     {
-        dy = box_height(&(mat[i]), psr);
-        /* mat[i].axis_x and _y were already set by draw_boxes */
+        dy = box_height(*m, psr);
+        /* m->axis_x and _y were already set by draw_boxes */
         if (psr->X.lineatzero)
         {
             ps_rgb(out, linecolors[psr->X.lineatzero]);
-            for (x = 0; (x < mat[i].nx); x++)
+            for (x = 0; (x < m->nx); x++)
             {
                 xx = xx00+(x+0.7)*psr->xboxsize;
                 /* draw lines whenever tick label almost zero (e.g. next trajectory) */
-                if (x != 0 && x < mat[i].nx-1 &&
-                    std::abs(mat[i].axis_x[x]) <
-                    0.1*std::abs(mat[i].axis_x[x+1]-mat[i].axis_x[x]) )
+                if (x != 0 && x < m->nx-1 &&
+                    std::abs(m->axis_x[x]) <
+                    0.1*std::abs(m->axis_x[x+1]-m->axis_x[x]) )
                 {
                     ps_line (out, xx, yy00, xx, yy00+dy+2);
                 }
@@ -548,38 +514,39 @@ static void draw_zerolines(t_psdata *out, real x0, real y0, real w,
         if (psr->Y.lineatzero)
         {
             ps_rgb(out, linecolors[psr->Y.lineatzero]);
-            for (y = 0; (y < mat[i].ny); y++)
+            for (y = 0; (y < m->ny); y++)
             {
                 yy = yy00+(y+0.7)*psr->yboxsize;
                 /* draw lines whenever tick label almost zero (e.g. next trajectory) */
-                if (y != 0 && y < mat[i].ny-1 &&
-                    std::abs(mat[i].axis_y[y]) <
-                    0.1*std::abs(mat[i].axis_y[y+1]-mat[i].axis_y[y]) )
+                if (y != 0 && y < m->ny-1 &&
+                    std::abs(m->axis_y[y]) <
+                    0.1*std::abs(m->axis_y[y+1]-m->axis_y[y]) )
                 {
                     ps_line (out, xx00, yy, xx00+w+2, yy);
                 }
             }
         }
-        yy00 += box_height(&(mat[i]), psr)+box_dh(psr)+box_dh_top(IS_ONCE, psr);
+        yy00 += box_height(*m, psr)+box_dh(psr)+box_dh_top(m+1 == mat.end(), psr);
     }
 }
 
-static void box_dim(int nmat, t_matrix mat[], t_matrix *mat2, t_psrec *psr,
+static void box_dim(gmx::ArrayRef<t_matrix> mat,
+                    gmx::ArrayRef<t_matrix> mat2, t_psrec *psr,
                     int elegend, gmx_bool bFrame,
                     real *w, real *h, real *dw, real *dh)
 {
-    int  i, maxytick;
+    int  maxytick;
     real ww, hh, dww, dhh;
 
     hh       = dww = dhh = 0;
     maxytick = 0;
 
     ww = 0;
-    for (i = 0; (i < nmat); i++)
+    for (const auto &m : mat)
     {
-        ww       = std::max(ww, mat[i].nx*psr->xboxsize);
-        hh      += box_height(&(mat[i]), psr);
-        maxytick = std::max(maxytick, mat[i].nx);
+        ww       = std::max(ww, m.nx*psr->xboxsize);
+        hh      += box_height(m, psr);
+        maxytick = std::max(maxytick, m.nx);
     }
     if (bFrame)
     {
@@ -602,9 +569,9 @@ static void box_dim(int nmat, t_matrix mat[], t_matrix *mat2, t_psrec *psr,
             dhh += psr->X.fontsize+2*DDD;
         }
         if ( /* fool emacs auto-indent */
-            (elegend == elBoth && (mat[0].legend[0] || (mat2 && mat2[0].legend[0]))) ||
+            (elegend == elBoth && (mat[0].legend[0] || (!mat2.empty() && mat2[0].legend[0]))) ||
             (elegend == elFirst && mat[0].legend[0]) ||
-            (elegend == elSecond && (mat2 && mat2[0].legend[0])) )
+            (elegend == elSecond && (!mat2.empty() && mat2[0].legend[0])) )
         {
             dhh += 2*(psr->legfontsize*FUDGE+2*DDD);
         }
@@ -621,11 +588,11 @@ static void box_dim(int nmat, t_matrix mat[], t_matrix *mat2, t_psrec *psr,
             dhh += psr->X.minorticklen;
         }
 
-        hh += (nmat-1)*box_dh(psr);
+        hh += (mat.size()-1)*box_dh(psr);
         hh += box_dh_top(TRUE, psr);
-        if (nmat > 1)
+        if (mat.size() > 1)
         {
-            hh += (nmat-1)*box_dh_top(FALSE, psr);
+            hh += (mat.size()-1)*box_dh_top(FALSE, psr);
         }
     }
     *w  = ww;
@@ -634,40 +601,36 @@ static void box_dim(int nmat, t_matrix mat[], t_matrix *mat2, t_psrec *psr,
     *dh = dhh;
 }
 
-static int add_maps(t_mapping **newmap,
-                    int nmap1, t_mapping map1[], int nmap2, t_mapping map2[])
+static std::vector<t_mapping> add_maps(gmx::ArrayRef<t_mapping> map1,
+                                       gmx::ArrayRef<t_mapping> map2)
 {
-    static char mapper[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789!@#$%^&*()-_=+{}|;:',<.>/?";
-    int         nsymbols;
-    int         nmap, j, k;
-    t_mapping  *map;
+    static char            mapper[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789!@#$%^&*()-_=+{}|;:',<.>/?";
+    std::vector<t_mapping> map(map1.size() + map2.size());
 
-    nsymbols = std::strlen(mapper);
-    nmap     = nmap1+nmap2;
-    if (nmap > nsymbols*nsymbols)
+    size_t                 nsymbols = std::strlen(mapper);
+    if (map.size() > nsymbols*nsymbols)
     {
         gmx_fatal(FARGS, "Not enough symbols to merge the two colormaps\n");
     }
-    printf("Combining colormaps of %d and %d elements into one of %d elements\n",
-           nmap1, nmap2, nmap);
-    snew(map, nmap);
-    for (j = 0; j < nmap1; j++)
+    printf("Combining colormaps of %zu and %zu elements into one of %zu elements\n",
+           map1.size(), map2.size(), map.size());
+    gmx::index k = 0;
+    for (gmx::index j = 0; j < gmx::ssize(map1) && k < gmx::ssize(map); ++j, ++k)
     {
-        map[j].code.c1 = mapper[j % nsymbols];
-        if (nmap > nsymbols)
+        map[k].code.c1 = mapper[k % nsymbols];
+        if (map.size() > nsymbols)
         {
-            map[j].code.c2 = mapper[j/nsymbols];
+            map[k].code.c2 = mapper[k/nsymbols];
         }
-        map[j].rgb.r = map1[j].rgb.r;
-        map[j].rgb.g = map1[j].rgb.g;
-        map[j].rgb.b = map1[j].rgb.b;
-        map[j].desc  = map1[j].desc;
+        map[k].rgb.r = map1[j].rgb.r;
+        map[k].rgb.g = map1[j].rgb.g;
+        map[k].rgb.b = map1[j].rgb.b;
+        map[k].desc  = map1[j].desc;
     }
-    for (j = 0; j < nmap2; j++)
+    for (gmx::index j = 0; j < gmx::ssize(map2) && k < gmx::ssize(map); ++j, ++k)
     {
-        k              = j+nmap1;
         map[k].code.c1 = mapper[k % nsymbols];
-        if (nmap > nsymbols)
+        if (map.size() > nsymbols)
         {
             map[k].code.c2 = mapper[k/nsymbols];
         }
@@ -677,61 +640,65 @@ static int add_maps(t_mapping **newmap,
         map[k].desc  = map2[j].desc;
     }
 
-    *newmap = map;
-    return nmap;
+    return map;
 }
 
-static void xpm_mat(const char *outf, int nmat, t_matrix *mat, t_matrix *mat2,
+static bool operator==(const t_mapping &lhs, const t_mapping &rhs)
+{
+    return (lhs.rgb.r == rhs.rgb.r && lhs.rgb.g == rhs.rgb.g && lhs.rgb.b == rhs.rgb.b);
+}
+
+static void xpm_mat(const char *outf,
+                    gmx::ArrayRef<t_matrix> mat,
+                    gmx::ArrayRef<t_matrix> mat2,
                     gmx_bool bDiag, gmx_bool bFirstDiag)
 {
     FILE      *out;
-    int        i, x, y, col;
-    int        nmap;
-    t_mapping *map = nullptr;
+    int        x, y, col;
 
     out = gmx_ffopen(outf, "w");
 
-    for (i = 0; i < nmat; i++)
+    GMX_RELEASE_ASSERT(mat.size() == mat2.size(), "Combined matrix write requires matrices of the same size");
+    for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
     {
-        if (!mat2 || !diff_maps(mat[i].nmap, mat[i].map, mat2[i].nmap, mat2[i].map))
+        // Color maps that differ only in RGB value are considered different
+        if (mat2.empty() || std::equal(mat[i].map.begin(), mat[i].map.end(), mat2[i].map.begin()))
         {
             write_xpm_m(out, mat[0]);
         }
         else
         {
-            nmap = add_maps(&map, mat[i].nmap, mat[i].map, mat2[i].nmap, mat2[i].map);
+            auto map = add_maps(mat[i].map, mat2[i].map);
             for (x = 0; (x < mat[i].nx); x++)
             {
                 for (y = 0; (y < mat[i].nx); y++)
                 {
                     if ((x < y) || ((x == y) && bFirstDiag)) /* upper left  -> map1 */
                     {
-                        col = mat[i].matrix[x][y];
+                        col = mat[i].matrix(x, y);
                     }
                     else /* lower right -> map2 */
                     {
-                        col = mat[i].nmap+mat[i].matrix[x][y];
+                        col = mat[i].map.size() + mat[i].matrix(x, y);
                     }
                     if ((bDiag) || (x != y))
                     {
-                        mat[i].matrix[x][y] = col;
+                        mat[i].matrix(x, y) = col;
                     }
                     else
                     {
-                        mat[i].matrix[x][y] = 0;
+                        mat[i].matrix(x, y) = 0;
                     }
                 }
             }
-            sfree(mat[i].map);
-            mat[i].nmap = nmap;
             mat[i].map  = map;
-            if (std::strcmp(mat[i].title, mat2[i].title) != 0)
+            if (mat[i].title != mat2[i].title)
             {
-                sprintf(mat[i].title+strlen(mat[i].title), " / %s", mat2[i].title);
+                mat[i].title += " / " + mat2[i].title;
             }
-            if (std::strcmp(mat[i].legend, mat2[i].legend) != 0)
+            if (mat[i].legend != mat2[i].legend)
             {
-                sprintf(mat[i].legend+strlen(mat[i].legend), " / %s", mat2[i].legend);
+                mat[i].legend += " / " + mat2[i].legend;
             }
             write_xpm_m(out, mat[i]);
         }
@@ -784,20 +751,19 @@ static void tick_spacing(int n, real axis[], real offset, char axisnm,
             axisnm, *major, *minor);
 }
 
-static void ps_mat(const char *outf, int nmat, t_matrix mat[], t_matrix mat2[],
+static void ps_mat(const char *outf,
+                   gmx::ArrayRef<t_matrix> mat,
+                   gmx::ArrayRef<t_matrix> mat2,
                    gmx_bool bFrame, gmx_bool bDiag, gmx_bool bFirstDiag,
                    gmx_bool bTitle, gmx_bool bTitleOnce, gmx_bool bYonce, int elegend,
                    real size, real boxx, real boxy, const char *m2p, const char *m2pout,
                    int mapoffset)
 {
-    char         *legend;
     t_psrec       psrec, *psr;
     int           W, H;
-    int           i, x, y, col, leg = 0;
+    int           x, y, col;
     real          x0, y0, xx;
     real          w, h, dw, dh;
-    int           nmap1 = 0, nmap2 = 0, leg_nmap;
-    t_mapping    *map1  = nullptr, *map2 = nullptr, *leg_map;
     gmx_bool      bMap1, bNextMap1, bDiscrete;
 
     try
@@ -811,7 +777,7 @@ static void ps_mat(const char *outf, int nmat, t_matrix mat[], t_matrix mat2[],
     if (psr->X.major <= 0)
     {
         tick_spacing((mat[0].flags & MAT_SPATIAL_X) ? mat[0].nx + 1 : mat[0].nx,
-                     mat[0].axis_x, psr->X.offset, 'X',
+                     mat[0].axis_x.data(), psr->X.offset, 'X',
                      &(psr->X.major), &(psr->X.minor) );
     }
     if (psr->X.minor <= 0)
@@ -821,7 +787,7 @@ static void ps_mat(const char *outf, int nmat, t_matrix mat[], t_matrix mat2[],
     if (psr->Y.major <= 0)
     {
         tick_spacing((mat[0].flags & MAT_SPATIAL_Y) ? mat[0].ny + 1 : mat[0].ny,
-                     mat[0].axis_y, psr->Y.offset, 'Y',
+                     mat[0].axis_y.data(), psr->Y.offset, 'Y',
                      &(psr->Y.major), &(psr->Y.minor) );
     }
     if (psr->Y.minor <= 0)
@@ -850,50 +816,49 @@ static void ps_mat(const char *outf, int nmat, t_matrix mat[], t_matrix mat2[],
         printf("Set the y-size of the box to %.3f\n", psr->yboxsize);
     }
 
-    nmap1 = 0;
-    for (i = 0; (i < nmat); i++)
+    gmx::ArrayRef<const t_mapping> map1;
+    int legendIndex = 0;
+    for (const auto &m : mat)
     {
-        if (mat[i].nmap > nmap1)
+        if (m.map.size() > map1.size())
         {
-            nmap1 = mat[i].nmap;
-            map1  = mat[i].map;
-            leg   = i+1;
+            if (map1.empty())
+            {
+                printf("Selected legend of matrix # %d for display\n", legendIndex);
+            }
+            map1 = m.map;
         }
+        ++legendIndex;
     }
-    if (leg != 1)
+    gmx::ArrayRef<const t_mapping> map2;
+    if (!mat2.empty())
     {
-        printf("Selected legend of matrix # %d for display\n", leg);
-    }
-    if (mat2)
-    {
-        nmap2 = 0;
-        for (i = 0; (i < nmat); i++)
+        for (const auto &m : mat2)
         {
-            if (mat2[i].nmap > nmap2)
+            if (m.map.size() > map2.size())
             {
-                nmap2 = mat2[i].nmap;
-                map2  = mat2[i].map;
-                leg   = i+1;
+                if (map2.empty())
+                {
+                    printf("Selected legend of matrix # %d for second display\n", legendIndex);
+                }
+                map2 = m.map;
             }
-        }
-        if (leg != 1)
-        {
-            printf("Selected legend of matrix # %d for second display\n", leg);
+            ++legendIndex;
         }
     }
-    if ( (mat[0].legend[0] == 0) && psr->legend)
+    if (mat[0].legend.empty() && psr->legend)
     {
-        std::strcpy(mat[0].legend, psr->leglabel);
+        mat[0].legend = psr->leglabel;
     }
 
-    bTitle          = bTitle     && (mat[nmat-1].title[0] != 0);
-    bTitleOnce      = bTitleOnce && (mat[nmat-1].title[0] != 0);
+    bTitle          = bTitle     && !mat.back().title.empty();
+    bTitleOnce      = bTitleOnce && !mat.back().title.empty();
     psr->bTitle     = bTitle;
     psr->bTitleOnce = bTitleOnce;
     psr->bYonce     = bYonce;
 
     /* Set up size of box for nice colors */
-    box_dim(nmat, mat, mat2, psr, elegend, bFrame, &w, &h, &dw, &dh);
+    box_dim(mat, mat2, psr, elegend, bFrame, &w, &h, &dw, &dh);
 
     /* Set up bounding box */
     W = static_cast<int>(w+dw);
@@ -918,29 +883,29 @@ static void ps_mat(const char *outf, int nmat, t_matrix mat[], t_matrix mat2[],
     if (bFrame)
     {
         ps_comment(&out, "Here starts the BOX drawing");
-        draw_boxes(&out, x0, y0, w, nmat, mat, psr);
+        draw_boxes(&out, x0, y0, w, mat, psr);
     }
 
-    for (i = 0; (i < nmat); i++)
+    for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
     {
-        if (bTitle || (bTitleOnce && i == nmat-1) )
+        if (bTitle || (bTitleOnce && i == gmx::ssize(mat)-1) )
         {
             /* Print title, if any */
             ps_rgb(&out, BLACK);
             ps_strfont(&out, psr->titfont, psr->titfontsize);
             std::string buf;
-            if (!mat2 || (std::strcmp(mat[i].title, mat2[i].title) == 0))
+            if (!mat2.empty() || mat[i].title == mat2[i].title)
             {
                 buf = mat[i].title;
             }
             else
             {
-                buf = gmx::formatString("%s / %s", mat[i].title, mat2[i].title);
+                buf = mat[i].title + " / " + mat2[i].title;
             }
-            ps_ctext(&out, x0+w/2, y0+box_height(&(mat[i]), psr)+psr->titfontsize,
+            ps_ctext(&out, x0+w/2, y0+box_height(mat[i], psr)+psr->titfontsize,
                      buf, eXCenter);
         }
-        ps_comment(&out, gmx::formatString("Here starts the filling of box #%d", i).c_str());
+        ps_comment(&out, gmx::formatString("Here starts the filling of box #%zd", i).c_str());
         for (x = 0; (x < mat[i].nx); x++)
         {
             int nexty;
@@ -949,10 +914,10 @@ static void ps_mat(const char *outf, int nmat, t_matrix mat[], t_matrix mat2[],
             xx = x0+x*psr->xboxsize;
             ps_moveto(&out, xx, y0);
             y     = 0;
-            bMap1 = ((mat2 == nullptr) || (x < y || (x == y && bFirstDiag)));
+            bMap1 = (mat2.empty() || (x < y || (x == y && bFirstDiag)));
             if ((bDiag) || (x != y))
             {
-                col = mat[i].matrix[x][y];
+                col = mat[i].matrix(x, y);
             }
             else
             {
@@ -960,7 +925,7 @@ static void ps_mat(const char *outf, int nmat, t_matrix mat[], t_matrix mat2[],
             }
             for (nexty = 1; (nexty <= mat[i].ny); nexty++)
             {
-                bNextMap1 = ((mat2 == nullptr) || (x < nexty || (x == nexty && bFirstDiag)));
+                bNextMap1 = (mat2.empty() || (x < nexty || (x == nexty && bFirstDiag)));
                 /* TRUE:  upper left  -> map1 */
                 /* FALSE: lower right -> map2 */
                 if ((nexty == mat[i].ny) || (!bDiag && (x == nexty)))
@@ -969,7 +934,7 @@ static void ps_mat(const char *outf, int nmat, t_matrix mat[], t_matrix mat2[],
                 }
                 else
                 {
-                    nextcol = mat[i].matrix[x][nexty];
+                    nextcol = mat[i].matrix(x, nexty);
                 }
                 if ( (nexty == mat[i].ny) || (col != nextcol) || (bMap1 != bNextMap1) )
                 {
@@ -981,7 +946,7 @@ static void ps_mat(const char *outf, int nmat, t_matrix mat[], t_matrix mat2[],
                         }
                         else
                         {
-                            assert(mat2);
+                            assert(!mat2.empty());
                             ps_rgb_nbox(&out, &(mat2[i].map[col].rgb), nexty-y);
                         }
                     }
@@ -995,7 +960,7 @@ static void ps_mat(const char *outf, int nmat, t_matrix mat[], t_matrix mat2[],
                 }
             }
         }
-        y0 += box_height(&(mat[i]), psr)+box_dh(psr)+box_dh_top(IS_ONCE, psr);
+        y0 += box_height(mat[i], psr)+box_dh(psr)+box_dh_top(i+1 == gmx::ssize(mat), psr);
     }
 
     if (psr->X.lineatzero || psr->Y.lineatzero)
@@ -1003,45 +968,45 @@ static void ps_mat(const char *outf, int nmat, t_matrix mat[], t_matrix mat2[],
         /* reset y0 for first box */
         y0 = dh;
         ps_comment(&out, "Here starts the zero lines drawing");
-        draw_zerolines(&out, x0, y0, w, nmat, mat, psr);
+        draw_zerolines(&out, x0, y0, w, mat, psr);
     }
 
     if (elegend != elNone)
     {
+        std::string                    legend;
+        gmx::ArrayRef<const t_mapping> leg_map;
         ps_comment(&out, "Now it's legend time!");
         ps_linewidth(&out, static_cast<int>(psr->linewidth));
-        if (mat2 == nullptr || elegend != elSecond)
+        if (mat2.empty() || elegend != elSecond)
         {
             bDiscrete = mat[0].bDiscrete;
             legend    = mat[0].legend;
-            leg_nmap  = nmap1;
             leg_map   = map1;
         }
         else
         {
             bDiscrete = mat2[0].bDiscrete;
             legend    = mat2[0].legend;
-            leg_nmap  = nmap2;
             leg_map   = map2;
         }
         if (bDiscrete)
         {
             leg_discrete(&out, psr->legfontsize, DDD, legend,
-                         psr->legfontsize, psr->legfont, leg_nmap, leg_map);
+                         psr->legfontsize, psr->legfont, leg_map);
         }
         else
         {
             if (elegend != elBoth)
             {
                 leg_continuous(&out, x0+w/2, w/2, DDD, legend,
-                               psr->legfontsize, psr->legfont, leg_nmap, leg_map,
+                               psr->legfontsize, psr->legfont, leg_map,
                                mapoffset);
             }
             else
             {
-                assert(mat2);
+                assert(!mat2.empty());
                 leg_bicontinuous(&out, x0+w/2, w, DDD, mat[0].legend, mat2[0].legend,
-                                 psr->legfontsize, psr->legfont, nmap1, map1, nmap2, map2);
+                                 psr->legfontsize, psr->legfont, map1, map2);
             }
         }
         ps_comment(&out, "Done processing");
@@ -1050,70 +1015,63 @@ static void ps_mat(const char *outf, int nmat, t_matrix mat[], t_matrix mat2[],
     ps_close(&out);
 }
 
-static void make_axis_labels(int nmat, t_matrix *mat)
+static void make_axis_labels(gmx::ArrayRef<t_matrix> mat1)
 {
-    int i, j;
-
-    for (i = 0; (i < nmat); i++)
+    for (auto &m : mat1)
     {
         /* Make labels for x axis */
-        if (mat[i].axis_x == nullptr)
+        if (m.axis_x.empty())
         {
-            snew(mat[i].axis_x, mat[i].nx);
-            for (j = 0; (j < mat[i].nx); j++)
-            {
-                mat[i].axis_x[j] = j;
-            }
+            m.axis_x.resize(m.nx);
+            std::iota(m.axis_x.begin(), m.axis_x.end(), 0);
         }
         /* Make labels for y axis */
-        if (mat[i].axis_y == nullptr)
+        if (m.axis_y.empty())
         {
-            snew(mat[i].axis_y, mat[i].ny);
-            for (j = 0; (j < mat[i].ny); j++)
-            {
-                mat[i].axis_y[j] = j;
-            }
+            m.axis_y.resize(m.ny);
+            std::iota(m.axis_y.begin(), m.axis_y.end(), 0);
         }
     }
 }
 
-static void prune_mat(int nmat, t_matrix *mat, t_matrix *mat2, int skip)
+static void prune_mat(gmx::ArrayRef<t_matrix> mat,
+                      gmx::ArrayRef<t_matrix> mat2,
+                      int                     skip)
 {
-    int i, x, y, xs, ys;
-
-    for (i = 0; i < nmat; i++)
+    GMX_RELEASE_ASSERT(mat.size() == mat2.size(), "Matrix pruning requires matrices of the same size");
+    for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
     {
         fprintf(stderr, "converting %dx%d matrix to %dx%d\n",
                 mat[i].nx, mat[i].ny,
                 (mat[i].nx+skip-1)/skip, (mat[i].ny+skip-1)/skip);
         /* walk through matrix */
-        xs = 0;
-        for (x = 0; (x < mat[i].nx); x++)
+        int xs = 0;
+        for (int x = 0; (x < mat[i].nx); x++)
         {
             if (x % skip == 0)
             {
                 mat[i].axis_x[xs] = mat[i].axis_x[x];
-                if (mat2)
+                if (!mat2.empty())
                 {
                     mat2[i].axis_x[xs] = mat2[i].axis_x[x];
                 }
-                ys = 0;
-                for (y = 0; (y < mat[i].ny); y++)
+                int ys = 0;
+                for (int y = 0; (y < mat[i].ny); y++)
                 {
                     if (x == 0)
                     {
                         mat[i].axis_y[ys] = mat[i].axis_y[y];
-                        if (mat2)
+                        if (!mat2.empty())
                         {
                             mat2[i].axis_y[ys] = mat2[i].axis_y[y];
                         }
                     }
                     if (y % skip == 0)
                     {
-                        mat[i].matrix[xs][ys] = mat[i].matrix[x][y];
-                        if (mat2)
+                        mat[i].matrix(xs, ys) = mat[i].matrix(x, y);
+                        if (!mat2.empty())
                         {
-                            mat2[i].matrix[xs][ys] = mat2[i].matrix[x][y];
+                            mat2[i].matrix(xs, ys) = mat2[i].matrix(x, y);
                         }
                         ys++;
                     }
@@ -1124,7 +1082,7 @@ static void prune_mat(int nmat, t_matrix *mat, t_matrix *mat2, int skip)
         /* adjust parameters */
         mat[i].nx = (mat[i].nx+skip-1)/skip;
         mat[i].ny = (mat[i].ny+skip-1)/skip;
-        if (mat2)
+        if (!mat2.empty())
         {
             mat2[i].nx = (mat2[i].nx+skip-1)/skip;
             mat2[i].ny = (mat2[i].ny+skip-1)/skip;
@@ -1132,15 +1090,15 @@ static void prune_mat(int nmat, t_matrix *mat, t_matrix *mat2, int skip)
     }
 }
 
-static void zero_lines(int nmat, t_matrix *mat, t_matrix *mat2)
+static void zero_lines(gmx::ArrayRef<t_matrix> mat,
+                       gmx::ArrayRef<t_matrix> mat2)
 {
-    int       i, x, y, m;
-    t_matrix *mats;
-
-    for (i = 0; i < nmat; i++)
+    GMX_RELEASE_ASSERT(mat.size() == mat2.size(), "zero_lines requires matrices of the same size");
+    for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
     {
-        for (m = 0; m < (mat2 ? 2 : 1); m++)
+        for (int m = 0; m < (!mat2.empty() ? 2 : 1); m++)
         {
+            gmx::ArrayRef<t_matrix> mats;
             if (m == 0)
             {
                 mats = mat;
@@ -1149,23 +1107,23 @@ static void zero_lines(int nmat, t_matrix *mat, t_matrix *mat2)
             {
                 mats = mat2;
             }
-            for (x = 0; x < mats[i].nx-1; x++)
+            for (int x = 0; x < mats[i].nx-1; x++)
             {
                 if (std::abs(mats[i].axis_x[x+1]) < 1e-5)
                 {
-                    for (y = 0; y < mats[i].ny; y++)
+                    for (int y = 0; y < mats[i].ny; y++)
                     {
-                        mats[i].matrix[x][y] = 0;
+                        mats[i].matrix(x, y) = 0;
                     }
                 }
             }
-            for (y = 0; y < mats[i].ny-1; y++)
+            for (int y = 0; y < mats[i].ny-1; y++)
             {
                 if (std::abs(mats[i].axis_y[y+1]) < 1e-5)
                 {
-                    for (x = 0; x < mats[i].nx; x++)
+                    for (int x = 0; x < mats[i].nx; x++)
                     {
-                        mats[i].matrix[x][y] = 0;
+                        mats[i].matrix(x, y) = 0;
                     }
                 }
             }
@@ -1174,20 +1132,21 @@ static void zero_lines(int nmat, t_matrix *mat, t_matrix *mat2)
 }
 
 static void write_combined_matrix(int ecombine, const char *fn,
-                                  int nmat, t_matrix *mat1, t_matrix *mat2,
+                                  gmx::ArrayRef<t_matrix> mat1,
+                                  gmx::ArrayRef<t_matrix> mat2,
                                   const real *cmin, const real *cmax)
 {
-    int        i, j, k, nlevels;
     FILE      *out;
     real     **rmat1, **rmat2;
     real       rhi, rlo;
 
     out = gmx_ffopen(fn, "w");
-    for (k = 0; k < nmat; k++)
+    GMX_RELEASE_ASSERT(mat1.size() == mat2.size(), "Combined matrix write requires matrices of the same size");
+    for (gmx::index k = 0; k != gmx::ssize(mat1); k++)
     {
         if (mat2[k].nx != mat1[k].nx || mat2[k].ny != mat1[k].ny)
         {
-            gmx_fatal(FARGS, "Size of frame %d in 1st (%dx%d) and 2nd matrix (%dx%d) do"
+            gmx_fatal(FARGS, "Size of frame %zd in 1st (%dx%d) and 2nd matrix (%dx%d) do"
                       " not match.\n", k, mat1[k].nx, mat1[k].ny, mat2[k].nx, mat2[k].ny);
         }
         printf("Combining two %dx%d matrices\n", mat1[k].nx, mat1[k].ny);
@@ -1201,9 +1160,9 @@ static void write_combined_matrix(int ecombine, const char *fn,
         }
         rlo = 1e38;
         rhi = -1e38;
-        for (j = 0; j < mat1[k].ny; j++)
+        for (int j = 0; j < mat1[k].ny; j++)
         {
-            for (i = 0; i < mat1[k].nx; i++)
+            for (int i = 0; i < mat1[k].nx; i++)
             {
                 switch (ecombine)
                 {
@@ -1226,124 +1185,107 @@ static void write_combined_matrix(int ecombine, const char *fn,
         {
             rhi = *cmax;
         }
-        nlevels = std::max(mat1[k].nmap, mat2[k].nmap);
+        int nlevels = static_cast<int>(std::max(mat1[k].map.size(), mat2[k].map.size()));
         if (rhi == rlo)
         {
             fprintf(stderr,
                     "combination results in uniform matrix (%g), no output\n", rhi);
         }
-        /*
-           else if (rlo>=0 || rhi<=0)
-           write_xpm(out, mat1[k].flags, mat1[k].title, mat1[k].legend,
-            mat1[k].label_x, mat1[k].label_y,
-            mat1[k].nx, mat1[k].ny, mat1[k].axis_x, mat1[k].axis_y,
-            rmat1, rlo, rhi, rhi<=0?red:white, rhi<=0?white:blue,
-            &nlevels);
-           else
-           write_xpm3(out, mat2[k].flags, mat1[k].title, mat1[k].legend,
-             mat1[k].label_x, mat1[k].label_y,
-             mat1[k].nx, mat1[k].ny, mat1[k].axis_x, mat1[k].axis_y,
-             rmat1, rlo, 0, rhi, red, white, blue, &nlevels);
-         */
         else
         {
             write_xpm(out, mat1[k].flags, mat1[k].title, mat1[k].legend,
                       mat1[k].label_x, mat1[k].label_y,
-                      mat1[k].nx, mat1[k].ny, mat1[k].axis_x, mat1[k].axis_y,
+                      mat1[k].nx, mat1[k].ny, mat1[k].axis_x.data(), mat1[k].axis_y.data(),
                       rmat1, rlo, rhi, white, black, &nlevels);
         }
     }
     gmx_ffclose(out);
 }
 
-static void do_mat(int nmat, t_matrix *mat, t_matrix *mat2,
+static void do_mat(gmx::ArrayRef<t_matrix> mat,
+                   gmx::ArrayRef<t_matrix> mat2,
                    gmx_bool bFrame, gmx_bool bZeroLine, gmx_bool bDiag, gmx_bool bFirstDiag, gmx_bool bTitle,
                    gmx_bool bTitleOnce, gmx_bool bYonce, int elegend,
                    real size, real boxx, real boxy,
                    const char *epsfile, const char *xpmfile, const char *m2p,
                    const char *m2pout, int skip, int mapoffset)
 {
-    int      i, j, k;
-
-    if (mat2)
+    GMX_RELEASE_ASSERT(mat.size() == mat2.size(), "Combined matrix write requires matrices of the same size");
+    if (!mat2.empty())
     {
-        for (k = 0; (k < nmat); k++)
+        for (gmx::index k = 0; k != gmx::ssize(mat); k++)
         {
             if ((mat2[k].nx != mat[k].nx) || (mat2[k].ny != mat[k].ny))
             {
-                gmx_fatal(FARGS, "WAKE UP!! Size of frame %d in 2nd matrix file (%dx%d) does not match size of 1st matrix (%dx%d) or the other way around.\n",
+                gmx_fatal(FARGS, "WAKE UP!! Size of frame %zd in 2nd matrix file (%dx%d) does not match size of 1st matrix (%dx%d) or the other way around.\n",
                           k, mat2[k].nx, mat2[k].ny, mat[k].nx, mat[k].ny);
             }
-            for (j = 0; (j < mat[k].ny); j++)
+            for (int j = 0; (j < mat[k].ny); j++)
             {
-                for (i = bFirstDiag ? j+1 : j; (i < mat[k].nx); i++)
+                for (int i = bFirstDiag ? j+1 : j; (i < mat[k].nx); i++)
                 {
-                    mat[k].matrix[i][j] = mat2[k].matrix[i][j];
+                    mat[k].matrix(i, j) = mat2[k].matrix(i, j);
                 }
             }
         }
     }
-    for (i = 0; (i < nmat); i++)
+    for (gmx::index i = 0; i != gmx::ssize(mat); i++)
     {
-        fprintf(stderr, "Matrix %d is %d x %d\n", i, mat[i].nx, mat[i].ny);
+        fprintf(stderr, "Matrix %zd is %d x %d\n", i, mat[i].nx, mat[i].ny);
     }
 
-    make_axis_labels(nmat, mat);
+    make_axis_labels(mat);
 
     if (skip > 1)
     {
-        prune_mat(nmat, mat, mat2, skip);
+        prune_mat(mat, mat2, skip);
     }
 
     if (bZeroLine)
     {
-        zero_lines(nmat, mat, mat);
+        zero_lines(mat, mat);
     }
 
     if (epsfile != nullptr)
     {
-        ps_mat(epsfile, nmat, mat, mat2, bFrame, bDiag, bFirstDiag,
+        ps_mat(epsfile, mat, mat2, bFrame, bDiag, bFirstDiag,
                bTitle, bTitleOnce, bYonce, elegend,
                size, boxx, boxy, m2p, m2pout, mapoffset);
     }
     if (xpmfile != nullptr)
     {
-        xpm_mat(xpmfile, nmat, mat, mat2, bDiag, bFirstDiag);
+        xpm_mat(xpmfile, mat, mat2, bDiag, bFirstDiag);
     }
 }
 
-static void gradient_map(const rvec grad, int nmap, t_mapping map[])
+static void gradient_map(const rvec grad, gmx::ArrayRef<t_mapping> map)
 {
-    int  i;
-    real c;
-
-    for (i = 0; i < nmap; i++)
+    int  i        = 0;
+    real fraction = 1.0 / (map.size() - 1.0);
+    for (auto &m : map)
     {
-        c            = i/(nmap-1.0);
-        map[i].rgb.r = 1-c*(1-grad[XX]);
-        map[i].rgb.g = 1-c*(1-grad[YY]);
-        map[i].rgb.b = 1-c*(1-grad[ZZ]);
+        real c    = i*fraction;
+        m.rgb.r = 1-c*(1-grad[XX]);
+        m.rgb.g = 1-c*(1-grad[YY]);
+        m.rgb.b = 1-c*(1-grad[ZZ]);
+        ++i;
     }
 }
 
-static void gradient_mat(rvec grad, int nmat, t_matrix mat[])
+static void gradient_mat(rvec grad, gmx::ArrayRef<t_matrix> mat)
 {
-    int m;
-
-    for (m = 0; m < nmat; m++)
+    for (auto &m : mat)
     {
-        gradient_map(grad, mat[m].nmap, mat[m].map);
+        gradient_map(grad, m.map);
     }
 }
 
-static void rainbow_map(gmx_bool bBlue, int nmap, t_mapping map[])
+static void rainbow_map(gmx_bool bBlue, gmx::ArrayRef<t_mapping> map)
 {
-    int  i;
-    real c, r, g, b;
-
-    for (i = 0; i < nmap; i++)
+    for (auto &m : map)
     {
-        c = (map[i].rgb.r + map[i].rgb.g + map[i].rgb.b)/3;
+        real c = (m.rgb.r + m.rgb.g + m.rgb.b)/3;
+        real r, g, b;
         if (c > 1)
         {
             c = 1;
@@ -1376,19 +1318,17 @@ static void rainbow_map(gmx_bool bBlue, int nmap, t_mapping map[])
             g = std::pow(4.0-4.0*c, 2.0/3.0);
             b = 0;
         }
-        map[i].rgb.r = r;
-        map[i].rgb.g = g;
-        map[i].rgb.b = b;
+        m.rgb.r = r;
+        m.rgb.g = g;
+        m.rgb.b = b;
     }
 }
 
-static void rainbow_mat(gmx_bool bBlue, int nmat, t_matrix mat[])
+static void rainbow_mat(gmx_bool bBlue, gmx::ArrayRef<t_matrix> mat)
 {
-    int m;
-
-    for (m = 0; m < nmat; m++)
+    for (auto &m : mat)
     {
-        rainbow_map(bBlue, mat[m].nmap, mat[m].map);
+        rainbow_map(bBlue, m.map);
     }
 }
 
@@ -1439,8 +1379,7 @@ int gmx_xpm2ps(int argc, char *argv[])
 
     gmx_output_env_t *oenv;
     const char       *fn, *epsfile = nullptr, *xpmfile = nullptr;
-    int               i, nmat, nmat2, etitle, elegend, ediag, erainbow, ecombine;
-    t_matrix         *mat = nullptr, *mat2 = nullptr;
+    int               i, etitle, elegend, ediag, erainbow, ecombine;
     gmx_bool          bTitle, bTitleOnce, bDiag, bFirstDiag, bGrad;
     static gmx_bool   bFrame = TRUE, bZeroLine = FALSE, bYonce = FALSE;
     static real       size   = 400, boxx = 0, boxy = 0, cmin = 0, cmax = 0;
@@ -1553,17 +1492,25 @@ int gmx_xpm2ps(int argc, char *argv[])
     bFirstDiag = ediag != edSecond;
 
     fn   = opt2fn("-f", NFILE, fnm);
-    nmat = read_xpm_matrix(fn, &mat);
-    fprintf(stderr, "There %s %d matri%s in %s\n", (nmat > 1) ? "are" : "is", nmat, (nmat > 1) ? "ces" : "x", fn);
+    std::vector<t_matrix> mat, mat2;
+    mat = read_xpm_matrix(fn);
+    fprintf(stderr, "There %s %zu matri%s in %s\n", (mat.size() > 1) ? "are" : "is", mat.size(), (mat.size() > 1) ? "ces" : "x", fn);
     fn = opt2fn_null("-f2", NFILE, fnm);
     if (fn)
     {
-        nmat2 = read_xpm_matrix(fn, &mat2);
-        fprintf(stderr, "There %s %d matri%s in %s\n", (nmat2 > 1) ? "are" : "is", nmat2, (nmat2 > 1) ? "ces" : "x", fn);
-        if (nmat != nmat2)
+        mat2 = read_xpm_matrix(fn);
+        fprintf(stderr, "There %s %zu matri%s in %s\n", (mat2.size() > 1) ? "are" : "is", mat2.size(), (mat2.size() > 1) ? "ces" : "x", fn);
+        if (mat.size() != mat2.size())
         {
             fprintf(stderr, "Different number of matrices, using the smallest number.\n");
-            nmat = nmat2 = std::min(nmat, nmat2);
+            if (mat.size() > mat2.size())
+            {
+                mat.resize(mat2.size());
+            }
+            else
+            {
+                mat2.resize(mat.size());
+            }
         }
     }
     else
@@ -1576,52 +1523,51 @@ int gmx_xpm2ps(int argc, char *argv[])
                     "         no matrix combination will be performed\n");
         }
         ecombine = 0;
-        nmat2    = 0;
     }
     bTitle     = etitle == etTop;
     bTitleOnce = etitle == etOnce;
     if (etitle == etYlabel)
     {
-        for (i = 0; (i < nmat); i++)
+        for (auto &m : mat)
         {
-            std::strcpy(mat[i].label_y, mat[i].title);
-            if (mat2)
-            {
-                std::strcpy(mat2[i].label_y, mat2[i].title);
-            }
+            m.label_y = m.title;
+        }
+        for (auto &m : mat2)
+        {
+            m.label_y = m.title;
         }
     }
     if (bGrad)
     {
-        gradient_mat(grad, nmat, mat);
-        if (mat2)
+        gradient_mat(grad, mat);
+        if (!mat2.empty())
         {
-            gradient_mat(grad, nmat2, mat2);
+            gradient_mat(grad, mat2);
         }
     }
     if (erainbow != erNo)
     {
-        rainbow_mat(erainbow == erBlue, nmat, mat);
-        if (mat2)
+        rainbow_mat(erainbow == erBlue, mat);
+        if (!mat2.empty())
         {
-            rainbow_mat(erainbow == erBlue, nmat2, mat2);
+            rainbow_mat(erainbow == erBlue, mat2);
         }
     }
 
-    if ((mat2 == nullptr) && (elegend != elNone))
+    if (mat2.empty() && (elegend != elNone))
     {
         elegend = elFirst;
     }
 
     if (ecombine && ecombine != ecHalves)
     {
-        write_combined_matrix(ecombine, xpmfile, nmat, mat, mat2,
+        write_combined_matrix(ecombine, xpmfile, mat, mat2,
                               opt2parg_bSet("-cmin", NPA, pa) ? &cmin : nullptr,
                               opt2parg_bSet("-cmax", NPA, pa) ? &cmax : nullptr);
     }
     else
     {
-        do_mat(nmat, mat, mat2, bFrame, bZeroLine, bDiag, bFirstDiag,
+        do_mat(mat, mat2, bFrame, bZeroLine, bDiag, bFirstDiag,
                bTitle, bTitleOnce, bYonce,
                elegend, size, boxx, boxy, epsfile, xpmfile,
                opt2fn_null("-di", NFILE, fnm), opt2fn_null("-do", NFILE, fnm), skip,
index 38d0d61d490e077b80b103d2a722e83c474a2af1..ad48b283ac7bbdeea90161c0daaacfb4a7d4057c 100644 (file)
 #include <algorithm>
 #include <numeric>
 
+#include "gromacs/math/multidimarray.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/utility/exceptions.h"
 
-#include "multidimarray.h"
-
 namespace gmx
 {
 
index eb11b74475f1ecae6ab903041f2cf4110190735c..15c6ebce6be827cfc5024e9ab3620f793a739dbc 100644 (file)
@@ -44,8 +44,7 @@
 #include "densityfittingforce.h"
 
 #include "gromacs/math/functions.h"
-
-#include "multidimarray.h"
+#include "gromacs/math/multidimarray.h"
 
 namespace gmx
 {
index 8674b5006970266ceecc5dc4b02271c9ecb3b95a..085dbf81ad42b1e7a0b0e80501e531508e3801d1 100644 (file)
 #include <array>
 
 #include "gromacs/math/functions.h"
+#include "gromacs/math/multidimarray.h"
 #include "gromacs/math/utilities.h"
 
-#include "multidimarray.h"
-
 namespace gmx
 {
 
index bc445ef085142b3ce44e892e5921dc847259e305..ef8546fde222c9c1d4c847b34f6f9164f4f05717 100644 (file)
@@ -45,6 +45,7 @@
 
 #include <vector>
 
+#include "gromacs/math/multidimarray.h"
 #include "gromacs/math/vectypes.h"
 #include "gromacs/mdspan/extensions.h"
 #include "gromacs/mdspan/mdspan.h"
@@ -52,8 +53,6 @@
 #include "gromacs/utility/classhelpers.h"
 #include "gromacs/utility/real.h"
 
-#include "multidimarray.h"
-
 namespace gmx
 {
 
index a9718ece68a40cb73c7b3f932fd103e8bd00c75a..9fcec3f3886f9dbd3df934a45e9deb635921e945 100644 (file)
 
 #include <array>
 
+#include "gromacs/math/multidimarray.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/utility/real.h"
 
-#include "multidimarray.h"
-
 namespace gmx
 {
 
index e1e83fa924a6915f06550ecd154dcaa351ef6c7d..9b6aee3f2f0b58ffb24c44959802a7b21c7d0f7f 100644 (file)
@@ -39,6 +39,7 @@
  * \author Christian Blau <cblau@gwdg.de>
  * \ingroup module_math
  * \ingroup module_mdspan
+ * \inlibraryapi
  */
 
 #ifndef GMX_MATH_MULTIDIMARRAY_H_