File: base.py

package info (click to toggle)
mdanalysis 2.10.0-1
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 116,696 kB
  • sloc: python: 92,135; ansic: 8,156; makefile: 215; sh: 138
file content (1135 lines) | stat: -rw-r--r-- 40,307 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
#
# MDAnalysis --- https://www.mdanalysis.org
# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors
# (see the file AUTHORS for the full list of names)
#
# Released under the Lesser GNU Public Licence, v2.1 or any higher version
#
# Please cite your use of MDAnalysis in published work:
#
# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,
# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.
# MDAnalysis: A Python package for the rapid analysis of molecular dynamics
# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th
# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.
# doi: 10.25080/majora-629e541a-00e
#
# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
#


"""
Auxiliary Readers --- :mod:`MDAnalysis.auxiliary.base`
======================================================

Base classes for deriving all auxiliary data readers. See the API in :mod:`MDAnalysis.auxiliary.__init__`.

.. autoclass:: AuxStep
   :members:

.. autoclass:: AuxReader
   :members:

.. autoclass:: AuxFileReader
   :members:

"""

import os
import numbers
import math
import warnings
from typing import Union, Optional, Dict
from collections import defaultdict
import pickle

import numpy as np

from ..lib.util import asiterable, anyopen
from .. import units

from . import _AUXREADERS
from .core import auxreader


class _AuxReaderMeta(type):
    # auto register on class creation
    def __init__(cls, name, bases, classdict):
        type.__init__(type, name, bases, classdict)
        try:
            fmt = asiterable(classdict["format"])
        except KeyError:
            pass
        else:
            for f in fmt:
                _AUXREADERS[f] = cls


class AuxStep(object):
    """Base class for auxiliary timesteps.

    Stores the auxiliary data for the current auxiliary step. On creation,
    ``step`` is set to -1.

    .. versionchanged:: 2.4.0
        Added memory_limit parameter to control raising memory usage warnings.

    Parameters
    ----------
    dt : float, optional
        Change in time between auxiliary steps (in ps). If not specified, will
        attempt to determine from auxiliary data; otherwise defaults to 1 ps.
        Ignored if ``constant_dt`` is False.
    initial_time : float, optional
        Time of first auxiliary step (in ps). If not specified, will attempt to
        determine from auxiliary data; otherwise defaults to 0 ps. Ignored if
        ``constant_dt`` is False.
    time_selector: optional
        Key to select 'time' value from the full set of data read for each
        step, if time selection is enabled; type will vary depending on the
        auxiliary data format (see individual AuxReader documentation). If
        ``None`` (default value), time is instead calculated as: ``time = step
        * dt + initial_time``
    data_selector: optional
        Key(s) to select auxiliary data values of interest from the full set of
        data read for each step, if data selection is enabled by the reader;
        type will vary depending on the auxiliary data format (see individual
        AuxReader documentation).
        If ``None`` (default value), the full set of data is returned.
    constant_dt : bool, optional
        (Default: True) Set to False if ``dt`` is not constant
        throughout the auxiliary data set, in which case a valid
        ``time_selector`` must be provided.
    memory_limit : float, optional
        Sets the threshold of memory usage by auxiliary data  (in GB) at which
        to issue a warning. Default: 1 GB.

    Attributes
    ----------
    step : int
        Number of the current auxiliary step (0-based).
    """

    def __init__(
        self,
        dt=1,
        initial_time=0,
        time_selector=None,
        data_selector=None,
        constant_dt=True,
        memory_limit=None,
    ):
        self.step = -1
        self._initial_time = initial_time
        self._dt = dt
        self._constant_dt = constant_dt
        # check for valid values when assigning _time/data_selector will fail
        # here as we don't have and _data yet, so set _time/data_selector_ directly;
        # if invalid, will catch later
        self._time_selector_ = time_selector
        self._data_selector_ = data_selector

    @property
    def time(self):
        """Time in ps of current auxiliary step (as float).

        Read from the set of auxiliary data read each step if time selection
        is enabled and a valid ``time_selector`` is specified; otherwise
        calculated as ``step * dt + initial_time``.
        """
        if self._time_selector is not None:
            return self._select_time(self._time_selector)
        elif self._constant_dt:
            # default to calculting time...
            return self.step * self._dt + self._initial_time
        else:
            raise ValueError(
                "If dt is not constant, must have a valid " "time selector"
            )

    @property
    def data(self):
        """Auxiliary values of interest for the current step (as ndarray).

        Read from the full set of data read for each step if data selection is
        enabled and a valid ``data_selector`` is specified; otherwise
        defaults to the full set of data.
        """
        if self._data_selector is not None:
            return self._select_data(self._data_selector)
        # default to full set of data...
        return self._data

    @property
    def _time_selector(self):
        """'Key' to select time from the full set of data read in each step.

        Will be passed to ``_select_time()``, defined separately for each
        auxiliary format, when returning the time of the current step.
        Format will depend on the auxiliary format. e.g. for the XVGReader,
        this is an index and ``_select_time()`` returns the value in that column
        of the current step data.

        Defaults to 'None' if time selection is not enabled.
        """
        try:
            self._select_time
        except AttributeError:
            warnings.warn(
                "{} does not support time selection. Reverting to "
                "default".format(self.__class__.__name__)
            )
            return None
        return self._time_selector_

    @_time_selector.setter
    def _time_selector(self, new):
        # check we have a select_time method
        try:
            select = self._select_time
        except AttributeError:
            warnings.warn(
                "{} does not support time selection".format(
                    self.__class__.__name__
                )
            )
        else:
            # check *new* is valid before setting; _select_time should raise
            # an error if not
            select(new)
            self._time_selector_ = new

    @property
    def _data_selector(self):
        """'Key' to select values of interest from full set of auxiliary data.
        These are the values that will be stored in ``data`` (and
        ``frame_data`` and ``frame_rep``).

        Will be passed to ``_select_data()``, defined separately for each
        auxiliary format, when returning the data of interest for the current
        step (``data``). Format will depend on the auxiliary format; e.g.
        for the XVGReader, this is a list of indices and `_select_data()` returns
        the value(s) in those columns of the current step data.

        Defaults to 'None' if data selection is not enabled.
        """
        try:
            self._select_data
        except AttributeError:
            warnings.warn(
                "{} does not support data selection. Reverting to "
                "default".format(self.__class__.__name__)
            )
            return None
        return self._data_selector_

    @_data_selector.setter
    def _data_selector(self, new):
        # check we have a select_data method
        try:
            select = self._select_data
        except AttributeError:
            warnings.warn(
                "{} does not support data selection".format(
                    self.__class__.__name__
                )
            )
        else:
            # check *new* is valid before setting; _select_data should raise an
            # error if not
            select(new)
            self._data_selector_ = new

    def _empty_data(self):
        """Create an 'empty' ``data``-like placeholder.

        Returns an ndarray in the format of ``data`` with all values set to
        np.nan; to use at the 'representative value' when no auxiliary steps
        are assigned to a trajectory timestep/within the cutoff.

        Default behaviour here works when ``data`` is a ndarray of floats. May
        need to overwrite in individual format's AuxSteps.

        .. versionchanged:: 2.4.0
            dtype changed to np.float64 from np.float_
        """
        return np.full_like(self.data, np.nan, dtype=np.float64)


class AuxReader(metaclass=_AuxReaderMeta):
    """Base class for auxiliary readers.

    Allows iteration over a set of data from a trajectory, additional
    ('auxiliary') to the regular positions/velocities/etc. This auxiliary
    data may be stored in e.g. an array or a separate file.

    See the :ref:`Auxiliary API` for more on use.

    .. versionchanged:: 2.4.0
        Behaviour of ``cutoff`` changed, default parameter which specifies
        not cutoff is set is now None, not -1.

    .. versionchanged:: 2.4.0
        :class:`AuxReader` instances now have a :meth:`copy` method which
        creates a deep copy of the instance.

    Parameters
    ----------
    auxname : str, optional
        Name for auxiliary data. When added to a trajectory, the representative
        auxiliary value(s) for the timestep may be accessed as ``ts.aux.auxname``
        or ``ts.aux['auxname']``.
    represent_ts_as : str
        Method to use to calculate representative value of auxiliary data for a
        trajectory timestep. See :func:`calc_representative` for valid options.
    cutoff : float, optional
        Auxiliary steps further from the trajectory timestep than *cutoff*
        (in ps) will be ignored when calculating representative values. If None
        (default), all auxiliary steps assigned to that timestep will be used.
    **kwargs
        Options to be passed to :class:`~AuxStep`


    Attributes
    ----------
    auxstep :
       :class:`~AuxStep` object containing data for current step.
    frame_data : dict
        Dictionary containing ``data`` from each auxiliary step assigned to the
        current trajectory timestep, indexed by the difference in time between
        the step and trajectory timestep (i.e. ``auxstep.time - ts.time``; in ps)
    frame_rep : ndarray
        Representative value(s) of auxiliary data for current trajectory timestep.


    Note
    ----
    Auxiliary data are assumed to be time ordered and contain no duplicates.
    """

    _Auxstep = AuxStep
    # update when add new options
    represent_options = ["closest", "average"]

    # list of attributes required to recreate the auxiliary
    required_attrs = [
        "represent_ts_as",
        "cutoff",
        "dt",
        "initial_time",
        "time_selector",
        "data_selector",
        "constant_dt",
        "auxname",
        "format",
        "_auxdata",
    ]

    def __init__(
        self, represent_ts_as="closest", auxname=None, cutoff=None, **kwargs
    ):
        # allow auxname to be optional for when using reader separate from
        # trajectory.
        self.auxname = auxname
        self.represent_ts_as = represent_ts_as
        # no cutoff was previously passed as -1. This check is to maintain this
        # behaviour: negative number is appropriately turned to None
        self.cutoff = cutoff if cutoff is not None and cutoff >= 0 else None
        self.frame_data = None
        self.frame_rep = None

        self.auxstep = self._Auxstep(**kwargs)
        self._read_next_step()

        # if dt is constant and auxiliary data includes time, calculate
        # initial time and dt
        if self.time_selector is not None and self.constant_dt:
            self.auxstep._initial_time = self.time
            self._read_next_step()
            self.auxstep._dt = self.time - self.initial_time
            self.rewind()

    def copy(self):
        """Returns a deep copy of the AuxReader"""
        orig_dict = pickle.dumps(self)
        new_reader = pickle.loads(orig_dict)
        return new_reader

    def __len__(self):
        """Number of steps in auxiliary data."""
        return self.n_steps

    def next(self):
        """Move to next step of auxiliary data."""
        return self._read_next_step()

    def __next__(self):
        """Move to next step of auxiliary data."""
        return self.next()

    def __iter__(self):
        """Iterate over all auxiliary steps."""
        self._restart()
        return self

    def _restart(self):
        """Reset back to start; calling next() should read first step."""
        # Overwrite as appropriate
        self.auxstep.step = -1

    def rewind(self):
        """Return to and read first step."""
        # Overwrite as appropriate
        #  could also use _go_to_step(0)
        self._restart()
        return self._read_next_step()

    def _read_next_step(self):
        """Move to next step and update auxstep.

        Should return the AuxStep instance corresponding to the next step.
        """
        # Define in each auxiliary reader
        raise NotImplementedError(
            "BUG: Override _read_next_step() in auxiliary reader!"
        )

    def update_ts(self, ts):
        """Read auxiliary steps corresponding to and update the trajectory
        timestep *ts*.

        Calls :meth:`read_ts`, then updates *ts* with the representative value.
        ``auxname`` must be set; the representative value will be accessible in
        *ts* as ``ts.aux.auxname`` or ``ts.aux['auxname']``.

        Parameters
        ----------
        ts : :class:`~MDAnalysis.coordinates.timestep.Timestep` object
            The trajectory timestep for which corresponding auxiliary data is
            to be read and updated.

        Returns
        -------
        :class:`~MDAnalysis.coordinates.timestep.Timestep`
            *ts* with the representative auxiliary
            value in ``ts.aux`` be updated appropriately.

        Raises
        ------
        ValueError
            If ``auxname`` is not set.

        See Also
        --------
        :meth:`read_ts`
        """
        if not self.auxname:
            raise ValueError(
                "Auxiliary name not set, cannot set representative "
                "value in timestep. Name auxiliary or use read_ts "
                "instead"
            )
        self.read_ts(ts)
        setattr(ts.aux, self.auxname, self.frame_rep)
        return ts

    def read_ts(self, ts):
        """Read auxiliary steps corresponding to the trajectory timestep *ts*.

        Read the auxiliary steps 'assigned' to *ts* (the steps that are within
        ``ts.dt/2`` of of the trajectory timestep/frame - ie. closer to *ts*
        than either the preceding or following frame). Then calculate a
        'representative value' for the timestep from the data in each of these
        auxiliary steps.

        To update *ts* with the representative value, use ``update_ts`` instead.

        Parameters
        ----------
        ts : :class:`~MDAnalysis.coordinates.timestep.Timestep` object
            The trajectory timestep for which corresponding auxiliary data is
            to be read.

        See Also
        --------
        :meth:`update_ts`

        Note
        ----
        The auxiliary reader will end up positioned at the last step assigned
        to the trajectory frame or, if the frame includes no auxiliary steps,
        (as when auxiliary data are less frequent), the most recent auxiliary
        step before the frame.
        """
        # Make sure our auxiliary step starts at the right point (just before
        # the frame being read): the current step should be assigned to a
        # previous frame, and the next step to either the frame being read or a
        # following frame. Move to right position if not.
        frame_for_step = self.step_to_frame(self.step, ts)
        frame_for_next_step = self.step_to_frame(self.step + 1, ts)
        if frame_for_step is not None:
            if frame_for_next_step is None:
                # self.step is the last auxiliary step in memory.
                if frame_for_step >= ts.frame:
                    self.move_to_ts(ts)
            elif not (frame_for_step < ts.frame <= frame_for_next_step):
                self.move_to_ts(ts)

        self._reset_frame_data()  # clear previous frame data
        # read through all the steps 'assigned' to ts.frame + add to frame_data
        while self.step_to_frame(self.step + 1, ts) == ts.frame:
            self._read_next_step()
            self._add_step_to_frame_data(ts.time)
        self.frame_rep = self.calc_representative()

    def attach_auxiliary(
        self,
        coord_parent,
        aux_spec: Optional[Union[str, Dict[str, str]]] = None,
        format: Optional[str] = None,
        **kwargs,
    ) -> None:
        """Attaches the data specified in `aux_spec` to the `coord_parent`

        This method is called from within
        :meth:`MDAnalysis.coordinates.base.ReaderBase.add_auxiliary()`.
        `add_auxiliary` should be agnostic to the type of AuxReader, so the
        method call leads here instead. First, some checks are done on
        the arguments to make sure the input is treated properly. Then,
        the AuxReader(s) with appropriate :attr:`data_selector` are
        associated with the `coord_parent` from which `add_auxiliary` was
        called.

        Parameters
        ----------
        coord_parent : MDAnalysis.coordinates.base.ReaderBase
            Reader object to which to attach the auxiliary data.

        aux_spec : str, Dict[str, str], None
            Specifies which data to add to `coord_parent`. String types are
            for :class:`MDAnalysis.auxiliary.XVG.XVGReader` only. Dictionaries
            are the standard way of providing `aux_spec` information (see also:
            :mod:`MDAnalysis.auxiliary.EDR`).
            Passing `None` causes all data to be added.

        Returns
        -------
        None

        Raises
        ------
        ValueError
            If trying to add data under an `aux_spec` key that is already
            assigned.

        """
        if "auxname" in kwargs:
            # This check is necessary for the tests in coordinates/base
            # `test_reload_auxiliaries_from_description`
            # because there, auxreaders are created from descriptions without
            # giving explicit `aux_spec`s
            aux_spec = {kwargs["auxname"]: None}

        elif aux_spec is None:
            # Add all terms found in the file if aux_spec is None
            aux_spec = {term: term for term in self.terms}

        elif isinstance(aux_spec, str):
            # This is to keep XVGReader functioning as-is, working with strings
            # for what used to be `auxname`. Often, no `auxterms` are specified
            # for XVG files. The below sets the data selector to None,
            # the default for XVGReader
            aux_spec = {aux_spec: None}

        for auxname in aux_spec:
            if auxname in coord_parent.aux_list:
                raise ValueError(
                    f"Auxiliary data with name {auxname} already " "exists"
                )
            if " " in auxname:
                warnings.warn(
                    f"Auxiliary name '{auxname}' contains a space. "
                    "Only dictionary style accession, not attribute "
                    f"style accession of '{auxname}' will work."
                )
            description = self.get_description()
            # make sure kwargs are also used
            description_kwargs = {**description, **kwargs}
            # Make a copy of the auxreader for every attribute to add
            # This is necessary so all attributes can be iterated over
            aux = auxreader(**description_kwargs)
            aux.auxname = auxname
            if aux.data_selector is None:
                # When calling ReaderBase.copy(), aux_spec information is lost
                # but data_selector is retained
                aux.data_selector = aux_spec[auxname]
            coord_parent._auxs[auxname] = aux
            coord_parent.ts = aux.update_ts(coord_parent.ts)

        aux_memory_usage = 0
        # Check if testing, which needs lower memory limit
        memory_limit = kwargs.get("memory_limit", 1e09)
        for reader in coord_parent._auxs.values():
            aux_memory_usage += reader._memory_usage()
        if aux_memory_usage > memory_limit:
            conv = 1e09  # convert to GB
            warnings.warn(
                "AuxReader: memory usage warning! "
                f"Auxiliary data takes up {aux_memory_usage/conv} "
                f"GB of memory (Warning limit: {memory_limit/conv} "
                "GB)."
            )

    def _memory_usage(self):
        raise NotImplementedError(
            "BUG: Override _memory_usage() " "in auxiliary reader!"
        )

    def step_to_frame(self, step, ts, return_time_diff=False):
        """Calculate closest trajectory frame for auxiliary step *step*.

        Calculated given dt, time and frame from *ts*::

            time_frame_0 = ts.time - ts.frame*ts.dt   # time at frame 0

            frame = floor((step_to_time(step) - time_frame_0 + ts.dt/2)/ts.dt))

        The difference in time between the step and the calculated frame can
        also optionally be returned with *return_time_diff*.

        Parameters
        ----------
        step : int
            Number of the auxiliary step to calculate closest trajectory frame
            for.
        ts : :class:`~MDAnalysis.coordinates.timestep.Timestep` object
            (Any) timestep from the trajectory the calculated frame number is to
            correspond to.
        return_time_diff : bool, optional
            (Default: False) Additionally return the time difference between
            *step* and returned frame.

        Returns
        -------
        frame_index : int or None
            Number of the trajectory frame closest (in time) to the given
            auxiliary step. If the step index is out of range for the auxiliary
            data, ``None`` is returned instead.
        time_diff : float (optional)
            Difference in time between *step* and *frame_index*.

        Note
        ----
        Assumes trajectory dt is consant.
        The returned frame number may be out of range for the trajectory.
        """
        if step >= self.n_steps or step < 0:
            return None
        time_frame_0 = ts.time - ts.frame * ts.dt  # assumes ts.dt is constant
        time_step = self.step_to_time(step)
        frame_index = int(
            math.floor((time_step - time_frame_0 + ts.dt / 2.0) / ts.dt)
        )
        if not return_time_diff:
            return frame_index
        else:
            time_frame = time_frame_0 + frame_index * ts.dt
            time_diff = abs(time_frame - time_step)
            return frame_index, time_diff

    def move_to_ts(self, ts):
        """Position auxiliary reader just before trajectory timestep *ts*.

        Calling ``next()`` should read the first auxiliary step 'assigned' to
        the trajectory timestep *ts* or, if no auxiliary steps are
        assigned to that timestep (as in the case of less frequent auxiliary
        data), the first auxiliary step after *ts*.

        Parameters
        ----------
        ts : :class:`~MDAnalysis.coordinates.timestep.Timestep` object
            The trajectory timestep before which the auxiliary reader is to
            be positioned.
        """
        # figure out what step we want to end up at
        if self.constant_dt:
            # if dt constant, calculate from dt/offset/etc
            step = int(
                math.floor((ts.time - ts.dt / 2 - self.initial_time) / self.dt)
            )
            # if we're out of range of the number of steps, reset back
            step = max(min(step, self.n_steps - 1), -1)
        else:
            # otherwise, go through steps till we find the right one
            for i in range(self.n_steps + 1):
                if self.step_to_frame(i) >= ts.frame:
                    break
            # we want the step before this
            step = i - 1
        if step == -1:
            self._restart()
        else:
            self._go_to_step(step)

    def next_nonempty_frame(self, ts):
        """Find the next trajectory frame for which a representative auxiliary
        value can be calculated.

        That is, the next trajectory frame to which one or more auxiliary steps
        are assigned and fall within the cutoff.

        Starts looking from the current step time. If the end of the auxiliary
        data is reached before a trajectory frame is found, None is returned.

        Parameters
        ----------
        ts : :class:`~MDAnalysis.coordinates.timestep.Timestep` object
            Any timestep from the trajectory for which the next 'non-empty'
            frame is to be found.

        Returns
        -------
        int
            Index of the next auxiliary-containing frame in the trajectory.

        Note
        ----
        The returned index may be out of range for the trajectory.
        """
        step = self.step
        while step < self.n_steps - 1:
            next_frame, time_diff = self.step_to_frame(
                self.step + 1, ts, return_time_diff=True
            )
            if self.cutoff is not None and time_diff > self.cutoff:
                # 'representative values' will be NaN; check next step
                step = step + 1
            else:
                return next_frame
        # we ran out of auxiliary steps...
        return None

    def __getitem__(self, i):
        """Return the AuxStep corresponding to the *i*-th auxiliary step(s)
        (0-based). Negative numbers are counted from the end.

        *i* may be an integer (in which case the corresponding AuxStep is
        returned) or a list of integers or slice (in which case an iterator is
        returned)::

             step_10 = aux_reader[10]

        will move to step 10 of the auxiliary and return the :class:`AuxStep`.
        By using a slice/list, we can iterated over specified steps in the
        auxiliary, e.g. when performing analysis ::

             for auxstep in aux_reader[100:200]:   # analyse only steps 100 to 200
                 run_analysis(auxstep)

             for auxstep in aux_reader[::10]       # analyse every 10th step
                 run_analysis(auxstep)
        """
        if isinstance(i, numbers.Integral):
            i = self._check_index(i)
            return self._go_to_step(i)

        elif isinstance(i, (list, np.ndarray)):
            return self._list_iter([self._check_index(x) for x in i])

        elif isinstance(i, slice):
            # default start to first frame (ie. 0)
            start = self._check_index(i.start) if i.start is not None else 0
            # default stop to after last frame (i.e. n_steps)
            # n_steps is a valid stop index but will fail _check_index;
            # deal with separately
            stop = (
                i.stop
                if i.stop == self.n_steps
                else (
                    self._check_index(i.stop)
                    if i.stop is not None
                    else self.n_steps
                )
            )
            step = i.step or 1
            if not isinstance(step, numbers.Integral) or step < 1:
                raise ValueError("Step must be positive integer")  # allow -ve?
            if start > stop:
                raise IndexError("Stop frame is lower than start frame")
            return self._slice_iter(slice(start, stop, step))
        else:
            raise TypeError("Index must be integer, list of integers or slice")

    def _check_index(self, i):
        if not isinstance(i, numbers.Integral):
            raise TypeError("Step indices must be integers")
        if i < 0:
            i = i + self.n_steps
        if i < 0 or i >= self.n_steps:
            raise IndexError(
                "{} is out of range of auxiliary (num. steps "
                "{})".format(i, self.n_steps)
            )
        return i

    def _list_iter(self, i):
        for j in i:
            yield self._go_to_step(j)

    def _slice_iter(self, i):
        for j in range(i.start, i.stop, i.step):
            yield self._go_to_step(j)

    def _go_to_step(self, i):
        """Move to and read i-th auxiliary step."""
        # Need to define in each auxiliary reader
        raise NotImplementedError(
            "BUG: Override _go_to_step() in auxiliary reader!"
        )

    def _reset_frame_data(self):
        self.frame_data = {}

    def _add_step_to_frame_data(self, ts_time):
        """Update ``frame_data`` with values for the current step.

        Parameters
        ----------
        ts_time : float
            the time of the timestep the current step is being 'added to'. Used
            to calculate difference in time between current step and timestep.
        """
        time_diff = self.time - ts_time
        self.frame_data[time_diff] = self.auxstep.data

    def calc_representative(self):
        """Calculate representative auxiliary value(s) from the data in
        *frame_data*.

        Currently implemented options for calculating representative value are:

          * `closest`: default; the value(s) from the step closest to in time
            to the trajectory timestep

          * `average`: average of the value(s) from steps 'assigned' to the
            trajectory timestep.

        Additionally, if ``cutoff`` is specified, only steps within this time
        of the trajectory timestep are considered in calculating the
        representative.

        If no auxiliary steps were assigned to the timestep, or none fall
        within the cutoff, representative values are set to ``np.nan``.

        Returns
        -------
        ndarray
            Array of auxiliary value(s) 'representative' for the timestep.
        """
        if self.cutoff is None:
            cutoff_data = self.frame_data
        else:
            cutoff_data = {
                key: val
                for key, val in self.frame_data.items()
                if abs(key) <= self.cutoff
            }

        if len(cutoff_data) == 0:
            # no steps are 'assigned' to this trajectory frame, so return
            # values of ``np.nan``
            value = self.auxstep._empty_data()
        elif self.represent_ts_as == "closest":
            min_diff = min([abs(i) for i in cutoff_data])
            # we don't know the original sign, and might have two equally-spaced
            # steps; check the earlier time first
            try:
                value = cutoff_data[-min_diff]
            except KeyError:
                value = cutoff_data[min_diff]
        elif self.represent_ts_as == "average":
            try:
                value = np.mean(
                    np.array([val for val in cutoff_data.values()]), axis=0
                )
            except TypeError:
                # for readers like EDRReader, the above does not work
                # because each step contains a dictionary of numpy arrays
                # as data
                value = defaultdict(float)
                for dataset in cutoff_data:
                    for term in self.terms:
                        value[term] += cutoff_data[dataset][term]
                for term in value:
                    value[term] = value[term] / len(cutoff_data)
        return value

    def __enter__(self):
        return self

    def __exit__(self, exc_type, exc_val, exc_tb):
        self.close()
        return False

    def close(self):
        # Overwrite as appropriate
        pass

    @property
    def n_steps(self):
        """Total number of steps in the auxiliary data."""
        try:
            return self._n_steps
        except AttributeError:
            self._n_steps = self._count_n_steps()
            return self._n_steps

    def step_to_time(self, i):
        """Return time of auxiliary step *i*.

        Calculated using ``dt`` and ``initial_time`` if ``constant_dt`` is True;
        otherwise from the list of times as read from the auxiliary data for
        each step.

        Parameters
        ----------
        i : int
            Index (0-based) of step to return time for

        Returns
        -------
        time : float
            Time (in ps) of step *i*

        Raises
        ------
        ValueError
            When *i* not in valid range
        """
        if i >= self.n_steps:
            raise ValueError(
                "{0} is not a valid step index (total number of "
                "steps is {1})".format(i, self.n_steps)
            )
        if self.constant_dt:
            return i * self.dt + self.initial_time
        else:
            try:
                return self._times[i]
            except AttributeError:
                self._times = self.read_all_times()
                return self._times[i]

    @property
    def represent_ts_as(self):
        """Method by which 'representative' timestep values of auxiliary data
        will be calculated.
        """
        return self._represent_ts_as

    @represent_ts_as.setter
    def represent_ts_as(self, new):
        if new not in self.represent_options:
            raise ValueError(
                "{0} is not a valid option for calculating "
                "representative value(s). Enabled options are: "
                "{1}".format(new, self.represent_options)
            )
        self._represent_ts_as = new

    def __del__(self):
        self.close()

    def get_description(self):
        """Get the values of the parameters necessary for replicating the
        AuxReader.

        An AuxReader can be duplicated using
        :func:`~MDAnalysis.auxiliary.core.auxreader`::

            description = original_aux.get_description()
            new_aux = MDAnalysis.auxiliary.auxreader(**description)

        The resulting dictionary may also be passed directly to
        :meth:`~MDAnalysis.coordinates.base.ProtoReader.add_auxiliary` to
        reload an auxiliary into a trajectory::

             trajectory.add_auxiliary(**description)

        Returns
        -------
        dict
            Key-word arguments and values that can be used to replicate the
            AuxReader.
        """
        description = {
            attr.strip("_"): getattr(self, attr)
            for attr in self.required_attrs
        }
        return description

    def __eq__(self, other):
        for attr in self.required_attrs:
            if getattr(self, attr) != getattr(other, attr):
                return False
        return True

    @property
    def step(self):
        """Number of the current auxiliary step (as stored in ``auxstep``;
        0-based)."""
        return self.auxstep.step

    @property
    def time(self):
        """Time of current auxiliary step (as stored in ``auxstep``; in ps)"""
        return self.auxstep.time

    @property
    def dt(self):
        """Change in time between auxiliary steps (as stored in ``auxstep``;
        in ps)"""
        return self.auxstep._dt

    @property
    def initial_time(self):
        """Time of first auxiliary step (as stored in ``auxstep``; in ps)"""
        return self.auxstep._initial_time

    @property
    def time_selector(self):
        """Key to select 'time' value from the full set of data read for each step.
        As stored in ``austep``.

        Type differs between auxiliary formats, depending how the data for each
        step is read in and stored; e.g. data from .xvg files is read in as a
        list and `time_selector` must be a valid index. If time selection is not
        enabled by the reader, ``time_selector`` will default to ``None``.

        See each individual auxiliary reader.
        """
        return self.auxstep._time_selector

    @time_selector.setter
    def time_selector(self, new):
        old = self.auxstep._time_selector
        self.auxstep._time_selector = new
        if old != new:
            # if constant_dt is False and so we're using a _times list, this will
            # now be made invalid
            try:
                del self._times
            except AttributeError:
                pass

    @property
    def data_selector(self):
        """Key(s) to select auxiliary data values of interest from the full set
        of data read for each step (as stored in ``auxstep``).

        Type differs between auxiliary formats, depending how the data for each
        step is read in and stored - e.g. data from .xvg files is read in as
        a list and `data_selector` must be a list of valid indicies. If data
        selection is not enabled by the reader, ``data_selector`` will default
        to ``None``.

        See each individual auxiliary reader.
        """
        return self.auxstep._data_selector

    @data_selector.setter
    def data_selector(self, new):
        self.auxstep._data_selector = new

    @property
    def constant_dt(self):
        """True if ``dt`` is constant throughout the auxiliary (as stored in
        ``auxstep``)"""
        return self.auxstep._constant_dt

    @constant_dt.setter
    def constant_dt(self, new):
        self.auxstep._constant_dt = new


class AuxFileReader(AuxReader):
    """Base class for auxiliary readers that read from file.

    Extends AuxReader with attributes and methods particular to reading
    auxiliary data from an open file, for use when auxiliary files may be too
    large to read in at once.

    Parameters
    ----------
    filename : str
       Location of the file containing the auxiliary data.
    **kwargs
       Other AuxReader options.

    See also
    --------
    :class:`AuxReader`

    Attributes
    ----------
    auxfile
       File object for the auxiliary file.

    """

    def __init__(self, filename, **kwargs):
        self.auxfile = anyopen(filename)
        self._auxdata = os.path.abspath(filename)
        super(AuxFileReader, self).__init__(**kwargs)

    def close(self):
        """Close *auxfile*."""
        if self.auxfile is None:
            return
        self.auxfile.close()
        self.auxfile = None

    def _restart(self):
        """Reposition to just before first step."""
        self.auxfile.seek(0)
        self.auxstep.step = -1

    def _reopen(self):
        """Close and then reopen *auxfile*."""
        if self.auxfile is not None:
            self.auxfile.close()
        self.auxfile = open(self._auxdata)
        self.auxstep.step = -1

    def _go_to_step(self, i):
        """Move to and read i-th auxiliary step.

        Parameters
        ----------
        i : int
            Step number (0-indexed) to move to

        Raises
        ------
        ValueError
            If step index not in valid range.

        Note
        ----
        Works by reading through all steps consecutively until correct step
        is reached. Overwrite if this can be done more efficiently.
        """
        ## could seek instead?
        if i >= self.n_steps:
            raise ValueError(
                "Step index {0} is not valid for auxiliary "
                "(num. steps {1}!".format(i, self.n_steps)
            )
        value = self.rewind()
        while self.step != i:
            value = self.next()
        return value