diff --git a/src/verifai/features/features.py b/src/verifai/features/features.py
index 8634ceb..2b08b58 100644
--- a/src/verifai/features/features.py
+++ b/src/verifai/features/features.py
@@ -9,6 +9,7 @@
 import random
 import itertools
 import functools
+from abc import ABC, abstractmethod
 from collections import OrderedDict, namedtuple
 
 import numpy as np
@@ -964,19 +965,26 @@ def distance(self, valueA, valueB):
         else:
             return self.distanceMetric(valueA, valueB)
 
-    @cached_property
-    def fixedDomains(self):
+    @staticmethod
+    def _timeExpandDomain(domain, timeBound):
+        return domain
+
+    def fixedDomains(self, timeBound):
         """Return the fixed-length Domains associated with this feature."""
+        timeExpandedDomain = self._timeExpandDomain(self.domain, timeBound) if timeBound is not None else self.domain
+
         if not self.lengthDomain:
-            return self.domain
-        domains = {}
-        for length in self.lengthDomain:
-            length = length[0]
-            domains[length] = Array(self.domain, (length,))
+            domains = timeExpandedDomain
+        else:
+            domains = {}
+            for length in self.lengthDomain:
+                length = length[0]
+                domains[length] = Array(timeExpandedDomain, (length,))
+        
         return domains
 
     def __repr__(self):
-        rep = f'Feature({self.domain}'
+        rep = f'{self.__class__.__name__}({self.domain}'
         if self.distribution is not None:
             rep += f', distribution={self.distribution}'
         if self.lengthDomain is not None:
@@ -987,11 +995,113 @@ def __repr__(self):
             rep += f', distanceMetric={self.distanceMetric}'
         return rep + ')'
 
-### Feature spaces
+class TimeSeriesFeature(Feature):
+    """A feature with a value at each timesetep of a simulation."""
+    @staticmethod
+    def _timeExpandDomain(domain, timeBound):
+        return Array(domain, (timeBound,))
+
+### Feature spaces and Samples
+
+class Sample(ABC):
+    """A sample from a feature space, containing static points and able to generate dynamic points.
+    
+    Args:
+        space (FeatureSpace): The feature space this Sample was sampled from.
+        dynamicSampleLengths (dict): A dictionary containing the lengths of each dynamic feature
+            with a length domain.
+    """
+    def __init__(self, space, dynamicSampleLengths):
+        self.space = space
+        self.dynamicSampleHistory = []
+        self.dynamicSampleLengths = dynamicSampleLengths
+
+    @property
+    @abstractmethod
+    def staticSample(self):
+        pass
+
+    @abstractmethod
+    def _getDynamicSample(self, info):
+        pass
+
+    def getDynamicSample(self, info=None):
+        sample = self._getDynamicSample(info)
+        self.dynamicSampleHistory.append(sample)
+        return sample
+
+    @abstractmethod
+    def update(self, rho):
+        pass
+
+    def __getattr__(self, attr):
+        space = super().__getattribute__("space")
+        if attr in space.staticFeatureNamed:
+            return getattr(self.staticSample, attr)
+        elif attr in space.dynamicFeatureNamed:
+            class DynamicFeatureHelper:
+                def __init__(self, dynamicSampleHistory, attr):
+                    self.dynamicSampleHistory = dynamicSampleHistory
+                    self.attr = attr
+
+                def __getitem__(self, i):
+                    if i > len(self.dynamicSampleHistory):
+                        raise ValueError("Attempting to access dynamic sample value that has not been sampled.")
+                    return getattr(self.dynamicSampleHistory[i], self.attr)
+
+            return DynamicFeatureHelper(self.dynamicSampleHistory, attr)
+        else:
+            return super().__getattr__(attr)
+
+
+class CompleteSample(Sample):
+    """A completed sample, which has fully computed static and dynamic points.
+    
+    Args:
+        space (FeatureSpace): The feature space this Sample was sampled from.
+        staticSample: The static point for this sample
+        dynamicSampleList: A list of the dynamic points for this sample.
+        updateCallback: A callback that is called with the value passed to update.
+        dynamicSampleLengths (dict): A dictionary containing the lengths of each dynamic feature
+            with a length domain.
+    """
+    def __init__(self, space, staticSample, dynamicSampleList, updateCallback, dynamicSampleLengths):
+        super().__init__(space, dynamicSampleLengths)
+        self._staticSample = staticSample
+        self._dynamicSampleList = dynamicSampleList
+        self._updateCallback = updateCallback
+        self._i = 0
+
+    @property
+    def staticSample(self):
+        return self._staticSample
+
+    def _getDynamicSample(self, info):
+        if self.space.timeBound == 0:
+            raise RuntimeError("Called `getDynamicSample` with `timeBound` of `FeatureSpace` set to 0")
+
+        if self._i >= self.space.timeBound:
+            raise RuntimeError("Exceeded `timeBound` of `FeatureSpace`")
+        
+        assert self._i < len(self._dynamicSampleList)
+
+        dynamic_sample = self._dynamicSampleList[self._i]
+        self._i += 1
+
+        return dynamic_sample
+
+    def update(self, rho):
+        return self._updateCallback(rho)
 
 class FeatureSpace:
     """A space consisting of named features.
 
+    Args:
+        features (Iterable): An iterable of the `Feature` objects in this space.
+        distanceMetric (function; optional): An optional distance metric to be used with this space.
+        timeBound (int; optional): An upper bound on the number of timesteps of a simulation using
+            this space.
+
     .. testcode::
 
         FeatureSpace({
@@ -1001,22 +1111,37 @@ class FeatureSpace:
         })
     """
 
-    def __init__(self, features, distanceMetric=None):
+    def __init__(self, features, distanceMetric=None, timeBound=0):
         self.namedFeatures = tuple(sorted(features.items(), key=lambda i: i[0]))
         self.featureNamed = OrderedDict(self.namedFeatures)
         self.features = tuple(self.featureNamed.values())
-        self.makePoint = namedtuple('SpacePoint', self.featureNamed.keys())
+
+        self.staticFeatureNamed = OrderedDict({name: feat for name, feat in self.featureNamed.items()
+                                               if not isinstance(feat, TimeSeriesFeature)})
+        self.dynamicFeatureNamed = OrderedDict({name: feat for name, feat in self.featureNamed.items()
+                                               if isinstance(feat, TimeSeriesFeature)})
+
+        self.hasTimeSeries = len(self.dynamicFeatureNamed) > 0
+
+        self.makeStaticPoint = namedtuple('StaticSpacePoint', self.staticFeatureNamed)
+        self.makeDynamicPoint = namedtuple('DynamicSpacePoint', self.dynamicFeatureNamed)
+
         self.distanceMetric = distanceMetric
+        self.timeBound = timeBound
+
+        if len(self.dynamicFeatureNamed) > 0 and self.timeBound == 0:
+            raise RuntimeError("FeatureSpace which includes TimeSeriesFeature has no timeBound.")
 
     @cached_property
     def domains(self):
-        """Return the domain or domains associated with this space.
+        """Return the expanded domain or domains associated with this space.
 
         Returns a pair consisting of the Domain of all lengths of feature
         lists, plus a dict mapping each (flattened) point in that Domain to the
         corresponding Domain of other features. If the FeatureSpace has no
         feature lists, then returns (None, dom) where dom is the fixed Domain
-        of all features.
+        of all features. If any Features are TimeSeriesFeatures then they are
+        expanded to a max of timeBound.
         """
         fixedDomains = {}
         lengthDomains = {}
@@ -1024,9 +1149,9 @@ def domains(self):
         for name, feature in self.namedFeatures:
             if feature.lengthDomain:
                 lengthDomains[name] = feature.lengthDomain
-                variableDomains[name] = feature.fixedDomains
+                variableDomains[name] = feature.fixedDomains(self.timeBound)
             else:
-                fixedDomains[name] = feature.domain
+                fixedDomains[name] = feature._timeExpandDomain(feature.domain, self.timeBound)
         if len(lengthDomains) == 0:
             return (None, Struct(fixedDomains))
         lengthDomain = Struct(lengthDomains)
@@ -1056,16 +1181,18 @@ def flatten(self, point, fixedDimension=False):
         """Flatten a point in this space. See Domain.flatten.
 
         If fixedDimension is True, the point is flattened out as if all feature
-        lists had their maximum lengths, with None as a placeholder. This means
-        that all points in the space will flatten to the same length.
+        lists had their maximum lengths and time steps, with None as a placeholder.
+        This means that all points in the space will flatten to the same length.
         """
+        assert isinstance(point, Sample)
+
         flattened = []
-        for feature, value in zip(self.features, point):
+        for feature, value in zip(self.staticFeatureNamed.values(), point.staticSample):
             domain = feature.domain
             if feature.lengthDomain:
                 length = len(value)
                 flattened.append(length)
-                fixedDomain = feature.fixedDomains[length]
+                fixedDomain = feature.fixedDomains(None)[length]
                 fixedDomain.flattenOnto(value, flattened)
                 if fixedDimension:      # add padding to maximum length
                     sizePerElt = domain.flattenedDimension
@@ -1074,6 +1201,43 @@ def flatten(self, point, fixedDimension=False):
                         flattened.append(None)
             else:
                 domain.flattenOnto(value, flattened)
+
+        if self.hasTimeSeries:
+            duration = len(point.dynamicSampleHistory)
+            flattened.append(duration)
+
+            for feature_i, f in enumerate(self.dynamicFeatureNamed.items()):
+                feature_name, feature = f
+                domain = feature.domain
+                sizePerElt = domain.flattenedDimension
+
+                if feature.lengthDomain:
+                    length = point.dynamicSampleLengths[feature_name]
+                else:
+                    length = None
+
+                flattened.append(length)
+
+                for dynamic_point in point.dynamicSampleHistory:
+                    value = dynamic_point[feature_i]
+
+                    if length is None:
+                        domain.flattenOnto(value, flattened)
+                    else:
+                        fixedDomain = feature.fixedDomains(None)[length]
+                        fixedDomain.flattenOnto(value, flattened)
+                        if fixedDimension:
+                            needed = (feature.maxLength - length) * sizePerElt
+                            for i in range(needed):
+                                flattened.append(None)
+
+                if fixedDimension:
+                    needed = (self.timeBound - len(point.dynamicSampleHistory)) * feature.maxLength * sizePerElt
+                    flattened += [None for _ in range(needed)]
+
+        flattened_point = tuple(flattened)
+        if fixedDimension:
+            assert len(flattened_point) == self.fixedFlattenedDimension
         return tuple(flattened)
 
     @cached_property
@@ -1083,13 +1247,16 @@ def fixedFlattenedDimension(self):
         Also an upper bound on the length of the vector returned by flatten
         by default, when fixedDimension=False."""
         dim = 0
+        if self.hasTimeSeries:
+            dim += 1 # Timesteps
         for feature in self.features:
             domain = feature.domain
+            timeMult = self.timeBound if isinstance(feature, TimeSeriesFeature) else 1
             if feature.lengthDomain:
                 dim += 1    # dimension storing length of the feature list
-                dim += feature.maxLength * domain.flattenedDimension
+                dim += timeMult * feature.maxLength * domain.flattenedDimension
             else:
-                dim += domain.flattenedDimension
+                dim += timeMult * domain.flattenedDimension
         return dim
 
     def meaningOfFlatCoordinate(self, index, pointName='point'):
@@ -1100,7 +1267,7 @@ def meaningOfFlatCoordinate(self, index, pointName='point'):
         have different meaning depending on the lengths of feature lists.
         """
         assert 0 <= index < self.fixedFlattenedDimension
-        for name, feature in self.namedFeatures:
+        for name, feature in self.staticFeatureNamed.items():
             domain = feature.domain
             if feature.lengthDomain:
                 if index == 0:
@@ -1120,6 +1287,35 @@ def meaningOfFlatCoordinate(self, index, pointName='point'):
                     return domain.meaningOfFlatCoordinate(index,
                                                           pointName=subPoint)
                 index -= domain.flattenedDimension
+
+        if index == 0:
+            return f'len({pointName}.dynamicSampleHistory)'
+        index -= 1
+
+        for name, feature in self.dynamicFeatureNamed.items():
+            domain = feature.domain
+            if feature.lengthDomain:
+                if index == 0:
+                    return f'{pointName}.dynamicSampleLengths["{name}"]'
+                else:
+                    index -= 1
+                    for time_i in range(self.timeBound):
+                        elem = index // domain.flattenedDimension
+                        if elem < feature.maxLength:
+                            subIndex = index % domain.flattenedDimension
+                            subPoint = f'{pointName}.{name}[{time_i}][{elem}]'
+                            return domain.meaningOfFlatCoordinate(subIndex,
+                                pointName=subPoint)
+                        index -= feature.maxLength * domain.flattenedDimension
+            else:
+                index -= 1
+                for time_i in range(self.timeBound):
+                    if index < domain.flattenedDimension:
+                        subPoint = f'{pointName}.{name}[{time_i}]'
+                        return domain.meaningOfFlatCoordinate(index,
+                                                            pointName=subPoint)
+                    index -= domain.flattenedDimension
+
         raise RuntimeError('impossible index arithmetic')
 
     def pandasIndexForFlatCoordinate(self, index):
@@ -1128,7 +1324,7 @@ def pandasIndexForFlatCoordinate(self, index):
         See meaningOfFlatCoordinate, and Domain.pandasIndexForFlatCoordinate.
         """
         assert 0 <= index < self.fixedFlattenedDimension
-        for name, feature in self.namedFeatures:
+        for name, feature in self.staticFeatureNamed.items():
             domain = feature.domain
             if feature.lengthDomain:
                 if index == 0:
@@ -1146,15 +1342,42 @@ def pandasIndexForFlatCoordinate(self, index):
                     panda = domain.pandasIndexForFlatCoordinate(index)
                     return (name,) + panda
                 index -= domain.flattenedDimension
+
+        if index == 0:
+            return ("dynamicSampleHistory", "length")
+        index -= 1
+
+        for name, feature in self.dynamicFeatureNamed.items():
+            domain = feature.domain
+            if feature.lengthDomain:
+                if index == 0:
+                    return (name, 'length')
+                else:
+                    index -= 1
+                    for time_i in range(self.timeBound):
+                        elem = index // domain.flattenedDimension
+                        if elem < feature.maxLength:
+                            subIndex = index % domain.flattenedDimension
+                            panda = domain.pandasIndexForFlatCoordinate(subIndex)
+                            return (name, elem) + panda
+                        index -= feature.maxLength * domain.flattenedDimension
+            else:
+                for time_i in range(self.timeBound):
+                    if index < domain.flattenedDimension:
+                        panda = domain.pandasIndexForFlatCoordinate(index)
+                        return (name,) + panda
+                    index -= domain.flattenedDimension
+
         raise RuntimeError('impossible index arithmetic')
 
+
     def coordinateIsNumerical(self, index):
         """Whether the value of a coordinate is intrinsically numerical.
 
         See meaningOfFlatCoordinate, and Domain.coordinateIsNumerical.
         """
         assert 0 <= index < self.fixedFlattenedDimension
-        for name, feature in self.namedFeatures:
+        for name, feature in self.staticFeatureNamed.items():
             domain = feature.domain
             if feature.lengthDomain:
                 if index == 0:
@@ -1170,29 +1393,99 @@ def coordinateIsNumerical(self, index):
                 if index < domain.flattenedDimension:
                     return domain.coordinateIsNumerical(index)
                 index -= domain.flattenedDimension
+
+        if index == 0:
+            return True
+        index -= 1
+
+        for name, feature in self.dynamicFeatureNamed.items():
+            domain = feature.domain
+            if feature.lengthDomain:
+                if index == 0:
+                    return True
+                else:
+                    index -= 1
+                    for time_i in range(self.timeBound):
+                        elem = index // domain.flattenedDimension
+                        if elem < feature.maxLength:
+                            subIndex = index % domain.flattenedDimension
+                            return domain.coordinateIsNumerical(subIndex)
+                        index -= feature.maxLength * domain.flattenedDimension
+            else:
+                for time_i in range(self.timeBound):
+                    if index < domain.flattenedDimension:
+                        return domain.coordinateIsNumerical(index)
+                    index -= domain.flattenedDimension
+
         raise RuntimeError('impossible index arithmetic')
 
     def unflatten(self, coords, fixedDimension=False):
         """Unflatten a tuple of coordinates to a point in this space."""
-        values = []
+        staticValues = []
         iterator = iter(coords)
-        for feature in self.features:
+
+        for feature in self.staticFeatureNamed.values():
             domain = feature.domain
             if feature.lengthDomain:
                 length = next(iterator)
-                fixedDomain = feature.fixedDomains[length]
-                values.append(fixedDomain.unflattenIterator(iterator))
+                fixedDomain = feature.fixedDomains(None)[length]
+                staticValues.append(fixedDomain.unflattenIterator(iterator))
                 if fixedDimension:      # consume padding
                     sizePerElt = domain.flattenedDimension
                     needed = (feature.maxLength - length) * sizePerElt
                     for i in range(needed):
                         next(iterator)
             else:
-                values.append(domain.unflattenIterator(iterator))
-        return self.makePoint(*values)
+                staticValues.append(domain.unflattenIterator(iterator))
+
+        staticSample = self.makeStaticPoint(*staticValues)
+
+        if self.hasTimeSeries:
+            duration = next(iterator)
+
+            dynamicValuesList = [[] for _ in range(duration)]
+            dynamicSampleLengths = {}
+
+            for feature_name, feature in self.dynamicFeatureNamed.items():
+                domain = feature.domain
+                sizePerElt = domain.flattenedDimension
+                length = next(iterator)
+                dynamicSampleLengths[feature_name] = length
+                
+                for time_i in range(duration):
+                    if length is None:
+                        dynamicValuesList[time_i].append(domain.unflattenIterator(iterator))
+                    else:
+                        fixedDomain = feature.fixedDomains(None)[length]
+                        dynamicValuesList[time_i].append(fixedDomain.unflattenIterator(iterator))
+                        if fixedDimension:      # consume padding
+                            needed = (feature.maxLength - length) * sizePerElt
+                            for _ in range(needed):
+                                next(iterator)
+
+                if fixedDimension:
+                    needed = (self.timeBound - duration) * length * sizePerElt
+                    for _ in range(needed):
+                        next(iterator)
+
+            dynamicSampleList = [self.makeDynamicPoint(*dynamicValues) for dynamicValues in dynamicValuesList]
+        else:
+            duration = 0
+            dynamicSampleList = []
+            dynamicSampleLengths = {}
+
+        updateCallback = lambda rho: None
+
+        sample = CompleteSample(space=self, staticSample=staticSample, dynamicSampleList=dynamicSampleList, updateCallback=updateCallback,dynamicSampleLengths=dynamicSampleLengths)
+
+        for _ in range(duration):
+            sample.getDynamicSample()
+
+        return sample
 
     def __repr__(self):
         rep = f'FeatureSpace({self.featureNamed}'
         if self.distanceMetric is not None:
             rep += f', distanceMetric={self.distanceMetric}'
         return rep + ')'
+    
diff --git a/src/verifai/samplers/domain_sampler.py b/src/verifai/samplers/domain_sampler.py
index 6ead563..a1ba7b2 100644
--- a/src/verifai/samplers/domain_sampler.py
+++ b/src/verifai/samplers/domain_sampler.py
@@ -47,22 +47,12 @@ def update(self, sample, info, rho):
         """
         pass
 
-    def nextSample(self, feedback=None):
-        """Generate the next sample, given feedback from the last sample.
-
-        This exists only for backwards-compatibility. It has been replaced by
-        the getSample and update APIs.
-        """
-        if self.last_sample is not None:
-            self.update(self.last_sample, self.last_info, feedback)
-        self.last_sample, self.last_info = self.getSample()
-        return self.last_sample
-
     def __iter__(self):
         try:
-            feedback = None
             while True:
-                feedback = yield self.nextSample(feedback)
+                sample, info = self.getSample()
+                rho = yield sample
+                self.update(sample, info, rho)
         except TerminationException:
             return
 
diff --git a/src/verifai/samplers/feature_sampler.py b/src/verifai/samplers/feature_sampler.py
index 5890505..c3e42ca 100644
--- a/src/verifai/samplers/feature_sampler.py
+++ b/src/verifai/samplers/feature_sampler.py
@@ -8,8 +8,10 @@
 import dill
 from dotmap import DotMap
 import numpy as np
+from abc import ABC, abstractmethod
+from contextlib import contextmanager
 
-from verifai.features import FilteredDomain
+from verifai.features import FilteredDomain, TimeSeriesFeature, Sample, CompleteSample
 from verifai.samplers.domain_sampler import SplitSampler, TerminationException
 from verifai.samplers.rejection import RejectionSampler
 from verifai.samplers.halton import HaltonSampler
@@ -23,7 +25,7 @@
 
 ### Samplers defined over FeatureSpaces
 
-class FeatureSampler:
+class FeatureSampler(ABC):
     """Abstract class for samplers over FeatureSpaces."""
 
     def __init__(self, space):
@@ -149,29 +151,14 @@ def makeDomainSampler(domain):
                 makeRandomSampler)
         return LateFeatureSampler(space, RandomSampler, makeDomainSampler)
 
-    def getSample(self):
-        """Generate a sample, along with any sampler-specific info.
-
-        Must return a pair consisting of the sample and arbitrary
-        sampler-specific info, which will be passed to the `update`
-        method after the sample is evaluated.
-        """
-        raise NotImplementedError('tried to use abstract FeatureSampler')
-    
-    def update(self, sample, info, rho):
+    def update(self, sample_id, rho):
         """Update the state of the sampler after evaluating a sample."""
         pass
 
-    def nextSample(self, feedback=None):
-        """Generate the next sample, given feedback from the last sample.
-
-        This function exists only for backwards compatibility. It has been
-        superceded by the `getSample` and `update` APIs.
-        """
-        if self.last_sample is not None:
-            self.update(self.last_sample, self.last_info, feedback)
-        self.last_sample, self.last_info = self.getSample()
-        return self.last_sample
+    @abstractmethod
+    def getSample(self):
+        """Returns a `Sample` object"""
+        pass
 
     def set_graph(self, graph):
         self.scenario.set_graph(graph)
@@ -194,15 +181,22 @@ def restoreFromFile(path):
 
     def __iter__(self):
         try:
-            feedback = None
             while True:
-                feedback = yield self.nextSample(feedback)
+                sample = self.getSample()
+                rho = yield sample
+                sample.update(rho)
         except TerminationException:
             return
 
+
 class LateFeatureSampler(FeatureSampler):
-    """FeatureSampler that works by first sampling only lengths of feature
-    lists, then sampling from the resulting fixed-dimensional Domain.
+    """ FeatureSampler that greedily samples a CompleteSample.
+    
+    FeatureSampler works as follows:
+        1. Sample lengths of feature lists.
+        2. Expand TimeSeriesFeatures into flattened features with of length
+           space.timeBound.
+        2. Sample from the resulting fixed-dimensional Domains.
 
     e.g. LateFeatureSampler(space, RandomSampler, HaltonSampler) creates a
     FeatureSampler which picks lengths uniformly at random and applies
@@ -211,10 +205,11 @@ class LateFeatureSampler(FeatureSampler):
 
     def __init__(self, space, makeLengthSampler, makeDomainSampler):
         super().__init__(space)
+
         lengthDomain, fixedDomains = space.domains
         if lengthDomain is None:    # space has no feature lists
             self.lengthSampler = None
-            self.domainSampler = makeDomainSampler(fixedDomains)
+            self.domainSamplers = {None: makeDomainSampler(fixedDomains)}
         else:
             self.lengthDomain = lengthDomain
             self.lengthSampler = makeLengthSampler(lengthDomain)
@@ -222,30 +217,64 @@ def __init__(self, space, makeLengthSampler, makeDomainSampler):
                 point: makeDomainSampler(domain)
                 for point, domain in fixedDomains.items()
             }
-            self.lastLength = None
+
+        self._id_metadata_dict = {}
+        self._last_id = 0
+
+    def _get_info_id(self, info, length, sample):
+        self._last_id += 1
+        self._id_metadata_dict[self._last_id] = (info, length, sample)
+        return self._last_id
 
     def getSample(self):
         if self.lengthSampler is None:
-            domainPoint, info = self.domainSampler.getSample()
+            length, info1 = None, None
         else:
             length, info1 = self.lengthSampler.getSample()
-            self.lastLength = length
-            domainPoint, info2 = self.domainSamplers[length].getSample()
-            info = (info1, info2)
-        return self.space.makePoint(*domainPoint), info
-    
-    def update(self, sample, info, rho):
+
+        domainPoint, info2 = self.domainSamplers[length].getSample()
+        info = (info1, info2)
+        
+        sample_id = self._get_info_id(info, length, domainPoint)
+        update_callback = lambda rho: self.update(sample_id, rho)
+
+        # Make static points and iterable over dynamic points
+        static_features = [v for v in domainPoint._asdict().items()
+                            if v[0] in self.space.staticFeatureNamed]
+        dynamic_features = [v for v in domainPoint._asdict().items()
+                            if v[0] not in self.space.staticFeatureNamed]
+        static_point = self.space.makeStaticPoint(*[v[1] for v in static_features])
+
+        dynamic_points = []
+        if any(isinstance(f, TimeSeriesFeature) for f in self.space.features):
+            for t in range(self.space.timeBound):
+                point_dict = {}
+
+                for f, val in dynamic_features:
+                    if not self.space.featureNamed[f].lengthDomain:
+                        point_dict[f] = val[t]
+                    else:
+                        point_dict[f] = tuple(v[t] for v in val)
+
+                dynamic_points.append(self.space.makeDynamicPoint(*point_dict.values()))
+
+
+        dynamicSampleLengths = ({feature_name: getattr(length, feature_name)[0]
+                                 for feature_name, feature in self.space.dynamicFeatureNamed.items()
+                                 if feature.lengthDomain}
+                                if self.lengthSampler else {})
+
+        return CompleteSample(self.space, static_point, dynamic_points, update_callback, dynamicSampleLengths)
+
+    def update(self, sample_id, rho):
+        info, lengthPoint, domainPoint = self._id_metadata_dict[sample_id]
+
         if self.lengthSampler is None:
-            self.domainSampler.update(sample, info, rho)
+            self.domainSamplers[None].update(domainPoint, info[1], rho)
         else:
-            self.lengthSampler.update(sample, info[0], rho)
-            lengths = []
-            for name, feature in self.space.namedFeatures:
-                if feature.lengthDomain:
-                    lengths.append((len(getattr(sample, name)),))
-            lengthPoint = self.lengthDomain.makePoint(*lengths)
-            self.domainSamplers[lengthPoint].update(sample, info[1], rho)
+            self.lengthSampler.update(domainPoint, info[0], rho)
 
+            self.domainSamplers[lengthPoint].update(domainPoint, info[1], rho)
 ### Utilities
 
 def makeRandomSampler(domain):
diff --git a/src/verifai/samplers/scenic_sampler.py b/src/verifai/samplers/scenic_sampler.py
index bcea7fa..d919182 100644
--- a/src/verifai/samplers/scenic_sampler.py
+++ b/src/verifai/samplers/scenic_sampler.py
@@ -16,7 +16,7 @@
 
 from verifai.features import (Constant, Categorical, Real, Box, Array, Struct,
                               Feature, FeatureSpace)
-from verifai.samplers.feature_sampler import FeatureSampler
+from verifai.samplers.feature_sampler import FeatureSampler, Sample
 from verifai.utils.frozendict import frozendict
 
 scenicMajorVersion = int(importlib.metadata.version('scenic').split('.')[0])
@@ -223,6 +223,22 @@ def spaceForScenario(scenario, ignoredProperties):
     })
     return space, quotedParams
 
+class ScenicSample(Sample):
+    def __init__(self, space, staticSample, updateCallback, dynamicSampleLengths):
+        super().__init__(space, dynamicSampleLengths)
+        self._staticSample = staticSample
+        self._updateCallback = updateCallback
+
+    @property
+    def staticSample(self):
+        return self._staticSample
+
+    def _getDynamicSample(self, info):
+        raise RuntimeError("ScenicSampler does not support dynamic sampling.")
+
+    def update(self, rho):
+        self._updateCallback(rho)
+
 class ScenicSampler(FeatureSampler):
     """Samples from the induced distribution of a Scenic scenario.
 
@@ -236,14 +252,16 @@ class ScenicSampler(FeatureSampler):
     def __init__(self, scenario, maxIterations=None, ignoredProperties=None):
         self.scenario = scenario
         self.maxIterations = 2000 if maxIterations is None else maxIterations
+        self._nextScene = None
         self.lastScene = None
+        self.lastFeedback = None
         if ignoredProperties is None:
             ignoredProperties = defaultIgnoredProperties
         space, self.quotedParams = spaceForScenario(scenario, ignoredProperties)
         super().__init__(space)
 
     @classmethod
-    def fromScenario(cls, path, maxIterations=None,
+    def fromScenario(cls, path, maxIterations=None, maxSteps=None,
                      ignoredProperties=None, **kwargs):
         """Create a sampler corresponding to a Scenic program.
 
@@ -262,25 +280,44 @@ def fromScenario(cls, path, maxIterations=None,
               e.g. ``params`` to override global parameters or ``model`` to set the
               :term:`world model`.
         """
+        if "params" not in kwargs:
+            kwargs["params"] = {}
+
+        kwargs["params"]["timeBound"] = maxSteps if maxSteps else 0
+
         scenario = scenic.scenarioFromFile(path, **kwargs)
         return cls(scenario, maxIterations=maxIterations,
                    ignoredProperties=ignoredProperties)
 
     @classmethod
-    def fromScenicCode(cls, code, maxIterations=None,
+    def fromScenicCode(cls, code, maxIterations=None, maxSteps=None,
                        ignoredProperties=None, **kwargs):
         """As above, but given a Scenic program as a string."""
+        if "params" not in kwargs:
+            kwargs["params"] = {}
+
+        kwargs["params"]["timeBound"] = maxSteps if maxSteps else 0
+
         scenario = scenic.scenarioFromString(code, **kwargs)
         return cls(scenario, maxIterations=maxIterations,
                    ignoredProperties=ignoredProperties)
 
-    def nextSample(self, feedback=None):
+    def getSample(self):
         ret = self.scenario.generate(
-            maxIterations=self.maxIterations, feedback=feedback, verbosity=0
+            maxIterations=self.maxIterations, feedback=self.lastFeedback, verbosity=0
         )
+
+        self.lastFeedback = None
         self.lastScene, _ = ret
+
         return self.pointForScene(self.lastScene)
 
+    def update(self, sample_id, rho):
+        assert sample_id == 0
+        if self.lastFeedback is not None:
+            raise RuntimeError("Called `update` twice in a row (ScenicSampler does not support non-sequential sampling)")
+        self.lastFeedback = rho
+
     def pointForScene(self, scene):
         """Convert a sampled Scenic :obj:`~scenic.core.scenarios.Scene` to a point in our feature space.
 
@@ -314,7 +351,12 @@ def pointForScene(self, scene):
             params[param] = pointForValue(subdom, scene.params[originalName])
         paramPoint = paramDomain.makePoint(**params)
 
-        return self.space.makePoint(objects=objPoint, params=paramPoint)
+        staticSample = self.space.makeStaticPoint(objects=objPoint, params=paramPoint)
+
+        updateCallback = lambda rho: self.update(0, rho)
+        dynamicSampleLengths = []
+
+        return ScenicSample(self.space, staticSample, updateCallback, dynamicSampleLengths)
 
     @staticmethod
     def nameForObject(i):
diff --git a/src/verifai/server.py b/src/verifai/server.py
index ef4b043..87a75fa 100644
--- a/src/verifai/server.py
+++ b/src/verifai/server.py
@@ -144,6 +144,8 @@ def __init__(self, sampling_data, monitor, options={}):
                 sampler_params=params
             )
 
+        if self.sample_space.hasTimeSeries:
+            raise ValueError("Sample space for `Server` cannot contain `TimeSeriesFeature`")
 
     def listen(self):
         client_socket, addr = self.socket.accept()
@@ -176,8 +178,8 @@ def terminate(self):
     def close_connection(self):
         self.client_socket.close()
 
-    def get_sample(self, feedback):
-        return self.sampler.nextSample(feedback)
+    def get_sample(self):
+        return self.sampler.getSample()
 
     def flatten_sample(self, sample):
         return self.sampler.space.flatten(sample)
@@ -193,9 +195,10 @@ def evaluate_sample(self, sample):
 
     def run_server(self):
         start = time.time()
-        sample = self.get_sample(self.lastValue)
+        sample = self.get_sample()
         after_sampling = time.time()
-        self.lastValue = self.evaluate_sample(sample)
+        self.lastValue = self.evaluate_sample(sample.staticSample)
+        sample.update(self.lastValue)
         after_simulation = time.time()
         timings = ServerTimings(sample_time=(after_sampling - start),
                                 simulate_time=(after_simulation - after_sampling))
diff --git a/tests/scenic/scenic_driving_behavior.scenic b/tests/scenic/scenic_driving_behavior.scenic
new file mode 100644
index 0000000..20c0e67
--- /dev/null
+++ b/tests/scenic/scenic_driving_behavior.scenic
@@ -0,0 +1,19 @@
+param map = localPath('Town01.xodr')
+param carla_map = 'Town01'
+
+model scenic.domains.driving.model
+
+foo = TimeSeries(VerifaiRange(0,0.01))
+
+behavior TestBehavior():
+    lastVal = None
+    while True:
+        newVal = foo.getSample()
+        assert newVal != lastVal, (newVal, lastVal)
+        lastVal = newVal
+        take SetThrottleAction(newVal)
+
+ego = new Car on road, with behavior TestBehavior()
+new Car behind ego by VerifaiRange(1,4)
+
+terminate after 5 seconds
diff --git a/tests/scenic/test_scenic.py b/tests/scenic/test_scenic.py
index 2b3cf78..c1fb3f7 100644
--- a/tests/scenic/test_scenic.py
+++ b/tests/scenic/test_scenic.py
@@ -15,7 +15,7 @@ def test_objects(new_Object):
         f'ego = {new_Object} at 4 @ 9',
         maxIterations=1
     )
-    sample = sampler.nextSample()
+    sample = sampler.getSample()
     objects = sample.objects
     assert len(objects) == 1
     pos = objects.object0.position
@@ -28,7 +28,7 @@ def test_params(new_Object):
         f'ego = {new_Object}',
         maxIterations=1
     )
-    sample = sampler.nextSample()
+    sample = sampler.getSample()
     x = sample.params.x
     assert type(x) is float
     assert 3 <= x <= 5
@@ -39,7 +39,7 @@ def test_quoted_param(new_Object):
         f'ego = {new_Object}',
         maxIterations=1
     )
-    sample = sampler.nextSample()
+    sample = sampler.getSample()
     v = sampler.paramDictForSample(sample)['x/y']
     assert type(v) is float
     assert 3 <= v <= 5
@@ -49,7 +49,7 @@ def test_lists(new_Object):
         f'ego = {new_Object} with foo [1, -1, 3.3]',
         maxIterations=1
     )
-    sample = sampler.nextSample()
+    sample = sampler.getSample()
     foo = sample.objects.object0.foo
     assert type(foo) is tuple
     assert foo == pytest.approx((1, -1, 3.3))
@@ -68,7 +68,7 @@ def test_object_order(new_Object):
         f'    {new_Object} at 2*i @ 0',
         maxIterations=1
     )
-    sample = sampler.nextSample()
+    sample = sampler.getSample()
     objects = sample.objects
     assert len(objects) == 11
     for i in range(len(objects)):
@@ -78,7 +78,7 @@ def test_object_order(new_Object):
 
     flat = sampler.space.flatten(sample)
     unflat = sampler.space.unflatten(flat)
-    assert unflat == sample
+    assert unflat.staticSample == sample.staticSample
 
 ## Active sampling
 
@@ -112,7 +112,7 @@ def test_active_save_restore(new_Object, tmpdir):
 
 def runSampler(sampler):
     for i in range(3):
-        sample = sampler.nextSample()
+        sample = sampler.getSample()
         print(f'Sample #{i}:')
         print(sample)
 
@@ -179,3 +179,56 @@ def test_driving_dynamic(pathToLocalFile):
                                   server_class=ScenicServer,
                                   server_options=server_options)
     falsifier.run_falsifier()
+
+def test_driving_dynamic_behavior(pathToLocalFile):
+    path = pathToLocalFile('scenic_driving_behavior.scenic')
+    sampler = ScenicSampler.fromScenario(
+        path,
+        model='scenic.simulators.newtonian.driving_model',
+        params=dict(render=False),
+        mode2D=True,
+        maxSteps=2
+    )
+    falsifier_params = DotMap(
+        n_iters=3,
+        save_error_table=False,
+        save_safe_table=False,
+    )
+    server_options = DotMap(maxSteps=2, verbosity=3)
+    falsifier = generic_falsifier(sampler=sampler,
+                                  falsifier_params=falsifier_params,
+                                  server_class=ScenicServer,
+                                  server_options=server_options)
+    falsifier.run_falsifier()
+
+double_access_scenario = """
+model scenic.simulators.newtonian.model
+foo = TimeSeries(VerifaiRange(0, 0.01))
+behavior TestBehavior():
+    while True:
+        foo.getSample()
+        foo.getSample()
+        wait
+ego = new Object with behavior TestBehavior()
+"""
+
+def test_double_time_series_access():
+    with pytest.raises(RuntimeError):
+        sampler = ScenicSampler.fromScenicCode(
+                double_access_scenario,
+                model='scenic.simulators.newtonian.model',
+                maxIterations=1,
+                maxSteps=2,
+                params=dict(render=False),
+            )
+        falsifier_params = DotMap(
+            n_iters=3,
+            save_error_table=False,
+            save_safe_table=False,
+        )
+        server_options = DotMap(maxSteps=2, verbosity=3)
+        falsifier = generic_falsifier(sampler=sampler,
+                                    falsifier_params=falsifier_params,
+                                    server_class=ScenicServer,
+                                    server_options=server_options)
+        falsifier.run_falsifier()
diff --git a/tests/test_examples.py b/tests/test_examples.py
index 764ae04..a03bdcd 100644
--- a/tests/test_examples.py
+++ b/tests/test_examples.py
@@ -15,9 +15,9 @@ def test_example():
     sampler = FeatureSampler.samplerFor(space)
 
     for i in range(3):
-        sample = sampler.nextSample()
+        sample = sampler.getSample()
         print(f'Sample #{i}:')
         print(sample)
         flat = space.flatten(sample)
         unflat = space.unflatten(flat)
-        assert sample == unflat
+        assert sample.staticSample == unflat.staticSample
diff --git a/tests/test_features.py b/tests/test_features.py
index 3722ecc..1497d79 100644
--- a/tests/test_features.py
+++ b/tests/test_features.py
@@ -10,12 +10,12 @@ def test_fs_flatten():
     })
     sampler = FeatureSampler.randomSamplerFor(space)
     for i in range(100):
-        point = sampler.nextSample()
+        point = sampler.getSample()
         flat = space.flatten(point)
         assert type(flat) is tuple
         assert len(flat) <= space.fixedFlattenedDimension
         unflat = space.unflatten(flat)
-        assert point == unflat
+        assert point.staticSample == unflat.staticSample
 
 def test_fs_flatten_fixed_dimension():
     space = FeatureSpace({
@@ -25,7 +25,7 @@ def test_fs_flatten_fixed_dimension():
     assert space.fixedFlattenedDimension == 4
     sampler = FeatureSampler.randomSamplerFor(space)
     for i in range(100):
-        point = sampler.nextSample()
+        point = sampler.getSample()
         flat = space.flatten(point, fixedDimension=True)
         assert type(flat) is tuple
         assert len(flat) == 4
@@ -42,7 +42,7 @@ def test_fs_flatten_fixed_dimension():
         else:
             assert eval(space.meaningOfFlatCoordinate(3)) == point.b[1][0]
         unflat = space.unflatten(flat, fixedDimension=True)
-        assert point == unflat
+        assert point.staticSample == unflat.staticSample
     assert space.pandasIndexForFlatCoordinate(0) == ('a', 0)
     assert space.pandasIndexForFlatCoordinate(1) == ('b', 'length')
     assert space.pandasIndexForFlatCoordinate(2) == ('b', 0, 0)
@@ -58,7 +58,7 @@ def test_fs_flatten_fixed_dimension2():
     assert space.fixedFlattenedDimension == 6
     sampler = FeatureSampler.randomSamplerFor(space)
     for i in range(100):
-        point = sampler.nextSample()
+        point = sampler.getSample()
         flat = space.flatten(point, fixedDimension=True)
         assert type(flat) is tuple
         assert len(flat) == 6
@@ -76,7 +76,7 @@ def test_fs_flatten_fixed_dimension2():
         else:
             assert eval(space.meaningOfFlatCoordinate(3)) == point.b[2]
         unflat = space.unflatten(flat, fixedDimension=True)
-        assert point == unflat
+        assert point.staticSample == unflat.staticSample
     assert space.pandasIndexForFlatCoordinate(0) == ('b', 'length')
     assert space.pandasIndexForFlatCoordinate(1) == ('b', 0)
     assert space.pandasIndexForFlatCoordinate(2) == ('b', 1)
@@ -87,12 +87,53 @@ def test_fs_flatten_fixed_dimension2():
     assert not any(space.coordinateIsNumerical(i) for i in range(1, 4))
     assert all(space.coordinateIsNumerical(i) for i in range(4, 6))
 
+def test_fs_flatten_fixed_dimension_dynamic():
+    space = FeatureSpace({
+        'a': Feature(DiscreteBox([0, 12])),
+        'b': TimeSeriesFeature(Box((0, 1)), lengthDomain=DiscreteBox((0, 2)))
+        },
+        timeBound=5
+    )
+    assert space.fixedFlattenedDimension == 13
+    sampler = FeatureSampler.randomSamplerFor(space)
+    for i in range(100):
+        point = sampler.getSample()
+        duration = random.randint(0, 5)
+        for _ in range(duration):
+            point.getDynamicSample()
+        flat = space.flatten(point, fixedDimension=True)
+        assert type(flat) is tuple
+        assert len(flat) == 13
+        assert eval(space.meaningOfFlatCoordinate(0)) == point.a[0]
+        assert flat[1] == duration
+        for t in range(duration):
+            offset = 2*t
+            bLen = len(point.dynamicSampleHistory[t].b)
+            assert flat[2] == bLen
+            assert eval(space.meaningOfFlatCoordinate(2)) == bLen
+            if bLen < 1:
+                assert flat[offset+3] is None
+            else:
+                assert eval(space.meaningOfFlatCoordinate(offset+3)) == point.dynamicSampleHistory[t].b[0][0]
+            if bLen < 2:
+                assert flat[offset+4] is None
+            else:
+                assert eval(space.meaningOfFlatCoordinate(offset+4)) == point.dynamicSampleHistory[t].b[1][0]
+        unflat = space.unflatten(flat, fixedDimension=True)
+        assert point.staticSample == unflat.staticSample
+        assert all(point.dynamicSampleHistory[t] == unflat.dynamicSampleHistory[t] for t in range(duration))
+    assert space.pandasIndexForFlatCoordinate(0) == ('a', 0)
+    assert space.pandasIndexForFlatCoordinate(2) == ('b', 'length')
+    assert space.pandasIndexForFlatCoordinate(3) == ('b', 0, 0)
+    assert space.pandasIndexForFlatCoordinate(4) == ('b', 1, 0)
+    assert all(space.coordinateIsNumerical(i) for i in range(4))
+
 def test_fs_distance():
     box = Box([0, 10])
     space = FeatureSpace({ 'a': Feature(box), 'b': Feature(box) })
     sampler = FeatureSampler.randomSamplerFor(space)
-    pointA = sampler.nextSample()
-    pointB = sampler.nextSample()
+    pointA = sampler.getSample().staticSample
+    pointB = sampler.getSample().staticSample
     assert pointA != pointB
     assert space.distance(pointA, pointA) == 0
     assert space.distance(pointB, pointB) == 0
diff --git a/tests/test_grid.py b/tests/test_grid.py
index de16685..3477efc 100644
--- a/tests/test_grid.py
+++ b/tests/test_grid.py
@@ -14,7 +14,9 @@ def test_grid():
     dict_samples = defaultdict(int)
     while True:
         try:
-            sample = sampler.nextSample()
+            sample = sampler.getSample()
+            sample.update(None)
+            sample = sample
             dict_samples[(sample.weather[0], sample.car_positions[0],
                           sample.car_positions[1])] = 0
         except TerminationException:
@@ -46,7 +48,9 @@ def f(sample):
     y_samples = []
 
     for i in range(21):
-        sample = sampler.nextSample()
+        sample = sampler.getSample()
+        sample.update(None)
+        sample = sample
         samples.append(sample)
         y_samples.append(f(sample))
 
@@ -59,6 +63,6 @@ def test_grid_non_standardizable():
         'b': Feature(FilteredDomain(Box([0,1]), lambda x: x[0] > 0.5))
     })
     sampler = FeatureSampler.gridSamplerFor(space)
-    samples = list(sampler)
+    samples = [s for s in sampler]
     assert len(samples) == 13
     assert all(sample.b[0] > 0.5 for sample in samples)
diff --git a/tests/test_halton.py b/tests/test_halton.py
index 4c3cb7f..6981fd8 100644
--- a/tests/test_halton.py
+++ b/tests/test_halton.py
@@ -24,7 +24,7 @@ def test_halton():
 
     for i in range(3):
         print(f'Sample #{i}:')
-        print(sampler.nextSample())
+        print(sampler.getSample().staticSample)
 
 def test_save_restore(tmpdir):
     space = FeatureSpace({
diff --git a/tests/test_samplers.py b/tests/test_samplers.py
index 4389e80..08b6c40 100644
--- a/tests/test_samplers.py
+++ b/tests/test_samplers.py
@@ -3,9 +3,41 @@
 import os.path
 
 from verifai.features import (Struct, Array, Box, DiscreteBox,
-                              Feature, FeatureSpace)
+                              Feature, TimeSeriesFeature, FeatureSpace)
 from verifai.samplers import RandomSampler, FeatureSampler
 
+def test_feature_sampling():
+    space = FeatureSpace({
+        'a': Feature(DiscreteBox([0, 12])),
+        'b': Feature(Box((0, 1)), lengthDomain=DiscreteBox((0, 2))),
+        'c': TimeSeriesFeature(Box((2,5))),
+        'd': TimeSeriesFeature(Box((5,6)), lengthDomain=DiscreteBox((0,2)))
+    },
+    timeBound=10)
+    sampler = FeatureSampler.randomSamplerFor(space)
+
+    sample = sampler.getSample()
+    static_point = sample
+
+    assert len(static_point.a) == 1
+    assert 0 <= static_point.a[0] <= 12
+    assert 0 <= len(static_point.b) <= 2
+    assert all(0 <= v[0] <= 1 for v in static_point.b)
+
+    for _ in range(space.timeBound):
+        dynamic_point = sample.getDynamicSample()
+
+        assert len(dynamic_point.c) == 1
+        assert 2 <= dynamic_point.c[0] <= 5
+        assert 0 <= len(dynamic_point.d) <= 2
+        assert all(5 <= v[0] <= 6 for v in dynamic_point.d)
+
+    for i in range(space.timeBound):
+        assert len(sample.c[i]) == 1
+        assert 2 <= sample.c[i][0] <= 5
+        assert 0 <= len(sample.d[i]) <= 2
+        assert all(5 <= v[0] <= 6 for v in sample.d[i])
+
 ## Random sampling
 
 def test_domain_random():
@@ -36,7 +68,7 @@ def check(samples):
         assert any(sample[0][0].position[0] < sample[1][1].position[0]
                    for sample in samples)
 
-    check([sampler.nextSample() for i in range(100)])
+    check([sampler.getSample()[0] for _ in range(100)])
     check(list(itertools.islice(sampler, 100)))
 
 def test_space_random():
@@ -48,7 +80,6 @@ def test_space_random():
 
     def check(samples):
         for sample in samples:
-            assert type(sample) is space.makePoint
             a = sample.a
             assert type(a) is tuple
             assert len(a) == 1
@@ -72,8 +103,8 @@ def check(samples):
         assert any(len(sample.b) == 1 for sample in samples)
         assert any(len(sample.b) == 2 for sample in samples)
 
-    check([sampler.nextSample() for i in range(100)])
-    check(list(itertools.islice(sampler, 100)))
+    check([sampler.getSample() for i in range(100)])
+    check(list(s for s in itertools.islice(sampler, 100)))
 
 def test_random_restore(tmpdir):
     space = FeatureSpace({
@@ -84,7 +115,7 @@ def test_random_restore(tmpdir):
 
     path = os.path.join(tmpdir, 'blah.dat')
     sampler.saveToFile(path)
-    sample1 = sampler.nextSample()
+    sample1 = sampler.getSample()
     sampler = FeatureSampler.restoreFromFile(path)
-    sample2 = sampler.nextSample()
-    assert sample1 == sample2
+    sample2 = sampler.getSample()
+    assert sample1.staticSample == sample2.staticSample
diff --git a/tests/utils.py b/tests/utils.py
index 5b5d76f..5a410aa 100644
--- a/tests/utils.py
+++ b/tests/utils.py
@@ -4,11 +4,11 @@
 from verifai.samplers import FeatureSampler
 
 def sampleWithFeedback(sampler, num_samples, f):
-    feedback = None
     samples = []
     for i in range(num_samples):
-        sample = sampler.nextSample(feedback)
+        sample = sampler.getSample()
         feedback = f(sample)
+        sample.update(feedback)
         print(f'Sample #{i}:')
         print(sample)
         samples.append(sample)
@@ -19,18 +19,21 @@ def checkSaveRestore(sampler, tmpdir, iterations=1):
     feedback = None
     for i in range(iterations):
         sampler.saveToFile(path)
-        sample1 = sampler.nextSample(feedback)
-        sample2 = sampler.nextSample(-1)
+        sample1 = sampler.getSample()
+        sample1.update(-1)
+        sample2 = sampler.getSample()
         sampler = FeatureSampler.restoreFromFile(path)
-        sample1b = sampler.nextSample(feedback)
-        sample2b = sampler.nextSample(-1)
-        assert sample1 != sample2
-        assert sample1 == sample1b
-        assert sample2 == sample2b
+        sample1b = sampler.getSample()
+        sample1b.update(-1)
+        sample2b = sampler.getSample()
+        sample2b.update(1)
+        assert sample1.staticSample != sample2.staticSample
+        assert sample1.staticSample == sample1b.staticSample
+        assert sample2.staticSample == sample2b.staticSample
         sampler.saveToFile(path)
-        sample3 = sampler.nextSample(1)
+        sample3 = sampler.getSample()
         sampler = FeatureSampler.restoreFromFile(path)
-        sample3b = sampler.nextSample(1)
-        assert sample3 not in (sample1, sample2)
-        assert sample3 == sample3b
-        feedback = 1
+        sample3b = sampler.getSample()
+        assert sample3.staticSample not in (sample1.staticSample, sample2.staticSample)
+        assert sample3.staticSample == sample3b.staticSample
+        sample3b.update(1)