Add unit tests for index groups in selections.
[alexxy/gromacs.git] / src / gromacs / selection / selelem.cpp
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2009,2010,2011,2012,2013, by the GROMACS development team, led by
5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
6  * others, as listed in the AUTHORS file in the top-level source
7  * directory and at http://www.gromacs.org.
8  *
9  * GROMACS is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public License
11  * as published by the Free Software Foundation; either version 2.1
12  * of the License, or (at your option) any later version.
13  *
14  * GROMACS is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17  * Lesser General Public License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with GROMACS; if not, see
21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
23  *
24  * If you want to redistribute modifications to GROMACS, please
25  * consider that scientific software is very special. Version
26  * control is crucial - bugs must be traceable. We will be happy to
27  * consider code for inclusion in the official distribution, but
28  * derived work must not be called official GROMACS. Details are found
29  * in the README & COPYING files - if they are missing, get the
30  * official version at http://www.gromacs.org.
31  *
32  * To help us fund GROMACS development, we humbly ask that you cite
33  * the research papers on the package. Check out http://www.gromacs.org.
34  */
35 /*! \internal \file
36  * \brief
37  * Implements functions in selelem.h.
38  *
39  * \author Teemu Murtola <teemu.murtola@gmail.com>
40  * \ingroup module_selection
41  */
42 #include <cstring>
43
44 #include "gromacs/legacyheaders/smalloc.h"
45
46 #include "gromacs/selection/indexutil.h"
47 #include "gromacs/selection/poscalc.h"
48 #include "gromacs/selection/position.h"
49 #include "gromacs/selection/selmethod.h"
50 #include "gromacs/utility/exceptions.h"
51 #include "gromacs/utility/gmxassert.h"
52
53 #include "keywords.h"
54 #include "mempool.h"
55 #include "selelem.h"
56
57 /*!
58  * \param[in] sel Selection for which the string is requested
59  * \returns   Pointer to a string that corresponds to \p sel->type.
60  *
61  * The return value points to a string constant and should not be \p free'd.
62  *
63  * The function returns NULL if \p sel->type is not one of the valid values.
64  */
65 const char *
66 _gmx_selelem_type_str(const gmx::SelectionTreeElement &sel)
67 {
68     switch (sel.type)
69     {
70         case SEL_CONST:      return "CONST";
71         case SEL_EXPRESSION: return "EXPR";
72         case SEL_BOOLEAN:    return "BOOL";
73         case SEL_ARITHMETIC: return "ARITH";
74         case SEL_ROOT:       return "ROOT";
75         case SEL_SUBEXPR:    return "SUBEXPR";
76         case SEL_SUBEXPRREF: return "REF";
77         case SEL_GROUPREF:   return "GROUPREF";
78         case SEL_MODIFIER:   return "MODIFIER";
79     }
80     return NULL;
81 }
82
83 /*!
84  * \param[in] val Value structore for which the string is requested.
85  * \returns   Pointer to a string that corresponds to \p val->type,
86  *   NULL if the type value is invalid.
87  *
88  * The return value points to a string constant and should not be \p free'd.
89  */
90 const char *
91 _gmx_sel_value_type_str(const gmx_ana_selvalue_t *val)
92 {
93     switch (val->type)
94     {
95         case NO_VALUE:       return "NONE";
96         case INT_VALUE:      return "INT";
97         case REAL_VALUE:     return "REAL";
98         case STR_VALUE:      return "STR";
99         case POS_VALUE:      return "VEC";
100         case GROUP_VALUE:    return "GROUP";
101     }
102     return NULL;
103 }
104
105 /*! \copydoc _gmx_selelem_type_str() */
106 const char *
107 _gmx_selelem_boolean_type_str(const gmx::SelectionTreeElement &sel)
108 {
109     switch (sel.u.boolt)
110     {
111         case BOOL_NOT:  return "NOT"; break;
112         case BOOL_AND:  return "AND"; break;
113         case BOOL_OR:   return "OR";  break;
114         case BOOL_XOR:  return "XOR"; break;
115     }
116     return NULL;
117 }
118
119
120 namespace gmx
121 {
122
123 SelectionTreeElement::SelectionTreeElement(e_selelem_t type)
124 {
125     this->type       = type;
126     this->flags      = (type != SEL_ROOT) ? SEL_ALLOCVAL : 0;
127     if (type == SEL_BOOLEAN)
128     {
129         this->v.type = GROUP_VALUE;
130         this->flags |= SEL_ALLOCDATA;
131     }
132     else
133     {
134         this->v.type = NO_VALUE;
135     }
136     _gmx_selvalue_clear(&this->v);
137     std::memset(&this->u, 0, sizeof(this->u));
138     this->evaluate   = NULL;
139     this->mempool    = NULL;
140     this->cdata      = NULL;
141 }
142
143 SelectionTreeElement::~SelectionTreeElement()
144 {
145     /* Free the children.
146      * Must be done before freeing other data, because the children may hold
147      * references to data in this element. */
148     child.reset();
149
150     freeValues();
151     freeExpressionData();
152     freeCompilerData();
153 }
154
155 void SelectionTreeElement::freeValues()
156 {
157     mempoolRelease();
158     if ((flags & SEL_ALLOCDATA) && v.u.ptr)
159     {
160         /* The number of position/group structures is constant, so the
161          * backup of using sel->v.nr should work for them.
162          * For strings, we report an error if we don't know the allocation
163          * size here. */
164         int n = (v.nalloc > 0) ? v.nalloc : v.nr;
165         switch (v.type)
166         {
167             case STR_VALUE:
168                 GMX_RELEASE_ASSERT(v.nalloc != 0,
169                                    "SEL_ALLOCDATA should only be set for allocated "
170                                    "STR_VALUE values");
171                 for (int i = 0; i < n; ++i)
172                 {
173                     sfree(v.u.s[i]);
174                 }
175                 break;
176             case POS_VALUE:
177                 for (int i = 0; i < n; ++i)
178                 {
179                     gmx_ana_pos_deinit(&v.u.p[i]);
180                 }
181                 break;
182             case GROUP_VALUE:
183                 for (int i = 0; i < n; ++i)
184                 {
185                     gmx_ana_index_deinit(&v.u.g[i]);
186                 }
187                 break;
188             default: /* No special handling for other types */
189                 break;
190         }
191     }
192     if (flags & SEL_ALLOCVAL)
193     {
194         sfree(v.u.ptr);
195     }
196     _gmx_selvalue_setstore(&v, NULL);
197     if (type == SEL_SUBEXPRREF && u.param)
198     {
199         u.param->val.u.ptr = NULL;
200     }
201 }
202
203 void
204 SelectionTreeElement::freeExpressionData()
205 {
206     if (type == SEL_EXPRESSION || type == SEL_MODIFIER)
207     {
208         _gmx_selelem_free_method(u.expr.method, u.expr.mdata);
209         u.expr.mdata  = NULL;
210         u.expr.method = NULL;
211         /* Free position data */
212         if (u.expr.pos)
213         {
214             gmx_ana_pos_free(u.expr.pos);
215             u.expr.pos = NULL;
216         }
217         /* Free position calculation data */
218         if (u.expr.pc)
219         {
220             gmx_ana_poscalc_free(u.expr.pc);
221             u.expr.pc = NULL;
222         }
223     }
224     if (type == SEL_ARITHMETIC)
225     {
226         sfree(u.arith.opstr);
227         u.arith.opstr = NULL;
228     }
229     if (type == SEL_SUBEXPR || type == SEL_ROOT
230         || (type == SEL_CONST && v.type == GROUP_VALUE))
231     {
232         gmx_ana_index_deinit(&u.cgrp);
233     }
234     if (type == SEL_GROUPREF)
235     {
236         sfree(u.gref.name);
237     }
238 }
239
240 void SelectionTreeElement::mempoolReserve(int count)
241 {
242     if (!mempool)
243     {
244         return;
245     }
246     switch (v.type)
247     {
248         case INT_VALUE:
249             v.u.i = static_cast<int *>(
250                     _gmx_sel_mempool_alloc(mempool, sizeof(*v.u.i)*count));
251             break;
252
253         case REAL_VALUE:
254             v.u.r = static_cast<real *>(
255                     _gmx_sel_mempool_alloc(mempool, sizeof(*v.u.r)*count));
256             break;
257
258         case GROUP_VALUE:
259             _gmx_sel_mempool_alloc_group(mempool, v.u.g, count);
260             break;
261
262         default:
263             GMX_THROW(gmx::InternalError("Memory pooling not implemented for requested type"));
264     }
265 }
266
267 void SelectionTreeElement::mempoolRelease()
268 {
269     if (!mempool)
270     {
271         return;
272     }
273     switch (v.type)
274     {
275         case INT_VALUE:
276         case REAL_VALUE:
277             _gmx_sel_mempool_free(mempool, v.u.ptr);
278             _gmx_selvalue_setstore(&v, NULL);
279             break;
280
281         case GROUP_VALUE:
282             if (v.u.g)
283             {
284                 _gmx_sel_mempool_free_group(mempool, v.u.g);
285             }
286             break;
287
288         default:
289             GMX_THROW(gmx::InternalError("Memory pooling not implemented for requested type"));
290     }
291 }
292
293 void SelectionTreeElement::fillNameIfMissing(const char *selectionText)
294 {
295     GMX_RELEASE_ASSERT(type == SEL_ROOT,
296                        "Should not be called for non-root elements");
297     if (name().empty())
298     {
299         // Check whether the actual selection given was from an external group,
300         // and if so, use the name of the external group.
301         SelectionTreeElementPointer child = this->child;
302         while (child->type == SEL_MODIFIER)
303         {
304             if (!child->child || child->child->type != SEL_SUBEXPRREF
305                 || !child->child->child)
306             {
307                 break;
308             }
309             child = child->child->child;
310         }
311         if (child->type == SEL_EXPRESSION
312             && child->child && child->child->type == SEL_SUBEXPRREF
313             && child->child->child)
314         {
315             if (child->child->child->type == SEL_CONST
316                 && child->child->child->v.type == GROUP_VALUE)
317             {
318                 setName(child->child->child->name());
319                 return;
320             }
321             // If the group reference is still unresolved, leave the name empty
322             // and fill it later.
323             if (child->child->child->type == SEL_GROUPREF)
324             {
325                 return;
326             }
327         }
328         // If there still is no name, use the selection string.
329         setName(selectionText);
330     }
331 }
332
333 } // namespace gmx
334
335 /*!
336  * \param[in,out] sel   Selection element to set the type for.
337  * \param[in]     vtype Value type for the selection element.
338  *
339  * If the new type is \ref GROUP_VALUE or \ref POS_VALUE, the
340  * \ref SEL_ALLOCDATA flag is also set.
341  *
342  * This function should only be called at most once for each element,
343  * preferably right after calling _gmx_selelem_create().
344  */
345 void
346 _gmx_selelem_set_vtype(const gmx::SelectionTreeElementPointer &sel,
347                        e_selvalue_t                            vtype)
348 {
349     GMX_RELEASE_ASSERT(sel->type != SEL_BOOLEAN || vtype == GROUP_VALUE,
350                        "Boolean elements must have a group value");
351     GMX_RELEASE_ASSERT(sel->v.type == NO_VALUE || vtype == sel->v.type,
352                        "_gmx_selelem_set_vtype() called more than once");
353     sel->v.type = vtype;
354     if (vtype == GROUP_VALUE || vtype == POS_VALUE)
355     {
356         sel->flags |= SEL_ALLOCDATA;
357     }
358 }
359
360 /*!
361  * \param[in] method Method to free.
362  * \param[in] mdata  Method data to free.
363  */
364 void
365 _gmx_selelem_free_method(gmx_ana_selmethod_t *method, void *mdata)
366 {
367     sel_freefunc free_func = NULL;
368
369     /* Save the pointer to the free function. */
370     if (method && method->free)
371     {
372         free_func = method->free;
373     }
374
375     /* Free the method itself.
376      * Has to be done before freeing the method data, because parameter
377      * values are typically stored in the method data, and here we may
378      * access them. */
379     if (method)
380     {
381         int  i, j;
382
383         /* Free the memory allocated for the parameters that are not managed
384          * by the selection method itself. */
385         for (i = 0; i < method->nparams; ++i)
386         {
387             gmx_ana_selparam_t *param = &method->param[i];
388
389             if (param->val.u.ptr)
390             {
391                 if (param->val.type == GROUP_VALUE)
392                 {
393                     for (j = 0; j < param->val.nr; ++j)
394                     {
395                         gmx_ana_index_deinit(&param->val.u.g[j]);
396                     }
397                 }
398                 else if (param->val.type == POS_VALUE)
399                 {
400                     for (j = 0; j < param->val.nr; ++j)
401                     {
402                         gmx_ana_pos_deinit(&param->val.u.p[j]);
403                     }
404                 }
405
406                 if (param->val.nalloc > 0)
407                 {
408                     sfree(param->val.u.ptr);
409                 }
410             }
411         }
412         sfree(method->param);
413         sfree(method);
414     }
415     /* Free method data. */
416     if (mdata)
417     {
418         if (free_func)
419         {
420             free_func(mdata);
421         }
422         else
423         {
424             sfree(mdata);
425         }
426     }
427 }
428
429 /*!
430  * \param[in] fp      File handle to receive the output.
431  * \param[in] sel     Root of the selection subtree to print.
432  * \param[in] bValues If true, the evaluated values of selection elements
433  *   are printed as well.
434  * \param[in] level   Indentation level, starting from zero.
435  */
436 void
437 _gmx_selelem_print_tree(FILE *fp, const gmx::SelectionTreeElement &sel,
438                         bool bValues, int level)
439 {
440     int          i;
441
442     fprintf(fp, "%*c %s %s", level*2+1, '*',
443             _gmx_selelem_type_str(sel), _gmx_sel_value_type_str(&sel.v));
444     if (!sel.name().empty())
445     {
446         fprintf(fp, " \"%s\"", sel.name().c_str());
447     }
448     fprintf(fp, " flg=");
449     if (sel.flags & SEL_FLAGSSET)
450     {
451         fprintf(fp, "s");
452     }
453     if (sel.flags & SEL_SINGLEVAL)
454     {
455         fprintf(fp, "S");
456     }
457     if (sel.flags & SEL_ATOMVAL)
458     {
459         fprintf(fp, "A");
460     }
461     if (sel.flags & SEL_VARNUMVAL)
462     {
463         fprintf(fp, "V");
464     }
465     if (sel.flags & SEL_DYNAMIC)
466     {
467         fprintf(fp, "D");
468     }
469     if (!(sel.flags & SEL_VALFLAGMASK))
470     {
471         fprintf(fp, "0");
472     }
473     if (sel.flags & SEL_ALLOCVAL)
474     {
475         fprintf(fp, "Av");
476     }
477     if (sel.flags & SEL_ALLOCDATA)
478     {
479         fprintf(fp, "Ad");
480     }
481     if (sel.mempool)
482     {
483         fprintf(fp, "P");
484     }
485     if (sel.type == SEL_CONST)
486     {
487         if (sel.v.type == INT_VALUE)
488         {
489             fprintf(fp, " %d", sel.v.u.i[0]);
490         }
491         else if (sel.v.type == REAL_VALUE)
492         {
493             fprintf(fp, " %f", sel.v.u.r[0]);
494         }
495         else if (sel.v.type == GROUP_VALUE)
496         {
497             const gmx_ana_index_t *g = sel.v.u.g;
498             if (!g || g->isize == 0)
499             {
500                 g = &sel.u.cgrp;
501             }
502             fprintf(fp, " (%d atoms)", g->isize);
503         }
504     }
505     else if (sel.type == SEL_BOOLEAN)
506     {
507         fprintf(fp, " %s", _gmx_selelem_boolean_type_str(sel));
508     }
509     else if (sel.type == SEL_EXPRESSION
510              && sel.u.expr.method->name == sm_compare.name)
511     {
512         _gmx_selelem_print_compare_info(fp, sel.u.expr.mdata);
513     }
514     if (sel.evaluate)
515     {
516         fprintf(fp, " eval=");
517         _gmx_sel_print_evalfunc_name(fp, sel.evaluate);
518     }
519     if (!(sel.flags & SEL_ALLOCVAL))
520     {
521         fprintf(fp, " (ext)");
522     }
523     fprintf(fp, "\n");
524
525     if ((sel.type == SEL_CONST && sel.v.type == GROUP_VALUE) || sel.type == SEL_ROOT)
526     {
527         const gmx_ana_index_t *g = sel.v.u.g;
528         if (!g || g->isize == 0 || sel.evaluate != NULL)
529         {
530             g = &sel.u.cgrp;
531         }
532         if (g->isize < 0)
533         {
534             fprintf(fp, "%*c group: (null)\n", level*2+1, ' ');
535         }
536         else if (g->isize > 0)
537         {
538             fprintf(fp, "%*c group:", level*2+1, ' ');
539             if (g->isize <= 20)
540             {
541                 for (i = 0; i < g->isize; ++i)
542                 {
543                     fprintf(fp, " %d", g->index[i] + 1);
544                 }
545             }
546             else
547             {
548                 fprintf(fp, " %d atoms", g->isize);
549             }
550             fprintf(fp, "\n");
551         }
552     }
553     else if (sel.type == SEL_EXPRESSION)
554     {
555         if (sel.u.expr.pc)
556         {
557             fprintf(fp, "%*c COM", level*2+3, '*');
558             fprintf(fp, "\n");
559         }
560     }
561
562     if (sel.cdata)
563     {
564         _gmx_selelem_print_compiler_info(fp, sel, level);
565     }
566
567     if (bValues && sel.type != SEL_CONST && sel.type != SEL_ROOT && sel.v.u.ptr)
568     {
569         fprintf(fp, "%*c value: ", level*2+1, ' ');
570         switch (sel.v.type)
571         {
572             case POS_VALUE:
573                 /* In normal use, the pointer should never be NULL, but it's
574                  * useful to have the check for debugging to avoid accidental
575                  * segfaults when printing the selection tree. */
576                 if (sel.v.u.p->x)
577                 {
578                     fprintf(fp, "(%f, %f, %f)",
579                             sel.v.u.p->x[0][XX], sel.v.u.p->x[0][YY],
580                             sel.v.u.p->x[0][ZZ]);
581                 }
582                 else
583                 {
584                     fprintf(fp, "(null)");
585                 }
586                 break;
587             case GROUP_VALUE:
588                 fprintf(fp, "%d atoms", sel.v.u.g->isize);
589                 if (sel.v.u.g->isize < 20)
590                 {
591                     if (sel.v.u.g->isize > 0)
592                     {
593                         fprintf(fp, ":");
594                     }
595                     for (i = 0; i < sel.v.u.g->isize; ++i)
596                     {
597                         fprintf(fp, " %d", sel.v.u.g->index[i] + 1);
598                     }
599                 }
600                 break;
601             default:
602                 fprintf(fp, "???");
603                 break;
604         }
605         fprintf(fp, "\n");
606     }
607
608     /* Print the subexpressions with one more level of indentation */
609     gmx::SelectionTreeElementPointer child = sel.child;
610     while (child)
611     {
612         if (!(sel.type == SEL_SUBEXPRREF && child->type == SEL_SUBEXPR))
613         {
614             _gmx_selelem_print_tree(fp, *child, bValues, level+1);
615         }
616         child = child->next;
617     }
618 }
619
620 /*!
621  * \param[in] root Root of the subtree to query.
622  * \returns true if \p root or any any of its elements require topology
623  *   information, false otherwise.
624  */
625 bool
626 _gmx_selelem_requires_top(const gmx::SelectionTreeElement &root)
627 {
628     if (root.type == SEL_EXPRESSION || root.type == SEL_MODIFIER)
629     {
630         if (root.u.expr.method && (root.u.expr.method->flags & SMETH_REQTOP))
631         {
632             return true;
633         }
634         if (root.u.expr.pc && gmx_ana_poscalc_requires_top(root.u.expr.pc))
635         {
636             return true;
637         }
638     }
639     gmx::SelectionTreeElementPointer child = root.child;
640     while (child)
641     {
642         if (_gmx_selelem_requires_top(*child))
643         {
644             return true;
645         }
646         child = child->next;
647     }
648     return false;
649 }