diff --git a/dimod/__init__.py b/dimod/__init__.py
index 6ea87623..e90ffac4 100644
--- a/dimod/__init__.py
+++ b/dimod/__init__.py
@@ -48,7 +48,7 @@ import dimod.higherorder
 
 from dimod.package_info import __version__, __author__, __authoremail__, __description__
 
-from dimod.sampleset import *
+from dimod.sampleset import as_samples, SampleSet, concatenate
 
 from dimod.response import *
 import dimod.response
diff --git a/dimod/sampleset.py b/dimod/sampleset.py
index db42668a..b6f5a4d7 100644
--- a/dimod/sampleset.py
+++ b/dimod/sampleset.py
@@ -27,6 +27,8 @@ from collections import namedtuple
 
 import numpy as np
 
+from numpy.lib import recfunctions
+
 from dimod.decorators import vartype_argument
 from dimod.utilities import resolve_label_conflict
 from dimod.vartypes import Vartype
@@ -147,6 +149,72 @@ def as_samples(samples_like, dtype=None):
     return samples_like, list(range(samples_like.shape[1]))
 
 
+def concatenate(samplesets, defaults=None):
+    """Combine SampleSets.
+
+    Args:
+        samplesets (iterable[:obj:`.SampleSet`):
+            An iterable of sample sets.
+
+        defaults (dict, optional):
+            Dictionary mapping data vector names to the corresponding default values.
+
+    Returns:
+        :obj:`.SampleSet`: A sample set with the same vartype and variable order as the first
+        given in `samplesets`.
+
+    Examples:
+        >>> a = dimod.SampleSet.from_samples(([-1, +1], 'ab'), dimod.SPIN, energy=-1)
+        >>> b = dimod.SampleSet.from_samples(([-1, +1], 'ba'), dimod.SPIN, energy=-1)
+        >>> ab = dimod.concatenate((a, b))
+        >>> ab.record.sample
+        array([[-1,  1],
+               [ 1, -1]], dtype=int8)
+
+    """
+
+    itertup = iter(samplesets)
+
+    try:
+        first = next(itertup)
+    except StopIteration:
+        raise ValueError("samplesets must contain at least one SampleSet")
+
+    vartype = first.vartype
+    variables = first.variables
+
+    records = [first.record]
+    records.extend(_iter_records(itertup, vartype, variables))
+
+    # dev note: I was able to get ~2x performance boost when trying to
+    # implement the same functionality here by hand (I didn't know that
+    # this function existed then). However I think it is better to use
+    # numpy's function and rely on their testing etc. If however this becomes
+    # a performance bottleneck in the future, it might be worth changing.
+    record = recfunctions.stack_arrays(records, defaults=defaults,
+                                       asrecarray=True, usemask=False)
+
+    return SampleSet(record, variables, {}, vartype)
+
+
+def _iter_records(samplesets, vartype, variables):
+    # coerce each record into the correct vartype and variable-order
+    for samples in samplesets:
+
+        # coerce vartype
+        if samples.vartype is not vartype:
+            samples = samples.change_vartype(vartype, inplace=False)
+
+        if samples.variables != variables:
+            new_record = samples.record.copy()
+            order = [samples.variables.index[v] for v in variables]
+            new_record.sample = samples.record.sample[:, order]
+            yield new_record
+        else:
+            # order matches so we're done
+            yield samples.record
+
+
 class SampleSet(Iterable, Sized):
     """Samples and any other data returned by dimod samplers.
 
@@ -174,8 +242,7 @@ class SampleSet(Iterable, Sized):
         >>> import dimod
         >>> import numpy as np
         ...
-        >>> dimod.SampleSet.from_samples(dimod.as_samples(np.ones(5, dtype='int8')),
-        ...                              'BINARY', 0)   # doctest: +SKIP
+        >>> dimod.SampleSet.from_samples(np.ones(5, dtype='int8'), 'BINARY', 0)   # doctest: +SKIP
         SampleSet(rec.array([([1, 1, 1, 1, 1], 0, 1)],
         ...       dtype=[('sample', 'i1', (5,)), ('energy', '<i4'), ('num_occurrences', '<i4')]),
         ...       [0, 1, 2, 3, 4], {}, 'BINARY')
@@ -222,7 +289,8 @@ class SampleSet(Iterable, Sized):
         self._vartype = vartype
 
     @classmethod
-    def from_samples(cls, samples_like, vartype, energy, info=None, num_occurrences=None, **vectors):
+    def from_samples(cls, samples_like, vartype, energy, info=None,
+                     num_occurrences=None, aggregate_samples=False, **vectors):
         """Build a :class:`SampleSet` from raw samples.
 
         Args:
@@ -243,7 +311,10 @@ class SampleSet(Iterable, Sized):
                 Information about the :class:`SampleSet` as a whole formatted as a dict.
 
             num_occurrences (array_like, optional):
-                Number of occurences for each sample. If not provided, defaults to a vector of 1s.
+                Number of occurrences for each sample. If not provided, defaults to a vector of 1s.
+
+            aggregate_samples (bool, optional, default=False):
+                If true, returned :obj:`.SampleSet` will have all unique samples.
 
             **vectors (array_like):
                 Other per-sample data.
@@ -265,6 +336,11 @@ class SampleSet(Iterable, Sized):
 
         .. _array_like:  https://docs.scipy.org/doc/numpy/user/basics.creation.html#converting-python-array-like-objects-to-numpy-arrays
         """
+        if aggregate_samples:
+            return cls.from_samples(samples_like, vartype, energy,
+                                    info=info, num_occurrences=num_occurrences,
+                                    aggregate_samples=False,
+                                    **vectors).aggregate()
 
         # get the samples, variable labels
         samples, variables = as_samples(samples_like)
@@ -299,6 +375,49 @@ class SampleSet(Iterable, Sized):
 
         return cls(record, variables, info, vartype)
 
+    @classmethod
+    def from_samples_bqm(cls, samples_like, bqm, **kwargs):
+        """Build a SampleSet from raw samples using a BinaryQuadraticModel to get energies and vartype.
+
+        Args:
+            samples_like:
+                A collection of raw samples. 'samples_like' is an extension of NumPy's array_like.
+                See :func:`.as_samples`.
+
+            bqm (:obj:`.BinaryQuadraticModel`):
+                A binary quadratic model. It is used to calculate the energies
+                and set the vartype.
+
+            info (dict, optional):
+                Information about the :class:`SampleSet` as a whole formatted as a dict.
+
+            num_occurrences (array_like, optional):
+                Number of occurrences for each sample. If not provided, defaults to a vector of 1s.
+
+            aggregate_samples (bool, optional, default=False):
+                If true, returned :obj:`.SampleSet` will have all unique samples.
+
+            **vectors (array_like):
+                Other per-sample data.
+
+        Returns:
+            :obj:`.SampleSet`
+
+        Examples:
+
+            >>> bqm = dimod.BinaryQuadraticModel.from_ising({}, {('a', 'b'): -1})
+            >>> samples = dimod.SampleSet.from_samples_bqm({'a': -1, 'b': 1}, bqm)
+            >>> samples =dimod.SampleSet.from_samples_bqm([[-1, 1], [1, -1]], bqm)
+
+        """
+        # more performant to do this once, here rather than again in bqm.energies
+        # and in cls.from_samples
+        samples_like = as_samples(samples_like)
+
+        energies = bqm.energies(samples_like)
+
+        return cls.from_samples(samples_like, energy=energies, vartype=bqm.vartype, **kwargs)
+
     @classmethod
     def from_future(cls, future, result_hook=None):
         """Construct a :class:`SampleSet` referencing the result of a future computation.
@@ -716,6 +835,33 @@ class SampleSet(Iterable, Sized):
         self._variables.relabel(mapping)
         return self
 
+    def aggregate(self):
+        """Create a new SampleSet with repeated samples aggregated.
+
+        Returns:
+            :obj:`.SampleSet`
+
+        Note:
+            :attr:`.SampleSet.record.num_occurrences` are accumulated but no
+            other fields are.
+
+        """
+
+        _, indices, inverse = np.unique(self.record.sample, axis=0,
+                                        return_index=True, return_inverse=True)
+
+        record = self.record[indices]
+
+        # fix the number of occurrences
+        record.num_occurrences = 0
+        for old_idx, new_idx in enumerate(inverse):
+            record[new_idx].num_occurrences += self.record[old_idx].num_occurrences
+
+        # dev note: we don't check the energies as they should be the same
+        # for individual samples
+
+        return self.__class__(record, self.variables, self.info, self.vartype)
+
     ###############################################################################################
     # Serialization
     ###############################################################################################
diff --git a/docs/reference/sampleset.rst b/docs/reference/sampleset.rst
index b6d2b212..74431514 100644
--- a/docs/reference/sampleset.rst
+++ b/docs/reference/sampleset.rst
@@ -57,6 +57,7 @@ Constructing a Response
    :toctree: generated/
 
    SampleSet.from_samples
+   SampleSet.from_samples_bqm
    SampleSet.from_future
 
 Transformations
@@ -76,6 +77,14 @@ Copy
 
    SampleSet.copy
 
+Aggregate
+~~~~~~~~~
+
+.. autosummary::
+   :toctree: generated/
+
+   SampleSet.aggregate
+
 Done
 ~~~~
 
@@ -83,3 +92,11 @@ Done
    :toctree: generated/
 
    SampleSet.done
+
+Utility Functions
+=================
+
+.. autosummary::
+   :toctree: generated/
+
+   concatenate

Ability to 'stack' multiple responses
**Application**
There are many cases when combining groups of samples (from different samplers say) is useful, however most of the naive approaches are inefficient. We also want to handle edge-cases like differently-ordered labels.

**Proposed Solution**
One of:

Pattern A:
```
combined = sampleset0 + sampleset1
```

Pattern B:
```
dimod.samples_stack(sampleset0, sampleset1)
```

@conta877 
diff --git a/tests/test_sampleset.py b/tests/test_sampleset.py
index 9ce1a7e5..08218875 100644
--- a/tests/test_sampleset.py
+++ b/tests/test_sampleset.py
@@ -123,6 +123,28 @@ class TestSampleSet(unittest.TestCase):
 
         self.assertEqual(len(dimod.SampleSet.from_samples(np.empty((0, 0)), dimod.SPIN, energy=[], a=1)), 0)
 
+    def test_from_samples_with_aggregation(self):
+        samples = dimod.SampleSet.from_samples(([[-1, 1], [-1, 1]], 'ab'), dimod.SPIN, energy=[0.0, 0.0],
+                                               aggregate_samples=True)
+        self.assertEqual(samples.aggregate(),
+                         dimod.SampleSet.from_samples(([-1, 1], 'ab'), dimod.SPIN, energy=0.0, num_occurrences=2))
+
+    def test_aggregate_simple(self):
+        samples = dimod.SampleSet.from_samples(([[-1, 1], [-1, 1]], 'ab'), dimod.SPIN, energy=[0.0, 0.0])
+
+        self.assertEqual(samples.aggregate(),
+                         dimod.SampleSet.from_samples(([-1, 1], 'ab'), dimod.SPIN, energy=0.0, num_occurrences=2))
+
+        # original should not be changed
+        self.assertEqual(samples,
+                         dimod.SampleSet.from_samples(([[-1, 1], [-1, 1]], 'ab'), dimod.SPIN, energy=[0.0, 0.0]))
+
+    def test_from_bqm_single_sample(self):
+        bqm = dimod.BinaryQuadraticModel.from_ising({}, {'ab': -1})
+        samples = dimod.SampleSet.from_samples_bqm({'a': -1, 'b': 1}, bqm)
+        self.assertEqual(samples,
+                         dimod.SampleSet.from_samples(([-1, 1], 'ab'), dimod.SPIN, energy=1))
+
 
 class TestSampleSetSerialization(unittest.TestCase):
 
@@ -178,3 +200,60 @@ class TestSampleSet_to_pandas_dataframe(unittest.TestCase):
                              columns=['a', 'b', 'c', 'energy', 'num_occurrences'])
 
         pd.testing.assert_frame_equal(df, other, check_dtype=False)
+
+
+class Test_concatenate(unittest.TestCase):
+    def test_simple(self):
+        ss0 = dimod.SampleSet.from_samples([-1, +1], dimod.SPIN, energy=-1)
+        ss1 = dimod.SampleSet.from_samples([+1, -1], dimod.SPIN, energy=-1)
+        ss2 = dimod.SampleSet.from_samples([[+1, +1], [-1, -1]], dimod.SPIN, energy=[1, 1])
+
+        comb = dimod.concatenate((ss0, ss1, ss2))
+
+        out = dimod.SampleSet.from_samples([[-1, +1], [+1, -1], [+1, +1], [-1, -1]], dimod.SPIN, energy=[-1, -1, 1, 1])
+
+        self.assertEqual(comb, out)
+        np.testing.assert_array_equal(comb.record.sample, out.record.sample)
+
+    def test_variables_order(self):
+        ss0 = dimod.SampleSet.from_samples(([-1, +1], 'ab'), dimod.SPIN, energy=-1)
+        ss1 = dimod.SampleSet.from_samples(([-1, +1], 'ba'), dimod.SPIN, energy=-1)
+        ss2 = dimod.SampleSet.from_samples(([[+1, +1], [-1, -1]], 'ab'), dimod.SPIN, energy=[1, 1])
+
+        comb = dimod.concatenate((ss0, ss1, ss2))
+
+        out = dimod.SampleSet.from_samples(([[-1, +1], [+1, -1], [+1, +1], [-1, -1]], 'ab'),
+                                           dimod.SPIN, energy=[-1, -1, 1, 1])
+
+        self.assertEqual(comb, out)
+        np.testing.assert_array_equal(comb.record.sample, out.record.sample)
+
+    def test_variables_order(self):
+        ss0 = dimod.SampleSet.from_samples(([-1, +1], 'ab'), dimod.SPIN, energy=-1)
+        ss1 = dimod.SampleSet.from_samples(([-1, +1], 'ba'), dimod.SPIN, energy=-1)
+        ss2 = dimod.SampleSet.from_samples(([[+1, +1], [-1, -1]], 'ab'), dimod.SPIN, energy=[1, 1])
+
+        comb = dimod.concatenate((ss0, ss1, ss2))
+
+        out = dimod.SampleSet.from_samples(([[-1, +1], [+1, -1], [+1, +1], [-1, -1]], 'ab'),
+                                           dimod.SPIN, energy=[-1, -1, 1, 1])
+
+        self.assertEqual(comb, out)
+        np.testing.assert_array_equal(comb.record.sample, out.record.sample)
+
+    def test_variables_order_and_vartype(self):
+        ss0 = dimod.SampleSet.from_samples(([-1, +1], 'ab'), dimod.SPIN, energy=-1)
+        ss1 = dimod.SampleSet.from_samples(([-1, +1], 'ba'), dimod.SPIN, energy=-1)
+        ss2 = dimod.SampleSet.from_samples(([[1, 1], [0, 0]], 'ab'), dimod.BINARY, energy=[1, 1])
+
+        comb = dimod.concatenate((ss0, ss1, ss2))
+
+        out = dimod.SampleSet.from_samples(([[-1, +1], [+1, -1], [+1, +1], [-1, -1]], 'ab'),
+                                           dimod.SPIN, energy=[-1, -1, 1, 1])
+
+        self.assertEqual(comb, out)
+        np.testing.assert_array_equal(comb.record.sample, out.record.sample)
+
+    def test_empty(self):
+        with self.assertRaises(ValueError):
+            dimod.concatenate([])

pip install -e .[all]
parse_log_pytest
requirements.txt
pytest
apt-get update apt-get install -y libboost-dev
3.6
requirements.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
attrs==22.2.0
certifi==2021.5.30
charset-normalizer==2.0.12
codecov==2.1.13
coverage==6.2
Cython==3.0.12
decorator==5.1.1
-e git+https://github.com/dwavesystems/dimod.git@5c8576958320c247750111e6936958412bad7acb#egg=dimod
idna==3.10
importlib-metadata==4.8.3
iniconfig==1.1.1
jsonschema==2.6.0
mock==2.0.0
networkx==2.0
numpy==1.15.0
packaging==21.3
pandas==0.22.0
pbr==6.1.1
pluggy==1.0.0
py==1.11.0
pymongo==3.7.1
pyparsing==3.1.4
pytest==7.0.1
python-dateutil==2.9.0.post0
pytz==2025.2
requests==2.27.1
six==1.11.0
tomli==1.2.3
typing_extensions==4.1.1
urllib3==1.26.20
zipp==3.6.0

name: dimod
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2021.5.30=py36h06a4308_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.3=he6710b0_2
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - ncurses=6.4=h6a678d5_0
  - openssl=1.1.1w=h7f8727e_0
  - pip=21.2.2=py36h06a4308_0
  - python=3.6.13=h12debd9_1
  - readline=8.2=h5eee18b_0
  - setuptools=58.0.4=py36h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - wheel=0.37.1=pyhd3eb1b0_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - attrs==22.2.0
      - charset-normalizer==2.0.12
      - codecov==2.1.13
      - coverage==6.2
      - cython==3.0.12
      - decorator==5.1.1
      - idna==3.10
      - importlib-metadata==4.8.3
      - iniconfig==1.1.1
      - jsonschema==2.6.0
      - mock==2.0.0
      - networkx==2.0
      - numpy==1.15.0
      - packaging==21.3
      - pandas==0.22.0
      - pbr==6.1.1
      - pluggy==1.0.0
      - py==1.11.0
      - pymongo==3.7.1
      - pyparsing==3.1.4
      - pytest==7.0.1
      - python-dateutil==2.9.0.post0
      - pytz==2025.2
      - requests==2.27.1
      - six==1.11.0
      - tomli==1.2.3
      - typing-extensions==4.1.1
      - urllib3==1.26.20
      - zipp==3.6.0
prefix: /opt/conda/envs/dimod

tests/test_sampleset.py::TestSampleSet::test_aggregate_simple tests/test_sampleset.py::TestSampleSet::test_from_bqm_single_sample tests/test_sampleset.py::TestSampleSet::test_from_samples_with_aggregation tests/test_sampleset.py::Test_concatenate::test_empty tests/test_sampleset.py::Test_concatenate::test_simple tests/test_sampleset.py::Test_concatenate::test_variables_order tests/test_sampleset.py::Test_concatenate::test_variables_order_and_vartype
tests/test_sampleset.py::TestSampleSet::test_eq_ordered tests/test_sampleset.py::TestSampleSet::test_from_samples tests/test_sampleset.py::TestSampleSet::test_from_samples_empty tests/test_sampleset.py::TestSampleSet::test_from_samples_fields_multiple tests/test_sampleset.py::TestSampleSet::test_from_samples_fields_single tests/test_sampleset.py::TestSampleSet::test_from_samples_iterator tests/test_sampleset.py::TestSampleSet::test_from_samples_single_sample tests/test_sampleset.py::TestSampleSet::test_from_samples_str_labels tests/test_sampleset.py::TestSampleSet::test_mismatched_shapes tests/test_sampleset.py::TestSampleSet::test_shorter_samples tests/test_sampleset.py::TestSampleSetSerialization::test_functional_json tests/test_sampleset.py::TestSampleSetSerialization::test_functional_simple_shapes tests/test_sampleset.py::TestSampleSetSerialization::test_functional_str tests/test_sampleset.py::TestSampleSet_to_pandas_dataframe::test_simple
docs/reference/sampleset.rst dimod/sampleset.py dimod/__init__.py
docs/reference/sampleset.rst dimod/sampleset.py dimod/__init__.py