Repository: dwavesystems/dimod

Issue:
extracting whether constraints introduced by make_quadratic are satisfied
under dimod.higher_order, we can add a function that takes the bqm object and the response and checks whether the constraints are satisfied for the samples. Can add that as a field to the solution object.

can have a bool field to remove/keep the extra penalty information after/during this step. 

Already have the function ready to do that - let me know if a PR is wanted.


Discussion:
arcondello: Looks good so far! Most of my comments are just nitpicks. Thanks @conta877 !
codecov-io: # [Codecov](https://codecov.io/gh/dwavesystems/dimod/pull/330?src=pr&el=h1) Report
> Merging [#330](https://codecov.io/gh/dwavesystems/dimod/pull/330?src=pr&el=desc) into [master](https://codecov.io/gh/dwavesystems/dimod/commit/96ef78f53b7762782491bdd635a2c2dd1d6c75d2?src=pr&el=desc) will **decrease** coverage by `2.25%`.
> The diff coverage is `14.08%`.

[![Impacted file tree graph](https://codecov.io/gh/dwavesystems/dimod/pull/330/graphs/tree.svg?width=650&token=KyxipO4fvT&height=150&src=pr)](https://codecov.io/gh/dwavesystems/dimod/pull/330?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master     #330      +/-   ##
==========================================
- Coverage    91.8%   89.55%   -2.26%     
==========================================
  Files          40       41       +1     
  Lines        2344     2413      +69     
==========================================
+ Hits         2152     2161       +9     
- Misses        192      252      +60
```


| [Impacted Files](https://codecov.io/gh/dwavesystems/dimod/pull/330?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [dimod/reference/composites/higherordercomposite.py](https://codecov.io/gh/dwavesystems/dimod/pull/330/diff?src=pr&el=tree#diff-ZGltb2QvcmVmZXJlbmNlL2NvbXBvc2l0ZXMvaGlnaGVyb3JkZXJjb21wb3NpdGUucHk=) | `0% <0%> (ø)` | |
| [dimod/higherorder.py](https://codecov.io/gh/dwavesystems/dimod/pull/330/diff?src=pr&el=tree#diff-ZGltb2QvaGlnaGVyb3JkZXIucHk=) | `63.49% <19.6%> (-30.02%)` | :arrow_down: |
| [dimod/io/json.py](https://codecov.io/gh/dwavesystems/dimod/pull/330/diff?src=pr&el=tree#diff-ZGltb2QvaW8vanNvbi5weQ==) | `96.26% <0%> (+0.93%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/dwavesystems/dimod/pull/330?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/dwavesystems/dimod/pull/330?src=pr&el=footer). Last update [96ef78f...4a926c9](https://codecov.io/gh/dwavesystems/dimod/pull/330?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


Patch:
diff --git a/.gitignore b/.gitignore
index 7c016dc1..bee3c8a5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -103,3 +103,7 @@ ENV/
 
 # cython intermediate files
 dimod/roof_duality/_fix_variables.cpp
+
+# benchmarking stuff
+benchmarks/html/
+benchmarks/results/
diff --git a/benchmarks/asv.conf.json b/benchmarks/asv.conf.json
new file mode 100644
index 00000000..c8c46ed6
--- /dev/null
+++ b/benchmarks/asv.conf.json
@@ -0,0 +1,160 @@
+{
+    // The version of the config file format.  Do not change, unless
+    // you know what you are doing.
+    "version": 1,
+
+    // The name of the project being benchmarked
+    "project": "dimod",
+
+    // The project's homepage
+    "project_url": "https://github.com/dwavesystems/dimod",
+
+    // The URL or local path of the source code repository for the
+    // project being benchmarked
+    "repo": "..",
+
+    // The Python project's subdirectory in your repo.  If missing or
+    // the empty string, the project is assumed to be located at the root
+    // of the repository.
+    // "repo_subdir": "",
+
+    // Customizable commands for building, installing, and
+    // uninstalling the project. See asv.conf.json documentation.
+    //
+    // "install_command": ["python -mpip install {wheel_file}"],
+    // "uninstall_command": ["return-code=any python -mpip uninstall -y {project}"],
+    "build_command": [
+        "python setup.py build_ext --inplace",
+        "PIP_NO_BUILD_ISOLATION=false python -mpip wheel --no-deps --no-index -w {build_cache_dir} {build_dir}"
+    ],
+
+    // List of branches to benchmark. If not provided, defaults to "master"
+    // (for git) or "default" (for mercurial).
+    "branches": ["master"], // for git
+    // "branches": ["default"],    // for mercurial
+
+    // The DVCS being used.  If not set, it will be automatically
+    // determined from "repo" by looking at the protocol in the URL
+    // (if remote), or by looking for special directories, such as
+    // ".git" (if local).
+    "dvcs": "git",
+
+    // The tool to use to create environments.  May be "conda",
+    // "virtualenv" or other value depending on the plugins in use.
+    // If missing or the empty string, the tool will be automatically
+    // determined by looking for tools on the PATH environment
+    // variable.
+    "environment_type": "virtualenv",
+
+    // timeout in seconds for installing any dependencies in environment
+    // defaults to 10 min
+    //"install_timeout": 600,
+
+    // the base URL to show a commit for the project.
+    "show_commit_url": "http://github.com/dwavesystems/dimod/commit/",
+
+    // The Pythons you'd like to test against.  If not provided, defaults
+    // to the current version of Python used to run `asv`.
+    "pythons": ["3.7"],
+
+    // The list of conda channel names to be searched for benchmark
+    // dependency packages in the specified order
+    // "conda_channels": ["conda-forge", "defaults"]
+
+    // The matrix of dependencies to test.  Each key is the name of a
+    // package (in PyPI) and the values are version numbers.  An empty
+    // list or empty string indicates to just test against the default
+    // (latest) version. null indicates that the package is to not be
+    // installed. If the package to be tested is only available from
+    // PyPi, and the 'environment_type' is conda, then you can preface
+    // the package name by 'pip+', and the package will be installed via
+    // pip (with all the conda available packages installed first,
+    // followed by the pip installed packages).
+    //
+    // "matrix": {
+    //     "numpy": ["1.6", "1.7"],
+    //     "six": ["", null],        // test with and without six installed
+    //     "pip+emcee": [""],   // emcee is only available for install with pip.
+    // },
+
+    // Combinations of libraries/python versions can be excluded/included
+    // from the set to test. Each entry is a dictionary containing additional
+    // key-value pairs to include/exclude.
+    //
+    // An exclude entry excludes entries where all values match. The
+    // values are regexps that should match the whole string.
+    //
+    // An include entry adds an environment. Only the packages listed
+    // are installed. The 'python' key is required. The exclude rules
+    // do not apply to includes.
+    //
+    // In addition to package names, the following keys are available:
+    //
+    // - python
+    //     Python version, as in the *pythons* variable above.
+    // - environment_type
+    //     Environment type, as above.
+    // - sys_platform
+    //     Platform, as in sys.platform. Possible values for the common
+    //     cases: 'linux2', 'win32', 'cygwin', 'darwin'.
+    //
+    // "exclude": [
+    //     {"python": "3.2", "sys_platform": "win32"}, // skip py3.2 on windows
+    //     {"environment_type": "conda", "six": null}, // don't run without six on conda
+    // ],
+    //
+    // "include": [
+    //     // additional env for python2.7
+    //     {"python": "2.7", "numpy": "1.8"},
+    //     // additional env if run on windows+conda
+    //     {"platform": "win32", "environment_type": "conda", "python": "2.7", "libpython": ""},
+    // ],
+
+    // The directory (relative to the current directory) that benchmarks are
+    // stored in.  If not provided, defaults to "benchmarks"
+    "benchmark_dir": "benchmarks",
+
+    // The directory (relative to the current directory) to cache the Python
+    // environments in.  If not provided, defaults to "env"
+    "env_dir": "env",
+
+    // The directory (relative to the current directory) that raw benchmark
+    // results are stored in.  If not provided, defaults to "results".
+    "results_dir": "results",
+
+    // The directory (relative to the current directory) that the html tree
+    // should be written to.  If not provided, defaults to "html".
+    "html_dir": "html",
+
+    // The number of characters to retain in the commit hashes.
+    "hash_length": 8,
+
+    // `asv` will cache results of the recent builds in each
+    // environment, making them faster to install next time.  This is
+    // the number of builds to keep, per environment.
+    // "build_cache_size": 2,
+
+    // The commits after which the regression search in `asv publish`
+    // should start looking for regressions. Dictionary whose keys are
+    // regexps matching to benchmark names, and values corresponding to
+    // the commit (exclusive) after which to start looking for
+    // regressions.  The default is to start from the first commit
+    // with results. If the commit is `null`, regression detection is
+    // skipped for the matching benchmark.
+    //
+    // "regressions_first_commits": {
+    //    "some_benchmark": "352cdf",  // Consider regressions only after this commit
+    //    "another_benchmark": null,   // Skip regression detection altogether
+    // },
+
+    // The thresholds for relative change in results, after which `asv
+    // publish` starts reporting regressions. Dictionary of the same
+    // form as in ``regressions_first_commits``, with values
+    // indicating the thresholds.  If multiple entries match, the
+    // maximum is taken. If no entry matches, the default is 5%.
+    //
+    // "regressions_thresholds": {
+    //    "some_benchmark": 0.01,     // Threshold of 1%
+    //    "another_benchmark": 0.5,   // Threshold of 50%
+    // },
+}
diff --git a/benchmarks/benchmarks/__init__.py b/benchmarks/benchmarks/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/benchmarks/benchmarks/bench_as_samples.py b/benchmarks/benchmarks/bench_as_samples.py
new file mode 100644
index 00000000..99d19bfc
--- /dev/null
+++ b/benchmarks/benchmarks/bench_as_samples.py
@@ -0,0 +1,82 @@
+# Copyright 2019 D-Wave Systems Inc.
+#
+#    Licensed under the Apache License, Version 2.0 (the "License");
+#    you may not use this file except in compliance with the License.
+#    You may obtain a copy of the License at
+#
+#        http://www.apache.org/licenses/LICENSE-2.0
+#
+#    Unless required by applicable law or agreed to in writing, software
+#    distributed under the License is distributed on an "AS IS" BASIS,
+#    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#    See the License for the specific language governing permissions and
+#    limitations under the License.
+#
+# =============================================================================
+import itertools
+
+import numpy as np
+import dimod
+
+
+from benchmarks.common import Benchmark
+
+
+class NumpyArray(Benchmark):
+    def setup(self):
+        self.empty = np.asarray([], dtype=np.int8)
+        self.arr100x100 = np.ones((100, 100), dtype=np.int8)
+        self.labels100 = list(range(100))
+
+    def time_empty(self):
+        dimod.as_samples(self.empty)
+
+    def time_100x100(self):
+        dimod.as_samples(self.arr100x100)
+
+    def time_empty_labelled(self):
+        dimod.as_samples((self.empty, []))
+
+    def time_100x100_labelled(self):
+        dimod.as_samples((self.arr100x100, self.labels100))
+
+
+class Dicts(Benchmark):
+    def setup(self):
+        # 1000 variables
+        variables = list(''.join(letters) for letters in itertools.product('abcdefghij', repeat=3))
+
+        self.dict1000 = {v: 1 for v in variables}
+        self.dict1x1000 = [{v: 1 for v in variables}]
+        self.dict1000x1000 = [{v: 1 for v in variables} for _ in range(1000)]
+
+    def time_empty(self):
+        dimod.as_samples({})
+
+    def time_dict1000(self):
+        dimod.as_samples(self.dict1000)
+
+    def time_dict1x1000(self):
+        dimod.as_samples(self.dict1x1000)
+
+    def time_dict1000x1000(self):
+        dimod.as_samples(self.dict1000x1000)
+
+
+class Lists(Benchmark):
+    def setup(self):
+        self.list1000 = [1]*1000
+        self.list1x1000 = [[1]*1000]
+        self.list1000x1000 = [[1]*1000 for _ in range(1000)]
+
+    def time_empty(self):
+        dimod.as_samples([])
+
+    def time_list1000(self):
+        dimod.as_samples(self.list1000)
+
+    def time_list1x1000(self):
+        dimod.as_samples(self.list1x1000)
+
+    def time_list1000x1000(self):
+        dimod.as_samples(self.list1000x1000)
diff --git a/benchmarks/benchmarks/bench_binary_quadratic_model.py b/benchmarks/benchmarks/bench_binary_quadratic_model.py
new file mode 100644
index 00000000..7602df69
--- /dev/null
+++ b/benchmarks/benchmarks/bench_binary_quadratic_model.py
@@ -0,0 +1,111 @@
+# Copyright 2019 D-Wave Systems Inc.
+#
+#    Licensed under the Apache License, Version 2.0 (the "License");
+#    you may not use this file except in compliance with the License.
+#    You may obtain a copy of the License at
+#
+#        http://www.apache.org/licenses/LICENSE-2.0
+#
+#    Unless required by applicable law or agreed to in writing, software
+#    distributed under the License is distributed on an "AS IS" BASIS,
+#    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#    See the License for the specific language governing permissions and
+#    limitations under the License.
+#
+# ================================================================================================
+import itertools
+import copy
+
+import numpy as np
+import dimod
+
+from benchmarks.common import Benchmark
+
+
+class Construction(Benchmark):
+    def setup(self):
+        self.h_K100 = {v: 0 for v in range(100)}
+        self.J_K100 = {edge: 0 for edge in itertools.combinations(range(100), 2)}
+
+    def time_k100(self):
+        dimod.BinaryQuadraticModel(self.h_K100, self.J_K100, 0.0, dimod.SPIN)
+
+    def mem_k100(self):
+        return dimod.BinaryQuadraticModel(self.h_K100, self.J_K100, 0.0, dimod.SPIN)
+
+    def time_empty(self):
+        dimod.BinaryQuadraticModel({}, {}, 0.0, dimod.SPIN)
+
+    def mem_empty(self):
+        return dimod.BinaryQuadraticModel({}, {}, 0.0, dimod.SPIN)
+
+
+class ConstructionVeryLarge(Benchmark):
+    def setup(self):
+        self.J_15000 = J = {}
+
+        # make a graph of approximately degree 15, 4394 nodes, 32617 edges
+        m = n = 13
+        t = 13
+        J.update((((i, j, 0, k0), (i, j, 1, k1)), 1)
+                 for i in range(n)
+                 for j in range(m)
+                 for k0 in range(t)
+                 for k1 in range(t))
+        J.update((((i, j, 1, k), (i, j+1, 1, k)), 1)
+                 for i in range(m)
+                 for j in range(n-1)
+                 for k in range(t))
+        J.update((((i, j, 0, k), (i+1, j, 0, k)), 1)
+                 for i in range(m-1)
+                 for j in range(n)
+                 for k in range(t))
+
+    def time_J15000(self):
+        dimod.BinaryQuadraticModel({}, self.J_15000, 0.0, dimod.SPIN)
+
+    def mem_J15000(self):
+        return dimod.BinaryQuadraticModel({}, self.J_15000, 0.0, dimod.SPIN)
+
+
+class Copy(Benchmark):
+    def setup(self):
+        h_K100 = {v: 0 for v in range(100)}
+        J_K100 = {edge: 0 for edge in itertools.combinations(range(100), 2)}
+        self.bqm_k100 = dimod.BinaryQuadraticModel(h_K100, J_K100, 0.0, dimod.SPIN)
+
+        self.bqm_empty = dimod.BinaryQuadraticModel({}, {}, 0.0, dimod.SPIN)
+
+    def time_k100(self):
+        self.bqm_k100.copy()
+
+    def time_empty(self):
+        self.bqm_empty.copy()
+
+
+class Deepcopy(Benchmark):
+    """How long it takes to do a deep copy"""
+    def setup(self):
+        h_K100 = {v: 0 for v in range(100)}
+        J_K100 = {edge: 0 for edge in itertools.combinations(range(100), 2)}
+        self.bqm_k100 = dimod.BinaryQuadraticModel(h_K100, J_K100, 0.0, dimod.SPIN)
+
+        self.bqm_empty = dimod.BinaryQuadraticModel({}, {}, 0.0, dimod.SPIN)
+
+    def time_k100(self):
+        copy.deepcopy(self.bqm_k100)
+
+    def time_empty(self):
+        copy.deepcopy(self.bqm_empty)
+
+
+class Energies(Benchmark):
+    def setup(self):
+        h_K100 = {v: 0 for v in range(100)}
+        J_K100 = {edge: 0 for edge in itertools.combinations(range(100), 2)}
+        self.bqm_k100 = dimod.BinaryQuadraticModel(h_K100, J_K100, 0.0, dimod.SPIN)
+
+        self.samples_1000x100 = np.ones((1000, 100), dtype=np.int8)
+
+    def time_k100(self):
+        self.bqm_k100.energies(self.samples_1000x100)
diff --git a/benchmarks/benchmarks/common.py b/benchmarks/benchmarks/common.py
new file mode 100644
index 00000000..073497d4
--- /dev/null
+++ b/benchmarks/benchmarks/common.py
@@ -0,0 +1,19 @@
+# Copyright 2019 D-Wave Systems Inc.
+#
+#    Licensed under the Apache License, Version 2.0 (the "License");
+#    you may not use this file except in compliance with the License.
+#    You may obtain a copy of the License at
+#
+#        http://www.apache.org/licenses/LICENSE-2.0
+#
+#    Unless required by applicable law or agreed to in writing, software
+#    distributed under the License is distributed on an "AS IS" BASIS,
+#    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#    See the License for the specific language governing permissions and
+#    limitations under the License.
+#
+# =============================================================================
+
+
+class Benchmark(object):
+    repeat = (1, 1000, .5)
diff --git a/benchmarks/requirements.txt b/benchmarks/requirements.txt
new file mode 100644
index 00000000..b3e29585
--- /dev/null
+++ b/benchmarks/requirements.txt
@@ -0,0 +1,4 @@
+-r ../requirements.txt
+
+asv
+virtualenv
diff --git a/dimod/__init__.py b/dimod/__init__.py
index e90ffac4..aeaa857a 100644
--- a/dimod/__init__.py
+++ b/dimod/__init__.py
@@ -43,12 +43,12 @@ import dimod.generators
 from dimod.exceptions import *
 import dimod.exceptions
 
-from dimod.higherorder import make_quadratic, poly_energy
+from dimod.higherorder import make_quadratic, poly_energy, poly_energies
 import dimod.higherorder
 
 from dimod.package_info import __version__, __author__, __authoremail__, __description__
 
-from dimod.sampleset import as_samples, SampleSet, concatenate
+from dimod.sampleset import as_samples, concatenate, SampleSet
 
 from dimod.response import *
 import dimod.response
diff --git a/dimod/binary_quadratic_model.py b/dimod/binary_quadratic_model.py
index ba05ea97..fa6a7a33 100644
--- a/dimod/binary_quadratic_model.py
+++ b/dimod/binary_quadratic_model.py
@@ -427,21 +427,22 @@ class BinaryQuadraticModel(abc.Sized, abc.Container, abc.Iterable):
         if vartype is not None and vartype is not self.vartype:
             if self.vartype is Vartype.SPIN and vartype is Vartype.BINARY:
                 # convert from binary to spin
-                bias /= 2.
+                bias /= 2
                 self.offset += bias
             elif self.vartype is Vartype.BINARY and vartype is Vartype.SPIN:
                 # convert from spin to binary
                 self.offset -= bias
-                bias *= 2.
+                bias *= 2
             else:
                 raise ValueError("unknown vartype")
 
-        # add the variable to linear and adj
-        linear = self.linear
-        if v in linear:
-            linear[v] += bias
+        # we used to do this using self.linear but working directly with _adj
+        # is much faster
+        _adj = self._adj
+        if v in _adj:
+            _adj[v][v] += bias
         else:
-            linear[v] = bias
+            _adj[v] = {v: bias}
 
         try:
             self._counterpart.add_variable(v, bias, vartype=self.vartype)
@@ -539,38 +540,46 @@ class BinaryQuadraticModel(abc.Sized, abc.Container, abc.Iterable):
         if u == v:
             raise ValueError("no self-loops allowed, therefore ({}, {}) is not an allowed interaction".format(u, v))
 
-        linear = self.linear
-        quadratic = self.quadratic
-        adj = self.adj
-
-        # so that they exist
-        self.add_variable(u, 0.0)
-        self.add_variable(v, 0.0)
-
         if vartype is not None and vartype is not self.vartype:
             if self.vartype is Vartype.SPIN and vartype is Vartype.BINARY:
                 # convert from binary to spin
                 bias /= 4
 
-                self.offset += bias
-                linear[u] += bias
-                linear[v] += bias
+                self.add_offset(bias)
+                self.add_variable(u, bias)
+                self.add_variable(v, bias)
 
             elif self.vartype is Vartype.BINARY and vartype is Vartype.SPIN:
                 # convert from spin to binary
 
-                self.offset += bias
-                linear[u] += -2 * bias
-                linear[v] += -2 * bias
+                self.add_offset(bias)
+                self.add_variable(u, -2 * bias)
+                self.add_variable(v, -2 * bias)
 
                 bias *= 4
             else:
                 raise ValueError("unknown vartype")
-
-        if (v, u) in quadratic:
-            quadratic[(v, u)] += bias
         else:
-            quadratic[(u, v)] = bias
+            # so that they exist.
+
+            # developer note: we could  try to match the type of the quadratic
+            # bias with type(bias)(0), but this sometimes fails. We could also
+            # put it in a try-except but I am not so sure that it will fail
+            # predictably. I think it is better to just use python int 0 as it
+            # is most likely to be compatible with other numeric types.
+
+            # We don't actually need to do the check whether is is present but
+            # it is slightly faster to do so.
+            if u not in self:
+                self.add_variable(u, 0)
+            if v not in self:
+                self.add_variable(v, 0)
+
+        _adj = self._adj
+        if u in _adj[v]:
+            _adj[u][v] = _adj[v][u] = _adj[u][v] + bias
+        else:
+            _adj[u][v] = _adj[v][u] = bias
 
         try:
             self._counterpart.add_interaction(u, v, bias, vartype=self.vartype)
@@ -1809,7 +1818,8 @@ class BinaryQuadraticModel(abc.Sized, abc.Container, abc.Iterable):
             ({0: 1, 1: -1, 2: 0.5}, {(0, 1): 0.5, (1, 2): 1.5}, 1.4)
 
         """
-        return self.spin.linear, self.spin.quadratic, self.spin.offset
+        # cast to a dict
+        return dict(self.spin.linear), dict(self.spin.quadratic), self.spin.offset
 
     @classmethod
     def from_ising(cls, h, J, offset=0.0):
@@ -1844,7 +1854,7 @@ class BinaryQuadraticModel(abc.Sized, abc.Container, abc.Iterable):
             BinaryQuadraticModel({1: 1, 2: 2, 3: 3, 4: 4}, {(1, 2): 12, (1, 3): 13, (1, 4): 14, (2, 3): 23, (3, 4): 34, (2, 4): 24}, 0.0, Vartype.SPIN)
 
         """
-        if isinstance(h, list):
+        if not isinstance(h, abc.Mapping):
             h = dict(enumerate(h))
 
         return cls(h, J, offset, Vartype.SPIN)
@@ -1874,14 +1884,8 @@ class BinaryQuadraticModel(abc.Sized, abc.Container, abc.Iterable):
             ({(0, 0): 1.0, (0, 1): 2.0, (1, 1): -6.0, (1, 2): 6.0, (2, 2): -2.0}, 2.9)
 
         """
-        qubo = {}
-
-        for v, bias in iteritems(self.binary.linear):
-            qubo[(v, v)] = bias
-
-        for edge, bias in iteritems(self.binary.quadratic):
-            qubo[edge] = bias
-
+        qubo = dict(self.binary.quadratic)
+        qubo.update(((v, v), bias) for v, bias in iteritems(self.binary.linear))
         return qubo, self.binary.offset
 
     @classmethod
diff --git a/dimod/decorators.py b/dimod/decorators.py
index 12fe1248..21aee399 100644
--- a/dimod/decorators.py
+++ b/dimod/decorators.py
@@ -17,7 +17,11 @@
 import inspect
 import itertools
 
-from collections import Sequence
+try:
+    import collections.abc as abc
+except ImportError:
+    import collections as abc
+
 from functools import wraps
 from numbers import Integral
 
@@ -266,7 +270,7 @@ def graph_argument(*arg_names):
                 # an integer, cast to a complete graph
                 kwargs[name] = (list(range(G)), list(itertools.combinations(range(G), 2)))
 
-            elif isinstance(G, Sequence) and len(G) == 2:
+            elif isinstance(G, abc.Sequence) and len(G) == 2:
                 # is a pair nodes/edges
                 if isinstance(G[0], integer_types):
                     # if nodes is an int
diff --git a/dimod/higherorder.py b/dimod/higherorder.py
index 24ff20c2..2f1121f3 100644
--- a/dimod/higherorder.py
+++ b/dimod/higherorder.py
@@ -15,13 +15,13 @@
 # ================================================================================================
 
 import itertools
-
 from collections import Counter
 
+import numpy as np
 from six import iteritems
 
 from dimod.binary_quadratic_model import BinaryQuadraticModel
-from dimod.decorators import vartype_argument
+from dimod.sampleset import as_samples
 from dimod.vartypes import Vartype
 
 __all__ = ['make_quadratic']
@@ -104,7 +104,7 @@ def make_quadratic(poly, strength, vartype=None, bqm=None):
     Examples:
 
         >>> poly = {(0,): -1, (1,): 1, (2,): 1.5, (0, 1): -1, (0, 1, 2): -2}
-        >>> bqm = make_quadratic(poly, 5.0, dimod.SPIN)
+        >>> bqm = dimod.make_quadratic(poly, 5.0, dimod.SPIN)
 
     """
 
@@ -133,6 +133,7 @@ def make_quadratic(poly, strength, vartype=None, bqm=None):
 
 
 def _reduce_degree(bqm, poly, vartype, scale):
+    """helper function for make_quadratic"""
 
     if all(len(term) <= 2 for term in poly):
         # termination criteria, we are already quadratic
@@ -193,6 +194,24 @@ def _reduce_degree(bqm, poly, vartype, scale):
     return _reduce_degree(bqm, new_poly, vartype, scale)
 
 
+def create_poly(linear, higherorder):
+    """Creates a polynomial dict
+
+    Args:
+        linear (dict): linear variables
+        higherorder (dict): quadratic and higher order variables
+
+    Returns
+        dict: a higher order problem dict that contains linear,
+              quadratic and n-order terms.
+
+    """
+
+    poly = {(k,): v for k, v in linear.items()}
+    poly.update(higherorder)
+    return poly
+
+
 def _prod(iterable):
     val = 1.
     for v in iterable:
@@ -200,21 +219,64 @@ def _prod(iterable):
     return val
 
 
-def poly_energy(sample, poly):
-    """Create a binary quadratic model from a higher order polynomial.
+def _prod_d(iterable, dim):
+    val = [1.] * dim
+    for v in iterable:
+        val *= v
+    return val
+
+
+def poly_energy(sample_like, poly):
+    """Calculates energy of a sample from a higher order polynomial.
 
     Args:
-        sample (dict):
-            Sample for which to calculate the energy, formatted as a dict where keys
-            are variables and values are the value associated with each variable.
+         sample (samples_like):
+            A raw sample. `samples_like` is an extension of NumPy's
+            array_like structure. See :func:`.as_samples`.
 
         poly (dict):
-            Polynomial as a dict of form {term: bias, ...}, where `term` is a tuple of
-            variables and `bias` the associated bias.
+            Polynomial as a dict of form {term: bias, ...}, where `term` is a
+            tuple of variables and `bias` the associated bias.
 
     Returns:
         float: The energy of the sample.
 
     """
-    return sum(_prod(sample[v] for v in variables) * bias
-               for variables, bias in poly.items())
+    sample, labels = as_samples(sample_like)
+    idx, label = zip(*enumerate(labels))
+    labeldict = dict(zip(label, idx))
+    num_samples = len(sample)
+
+    if num_samples > 1:
+        raise ValueError('poly_energy accepts a single sample. For multiple '
+                         'samples use poly_energies')
+    else:
+        return sum(_prod(sample[0][labeldict[v]] for v in variables) * bias
+                   for variables, bias in poly.items())
+
+
+def poly_energies(samples_like, poly):
+    """Calculates energy of samples from a higher order polynomial.
+
+    Args:
+        sample (samples_like):
+            A collection of raw samples. `samples_like` is an extension of
+            NumPy's array_like structure. See :func:`.as_samples`.
+
+        poly (dict):
+            Polynomial as a dict of form {term: bias, ...}, where `term` is a 
+            tuple of variables and `bias` the associated bias. Variable
+            labeling/indexing of terms in poly dict must match that of the
+            sample(s).
+
+    Returns:
+        list/:obj:`numpy.ndarray`: The energy of the sample(s).
+
+    """
+    sample, labels = as_samples(samples_like)
+    idx, label = zip(*enumerate(labels))
+    labeldict = dict(zip(label, idx))
+
+    num_samples = len(sample)
+    return sum(_prod_d([sample[:, labeldict[v]] for v in variables],
+                       num_samples) * bias for variables, bias in poly.items())
diff --git a/dimod/package_info.py b/dimod/package_info.py
index ad0734e3..27360166 100644
--- a/dimod/package_info.py
+++ b/dimod/package_info.py
@@ -14,7 +14,7 @@
 #
 # ================================================================================================
 
-__version__ = '0.8.1'
+__version__ = '0.8.2'
 __author__ = 'D-Wave Systems Inc.'
 __authoremail__ = 'acondello@dwavesys.com'
 __description__ = 'A shared API for binary quadratic model samplers.'
diff --git a/dimod/reference/composites/__init__.py b/dimod/reference/composites/__init__.py
index 1d1f4794..c6e934cd 100644
--- a/dimod/reference/composites/__init__.py
+++ b/dimod/reference/composites/__init__.py
@@ -17,3 +17,4 @@
 from dimod.reference.composites.spin_transform import *
 from dimod.reference.composites.structure import *
 from dimod.reference.composites.fixedvariable import *
+from dimod.reference.composites.higherordercomposite import *
diff --git a/dimod/reference/composites/higherordercomposite.py b/dimod/reference/composites/higherordercomposite.py
new file mode 100644
index 00000000..19a68d66
--- /dev/null
+++ b/dimod/reference/composites/higherordercomposite.py
@@ -0,0 +1,261 @@
+# Copyright 2018 D-Wave Systems Inc.
+#
+#    Licensed under the Apache License, Version 2.0 (the "License");
+#    you may not use this file except in compliance with the License.
+#    You may obtain a copy of the License at
+#
+#        http://www.apache.org/licenses/LICENSE-2.0
+#
+#    Unless required by applicable law or agreed to in writing, software
+#    distributed under the License is distributed on an "AS IS" BASIS,
+#    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#    See the License for the specific language governing permissions and
+#    limitations under the License.
+#
+# ================================================================================================
+
+"""
+A composite that convert a higher order polynomial problem into a bqm by
+introducing penalties before sending the bqm to its child sampler.
+"""
+import numpy as np
+
+from dimod.binary_quadratic_model import BinaryQuadraticModel
+from dimod.core.composite import ComposedSampler
+from dimod.higherorder import make_quadratic, poly_energies
+from dimod.response import SampleSet
+
+__all__ = ['HigherOrderComposite']
+
+
+class HigherOrderComposite(ComposedSampler):
+    """Reduces a HUBO to bqm by introducing penalties.
+
+    Energies of the returned samples do not include the penalties.
+
+   Args:
+       sampler (:obj:`dimod.Sampler`):
+            A dimod sampler
+
+   Example:
+       This example uses :class:`.HigherOrderComposite` to instantiate a
+       composed sampler that submits a simple Ising problem to a sampler.
+       The composed sampler creates a bqm from a higher order problem.
+
+       >>> sampler = dimod.HigherOrderComposite(dimod.ExactSolver())
+       >>> linear = {0: -0.5, 1: -0.3, 2: -0.8}
+       >>> quadratic = {(0, 1, 2): -1.7}
+       >>> response = sampler.sample_ising(linear,quadratic,penalty_strength=2,
+       ...                              keep_penalty_variables=False,
+       ...                              discard_unsatisfied=True)
+       >>> print(response.first)  # doctest: +SKIP
+       Sample(sample={0: 1, 1: 1, 2: 1}, energy=-3.3, num_occurrences=1,
+              penalty_satisfaction=True)
+
+        """
+
+    def __init__(self, child_sampler):
+        self._children = [child_sampler]
+
+    @property
+    def children(self):
+        return self._children
+
+    @property
+    def parameters(self):
+        param = self.child.parameters.copy()
+        param['penalty_strength'] = []
+        return param
+
+    @property
+    def properties(self):
+        return {'child_properties': self.child.properties.copy()}
+
+    def sample_ising(self, h, J, offset=0, penalty_strength=1.0,
+                     keep_penalty_variables=False,
+                     discard_unsatisfied=False, **parameters):
+        """ Sample from the problem provided by h, J, offset.
+
+        Takes in linear variables in h and quadratic and higher order
+        terms in J. Introducing penalties, reduces the higher-order problem
+        into a quadratic problem and send it to its child sampler.
+
+        Args:
+            h (dict): linear biases corresponding to the HUBO form
+
+            J (dict): higher order biases corresponding to the HUBO form
+
+            offset (float, optional): constant energy offset
+
+            penalty_strength (float, optional): Strength of the reduction constraint.
+                Insufficient strength can result in the binary quadratic model
+                not having the same minimization as the polynomial.
+
+            keep_penalty_variables (bool, optional): default is True. if False
+                will remove the variables used for penalty from the samples
+
+            discard_unsatisfied (bool, optional): default is False. If True
+                will discard samples that do not satisfy the penalty conditions.
+
+            **parameters: Parameters for the sampling method, specified by
+            the child sampler.
+
+        Returns:
+            :obj:`dimod.SampleSet`
+        """
+
+        # solve the problem on the child system
+
+        bqm = BinaryQuadraticModel.from_ising(h, {})
+        bqm = make_quadratic(J, penalty_strength, bqm=bqm)
+        response = self.child.sample(bqm, **parameters)
+
+        return polymorph_response(response, h, J, bqm,
+                                  offset=offset,
+                                  penalty_strength=penalty_strength,
+                                  keep_penalty_variables=keep_penalty_variables,
+                                  discard_unsatisfied=discard_unsatisfied)
+
+
+def penalty_satisfaction(response, bqm):
+    """ Creates a penalty satisfaction list
+
+    Given a sampleSet and a bqm object, will create a binary list informing
+    whether the penalties introduced during degree reduction are satisfied for
+    each sample in sampleSet
+
+    Args:
+        response (:obj:`.SampleSet`): Samples corresponding to provided bqm
+
+        bqm (:obj:`.BinaryQuadraticModel`): a bqm object that contains
+            its reduction info.
+
+    Returns:
+        :obj:`numpy.ndarray`: a binary array of penalty satisfaction information
+
+    """
+    record = response.record
+    label_dict = response.variables.index
+
+    if len(bqm.info['reduction']) == 0:
+        return np.array([1] * len(record.sample))
+
+    penalty_vector = np.prod([record.sample[:, label_dict[qi]] *
+                              record.sample[:, label_dict[qj]]
+                              == record.sample[:,
+                                 label_dict[valdict['product']]]
+                              for (qi, qj), valdict in
+                              bqm.info['reduction'].items()], axis=0)
+    return penalty_vector
+
+
+def _relabeled_poly(h, j, label_dict):
+    """ Creates relabeled polynomial dict.
+
+    given h, j and a relabeling dictionary, will create a polynomial
+    with the new labels
+
+    Args:
+        h (dict): a dict of linear variables
+
+        j (dict): a dict of quadratic and higher order variables
+
+        label_dict (dict): a dict for relabeling e.g. {old_label:new_label}
+
+    Returns:
+        dict: a higher order problem dict that contains linear,
+              quadratic and n-order terms written in a new labeling
+              scheme provided by label_dict
+
+    """
+
+    poly = {}
+    for k, v in h.items():
+        new_tup = (label_dict[k],)
+        poly[new_tup] = v
+    for k, v in j.items():
+        new_tup = tuple((label_dict[vidx] for vidx in list(k)))
+        poly[new_tup] = v
+    return poly
+
+
+def polymorph_response(response, h, J, bqm, offset=0,
+                       penalty_strength=None,
+                       keep_penalty_variables=True,
+                       discard_unsatisfied=False):
+    """ Transforms the sampleset for the higher order problem.
+
+    Given a response of a penalized HUBO, this function creates a new sampleset
+    object, taking into account penalty information and calculates the
+    energies of samples for the higherorder problem.
+
+    Args:
+        response (:obj:`.SampleSet`): response for a penalized hubo.
+
+        h (dict): linear biases corresponding to the HUBO form
+
+        J (dict): higher order biases corresponding to the HUBO form
+
+        bqm (:obj:`dimod.BinaryQuadraticModel`): Binary quadratic model of the
+            reduced problem.
+
+        offset (float, optional): constant energy offset
+
+        penalty_strength (float, optional): default is None, if provided,
+            will be added to the info field of the returned sampleSet object.
+
+        keep_penalty_variables (bool, optional): default is True. if False
+            will remove the variables used for penalty from the samples
+
+        discard_unsatisfied (bool, optional): default is False. If True
+            will discard samples that do not satisfy the penalty conditions.
+
+    Returns:
+        (:obj:`.SampleSet'): A sampleSet object that has additional penalty
+            information. The energies of samples are calculated for the HUBO
+            ignoring the penalty variables.
+
+    """
+    record = response.record
+    penalty_vector = penalty_satisfaction(response, bqm)
+    original_variables = bqm.variables
+
+    if discard_unsatisfied:
+        samples_to_keep = list(map(bool, list(penalty_vector)))
+        penalty_vector = np.array([True] * np.sum(samples_to_keep))
+    else:
+        samples_to_keep = list(map(bool, [1] * len(record.sample)))
+
+    poly = _relabeled_poly(h, J, response.variables.index)
+    samples = record.sample[samples_to_keep]
+    energy_vector = np.add(poly_energies(samples, poly), offset)
+
+    if not keep_penalty_variables:
+        original_variables = set(h).union(*J)
+        idxs = [response.variables.index[v] for v in original_variables]
+        samples = np.asarray(samples[:, idxs])
+
+    num_samples, num_variables = np.shape(samples)
+
+    datatypes = [('sample', np.dtype(np.int8), (num_variables,)),
+                 ('energy', energy_vector.dtype),
+                 ('penalty_satisfaction',
+                  penalty_vector.dtype)]
+    datatypes.extend((name, record[name].dtype, record[name].shape[1:])
+                     for name in record.dtype.names if
+                     name not in {'sample',
+                                  'energy'})
+
+    data = np.rec.array(np.empty(num_samples, dtype=datatypes))
+    data.sample = samples
+    data.energy = energy_vector
+    for name in record.dtype.names:
+        if name not in {'sample', 'energy'}:
+            data[name] = record[name][samples_to_keep]
+
+    data['penalty_satisfaction'] = penalty_vector
+    response.info['reduction'] = bqm.info['reduction']
+    if penalty_strength is not None:
+        response.info['penalty_strength'] = penalty_strength
+    return SampleSet(data, original_variables, response.info,
+                     response.vartype)
diff --git a/dimod/response.py b/dimod/response.py
index 805db27e..1f0aa936 100644
--- a/dimod/response.py
+++ b/dimod/response.py
@@ -17,7 +17,6 @@ from __future__ import division
 
 import itertools
 
-from collections import ItemsView, Iterable, Mapping, Sized, ValuesView
 from collections import namedtuple
 
 import numpy as np
diff --git a/dimod/sampleset.py b/dimod/sampleset.py
index b6f5a4d7..d19f6142 100644
--- a/dimod/sampleset.py
+++ b/dimod/sampleset.py
@@ -13,16 +13,14 @@
 #    limitations under the License.
 #
 # ================================================================================================
-"""
-Samples can be provided in many forms, generally referred to as samples_like, an extension
-of NumPy's array_like_ to allow for arbitrary variable labels.
-
-.. _array_like: https://docs.scipy.org/doc/numpy/user/basics.creation.html#converting-python-array-like-objects-to-numpy-arrays
-
-"""
 import itertools
+import numbers
+
+try:
+    import collections.abc as abc
+except ImportError:
+    import collections as abc
 
-from collections import Iterable, Sized, Mapping, Iterator
 from collections import namedtuple
 
 import numpy as np
@@ -30,32 +28,39 @@ 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
 from dimod.views import SampleView
 from dimod.variables import Variables
 
-__all__ = 'SampleSet', 'as_samples'
+__all__ = 'as_samples', 'concatenate', 'SampleSet'
 
 
-def as_samples(samples_like, dtype=None):
+def as_samples(samples_like, dtype=None, copy=False, order='C'):
     """Convert a samples_like object to a NumPy array and list of labels.
 
     Args:
         samples_like (samples_like):
-            A collection of raw samples. `samples_like` is an extension of NumPy's array_like_
-            structure. See examples below.
+            A collection of raw samples. `samples_like` is an extension of
+            NumPy's array_like_ structure. See examples below.
 
         dtype (data-type, optional):
-            dtype for the returned samples array. If not provided, it is either derived from
-            `samples_like`, if that object has a dtype, or set to :class:`numpy.int8`.
+            dtype for the returned samples array. If not provided, it is either
+            derived from `samples_like`, if that object has a dtype, or set to
+            :class:`numpy.int8`.
+
+        copy (bool, optional, default=False):
+            If true, then samples_like is guaranteed to be copied, otherwise
+            it is only copied if necessary.
+
+        order ({'K', 'A', 'C', 'F'}, optional, default='C'):
+            Specify the memory layout of the array. See :func:`numpy.array`.
 
     Returns:
         tuple: A 2-tuple containing:
 
             :obj:`numpy.ndarray`: Samples.
 
-            list: Variable labels as a list.
+            list: Variable labels
 
     Examples:
         The following examples convert a variety of samples_like objects:
@@ -99,54 +104,74 @@ def as_samples(samples_like, dtype=None):
                 [0, 0],
                 [0, 0]], dtype=int8), ['in', 'out'])
 
-    .. _array_like:  https://docs.scipy.org/doc/numpy/user/basics.creation.html#converting-python-array-like-objects-to-numpy-arrays
+    .. _array_like: https://docs.scipy.org/doc/numpy/user/basics.creation.html
 
     """
-
-    if isinstance(samples_like, Iterator):
-        raise TypeError('samples_like cannot be an iterator')
-
     if isinstance(samples_like, tuple) and len(samples_like) == 2:
-        # (samples_like, labels)
         samples_like, labels = samples_like
 
-        samples, __ = as_samples(samples_like)
+        if not isinstance(labels, list) and labels is not None:
+            labels = list(labels)
+    else:
+        labels = None
 
-        labels = list(labels)  # coerce and/or shallow copy
+    if isinstance(samples_like, abc.Iterator):
+        # if we don't check this case we can get unexpected behaviour where an
+        # iterator can be depleted
+        raise TypeError('samples_like cannot be an iterator')
 
-        if len(labels) != samples.shape[1]:
-            raise ValueError("labels and samples_like dimensions do not match")
+    if isinstance(samples_like, abc.Mapping):
+        return as_samples(([samples_like], labels), dtype=dtype)
 
-        return samples, labels
+    if (isinstance(samples_like, list) and samples_like and
+            isinstance(samples_like[0], numbers.Number)):
+        # this is not actually necessary but it speeds up the
+        # samples_like = [1, 0, 1,...] case significantly
+        return as_samples(([samples_like], labels), dtype=dtype)
 
-    # if no ` is specified and the array_like doesn't already have a dtype, we default to int8
-    if dtype is None and not hasattr(samples_like, 'dtype'):
-        dtype = np.int8
+    if not isinstance(samples_like, np.ndarray):
+        if any(isinstance(sample, abc.Mapping) for sample in samples_like):
+            # go through samples-like, turning the dicts into lists
+            samples_like, old = list(samples_like), samples_like
 
-    if len(samples_like) == 0:
-        return np.empty((0, 0), dtype=dtype), []
-
-    if isinstance(samples_like, Mapping):
-        return as_samples([samples_like], dtype=dtype)
+            if labels is None:
+                first = samples_like[0]
+                if isinstance(first, abc.Mapping):
+                    labels = list(first)
+                else:
+                    labels = list(range(len(first)))
 
-    if isinstance(samples_like, Iterable) and all(isinstance(sample, Mapping) for sample in samples_like):
-        # list of dicts
-        return _samples_dicts_to_array(samples_like)
+            for idx, sample in enumerate(old):
+                if isinstance(sample, abc.Mapping):
+                    try:
+                        samples_like[idx] = [sample[v] for v in labels]
+                    except KeyError:
+                        raise ValueError("samples_like and labels do not match")
 
-    # anything else should be array_like, which covers ndarrays, lists of lists, etc
+    if dtype is None and not hasattr(samples_like, 'dtype'):
+        dtype = np.int8
 
-    try:
-        samples_like = np.asarray(samples_like, dtype=dtype)
-    except (ValueError, TypeError):
-        raise TypeError("unknown format for samples_like")
+    # samples-like should now be array-like
+    arr = np.array(samples_like, dtype=dtype, copy=copy, order=order)
 
-    # want 2D array
-    if samples_like.ndim == 1:
-        samples_like = np.expand_dims(samples_like, 0)
-    elif samples_like.ndim > 2:
-        ValueError("expected sample_like to be <= 2 dimensions")
+    if arr.ndim > 2:
+        raise ValueError("expected samples_like to be <= 2 dimensions")
+    if arr.ndim < 2:
+        if arr.size:
+            arr = np.atleast_2d(arr)
+        elif labels:  # is not None and len > 0
+            arr = arr.reshape((0, len(labels)))
+        else:
+            arr = arr.reshape((0, 0))
 
-    return samples_like, list(range(samples_like.shape[1]))
+    # ok we're basically done, just need to check against the labels
+    if labels is None:
+        return arr, list(range(arr.shape[1]))
+    elif len(labels) != arr.shape[1]:
+        print(arr, arr.shape, samples_like, labels, len(labels))
+        raise ValueError("samples_like and labels dimensions do not match")
+    else:
+        return arr, labels
 
 
 def concatenate(samplesets, defaults=None):
@@ -215,7 +240,7 @@ def _iter_records(samplesets, vartype, variables):
             yield samples.record
 
 
-class SampleSet(Iterable, Sized):
+class SampleSet(abc.Iterable, abc.Sized):
     """Samples and any other data returned by dimod samplers.
 
     Args:
@@ -946,30 +971,3 @@ class SampleSet(Iterable, Sized):
             df.loc[:, field] = self.record[field]
 
         return df
-
-
-def _samples_dicts_to_array(samples_dicts):
-    """Convert an iterable of samples where each sample is a dict to a numpy 2d array. Also
-    determines the labels is they are None.
-    """
-    itersamples = iter(samples_dicts)
-
-    first_sample = next(itersamples)
-
-    labels = list(first_sample)
-
-    num_variables = len(labels)
-
-    def _iter_samples():
-        yield np.fromiter((first_sample[v] for v in labels),
-                          count=num_variables, dtype=np.int8)
-
-        try:
-            for sample in itersamples:
-                yield np.fromiter((sample[v] for v in labels),
-                                  count=num_variables, dtype=np.int8)
-        except KeyError:
-            msg = ("Each dict in 'samples' must have the same keys.")
-            raise ValueError(msg)
-
-    return np.stack(list(_iter_samples())), labels
diff --git a/docs/reference/bqm_functions.rst b/docs/reference/bqm_functions.rst
index 38156965..29f93774 100644
--- a/docs/reference/bqm_functions.rst
+++ b/docs/reference/bqm_functions.rst
@@ -27,6 +27,8 @@ Higher-order Polynomials
 
    make_quadratic
    poly_energy
+   poly_energies
+   create_poly
 
 
 Roof-duality
diff --git a/docs/reference/ref_samplers_composites.rst b/docs/reference/ref_samplers_composites.rst
index fae35f45..f958f5cd 100644
--- a/docs/reference/ref_samplers_composites.rst
+++ b/docs/reference/ref_samplers_composites.rst
@@ -114,7 +114,24 @@ Methods
 
    StructureComposite.sample
 
+Higher Order Composite
+----------------------
 
+.. automodule:: dimod.reference.composites.higherordercomposite
+
+Class
+~~~~~
+
+.. currentmodule:: dimod.reference.composites
+.. autoclass:: HigherOrderComposite
+
+Methods
+~~~~~~~
+
+.. autosummary::
+   :toctree: generated/
+
+   HigherOrderComposite.sample_ising
 
 Fixed Variable Composite
 ------------------------
@@ -133,4 +150,4 @@ Methods
 .. autosummary::
    :toctree: generated/
 
-   FixedVariableComposite.sample
\ No newline at end of file
+   FixedVariableComposite.sample
diff --git a/docs/reference/sampleset.rst b/docs/reference/sampleset.rst
index 74431514..d9705d80 100644
--- a/docs/reference/sampleset.rst
+++ b/docs/reference/sampleset.rst
@@ -4,12 +4,10 @@
 Samples
 =======
 
-This module is under development as a possible replacement for :obj:`.Response`.
-
 sample_like Objects
 ===================
 
-.. automodule:: dimod.sampleset
+.. currentmodule:: dimod
 
 .. autosummary::
    :toctree: generated/
@@ -22,7 +20,6 @@ SampleSet
 Class
 -----
 
-.. currentmodule:: dimod
 .. autoclass:: SampleSet
 
 Properties


Test patch:
diff --git a/dimod/testing/asserts.py b/dimod/testing/asserts.py
index 5071f147..b28d929f 100644
--- a/dimod/testing/asserts.py
+++ b/dimod/testing/asserts.py
@@ -18,7 +18,10 @@ This module contains asserts that can be used to test the correctness of dimod s
 composites and responses. This is useful for checking that a created sampler correctly fulfills
 the dimod API.
 """
-from collections import Mapping, Sequence, Set
+try:
+    import collections.abc as abc
+except ImportError:
+    import collections as abc
 
 import dimod
 
@@ -58,12 +61,12 @@ def assert_sampler_api(sampler):
     msg = "instantiated sampler must have a 'parameters' property, set to a Mapping"
     assert hasattr(sampler, 'parameters'), msg
     assert not callable(sampler.parameters), msg
-    assert isinstance(sampler.parameters, Mapping), msg
+    assert isinstance(sampler.parameters, abc.Mapping), msg
 
     msg = "instantiated sampler must have a 'properties' property, set to a Mapping"
     assert hasattr(sampler, 'properties'), msg
     assert not callable(sampler.properties), msg
-    assert isinstance(sampler.properties, Mapping), msg
+    assert isinstance(sampler.properties, abc.Mapping), msg
 
 
 def assert_composite_api(composed_sampler):
@@ -90,7 +93,7 @@ def assert_composite_api(composed_sampler):
 
     msg = "instantiated composed sampler must have a 'children' property, set to a list (or Sequence)"
     assert hasattr(composed_sampler, 'children'), msg
-    assert isinstance(composed_sampler.children, Sequence), msg
+    assert isinstance(composed_sampler.children, abc.Sequence), msg
 
     msg = "instantiated composed sampler must have a 'child' property, set to one of sampler.children"
     assert hasattr(composed_sampler, 'child'), msg
@@ -123,24 +126,24 @@ def assert_structured_api(sampler):
     msg = ("instantiated structured sampler must have an 'adjacency' property formatted as a dict "
            "where the keys are the nodes and the values are sets of all node adjacency to the key")
     assert hasattr(sampler, 'adjacency'), msg
-    assert isinstance(sampler.adjacency, Mapping), msg
+    assert isinstance(sampler.adjacency, abc.Mapping), msg
     for u, neighborhood in sampler.adjacency.items():
-        assert isinstance(neighborhood, Set), msg
+        assert isinstance(neighborhood, abc.Set), msg
         for v in neighborhood:
             assert v in sampler.adjacency, msg
             assert u in sampler.adjacency[v], msg
 
     msg = "instantiated structured sampler must have a 'nodelist' property, set to a list"
     assert hasattr(sampler, 'nodelist'), msg
-    assert isinstance(sampler.nodelist, Sequence), msg
+    assert isinstance(sampler.nodelist, abc.Sequence), msg
     for v in sampler.nodelist:
         assert v in sampler.adjacency, msg
 
     msg = "instantiated structured sampler must have a 'edge' property, set to a list of 2-lists/tuples"
     assert hasattr(sampler, 'edgelist'), msg
-    assert isinstance(sampler.edgelist, Sequence), msg
+    assert isinstance(sampler.edgelist, abc.Sequence), msg
     for edge in sampler.edgelist:
-        assert isinstance(edge, Sequence), msg
+        assert isinstance(edge, abc.Sequence), msg
         assert len(edge) == 2, msg
 
         u, v = edge
@@ -181,7 +184,7 @@ def assert_response_energies(response, bqm, precision=7):
     assert isinstance(response, dimod.SampleSet), "expected response to be a dimod SampleSet object"
 
     for sample, energy in response.data(['sample', 'energy']):
-        assert isinstance(sample, Mapping), "'for sample in response', each sample should be a Mapping"
+        assert isinstance(sample, abc.Mapping), "'for sample in response', each sample should be a Mapping"
 
         for v, value in sample.items():
             assert v in bqm.linear, 'sample contains a variable not in the given bqm'
diff --git a/tests/test_higherorder.py b/tests/test_higherorder.py
index 47c849eb..1dbcb381 100644
--- a/tests/test_higherorder.py
+++ b/tests/test_higherorder.py
@@ -13,12 +13,12 @@
 #    limitations under the License.
 #
 # ================================================================================================
-import unittest
 import itertools
+import unittest
 
 import dimod
+from dimod import make_quadratic, poly_energy, poly_energies
 
-from dimod import make_quadratic, poly_energy
 
 class TestMakeQuadratic(unittest.TestCase):
 
@@ -30,7 +30,8 @@ class TestMakeQuadratic(unittest.TestCase):
             self.assertIn(v, bqm)
         self.assertEqual(len(bqm), 4)  # has an auxiliary variable
 
-        for sample, energy in dimod.ExactSolver().sample(bqm).data(['sample', 'energy']):
+        for sample, energy in dimod.ExactSolver().sample(bqm).data(
+                ['sample', 'energy']):
             if energy == 0:
                 self.assertEqual(sample['a'] * sample['b'], sample['p'])
             if sample['a'] * sample['b'] != sample['p']:
@@ -47,7 +48,8 @@ class TestMakeQuadratic(unittest.TestCase):
         self.assertEqual(len(bqm), 3)
 
         seen_configs = set()
-        for sample, energy in dimod.ExactSolver().sample(bqm).data(['sample', 'energy']):
+        for sample, energy in dimod.ExactSolver().sample(bqm).data(
+                ['sample', 'energy']):
             if energy == 0:
                 self.assertEqual(sample['a'] * sample['b'], sample['p'])
                 seen_configs.add(tuple(sample[v] for v in variables))
@@ -70,7 +72,8 @@ class TestMakeQuadratic(unittest.TestCase):
             energy = poly_energy(sample, poly)
 
             reduced_energies = []
-            for aux_config in itertools.product((-1, 1), repeat=len(aux_variables)):
+            for aux_config in itertools.product((-1, 1),
+                                                repeat=len(aux_variables)):
                 aux_sample = dict(zip(aux_variables, aux_config))
                 aux_sample.update(sample)
                 reduced_energies.append(bqm.energy(aux_sample))
@@ -92,7 +95,8 @@ class TestMakeQuadratic(unittest.TestCase):
             energy = poly_energy(sample, poly)
 
             reduced_energies = []
-            for aux_config in itertools.product((-1, 1), repeat=len(aux_variables)):
+            for aux_config in itertools.product((-1, 1),
+                                                repeat=len(aux_variables)):
                 aux_sample = dict(zip(aux_variables, aux_config))
                 aux_sample.update(sample)
                 reduced_energies.append(bqm.energy(aux_sample))
@@ -104,7 +108,8 @@ class TestMakeQuadratic(unittest.TestCase):
                 (0,): .4,
                 (): .5}
 
-        bqm = make_quadratic(poly, 5.0, bqm=dimod.BinaryQuadraticModel.empty(dimod.SPIN))
+        bqm = make_quadratic(poly, 5.0,
+                             bqm=dimod.BinaryQuadraticModel.empty(dimod.SPIN))
 
         variables = set().union(*poly)
         aux_variables = tuple(set(bqm.linear) - variables)
@@ -116,7 +121,8 @@ class TestMakeQuadratic(unittest.TestCase):
             energy = poly_energy(sample, poly)
 
             reduced_energies = []
-            for aux_config in itertools.product((-1, 1), repeat=len(aux_variables)):
+            for aux_config in itertools.product((-1, 1),
+                                                repeat=len(aux_variables)):
                 aux_sample = dict(zip(aux_variables, aux_config))
                 aux_sample.update(sample)
                 reduced_energies.append(bqm.energy(aux_sample))
@@ -132,7 +138,9 @@ class TestMakeQuadratic(unittest.TestCase):
         poly.update({(v,): bias for v, bias in h.items()})
         poly[()] = off
 
-        bqm = make_quadratic(J, 5.0, bqm=dimod.BinaryQuadraticModel.from_ising(h, {}, off))
+        bqm = make_quadratic(J, 5.0,
+                             bqm=dimod.BinaryQuadraticModel.from_ising(h, {},
+                                                                       off))
 
         variables = set(h).union(*J)
         aux_variables = tuple(set(bqm.linear) - variables)
@@ -144,7 +152,8 @@ class TestMakeQuadratic(unittest.TestCase):
             energy = poly_energy(sample, poly)
 
             reduced_energies = []
-            for aux_config in itertools.product((-1, 1), repeat=len(aux_variables)):
+            for aux_config in itertools.product((-1, 1),
+                                                repeat=len(aux_variables)):
                 aux_sample = dict(zip(aux_variables, aux_config))
                 aux_sample.update(sample)
                 reduced_energies.append(bqm.energy(aux_sample))
@@ -160,7 +169,9 @@ class TestMakeQuadratic(unittest.TestCase):
         poly.update({(v,): bias for v, bias in h.items()})
         poly[()] = off
 
-        bqm = make_quadratic(J, 10.0, bqm=dimod.BinaryQuadraticModel.from_ising(h, {}, off))
+        bqm = make_quadratic(J, 10.0,
+                             bqm=dimod.BinaryQuadraticModel.from_ising(h, {},
+                                                                       off))
 
         variables = set(h).union(*J)
         aux_variables = tuple(set(bqm.linear) - variables)
@@ -172,7 +183,8 @@ class TestMakeQuadratic(unittest.TestCase):
             energy = poly_energy(sample, poly)
 
             reduced_energies = []
-            for aux_config in itertools.product((-1, 1), repeat=len(aux_variables)):
+            for aux_config in itertools.product((-1, 1),
+                                                repeat=len(aux_variables)):
                 aux_sample = dict(zip(aux_variables, aux_config))
                 aux_sample.update(sample)
                 reduced_energies.append(bqm.energy(aux_sample))
@@ -185,7 +197,8 @@ class TestMakeQuadratic(unittest.TestCase):
 
         bqm = dimod.make_quadratic(poly, 10.0, dimod.BINARY)
 
-        self.assertEqual(bqm, dimod.BinaryQuadraticModel({0: .5}, {}, 1.3, dimod.BINARY))
+        self.assertEqual(bqm, dimod.BinaryQuadraticModel({0: .5}, {}, 1.3,
+                                                         dimod.BINARY))
 
     def test_binary_polynomial(self):
 
@@ -203,9 +216,41 @@ class TestMakeQuadratic(unittest.TestCase):
             energy = poly_energy(sample, HUBO)
 
             reduced_energies = []
-            for aux_config in itertools.product((0, 1), repeat=len(aux_variables)):
+            for aux_config in itertools.product((0, 1),
+                                                repeat=len(aux_variables)):
                 aux_sample = dict(zip(aux_variables, aux_config))
                 aux_sample.update(sample)
                 reduced_energies.append(bqm.energy(aux_sample))
 
             self.assertAlmostEqual(energy, min(reduced_energies))
+
+    def test_poly_energies(self):
+        from dimod.higherorder import create_poly
+
+        linear = {0: 1.0, 1: 1.0}
+        j = {(0, 1, 2): 0.5}
+        poly = create_poly(linear, j)
+        samples = [[1, 1, -1], [1, -1, 1], [1, 1, 1], [-1, 1, -1]]
+
+        en = poly_energies(samples, poly)
+        self.assertListEqual(list(en), [1.5, -0.5, 2.5, 0.5])
+
+        en = poly_energy(samples[0], poly)
+        self.assertAlmostEqual(en, 1.5)
+
+        with self.assertRaises(ValueError):
+            poly_energy(samples, poly)
+
+        poly = {('a',): 1.0, ('b',): 1.0, ('a', 'b', 'c'): 0.5}
+        samples = [{'a': 1, 'b': 1, 'c': -1},
+                   {'a': 1, 'b': -1, 'c': 1},
+                   {'a': 1, 'b': 1, 'c': 1},
+                   {'a': -1, 'b': 1, 'c': -1}]
+        en = poly_energies(samples, poly)
+        self.assertListEqual(list(en), [1.5, -0.5, 2.5, 0.5])
+
+        en = poly_energy(samples[0], poly)
+        self.assertAlmostEqual(en, 1.5)
+
+        with self.assertRaises(ValueError):
+            poly_energy(samples, poly)
diff --git a/tests/test_higherordercomposite.py b/tests/test_higherordercomposite.py
new file mode 100644
index 00000000..1b94c998
--- /dev/null
+++ b/tests/test_higherordercomposite.py
@@ -0,0 +1,94 @@
+# Copyright 2018 D-Wave Systems Inc.
+#
+#    Licensed under the Apache License, Version 2.0 (the "License");
+#    you may not use this file except in compliance with the License.
+#    You may obtain a copy of the License at
+#
+#        http://www.apache.org/licenses/LICENSE-2.0
+#
+#    Unless required by applicable law or agreed to in writing, software
+#    distributed under the License is distributed on an "AS IS" BASIS,
+#    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#    See the License for the specific language governing permissions and
+#    limitations under the License.
+#
+# ================================================================================================
+
+import numpy as np
+import unittest
+
+import dimod.testing as dtest
+
+from dimod import HigherOrderComposite, ExactSolver
+
+
+class TestFixedVariableComposite(unittest.TestCase):
+
+    def test_instantiation_smoketest(self):
+        sampler = HigherOrderComposite(ExactSolver())
+
+        dtest.assert_sampler_api(sampler)
+
+    def test_sample(self):
+        linear = {0: -0.5, 1: -0.3, 2: -0.8}
+        quadratic = {(0, 1, 2): -1.7}
+
+        sampler = HigherOrderComposite(ExactSolver())
+        response = sampler.sample_ising(linear, quadratic, penalty_strength=10,
+                                        keep_penalty_variables=False,
+                                        discard_unsatisfied=False)
+
+        self.assertEqual(response.first.sample, {0: 1, 1: 1, 2: 1})
+        self.assertAlmostEqual(response.first.energy, -3.3)
+        self.assertFalse(np.prod(response.record.penalty_satisfaction))
+
+    def test_discard(self):
+        linear = {0: -0.5, 1: -0.3, 2: -0.8}
+        quadratic = {(0, 1, 2): -1.7}
+
+        sampler = HigherOrderComposite(ExactSolver())
+        response = sampler.sample_ising(linear, quadratic, penalty_strength=10,
+                                        discard_unsatisfied=True)
+
+        self.assertEqual(response.first.sample, {0: 1, 1: 1, 2: 1})
+        self.assertAlmostEqual(response.first.energy, -3.3)
+        self.assertTrue(np.prod(response.record.penalty_satisfaction))
+
+    def test_penalty_variables(self):
+        linear = {0: -0.5, 1: -0.3, 2: -0.8}
+        quadratic = {(0, 1, 2): -1.7}
+
+        sampler = HigherOrderComposite(ExactSolver())
+        response = sampler.sample_ising(linear, quadratic, penalty_strength=10,
+                                        keep_penalty_variables=True,
+                                        discard_unsatisfied=True)
+
+        self.assertEqual(len(response.first.sample),5)
+        self.assertAlmostEqual(response.first.energy, -3.3)
+        self.assertTrue(response.first.penalty_satisfaction)
+
+    def test_already_qubo(self):
+        linear = {0: -0.5, 1: -0.3}
+        quadratic = {(0, 1): -1.7}
+
+        sampler = HigherOrderComposite(ExactSolver())
+        response = sampler.sample_ising(linear, quadratic,
+                                        keep_penalty_variables=True,
+                                        discard_unsatisfied=False)
+
+        self.assertEqual(response.first.sample, {0: 1, 1: 1})
+        self.assertAlmostEqual(response.first.energy, -2.5)
+        self.assertTrue(response.first.penalty_satisfaction)
+
+    def test_already_qubo_2(self):
+        linear = {0: -0.5, 1: -0.3}
+        quadratic = {(0, 1): -1.7}
+
+        sampler = HigherOrderComposite(ExactSolver())
+        response = sampler.sample_ising(linear, quadratic, penalty_strength=10,
+                                        keep_penalty_variables=True,
+                                        discard_unsatisfied=True)
+
+        self.assertEqual(response.first.sample, {0: 1, 1: 1})
+        self.assertAlmostEqual(response.first.energy, -2.5)
+        self.assertTrue(response.first.penalty_satisfaction)
diff --git a/tests/test_sampleset.py b/tests/test_sampleset.py
index 08218875..eefa64fd 100644
--- a/tests/test_sampleset.py
+++ b/tests/test_sampleset.py
@@ -16,6 +16,8 @@
 import unittest
 import json
 
+from collections import OrderedDict
+
 import numpy as np
 
 import dimod
@@ -28,6 +30,102 @@ except ImportError:
     _pandas = False
 
 
+class Test_as_samples(unittest.TestCase):
+
+    def test_copy_false(self):
+        samples_like = np.ones((5, 5))
+        labels = list('abcde')
+        arr, lab = dimod.as_samples((samples_like, labels))
+        np.testing.assert_array_equal(arr, np.ones((5, 5)))
+        self.assertEqual(lab, list('abcde'))
+        self.assertIs(labels, lab)
+        self.assertTrue(np.shares_memory(arr, samples_like))
+
+    def test_dict_with_inconsistent_labels(self):
+        with self.assertRaises(ValueError):
+            dimod.as_samples(({'a': -1}, 'b'))
+
+    def test_dict_with_labels(self):
+        arr, labels = dimod.as_samples(({'a': -1}, 'a'))
+        np.testing.assert_array_equal(arr, [[-1]])
+        self.assertEqual(labels, ['a'])
+
+    def test_empty_dict(self):
+        # one sample, no labels
+        arr, labels = dimod.as_samples({})
+        np.testing.assert_array_equal(arr, np.zeros((1, 0)))
+        self.assertEqual(labels, [])
+
+    def test_empty_list(self):
+        # no samples, no labels
+        arr, labels = dimod.as_samples([])
+        np.testing.assert_array_equal(arr, np.zeros((0, 0)))
+        self.assertEqual(labels, [])
+
+    def test_empty_list_labelled(self):
+        # no samples, no labels
+        arr, labels = dimod.as_samples(([], []))
+        np.testing.assert_array_equal(arr, np.zeros((0, 0)))
+        self.assertEqual(labels, [])
+
+        # no samples, 1 label
+        arr, labels = dimod.as_samples(([], ['a']))
+        np.testing.assert_array_equal(arr, np.zeros((0, 1)))
+        self.assertEqual(labels, ['a'])
+
+        # no samples, 2 labels
+        arr, labels = dimod.as_samples(([], ['a', 'b']))
+        np.testing.assert_array_equal(arr, np.zeros((0, 2)))
+        self.assertEqual(labels, ['a', 'b'])
+
+    def test_empty_ndarray(self):
+        arr, labels = dimod.as_samples(np.ones(0))
+        np.testing.assert_array_equal(arr, np.zeros((0, 0)))
+        self.assertEqual(labels, [])
+
+    def test_iterator(self):
+        with self.assertRaises(TypeError):
+            dimod.as_samples(([-1] for _ in range(10)))
+
+    def test_iterator_labelled(self):
+        with self.assertRaises(TypeError):
+            dimod.as_samples(([-1] for _ in range(10)), 'a')
+
+    def test_list_of_empty(self):
+        arr, labels = dimod.as_samples([[], [], []])
+        np.testing.assert_array_equal(arr, np.empty((3, 0)))
+        self.assertEqual(labels, [])
+
+        arr, labels = dimod.as_samples([{}, {}, {}])
+        np.testing.assert_array_equal(arr, np.empty((3, 0)))
+        self.assertEqual(labels, [])
+
+        arr, labels = dimod.as_samples(np.empty((3, 0)))
+        np.testing.assert_array_equal(arr, np.empty((3, 0)))
+        self.assertEqual(labels, [])
+
+    def test_mixed_sampletype(self):
+        s1 = [0, 1]
+        s2 = OrderedDict([(1, 0), (0, 1)])
+        s3 = OrderedDict([(0, 1), (1, 0)])
+
+        arr, labels = dimod.as_samples([s1, s2, s3])
+        np.testing.assert_array_equal(arr, [[0, 1], [1, 0], [1, 0]])
+        self.assertEqual(labels, [0, 1])
+
+    def test_ndarray(self):
+        arr, labels = dimod.as_samples(np.ones(5, dtype=np.int32))
+        np.testing.assert_array_equal(arr, np.ones((1, 5)))
+        self.assertEqual(labels, list(range(5)))
+        self.assertEqual(arr.dtype, np.int32)
+
+    def test_ndarray_labelled(self):
+        arr, labels = dimod.as_samples((np.ones(5, dtype=np.int32), 'abcde'))
+        np.testing.assert_array_equal(arr, np.ones((1, 5)))
+        self.assertEqual(labels, ['a', 'b', 'c', 'd', 'e'])
+        self.assertEqual(arr.dtype, np.int32)
+
+
 class TestSampleSet(unittest.TestCase):
     def test_from_samples(self):
 
@@ -119,8 +217,6 @@ class TestSampleSet(unittest.TestCase):
 
         self.assertEqual(len(dimod.SampleSet.from_samples([], dimod.SPIN, energy=[], a=1)), 0)
 
-        self.assertEqual(len(dimod.SampleSet.from_samples({}, dimod.SPIN, energy=[], a=1)), 0)
-
         self.assertEqual(len(dimod.SampleSet.from_samples(np.empty((0, 0)), dimod.SPIN, energy=[], a=1)), 0)
 
     def test_from_samples_with_aggregation(self):