Repository: dwavesystems/dimod

Issue:
Remove unnecessary type check in BinaryQuadraticModel.scale
https://github.com/dwavesystems/dimod/blob/7fda1ba0f92a44358d12c0e8f9831dc68c7768e7/dimod/binary_quadratic_model.py#L861

which will allow scaling by other number-like object in python (like sympy symbols).

Discussion:
codecov-io: # [Codecov](https://codecov.io/gh/dwavesystems/dimod/pull/367?src=pr&el=h1) Report
> Merging [#367](https://codecov.io/gh/dwavesystems/dimod/pull/367?src=pr&el=desc) into [master](https://codecov.io/gh/dwavesystems/dimod/commit/7fda1ba0f92a44358d12c0e8f9831dc68c7768e7?src=pr&el=desc) will **increase** coverage by `0.01%`.
> The diff coverage is `100%`.

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

```diff
@@            Coverage Diff             @@
##           master     #367      +/-   ##
==========================================
+ Coverage   91.55%   91.56%   +0.01%     
==========================================
  Files          37       37              
  Lines        2154     2158       +4     
==========================================
+ Hits         1972     1976       +4     
  Misses        182      182
```


| [Impacted Files](https://codecov.io/gh/dwavesystems/dimod/pull/367?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [dimod/binary\_quadratic\_model.py](https://codecov.io/gh/dwavesystems/dimod/pull/367/diff?src=pr&el=tree#diff-ZGltb2QvYmluYXJ5X3F1YWRyYXRpY19tb2RlbC5weQ==) | `95.36% <100%> (ø)` | :arrow_up: |
| [dimod/sampleset.py](https://codecov.io/gh/dwavesystems/dimod/pull/367/diff?src=pr&el=tree#diff-ZGltb2Qvc2FtcGxlc2V0LnB5) | `94.18% <100%> (+0.06%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/dwavesystems/dimod/pull/367?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/367?src=pr&el=footer). Last update [7fda1ba...2abffa3](https://codecov.io/gh/dwavesystems/dimod/pull/367?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


Patch:
diff --git a/dimod/__init__.py b/dimod/__init__.py
index 4fd70e1b..b0b5000f 100644
--- a/dimod/__init__.py
+++ b/dimod/__init__.py
@@ -33,7 +33,7 @@ else:
 
 import dimod.testing
 
-from dimod.binary_quadratic_model import BinaryQuadraticModel
+from dimod.binary_quadratic_model import BinaryQuadraticModel, BQM
 import dimod.binary_quadratic_model
 
 import dimod.decorators
@@ -43,7 +43,7 @@ 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__
diff --git a/dimod/binary_quadratic_model.py b/dimod/binary_quadratic_model.py
index fa6a7a33..77f60957 100644
--- a/dimod/binary_quadratic_model.py
+++ b/dimod/binary_quadratic_model.py
@@ -66,6 +66,8 @@ from dimod.utilities import resolve_label_conflict
 from dimod.views import LinearView, QuadraticView, AdjacencyView, SampleView
 from dimod.vartypes import Vartype
 
+__all__ = 'BinaryQuadraticModel', 'BQM'
+
 
 class BinaryQuadraticModel(abc.Sized, abc.Container, abc.Iterable):
     """Encodes a binary quadratic model.
@@ -858,8 +860,6 @@ class BinaryQuadraticModel(abc.Sized, abc.Container, abc.Iterable):
             0.5
 
         """
-        if not isinstance(scalar, Number):
-            raise TypeError("expected scalar to be a Number")
 
         if ignored_variables is None:
             ignored_variables = set()
@@ -1854,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 not isinstance(h, abc.Mapping):
+        if isinstance(h, abc.Sequence):
             h = dict(enumerate(h))
 
         return cls(h, J, offset, Vartype.SPIN)
@@ -2342,3 +2342,7 @@ class BinaryQuadraticModel(abc.Sized, abc.Container, abc.Iterable):
             bqm.add_interaction(u, v, 0.0)
 
         return bqm
+
+
+BQM = BinaryQuadraticModel
+"""Alias for :obj:`.BinaryQuadraticModel`"""
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/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/sampleset.py b/dimod/sampleset.py
index d19f6142..b8636fed 100644
--- a/dimod/sampleset.py
+++ b/dimod/sampleset.py
@@ -559,12 +559,20 @@ class SampleSet(abc.Iterable, abc.Sized):
     def first(self):
         """Sample with the lowest-energy.
 
+        Raises:
+            ValueError: If empty.
+
         Example:
 
-            >>> sampleset.first   # doctest: +SKIP
+            >>> sampleset = dimod.ExactSolver().sample_ising({'a': 1}, {('a', 'b'): 1})
+            >>> sampleset.first
             Sample(sample={'a': -1, 'b': -1}, energy=-1.0, num_occurrences=1)
+
         """
-        return next(self.data(sorted_by='energy', name='Sample'))
+        try:
+            return next(self.data(sorted_by='energy', name='Sample'))
+        except StopIteration:
+            raise ValueError('{} is empty'.format(self.__class__.__name__))
 
     @property
     def info(self):
@@ -683,7 +691,7 @@ class SampleSet(abc.Iterable, abc.Sized):
             for sample in itertools.islice(self.samples(n=None, sorted_by=sorted_by), n):
                 yield sample
 
-    def data(self, fields=None, sorted_by='energy', name='Sample'):
+    def data(self, fields=None, sorted_by='energy', name='Sample', reverse=False):
         """Iterate over the data in the :class:`SampleSet`.
 
         Args:
@@ -698,6 +706,9 @@ class SampleSet(abc.Iterable, abc.Sized):
             name (str/None, optional, default='Sample'):
                 Name of the yielded namedtuples or None to yield regular tuples.
 
+            reverse (bool, optional, default=False):
+                If True, yield in reverse order.
+
         Yields:
             namedtuple/tuple: The data in the :class:`SampleSet`, in the order specified by the input
             `fields`.
@@ -736,6 +747,9 @@ class SampleSet(abc.Iterable, abc.Sized):
         else:
             order = np.argsort(record[sorted_by])
 
+        if reverse:
+            order = np.flip(order)
+
         if name is None:
             # yielding a tuple
             def _pack(values):
diff --git a/docs/reference/binary_quadratic_model.rst b/docs/reference/binary_quadratic_model.rst
index 1923308c..f00d567c 100644
--- a/docs/reference/binary_quadratic_model.rst
+++ b/docs/reference/binary_quadratic_model.rst
@@ -110,3 +110,7 @@ Converting to other types
    BinaryQuadraticModel.to_qubo
    BinaryQuadraticModel.to_pandas_dataframe
    BinaryQuadraticModel.to_serializable
+
+Alias
+=====
+.. autoclass:: BQM
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


Test patch:
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 eefa64fd..a2b34bc7 100644
--- a/tests/test_sampleset.py
+++ b/tests/test_sampleset.py
@@ -242,6 +242,16 @@ class TestSampleSet(unittest.TestCase):
                          dimod.SampleSet.from_samples(([-1, 1], 'ab'), dimod.SPIN, energy=1))
 
 
+class TestSampleSetIteration(unittest.TestCase):
+    def test_data_reverse(self):
+        bqm = dimod.BinaryQuadraticModel.from_ising({}, {'ab': -1})
+        sampleset = dimod.SampleSet.from_samples_bqm([{'a': -1, 'b': 1}, {'a': 1, 'b': 1}], bqm)
+
+        samples = list(sampleset.data())
+        reversed_samples = list(sampleset.data(reverse=True))
+        self.assertEqual(samples, list(reversed(reversed_samples)))
+
+
 class TestSampleSetSerialization(unittest.TestCase):
 
     def test_functional_simple_shapes(self):
@@ -298,6 +308,12 @@ class TestSampleSet_to_pandas_dataframe(unittest.TestCase):
         pd.testing.assert_frame_equal(df, other, check_dtype=False)
 
 
+class TestSampleSet_first(unittest.TestCase):
+    def test_empty(self):
+        with self.assertRaises(ValueError):
+            dimod.SampleSet.from_samples([], dimod.SPIN, energy=[]).first
+
+
 class Test_concatenate(unittest.TestCase):
     def test_simple(self):
         ss0 = dimod.SampleSet.from_samples([-1, +1], dimod.SPIN, energy=-1)