src/gmxlib/selection/sm_compare.c

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