src/gmxlib/selection/sm_compare.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2009, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /*! \internal \file
  39  * \brief
  40  * Implementation of internal selection method for comparison expressions.
  41  */
  42 #ifdef HAVE_CONFIG_H
  43 #include <config.h>
  44 #endif
  45
  46 #include <maths.h>
  47 #include <macros.h>
  48 #include <smalloc.h>
  49 #include <gmx_fatal.h>
  50
  51 #include <selmethod.h>
  52
  53 /** Defines the comparison operator for comparison expressions. */
  54 typedef enum
  55 {
  56     CMP_INVALID,        /**< Indicates an error */
  57     CMP_LESS,           /**< '<' */
  58     CMP_LEQ,            /**< '<=' */
  59     CMP_GTR,            /**< '>' */
  60     CMP_GEQ,            /**< '>=' */
  61     CMP_EQUAL,          /**< '==' */
  62     CMP_NEQ             /**< '!=' */
  63 } e_comparison_t;
  64
  65 /** The operand has a single value. */
  66 #define CMP_SINGLEVAL  1
  67 /** The operand value is dynamic. */
  68 #define CMP_DYNAMICVAL 2
  69 /** The value is real. */
  70 #define CMP_REALVAL    4
  71 /** The integer array is allocated. */
  72 #define CMP_ALLOCINT   16
  73 /** The real array is allocated. */
  74 #define CMP_ALLOCREAL  32
  75
  76 /*! \internal \brief
  77  * Data structure for comparison expression operand values.
  78  */
  79 typedef struct
  80 {
  81     /** Flags that describe the type of the operand. */
  82     int             flags;
  83     /** (Array of) integer value(s). */
  84     int        *i;
  85     /** (Array of) real value(s). */
  86     real       *r;
  87 } t_compare_value;
  88
  89 /*! \internal \brief
  90  * Data structure for comparison expression evaluation.
  91  */
  92 typedef struct
  93 {
  94     /** Comparison operator as a string. */
  95     char            *cmpop;
  96     /** Comparison operator type. */
  97     e_comparison_t   cmpt;
  98     /** Left value. */
  99     t_compare_value  left;
 100     /** Right value. */
 101     t_compare_value  right;
 102 } t_methoddata_compare;
 103
 104 /** Allocates data for comparison expression evaluation. */
 105 static void *
 106 init_data_compare(int npar, gmx_ana_selparam_t *param);
 107 /** Initializes data for comparison expression evaluation. */
 108 static int
 109 init_compare(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
 110 /** Frees the memory allocated for comparison expression evaluation. */
 111 static void
 112 free_data_compare(void *data);
 113 /** Evaluates comparison expressions. */
 114 static int
 115 evaluate_compare(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 116                  gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
 117
 118 /** Parameters for comparison expression evaluation. */
 119 static gmx_ana_selparam_t smparams_compare[] = {
 120     {"int1",  {INT_VALUE,  -1, {NULL}}, NULL,
 121      SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
 122     {"real1", {REAL_VALUE, -1, {NULL}}, NULL,
 123      SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
 124     {"op",    {STR_VALUE,   1, {NULL}}, NULL, 0},
 125     {"int2",  {INT_VALUE,  -1, {NULL}}, NULL,
 126      SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
 127     {"real2", {REAL_VALUE, -1, {NULL}}, NULL,
 128      SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
 129 };
 130
 131 /** \internal Selection method data for comparison expression evaluation. */
 132 gmx_ana_selmethod_t sm_compare = {
 133     "cmp", GROUP_VALUE, SMETH_SINGLEVAL,
 134     asize(smparams_compare), smparams_compare,
 135     &init_data_compare,
 136     NULL,
 137     &init_compare,
 138     NULL,
 139     &free_data_compare,
 140     NULL,
 141     &evaluate_compare,
 142     NULL,
 143     {NULL, 0, NULL},
 144 };
 145
 146 /*! \brief
 147  * Returns a \c e_comparison_t value corresponding to an operator.
 148  *
 149  * \param[in] str  String to process.
 150  * \returns   The comparison type corresponding to the first one or two
 151  *   characters of \p str.
 152  *
 153  * \p str can contain any number of characters; only the first two
 154  * are used.
 155  * If the beginning of \p str does not match any of the recognized types,
 156  * \ref CMP_INVALID is returned.
 157  */
 158 static e_comparison_t
 159 comparison_type(char *str)
 160 {
 161     switch (str[0])
 162     {
 163         case '<': return (str[1] == '=') ? CMP_LEQ   : CMP_LESS;
 164         case '>': return (str[1] == '=') ? CMP_GEQ   : CMP_GTR;
 165         case '=': return (str[1] == '=') ? CMP_EQUAL : CMP_INVALID;
 166         case '!': return (str[1] == '=') ? CMP_NEQ   : CMP_INVALID;
 167     }
 168     return CMP_INVALID;
 169 }
 170
 171 /*! \brief
 172  * Returns a string corresponding to a \c e_comparison_t value.
 173  *
 174  * \param[in] cmpt  Comparison type to convert.
 175  * \returns   Pointer to a string that corresponds to \p cmpt.
 176  *
 177  * The return value points to a string constant and should not be \p free'd.
 178  *
 179  * The function returns NULL if \p cmpt is not one of the valid values.
 180  */
 181 static const char *
 182 comparison_type_str(e_comparison_t cmpt)
 183 {
 184     switch (cmpt)
 185     {
 186         case CMP_INVALID: return "INVALID"; break;
 187         case CMP_LESS:    return "<";  break;
 188         case CMP_LEQ:     return "<="; break;
 189         case CMP_GTR:     return ">";  break;
 190         case CMP_GEQ:     return ">="; break;
 191         case CMP_EQUAL:   return "=="; break;
 192         case CMP_NEQ:     return "!="; break;
 193     }
 194     return NULL;
 195 }
 196
 197 /*!
 198  * \param[in] fp    File to receive the output.
 199  * \param[in] data  Should point to a \c t_methoddata_compare.
 200  */
 201 void
 202 _gmx_selelem_print_compare_info(FILE *fp, void *data)
 203 {
 204     t_methoddata_compare *d = (t_methoddata_compare *)data;
 205
 206     fprintf(fp, " \"");
 207     /* Print the left value */
 208     if ((d->left.flags & CMP_SINGLEVAL) && !(d->left.flags & CMP_DYNAMICVAL))
 209     {
 210         if (d->left.flags & CMP_REALVAL)
 211         {
 212             fprintf(fp, "%f ", d->left.r[0]);
 213         }
 214         else
 215         {
 216             fprintf(fp, "%d ", d->left.i[0]);
 217         }
 218     }
 219     /* Print the operator */
 220     if (d->cmpt != CMP_INVALID)
 221     {
 222         fprintf(fp, "%s", comparison_type_str(d->cmpt));
 223     }
 224     else
 225     {
 226         fprintf(fp, "%s", d->cmpop);
 227     }
 228     /* Print the right value */
 229     if ((d->right.flags & CMP_SINGLEVAL) && !(d->right.flags & CMP_DYNAMICVAL))
 230     {
 231         if (d->right.flags & CMP_REALVAL)
 232         {
 233             fprintf(fp, " %f", d->right.r[0]);
 234         }
 235         else
 236         {
 237             fprintf(fp, " %d", d->right.i[0]);
 238         }
 239     }
 240     fprintf(fp, "\"");
 241 }
 242
 243 /*!
 244  * \param[in]     npar  Not used (should be 5).
 245  * \param[in,out] param Method parameters (should point to a copy of
 246  *   \ref smparams_compare).
 247  * \returns       Pointer to the allocated data (\c t_methoddata_compare).
 248  *
 249  * Allocates memory for a \c t_methoddata_compare structure.
 250  */
 251 static void *
 252 init_data_compare(int npar, gmx_ana_selparam_t *param)
 253 {
 254     t_methoddata_compare *data;
 255
 256     snew(data, 1);
 257     param[2].val.u.s = &data->cmpop;
 258     return data;
 259 }
 260
 261 /* \brief
 262  * Reverses a comparison operator.
 263  *
 264  * \param[in] type  Comparison operator to reverse.
 265  * \returns   The correct comparison operator that equals \p type when the
 266  *   left and right sides are interchanged.
 267  */
 268 static e_comparison_t
 269 reverse_comparison_type(e_comparison_t type)
 270 {
 271     switch (type)
 272     {
 273         case CMP_LESS: return CMP_GTR;
 274         case CMP_LEQ:  return CMP_GEQ;
 275         case CMP_GTR:  return CMP_LESS;
 276         case CMP_GEQ:  return CMP_LEQ;
 277         default:       break;
 278     }
 279     return type;
 280 }
 281
 282 /*! \brief
 283  * Initializes the value storage for comparison expression.
 284  *
 285  * \param[out] val   Value structure to initialize.
 286  * \param[in]  param Parameters to use for initialization.
 287  * \returns    The number of values provided for the value, 0 on error.
 288  */
 289 static int
 290 init_comparison_value(t_compare_value *val, gmx_ana_selparam_t param[2])
 291 {
 292     int  n;
 293
 294     val->flags = 0;
 295     if (param[0].flags & SPAR_SET)
 296     {
 297         val->flags |=  (param[0].flags & SPAR_DYNAMIC) ? CMP_DYNAMICVAL : 0;
 298         val->flags |= !(param[0].flags & SPAR_ATOMVAL) ? CMP_SINGLEVAL  : 0;
 299         n           = param[0].val.nr;
 300         val->i      = param[0].val.u.i;
 301     }
 302     else if (param[1].flags & SPAR_SET)
 303     {
 304         val->flags |=  (param[1].flags & SPAR_DYNAMIC) ? CMP_DYNAMICVAL : 0;
 305         val->flags |= !(param[1].flags & SPAR_ATOMVAL) ? CMP_SINGLEVAL  : 0;
 306         val->flags |= CMP_REALVAL;
 307         n           = param[1].val.nr;
 308         val->r      = param[1].val.u.r;
 309     }
 310     else
 311     {
 312         n           = 0;
 313         val->i      = NULL;
 314         val->r      = NULL;
 315     }
 316     return n;
 317 }
 318
 319 /* \brief
 320  * Converts an integer value to floating point.
 321  *
 322  * \param[in]     n   Number of values in the \p val->u array.
 323  * \param[in,out] val Value to convert.
 324  */
 325 static void
 326 convert_int_real(int n, t_compare_value *val)
 327 {
 328     int   i;
 329     real *rv;
 330
 331     snew(rv, n);
 332     for (i = 0; i < n; ++i)
 333     {
 334         rv[i] = (real)val->i[i];
 335     }
 336     /* Free the previous value if one is present. */
 337     sfree(val->r);
 338     val->r      = rv;
 339     val->flags |= CMP_REALVAL | CMP_ALLOCREAL;
 340 }
 341
 342 /* \brief
 343  * Converts a floating point value to integer.
 344  *
 345  * \param[in]     n      Number of values in the \p val->u array.
 346  * \param[in,out] val    Value to convert.
 347  * \param[in]     cmpt   Comparison operator type.
 348  * \param[in]     bRight TRUE if \p val appears on the right hand size of
 349  *   \p cmpt.
 350  * \returns       0 on success, EINVAL on error.
 351  *
 352  * The values are rounded such that the same comparison operator can be used.
 353  */
 354 static int
 355 convert_real_int(int n, t_compare_value *val, e_comparison_t cmpt, gmx_bool bRight)
 356 {
 357     int   i;
 358     int  *iv;
 359
 360     if (!bRight)
 361     {
 362         cmpt = reverse_comparison_type(cmpt);
 363     }
 364     snew(iv, n);
 365     /* Round according to the comparison type */
 366     for (i = 0; i < n; ++i)
 367     {
 368         switch (cmpt)
 369         {
 370             case CMP_LESS:
 371             case CMP_GEQ:
 372                 iv[i] = (int)ceil(val->r[i]);
 373                 break;
 374             case CMP_GTR:
 375             case CMP_LEQ:
 376                 iv[i] = (int)floor(val->r[i]);
 377                 break;
 378             case CMP_EQUAL:
 379             case CMP_NEQ:
 380                 fprintf(stderr, "comparing equality an integer expression and a real value\n");
 381                 sfree(iv);
 382                 return EINVAL;
 383             case CMP_INVALID: /* Should not be reached */
 384                 gmx_bug("internal error");
 385                 sfree(iv);
 386                 return EINVAL;
 387         }
 388     }
 389     /* Free the previous value if one is present. */
 390     sfree(val->i);
 391     val->i      = iv;
 392     val->flags &= ~CMP_REALVAL;
 393     val->flags |= CMP_ALLOCINT;
 394     return 0;
 395 }
 396
 397 /*!
 398  * \param[in] top   Not used.
 399  * \param[in] npar  Not used (should be 5).
 400  * \param[in] param Method parameters (should point to \ref smparams_compare).
 401  * \param[in] data  Should point to a \c t_methoddata_compare.
 402  * \returns   0 if the input data is valid, -1 on error.
 403  */
 404 static int
 405 init_compare(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data)
 406 {
 407     t_methoddata_compare *d = (t_methoddata_compare *)data;
 408     int                   n1, n2;
 409
 410     /* Store the values */
 411     n1 = init_comparison_value(&d->left, &param[0]);
 412     n2 = init_comparison_value(&d->right, &param[3]);
 413     if (n1 == 0 || n2 == 0)
 414     {
 415         gmx_bug("one of the values for comparison missing");
 416         return -1;
 417     }
 418     /* Store the comparison type */
 419     d->cmpt = comparison_type(d->cmpop);
 420     if (d->cmpt == CMP_INVALID)
 421     {
 422         gmx_bug("invalid comparison type");
 423         return -1;
 424     }
 425     /* Convert the values to the same type */
 426     if ((d->left.flags & CMP_REALVAL) && !(d->right.flags & CMP_REALVAL))
 427     {
 428         if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
 429         {
 430             /* Nothing can be done */
 431         }
 432         else if (!(d->right.flags & CMP_DYNAMICVAL))
 433         {
 434             convert_int_real(n2, &d->right);
 435         }
 436         else /* d->left is static */
 437         {
 438             if (convert_real_int(n1, &d->left, d->cmpt, FALSE))
 439             {
 440                 return -1;
 441             }
 442         }
 443     }
 444     else if (!(d->left.flags & CMP_REALVAL) && (d->right.flags & CMP_REALVAL))
 445     {
 446         if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
 447         {
 448             /* Reverse the sides to place the integer on the right */
 449             int    flags;
 450             d->left.r  = d->right.r;
 451             d->right.r = NULL;
 452             d->right.i = d->left.i;
 453             d->left.i  = NULL;
 454             flags          = d->left.flags;
 455             d->left.flags  = d->right.flags;
 456             d->right.flags = flags;
 457             d->cmpt = reverse_comparison_type(d->cmpt);
 458         }
 459         else if (!(d->left.flags & CMP_DYNAMICVAL))
 460         {
 461             convert_int_real(n1, &d->left);
 462         }
 463         else /* d->right is static */
 464         {
 465             if (convert_real_int(n2, &d->right, d->cmpt, TRUE))
 466             {
 467                 return -1;
 468             }
 469         }
 470     }
 471     return 0;
 472 }
 473
 474 /*!
 475  * \param data Data to free (should point to a \c t_methoddata_compare).
 476  *
 477  * Frees the memory allocated for \c t_methoddata_compare.
 478  */
 479 static void
 480 free_data_compare(void *data)
 481 {
 482     t_methoddata_compare *d = (t_methoddata_compare *)data;
 483
 484     sfree(d->cmpop);
 485     if (d->left.flags & CMP_ALLOCINT)
 486     {
 487         sfree(d->left.i);
 488     }
 489     if (d->left.flags & CMP_ALLOCREAL)
 490     {
 491         sfree(d->left.r);
 492     }
 493     if (d->right.flags & CMP_ALLOCINT)
 494     {
 495         sfree(d->right.i);
 496     }
 497     if (d->right.flags & CMP_ALLOCREAL)
 498     {
 499         sfree(d->right.r);
 500     }
 501 }
 502
 503 /*!
 504  * \param[in]  top   Not used.
 505  * \param[in]  fr    Not used.
 506  * \param[in]  pbc   Not used.
 507  * \param[in]  g     Evaluation index group.
 508  * \param[out] out   Output data structure (\p out->u.g is used).
 509  * \param[in]  data  Should point to a \c t_methoddata_compare.
 510  * \returns    0 for success.
 511  */
 512 static int
 513 evaluate_compare_int(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 514                      gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
 515 {
 516     t_methoddata_compare *d = (t_methoddata_compare *)data;
 517     int                   i, i1, i2, ig;
 518     int                   a, b;
 519     gmx_bool                  bAccept;
 520
 521     for (i = i1 = i2 = ig = 0; i < g->isize; ++i)
 522     {
 523         a = d->left.i[i1];
 524         b = d->right.i[i2];
 525         bAccept = FALSE;
 526         switch (d->cmpt)
 527         {
 528             case CMP_INVALID: break;
 529             case CMP_LESS:    bAccept = a <  b; break;
 530             case CMP_LEQ:     bAccept = a <= b; break;
 531             case CMP_GTR:     bAccept = a >  b; break;
 532             case CMP_GEQ:     bAccept = a >= b; break;
 533             case CMP_EQUAL:   bAccept = a == b; break;
 534             case CMP_NEQ:     bAccept = a != b; break;
 535         }
 536         if (bAccept)
 537         {
 538             out->u.g->index[ig++] = g->index[i];
 539         }
 540         if (!(d->left.flags & CMP_SINGLEVAL))
 541         {
 542             ++i1;
 543         }
 544         if (!(d->right.flags & CMP_SINGLEVAL))
 545         {
 546             ++i2;
 547         }
 548     }
 549     out->u.g->isize = ig;
 550     return 0;
 551 }
 552
 553 /*!
 554  * \param[in]  top   Not used.
 555  * \param[in]  fr    Not used.
 556  * \param[in]  pbc   Not used.
 557  * \param[in]  g     Evaluation index group.
 558  * \param[out] out   Output data structure (\p out->u.g is used).
 559  * \param[in]  data  Should point to a \c t_methoddata_compare.
 560  * \returns    0 for success.
 561  */
 562 static int
 563 evaluate_compare_real(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 564                       gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
 565 {
 566     t_methoddata_compare *d = (t_methoddata_compare *)data;
 567     int                   i, i1, i2, ig;
 568     real                  a, b;
 569     gmx_bool                  bAccept;
 570
 571     for (i = i1 = i2 = ig = 0; i < g->isize; ++i)
 572     {
 573         a = d->left.r[i1];
 574         b = (d->right.flags & CMP_REALVAL) ? d->right.r[i2] : d->right.i[i2];
 575         bAccept = FALSE;
 576         switch (d->cmpt)
 577         {
 578             case CMP_INVALID: break;
 579             case CMP_LESS:    bAccept = a <  b; break;
 580             case CMP_LEQ:     bAccept = a <= b; break;
 581             case CMP_GTR:     bAccept = a >  b; break;
 582             case CMP_GEQ:     bAccept = a >= b; break;
 583             case CMP_EQUAL:   bAccept =  gmx_within_tol(a, b, GMX_REAL_EPS); break;
 584             case CMP_NEQ:     bAccept = !gmx_within_tol(a, b, GMX_REAL_EPS); break;
 585         }
 586         if (bAccept)
 587         {
 588             out->u.g->index[ig++] = g->index[i];
 589         }
 590         if (!(d->left.flags & CMP_SINGLEVAL))
 591         {
 592             ++i1;
 593         }
 594         if (!(d->right.flags & CMP_SINGLEVAL))
 595         {
 596             ++i2;
 597         }
 598     }
 599     out->u.g->isize = ig;
 600     return 0;
 601 }
 602
 603 /*!
 604  * \param[in]  top   Not used.
 605  * \param[in]  fr    Not used.
 606  * \param[in]  pbc   Not used.
 607  * \param[in]  g     Evaluation index group.
 608  * \param[out] out   Output data structure (\p out->u.g is used).
 609  * \param[in]  data  Should point to a \c t_methoddata_compare.
 610  * \returns    0 for success.
 611  */
 612 static int
 613 evaluate_compare(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 614                  gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
 615 {
 616     t_methoddata_compare *d = (t_methoddata_compare *)data;
 617
 618     if (!((d->left.flags | d->right.flags) & CMP_REALVAL))
 619     {
 620         return evaluate_compare_int(top, fr, pbc, g, out, data);
 621     }
 622     else
 623     {
 624         return evaluate_compare_real(top, fr, pbc, g, out, data);
 625     }
 626     return 0;
 627 }