Repository: dwavesystems/dimod

Issue:
Support higher order polynomials
Right now in dimod there is 'dimod/higherorder.py` which supports higher-order Ising problems. But as this is proving to be a more common use case, we should consider creating a proper `HUBO` (or similar)
class.

We could also define a Sampler ABC similar to `Structured`, to signal samplers/composites that support HUBOs. These would then implement a `sample_hubo` method.

Discussion:
arcondello: One thing that occurs to me is that using frozenset under the hood rules out `a^2` for instance.
codecov-io: # [Codecov](https://codecov.io/gh/dwavesystems/dimod/pull/396?src=pr&el=h1) Report
> Merging [#396](https://codecov.io/gh/dwavesystems/dimod/pull/396?src=pr&el=desc) into [master](https://codecov.io/gh/dwavesystems/dimod/commit/eafc63e2840417d278dde055c252d93b81d452d1?src=pr&el=desc) will **decrease** coverage by `2.27%`.
> The diff coverage is `78.27%`.

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

```diff
@@            Coverage Diff             @@
##           master     #396      +/-   ##
==========================================
- Coverage    92.9%   90.62%   -2.28%     
==========================================
  Files          40       43       +3     
  Lines        2494     2721     +227     
==========================================
+ Hits         2317     2466     +149     
- Misses        177      255      +78
```


| [Impacted Files](https://codecov.io/gh/dwavesystems/dimod/pull/396?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [dimod/higherorder/\_\_init\_\_.py](https://codecov.io/gh/dwavesystems/dimod/pull/396/diff?src=pr&el=tree#diff-ZGltb2QvaGlnaGVyb3JkZXIvX19pbml0X18ucHk=) | `100% <100%> (ø)` | |
| [dimod/core/\_\_init\_\_.py](https://codecov.io/gh/dwavesystems/dimod/pull/396/diff?src=pr&el=tree#diff-ZGltb2QvY29yZS9fX2luaXRfXy5weQ==) | `100% <100%> (ø)` | :arrow_up: |
| [dimod/higherorder/utils.py](https://codecov.io/gh/dwavesystems/dimod/pull/396/diff?src=pr&el=tree#diff-ZGltb2QvaGlnaGVyb3JkZXIvdXRpbHMucHk=) | `80.18% <100%> (ø)` | |
| [dimod/reference/composites/\_\_init\_\_.py](https://codecov.io/gh/dwavesystems/dimod/pull/396/diff?src=pr&el=tree#diff-ZGltb2QvcmVmZXJlbmNlL2NvbXBvc2l0ZXMvX19pbml0X18ucHk=) | `100% <100%> (ø)` | :arrow_up: |
| [dimod/reference/composites/scalecomposite.py](https://codecov.io/gh/dwavesystems/dimod/pull/396/diff?src=pr&el=tree#diff-ZGltb2QvcmVmZXJlbmNlL2NvbXBvc2l0ZXMvc2NhbGVjb21wb3NpdGUucHk=) | `94.69% <100%> (+0.35%)` | :arrow_up: |
| [dimod/core/polysampler.py](https://codecov.io/gh/dwavesystems/dimod/pull/396/diff?src=pr&el=tree#diff-ZGltb2QvY29yZS9wb2x5c2FtcGxlci5weQ==) | `67.74% <67.74%> (ø)` | |
| [dimod/higherorder/polynomial.py](https://codecov.io/gh/dwavesystems/dimod/pull/396/diff?src=pr&el=tree#diff-ZGltb2QvaGlnaGVyb3JkZXIvcG9seW5vbWlhbC5weQ==) | `74.35% <74.35%> (ø)` | |
| [...imod/reference/composites/higherordercomposites.py](https://codecov.io/gh/dwavesystems/dimod/pull/396/diff?src=pr&el=tree#diff-ZGltb2QvcmVmZXJlbmNlL2NvbXBvc2l0ZXMvaGlnaGVyb3JkZXJjb21wb3NpdGVzLnB5) | `87.09% <82.97%> (ø)` | |
| ... and [3 more](https://codecov.io/gh/dwavesystems/dimod/pull/396/diff?src=pr&el=tree-more) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/dwavesystems/dimod/pull/396?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/396?src=pr&el=footer). Last update [eafc63e...dc1e1bb](https://codecov.io/gh/dwavesystems/dimod/pull/396?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 28351d42..6f3e3e87 100644
--- a/dimod/__init__.py
+++ b/dimod/__init__.py
@@ -43,7 +43,7 @@ import dimod.generators
 from dimod.exceptions import *
 import dimod.exceptions
 
-from dimod.higherorder import make_quadratic, poly_energy, poly_energies
+from dimod.higherorder import make_quadratic, poly_energy, poly_energies, BinaryPolynomial
 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 61d713bb..252f0243 100644
--- a/dimod/binary_quadratic_model.py
+++ b/dimod/binary_quadratic_model.py
@@ -960,7 +960,7 @@ class BinaryQuadraticModel(abc.Sized, abc.Container, abc.Iterable):
                                                   quadratic_range))
 
         lin_min, lin_max = min_and_max([v for k, v in self.linear.items()
-                                    if k not in ignored_variables])
+                                     if k not in ignored_variables])
         quad_min, quad_max = min_and_max([v for (a,b),
                                                 v in self.quadratic.items()
                                           if ((a, b) not in ignored_interactions
diff --git a/dimod/core/__init__.py b/dimod/core/__init__.py
index c3eb98f8..e1f5d0f0 100644
--- a/dimod/core/__init__.py
+++ b/dimod/core/__init__.py
@@ -14,6 +14,7 @@
 #
 # ================================================================================================
 
-from dimod.core.composite import *
-from dimod.core.sampler import *
-from dimod.core.structured import *
+from dimod.core.composite import Composite, ComposedSampler
+from dimod.core.polysampler import PolySampler, ComposedPolySampler
+from dimod.core.sampler import Sampler
+from dimod.core.structured import Structured
diff --git a/dimod/core/polysampler.py b/dimod/core/polysampler.py
new file mode 100644
index 00000000..ed676d1a
--- /dev/null
+++ b/dimod/core/polysampler.py
@@ -0,0 +1,128 @@
+# 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.
+#
+# ============================================================================
+"""
+It is possible to construct samplers that handle binary polynomials - problems
+that have binary variables but they are not constrained to quadratic
+interactions.
+
+"""
+import abc
+import warnings
+
+from six import add_metaclass
+
+from dimod.core.composite import Composite
+from dimod.higherorder.polynomial import BinaryPolynomial
+
+__all__ = 'PolySampler', 'ComposedPolySampler'
+
+
+@add_metaclass(abc.ABCMeta)
+class PolySampler:
+    """Sampler supports binary polynomials.
+
+    Binary polynomials are an extension of binary quadratic models that allow
+    higher-order interactions.
+
+    """
+    @abc.abstractproperty  # for python2 compatibility
+    def parameters(self):
+        """dict: A dict where keys are the keyword parameters accepted by the sampler
+        methods and values are lists of the properties relevant to each parameter.
+        """
+        pass
+
+    @abc.abstractproperty  # for python2 compatibility
+    def properties(self):
+        """dict: A dict containing any additional information about the sampler.
+        """
+        pass
+
+    @abc.abstractmethod
+    def sample_poly(self, polynomial, **kwargs):
+        """Sample from a higher-order polynomial."""
+        pass
+
+    def sample_hising(self, h, J, **kwargs):
+        """Sample from a higher-order Ising model.
+
+        Convert the given higher-order Ising model to a :obj:`.BinaryPolynomial`
+        and invoke `sample_poly`.
+
+        Args:
+            h (dict):
+                The variable biases of the Ising problem.
+
+            J (dict):
+                The interaction biases of the Ising problem.
+
+            **kwargs:
+                See `sample_poly` for additional keyword definitions.
+
+        Returns:
+            :obj:`.SampleSet`
+
+        """
+        return self.sample_poly(BinaryPolynomial.from_hising(h, J), **kwargs)
+
+    def sample_hubo(self, H, **kwargs):
+        """Sample from a higher-order unconstrained binary optimization problem.
+
+        Convert the given higher-order unconstrained binary optimization
+        problem to a :obj:`.BinaryPolynomial` and invoke `sample_poly`.
+
+        Args:
+            H (dict):
+                The coefficients of the HUBO.
+
+            **kwargs:
+                See `sample_poly` for additional keyword definitions.
+
+        Returns:
+            :obj:`.SampleSet`
+
+        """
+        return self.sample_poly(BinaryPolynomial.from_hubo(H), **kwargs)
+
+    def sample(self, bqm, *args, **kwargs):
+        msg = ("PolySampler.sample is deprecated and will be removed in dimod "
+               "0.9.0. In the future, when using PolySamplers, you should use "
+               ".sample_poly")
+        warnings.warn(msg, DeprecationWarning)
+        poly = BinaryPolynomial(bqm.quadratic)
+        poly.update(((v,), bias) for v, bias in bqm.linear.items())
+        return self.sample_poly(poly, *args, **kwargs)
+
+    def sample_ising(self, *args, **kwargs):
+        msg = ("PolySampler.sample_ising is deprecated and will be removed in dimod "
+               "0.9.0. In the future, when using PolySamplers, you should use "
+               ".sample_hising")
+        warnings.warn(msg, DeprecationWarning)
+        return self.sample_hising(*args, **kwargs)
+
+    def sample_qubo(self, *args, **kwargs):
+        msg = ("PolySampler.sample_qubo is deprecated and will be removed in dimod "
+               "0.9.0. In the future, when using PolySamplers, you should use "
+               ".sample_hubo")
+        warnings.warn(msg, DeprecationWarning)
+        return self.sample_hubo(*args, **kwargs)
+
+
+class ComposedPolySampler(PolySampler, Composite):
+    """Abstract base class for dimod composed polynomial samplers.
+
+    Inherits from :class:`.PolySampler` and :class:`.Composite`."""
+    pass
diff --git a/dimod/higherorder/__init__.py b/dimod/higherorder/__init__.py
new file mode 100644
index 00000000..ba0b7494
--- /dev/null
+++ b/dimod/higherorder/__init__.py
@@ -0,0 +1,18 @@
+# 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.
+#
+# ============================================================================
+
+from dimod.higherorder.polynomial import BinaryPolynomial
+from dimod.higherorder.utils import make_quadratic, poly_energy, poly_energies
diff --git a/dimod/higherorder/polynomial.py b/dimod/higherorder/polynomial.py
new file mode 100644
index 00000000..4b2b47ec
--- /dev/null
+++ b/dimod/higherorder/polynomial.py
@@ -0,0 +1,528 @@
+# 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.
+#
+# ============================================================================
+from __future__ import division
+
+import itertools
+
+try:
+    import collections.abc as abc
+except ImportError:
+    import collections as abc
+
+from numbers import Number
+
+import numpy as np
+
+from dimod.decorators import vartype_argument
+from dimod.sampleset import as_samples
+from dimod.utilities import resolve_label_conflict
+from dimod.vartypes import Vartype
+
+__all__ = 'BinaryPolynomial',
+
+
+def asfrozenset(term):
+    """Convert to frozenset if it is not already"""
+    return term if isinstance(term, frozenset) else frozenset(term)
+
+
+class BinaryPolynomial(abc.MutableMapping):
+    """A polynomial with binary variables and real-valued coefficients.
+
+    Args:
+        poly (mapping/iterable):
+            Polynomial as a mapping of form {term: bias, ...}, where `term` is
+            a collection of variables and `bias` the associated bias. It can also
+            be an iterable of 2-tuples (term, bias).
+
+        vartype (:class:`.Vartype`/str/set):
+            Variable type for the binary quadratic model. Accepted input values:
+
+            * :class:`.Vartype.SPIN`, ``'SPIN'``, ``{-1, 1}``
+            * :class:`.Vartype.BINARY`, ``'BINARY'``, ``{0, 1}``
+
+    Attributes:
+        degree (int):
+            The degree of the polynomial.
+
+        variables (set):
+            The variables.
+
+        vartype (:class:`.Vartype`):
+            One of :class:`.Vartype.SPIN` or :class:`.Vartype.BINARY`.
+
+    Examples:
+
+        Binary polynomials can be constructed in many different ways. The
+        following are all equivalent
+
+        >>> poly = dimod.BinaryPolynomial({'a': -1, 'ab': 1}, dimod.SPIN)
+        >>> poly = dimod.BinaryPolynomial({('a',): -1, ('a', 'b'): 1}, dimod.SPIN)
+        >>> poly = dimod.BinaryPolynomial([('a', -1), (('a', 'b'), 1)], dimod.SPIN)
+        >>> poly = dimod.BinaryPolynomial({'a': -1, 'ab': .5, 'ba': .5}, dimod.SPIN)
+
+        Binary polynomials act a mutable mappings but the terms can be accessed with
+        any sequence.
+
+        >>> poly = dimod.BinaryPolynomial({'a': -1, 'ab': 1}, dimod.BINARY)
+        >>> poly['ab']
+        1
+        >>> poly['ba']
+        1
+        >>> poly[{'a', 'b'}]
+        1
+        >>> poly[('a', 'b')]
+        1
+        >>> poly['cd'] = 4
+        >>> poly['dc']
+        4
+
+    """
+    @vartype_argument('vartype')
+    def __init__(self, poly, vartype):
+        if isinstance(poly, abc.Mapping):
+            poly = poly.items()
+
+        # we need to aggregate the repeated terms
+        self._terms = terms = {}
+        for term, bias in poly:
+
+            fsterm = asfrozenset(term)
+
+            # when SPIN-valued, s^2 == 1, so we need to handle that case
+            # in BINARY, x^2 == x
+            if len(fsterm) < len(term) and vartype is Vartype.SPIN:
+                new = set()
+                term = tuple(term)  # make sure it has .count
+                for v in fsterm:
+                    if term.count(v) % 2:
+                        new.add(v)
+                fsterm = frozenset(new)
+
+            if fsterm in terms:
+                terms[fsterm] += bias
+            else:
+                terms[fsterm] = bias
+
+        self.vartype = vartype
+
+    def __contains__(self, term):
+        return asfrozenset(term) in self._terms
+
+    def __delitem__(self, term):
+        del self._terms[asfrozenset(term)]
+
+    def __eq__(self, other):
+        if not isinstance(other, BinaryPolynomial):
+            try:
+                other = type(self)(other, self.vartype)
+            except Exception:
+                # not a polynomial
+                return False
+
+        if self.vartype != other.vartype:
+            return False
+
+        for term, bias in self.items():
+            if bias and other[term] != bias:
+                return False
+        for term, bias in other.items():
+            if bias and self[term] != bias:
+                return False
+        return True
+
+    def __ne__(self, other):
+        return not (self == other)
+
+    def __getitem__(self, term):
+        return self._terms[asfrozenset(term)]
+
+    def __iter__(self):
+        return iter(self._terms)
+
+    def __len__(self):
+        return len(self._terms)
+
+    def __setitem__(self, term, bias):
+        self._terms[asfrozenset(term)] = bias
+
+    def __repr__(self):
+        return '{!s}({!r}, {!r})'.format(self.__class__.__name__,
+                                         self._terms, self.vartype.name)
+
+    @property
+    def variables(self):
+        """Variables of the polynomial."""
+        return set().union(*self._terms)
+
+    @property
+    def degree(self):
+        """Degree of the polynomial."""
+        if len(self) == 0:
+            return 0
+        return max(map(len, self._terms))
+
+    def copy(self):
+        """Create a shallow copy."""
+        return type(self)(self, self.vartype)
+
+    def energy(self, sample_like, dtype=np.float):
+        """The energy of the given sample.
+
+        Args:
+            sample_like (samples_like):
+                A raw sample. `sample_like` is an extension of
+                NumPy's array_like structure. See :func:`.as_samples`.
+
+            dtype (:class:`numpy.dtype`, optional):
+                The data type of the returned energies. Defaults to float.
+
+        Returns:
+            The energy.
+
+        """
+        energy, = self.energies(sample_like, dtype=dtype)
+        return energy
+
+    def energies(self, samples_like, dtype=np.float):
+        """The energies of the given samples.
+
+        Args:
+            samples_like (samples_like):
+                A collection of raw samples. `samples_like` is an extension of
+                NumPy's array_like structure. See :func:`.as_samples`.
+
+            dtype (:class:`numpy.dtype`, optional):
+                The data type of the returned energies. Defaults to float.
+
+        Returns:
+            :obj:`numpy.ndarray`: The energies.
+
+        """
+        samples, labels = as_samples(samples_like)
+        if labels:
+            idx, label = zip(*enumerate(labels))
+            labeldict = dict(zip(label, idx))
+        else:
+            labeldict = {}
+
+        num_samples = samples.shape[0]
+
+        energies = np.zeros(num_samples, dtype=dtype)
+        for term, bias in self.items():
+            if len(term) == 0:
+                energies += bias
+            else:
+                energies += np.prod([samples[:, labeldict[v]] for v in term], axis=0) * bias
+
+        return energies
+
+    def relabel_variables(self, mapping, inplace=True):
+        """Relabel variables of a binary polynomial as specified by mapping.
+
+        Args:
+            mapping (dict):
+                Dict mapping current variable labels to new ones. If an
+                incomplete mapping is provided, unmapped variables retain their
+                current labels.
+
+            inplace (bool, optional, default=True):
+                If True, the binary polynomial is updated in-place; otherwise, a
+                new binary polynomial is returned.
+
+        Returns:
+            :class:`.BinaryPolynomial`: A binary polynomial with the variables
+            relabeled. If `inplace` is set to True, returns itself.
+
+        """
+        if not inplace:
+            return self.copy().relabel_variables(mapping, inplace=True)
+
+        try:
+            old_labels = set(mapping)
+            new_labels = set(mapping.values())
+        except TypeError:
+            raise ValueError("mapping targets must be hashable objects")
+
+        variables = self.variables
+        for v in new_labels:
+            if v in variables and v not in old_labels:
+                raise ValueError(('A variable cannot be relabeled "{}" without also relabeling '
+                                  "the existing variable of the same name").format(v))
+
+        shared = old_labels & new_labels
+        if shared:
+            old_to_intermediate, intermediate_to_new = resolve_label_conflict(mapping, old_labels, new_labels)
+
+            self.relabel_variables(old_to_intermediate, inplace=True)
+            self.relabel_variables(intermediate_to_new, inplace=True)
+            return self
+
+        for oldterm, bias in list(self.items()):
+            newterm = frozenset((mapping.get(v, v) for v in oldterm))
+
+            if newterm != oldterm:
+                self[newterm] = bias
+                del self[oldterm]
+
+        return self
+
+    def normalize(self, bias_range=1, poly_range=None, ignored_terms=None):
+        """Normalizes the biases of the binary polynomial such that they fall in
+        the provided range(s).
+
+        If `poly_range` is provided, then `bias_range` will be treated as
+        the range for the linear biases and `poly_range` will be used for
+        the range of the other biases.
+
+        Args:
+            bias_range (number/pair):
+                Value/range by which to normalize the all the biases, or if
+                `poly_range` is provided, just the linear biases.
+
+            poly_range (number/pair, optional):
+                Value/range by which to normalize the higher order biases.
+
+            ignored_terms (iterable, optional):
+                Biases associated with these terms are not scaled.
+
+        """
+
+        def parse_range(r):
+            if isinstance(r, Number):
+                return -abs(r), abs(r)
+            return r
+
+        if ignored_terms is None:
+            ignored_terms = set()
+        else:
+            ignored_terms = {asfrozenset(term) for term in ignored_terms}
+
+        if poly_range is None:
+            linear_range, poly_range = bias_range, bias_range
+        else:
+            linear_range = bias_range
+
+        lin_range, poly_range = map(parse_range, (linear_range, poly_range))
+
+        # determine the current ranges for linear, higherorder
+        lmin = lmax = 0
+        pmin = pmax = 0
+        for term, bias in self.items():
+
+            if term in ignored_terms:
+                # we don't use the ignored terms to calculate the scaling
+                continue
+
+            if len(term) == 1:
+                lmin = min(bias, lmin)
+                lmax = max(bias, lmax)
+            elif len(term) > 1:
+                pmin = min(bias, pmin)
+                pmax = max(bias, pmax)
+
+        inv_scalar = max(lmin / lin_range[0], lmax / lin_range[1],
+                         pmin / poly_range[0], pmax / poly_range[1])
+
+        if inv_scalar != 0:
+            self.scale(1 / inv_scalar, ignored_terms=ignored_terms)
+
+    def scale(self, scalar, ignored_terms=None):
+        """Multiply the polynomial by the given scalar.
+
+        Args:
+            scalar (number):
+                Value to multiply the polynomial by.
+
+            ignored_terms (iterable, optional):
+                Biases associated with these terms are not scaled.
+
+        """
+
+        if ignored_terms is None:
+            ignored_terms = set()
+        else:
+            ignored_terms = {asfrozenset(term) for term in ignored_terms}
+
+        for term in self:
+            if term not in ignored_terms:
+                self[term] *= scalar
+
+    @classmethod
+    def from_hising(cls, h, J, offset=None):
+        """Construct a binary polynomial from a higher-order Ising problem.
+
+        Args:
+            h (dict):
+                The linear biases.
+
+            J (dict):
+                The higher-order biases.
+
+            offset (optional, default=0.0):
+                Constant offset applied to the model.
+
+        Returns:
+            :obj:`.BinaryPolynomial`
+
+        Examples:
+            >>> poly = dimod.BinaryPolynomial.from_hising({'a': 2}, {'ab': -1}, 0)
+
+        """
+        poly = {(k,): v for k, v in h.items()}
+        poly.update(J)
+        if offset is not None:
+            poly[frozenset([])] = offset
+        return cls(poly, Vartype.SPIN)
+
+    def to_hising(self):
+        """Construct a higher-order Ising problem from a binary polynomial.
+
+        Returns:
+            tuple: A 3-tuple of the form (`h`, `J`, `offset`) where `h` includes
+            the linear biases, `J` has the higher-order biases and `offset` is
+            the linear offset.
+
+        Examples:
+            >>> poly = dimod.BinaryPolynomial({'a': -1, 'ab': 1, 'abc': -1})
+            >>> h, J, off = poly.to_hising()
+            >>> h
+            {'a': -1}
+
+        """
+        if self.vartype is Vartype.BINARY:
+            return self.to_spin().to_hising()
+
+        h = {}
+        J = {}
+        offset = 0
+        for term, bias in self.items():
+            if len(term) == 0:
+                offset += bias
+            elif len(term) == 1:
+                v, = term
+                h[v] = bias
+            else:
+                J[tuple(term)] = bias
+
+        return h, J, offset
+
+    @classmethod
+    def from_hubo(cls, H, offset=None):
+        """Construct a binary polynomial from a higher-order unconstrained
+        binary optimization (HUBO) problem.
+
+        Args:
+            H (dict):
+                Coefficients of a higher-order unconstrained binary optimization
+                (HUBO) model.
+
+        Returns:
+            :obj:`.BinaryPolynomial`
+
+        Examples:
+            >>> poly = dimod.BinaryPolynomial.from_hubo({('a', 'b', 'c'): -1})
+
+        """
+        poly = cls(H, Vartype.BINARY)
+        if offset is not None:
+            poly[()] = poly.get((), 0) + offset
+        return poly
+
+    def to_hubo(self):
+        """Construct a higher-order unconstrained binary optimization (HUBO)
+        problem from a binary polynomial.
+
+        Returns:
+            tuple: A 2-tuple of the form (`H`, `offset`) where `H` is the HUBO
+            and `offset` is the linear offset.
+
+        """
+        if self.vartype is Vartype.SPIN:
+            return self.to_binary().to_hubo()
+
+        H = {tuple(term): bias for term, bias in self.items() if term}
+        offset = self[tuple()] if tuple() in self else 0
+        return H, offset
+
+    def to_binary(self, copy=False):
+        """Return a binary polynomial over `{0, 1}` variables.
+
+        Args:
+            copy (optional, default=False):
+                If True, the returned polynomial is always a copy. Otherwise,
+                if the polynomial is binary-valued already it returns itself.
+
+        Returns:
+            :obj:`.BinaryPolynomial`
+
+        """
+        if self.vartype is Vartype.BINARY:
+            if copy:
+                return self.copy()
+            else:
+                return self
+
+        new = BinaryPolynomial({}, Vartype.BINARY)
+
+        # s = 2x - 1
+        for term, bias in self.items():
+            for t in map(frozenset, powerset(term)):
+                newbias = bias * 2**len(t) * (-1)**(len(term) - len(t))
+
+                if t in new:
+                    new[t] += newbias
+                else:
+                    new[t] = newbias
+
+        return new
+
+    def to_spin(self, copy=False):
+        """Return a binary polynomial over `{-1, +1}` variables.
+
+        Args:
+            copy (optional, default=False):
+                If True, the returned polynomial is always a copy. Otherwise,
+                if the polynomial is spin-valued already it returns itself.
+
+        Returns:
+            :obj:`.BinaryPolynomial`
+
+        """
+        if self.vartype is Vartype.SPIN:
+            if copy:
+                return self.copy()
+            else:
+                return self
+
+        new = BinaryPolynomial({}, Vartype.SPIN)
+
+        # x = (s + 1) / 2
+        for term, bias in self.items():
+            newbias = bias / (2**len(term))
+
+            for t in map(frozenset, powerset(term)):
+                if t in new:
+                    new[t] += newbias
+                else:
+                    new[t] = newbias
+
+        return new
+
+
+def powerset(iterable):
+    return itertools.chain.from_iterable(itertools.combinations(iterable, r)
+                                         for r in range(len(iterable)+1))
diff --git a/dimod/higherorder.py b/dimod/higherorder/utils.py
similarity index 76%
rename from dimod/higherorder.py
rename to dimod/higherorder/utils.py
index b8909a41..5b46d050 100644
--- a/dimod/higherorder.py
+++ b/dimod/higherorder/utils.py
@@ -12,15 +12,19 @@
 #    See the License for the specific language governing permissions and
 #    limitations under the License.
 #
-# ================================================================================================
+# =============================================================================
 
 import itertools
+import warnings
+
 from collections import Counter
 
 import numpy as np
+
 from six import iteritems
 
 from dimod.binary_quadratic_model import BinaryQuadraticModel
+from dimod.higherorder.polynomial import BinaryPolynomial
 from dimod.sampleset import as_samples
 from dimod.vartypes import Vartype
 
@@ -194,68 +198,6 @@ def _reduce_degree(bqm, poly, vartype, scale):
     return _reduce_degree(bqm, new_poly, vartype, scale)
 
 
-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 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:
-        val *= v
-    return val
-
-
-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.
 
@@ -272,17 +214,13 @@ def poly_energy(sample_like, poly):
         float: The energy of the sample.
 
     """
-    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())
+
+    msg = ("poly_energy is deprecated and will be removed in dimod 0.9.0."
+           "In the future, use BinaryPolynomial.energy")
+    warnings.warn(msg, DeprecationWarning)
+    # dev note the vartype is not used in the energy calculation and this will
+    # be deprecated in the future
+    return BinaryPolynomial(poly, 'SPIN').energy(sample_like)
 
 
 def poly_energies(samples_like, poly):
@@ -294,7 +232,7 @@ def poly_energies(samples_like, poly):
             NumPy's array_like structure. See :func:`.as_samples`.
 
         poly (dict):
-            Polynomial as a dict of form {term: bias, ...}, where `term` is a 
+            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).
@@ -303,15 +241,9 @@ def poly_energies(samples_like, poly):
         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 np.sum([_prod_d([sample[:, labeldict[v]] for v in variables],
-                       num_samples) * bias for variables, bias in poly.items()],
-                  axis=0)
-
-
-def check_isin(key, key_list):
-    return sum(set(key) == set(key_tmp) for key_tmp in key_list)
+    msg = ("poly_energies is deprecated and will be removed in dimod 0.9.0."
+           "In the future, use BinaryPolynomial.energies")
+    warnings.warn(msg, DeprecationWarning)
+    # dev note the vartype is not used in the energy calculation and this will
+    # be deprecated in the future
+    return BinaryPolynomial(poly, 'SPIN').energies(samples_like)
diff --git a/dimod/reference/composites/__init__.py b/dimod/reference/composites/__init__.py
index c181a76a..1af76669 100644
--- a/dimod/reference/composites/__init__.py
+++ b/dimod/reference/composites/__init__.py
@@ -15,7 +15,7 @@
 # ================================================================================================
 
 from dimod.reference.composites.fixedvariable import *
-from dimod.reference.composites.higherordercomposite import *
+from dimod.reference.composites.higherordercomposites import *
 from dimod.reference.composites.scalecomposite import *
 from dimod.reference.composites.spin_transform import *
 from dimod.reference.composites.structure import *
diff --git a/dimod/reference/composites/higherordercomposite.py b/dimod/reference/composites/higherordercomposite.py
deleted file mode 100644
index 7fff7ca4..00000000
--- a/dimod/reference/composites/higherordercomposite.py
+++ /dev/null
@@ -1,231 +0,0 @@
-# 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, _relabeled_poly
-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 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/reference/composites/higherordercomposites.py b/dimod/reference/composites/higherordercomposites.py
new file mode 100644
index 00000000..bbc44a2c
--- /dev/null
+++ b/dimod/reference/composites/higherordercomposites.py
@@ -0,0 +1,347 @@
+# 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.
+#
+# =============================================================================
+"""
+Composites that convert binary quadratic model samplers into polynomial samplers
+or that work with binary polynomials.
+
+Higher-order composites implement three sampling methods (similar to
+:class:`.Sampler`):
+
+* :meth:`.PolySampler.sample_poly`
+* :meth:`.PolySampler.sample_hising`
+* :meth:`.PolySampler.sample_hubo`
+
+"""
+from __future__ import division
+
+import numpy as np
+
+from dimod.core.polysampler import ComposedPolySampler, PolySampler
+from dimod.higherorder.polynomial import BinaryPolynomial
+from dimod.higherorder.utils import make_quadratic, poly_energies
+from dimod.response import SampleSet
+
+__all__ = 'HigherOrderComposite', 'PolyScaleComposite'
+
+
+class HigherOrderComposite(ComposedPolySampler):
+    """Convert a binary quadratic model sampler to a binary polynomial sampler.
+
+    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())
+        >>> h = {0: -0.5, 1: -0.3, 2: -0.8}
+        >>> J = {(0, 1, 2): -1.7}
+        >>> sampleset = sampler.sample_hising(h, J, discard_unsatisfied=True)
+        >>> sampleset.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):
+        """A list containing the wrapped sampler."""
+        return self._children
+
+    @property
+    def parameters(self):
+        param = self.child.parameters.copy()
+        param['penalty_strength'] = []
+        param['discard_unsatisfied'] = []
+        param['keep_penalty_variables'] = []
+        return param
+
+    @property
+    def properties(self):
+        return {'child_properties': self.child.properties.copy()}
+
+    def sample_ising(self, h, J, offset=0, *args, **kwargs):
+        # need to handle offset input for backwards compatibility
+        if offset:
+            J[()] = offset
+        return ComposedPolySampler.sample_ising(self, h, J, *args, **kwargs)
+
+    def sample_poly(self, poly, penalty_strength=1.0,
+                    keep_penalty_variables=False,
+                    discard_unsatisfied=False, **parameters):
+        """Sample from the given binary polynomial.
+
+        Takes the given binary polynomial, introduces penalties, reduces the
+        higher-order problem into a quadratic problem and sends it to its child
+        sampler.
+
+        Args:
+            poly (:obj:`.BinaryPolynomial`):
+                A binary polynomial.
+
+            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`
+
+        """
+
+        bqm = make_quadratic(poly, penalty_strength, vartype=poly.vartype)
+        response = self.child.sample(bqm, **parameters)
+
+        return polymorph_response(response, poly, bqm,
+                                  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 polymorph_response(response, poly, bqm,
+                       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.
+
+        poly (:obj:`.BinaryPolynomial`):
+            A binary polynomial.
+
+        bqm (:obj:`dimod.BinaryQuadraticModel`): Binary quadratic model of the
+            reduced problem.
+
+        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)))
+
+    samples = record.sample[samples_to_keep]
+    energy_vector = poly.energies((samples, response.variables))
+
+    if not keep_penalty_variables:
+        original_variables = poly.variables
+        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)
+
+
+class PolyScaleComposite(ComposedPolySampler):
+    """Composite to scale biases of a binary polynomial.
+
+    Args:
+        child (:obj:`.PolySampler`):
+            A binary polynomial sampler.
+
+    Examples:
+
+       >>> linear = {'a': -4.0, 'b': -4.0}
+       >>> quadratic = {('a', 'b'): 3.2, ('a', 'b', 'c'): 1}
+       >>> sampler = dimod.PolyScaleComposite(dimod.HigherOrderComposite(dimod.ExactSolver()))
+       >>> response = sampler.sample_hising(linear, quadratic, scalar=0.5,
+       ...                ignored_terms=[('a','b')])
+
+    """
+
+    def __init__(self, child):
+        if not isinstance(child, PolySampler):
+            raise TypeError("Child sampler must be a PolySampler")
+        self._children = [child]
+
+    @property
+    def children(self):
+        """The child sampler in a list"""
+        return self._children
+
+    @property
+    def parameters(self):
+        param = self.child.parameters.copy()
+        param.update({'scalar': [],
+                      'bias_range': [],
+                      'poly_range': [],
+                      'ignored_terms': [],
+                      })
+        return param
+
+    @property
+    def properties(self):
+        return {'child_properties': self.child.properties.copy()}
+
+    def sample_poly(self, poly, scalar=None, bias_range=1, poly_range=None,
+                    ignored_terms=None, **parameters):
+        """Scale and sample from the given binary polynomial.
+
+        If scalar is not given, problem is scaled based on bias and polynomial
+        ranges. See :meth:`.BinaryPolynomial.scale` and
+        :meth:`.BinaryPolynomial.normalize`
+
+        Args:
+            poly (obj:`.BinaryPolynomial`): A binary polynomial.
+
+            scalar (number, optional):
+                Value by which to scale the energy range of the binary polynomial.
+
+            bias_range (number/pair, optional, default=1):
+                Value/range by which to normalize the all the biases, or if
+                `poly_range` is provided, just the linear biases.
+
+            poly_range (number/pair, optional):
+                Value/range by which to normalize the higher order biases.
+
+            ignored_terms (iterable, optional):
+                Biases associated with these terms are not scaled.
+
+            **parameters:
+                Other parameters for the sampling method, specified by
+                the child sampler.
+
+        """
+
+        if ignored_terms is None:
+            ignored_terms = set()
+        else:
+            ignored_terms = {frozenset(term) for term in ignored_terms}
+
+        # scale and normalize happen in-place so we need to make a copy
+        original, poly = poly, poly.copy()
+
+        if scalar is not None:
+            poly.scale(scalar, ignored_terms=ignored_terms)
+        else:
+            poly.normalize(bias_range=bias_range, poly_range=poly_range,
+                           ignored_terms=ignored_terms)
+
+            # we need to know how much we scaled by, which we can do by looking
+            # at the biases
+            try:
+                v = next(v for v, bias in original.items()
+                         if bias and v not in ignored_terms)
+            except StopIteration:
+                # nothing to scale
+                scalar = 1
+            else:
+                scalar = poly[v] / original[v]
+
+        sampleset = self.child.sample_poly(poly, **parameters)
+
+        if ignored_terms:
+            # we need to recalculate the energy
+            sampleset.record.energy = original.energies((sampleset.record.sample,
+                                                         sampleset.variables))
+        else:
+            sampleset.record.energy /= scalar
+
+        return sampleset
diff --git a/dimod/reference/composites/scalecomposite.py b/dimod/reference/composites/scalecomposite.py
index a2f49e28..e86a4d47 100644
--- a/dimod/reference/composites/scalecomposite.py
+++ b/dimod/reference/composites/scalecomposite.py
@@ -31,9 +31,8 @@ import numpy as np
 
 from dimod.binary_quadratic_model import BinaryQuadraticModel
 from dimod.core.composite import ComposedSampler
-from dimod.higherorder import poly_energies, _relabeled_poly, check_isin
 
-__all__ = ['ScaleComposite']
+__all__ = 'ScaleComposite',
 
 
 class ScaleComposite(ComposedSampler):
@@ -179,35 +178,40 @@ class ScaleComposite(ComposedSampler):
 
         """
 
-        # if quadratic, create a bqm and send to sample
-
-        if max(map(len, J.keys())) == 2:
-            bqm = BinaryQuadraticModel.from_ising(h, J, offset=offset)
-            return self.sample(bqm, scalar=scalar,
-                               bias_range=bias_range,
-                               quadratic_range=quadratic_range,
-                               ignored_variables=ignored_variables,
-                               ignored_interactions=ignored_interactions,
-                               ignore_offset=ignore_offset, **parameters)
-
-        # handle HUBO
-
-        ignored_variables, ignored_interactions = _check_params(
-            ignored_variables,
-            ignored_interactions)
-
-        child = self.child
-        h_sc, J_sc, offset_sc = _scaled_hubo(h, J, offset, scalar, bias_range,
-                                             quadratic_range, ignored_variables,
-                                             ignored_interactions,
-                                             ignore_offset)
-        response = child.sample_ising(h_sc, J_sc, offset=offset_sc,
-                                      **parameters)
-
-        poly = _relabeled_poly(h, J, response.variables.index)
-        response.record.energy = np.add(poly_energies(response.record.sample,
-                                                      poly), offset)
-        return response
+        if any(len(inter) > 2 for inter in J):
+            # handle HUBO
+            import warnings
+            msg = ("Support for higher order Ising models in ScaleComposite is "
+                   "deprecated and will be removed in dimod 0.9.0. Please use "
+                   "PolyScaleComposite.sample_hising instead.")
+            warnings.warn(msg, DeprecationWarning)
+
+            from dimod.reference.composites.higherordercomposites import PolyScaleComposite
+            from dimod.higherorder.polynomial import BinaryPolynomial
+
+            poly = BinaryPolynomial.from_hising(h, J, offset=offset)
+
+            ignored_terms = set()
+            if ignored_variables is not None:
+                ignored_terms.update(frozenset(v) for v in ignored_variables)
+            if ignored_interactions is not None:
+                ignored_terms.update(frozenset(inter) for inter in ignored_interactions)
+            if ignore_offset:
+                ignored_terms.add(frozenset())
+
+            return PolyScaleComposite(self.child).sample_poly(poly, scalar=scalar,
+                                                              bias_range=bias_range,
+                                                              poly_range=quadratic_range,
+                                                              ignored_terms=ignored_terms,
+                                                              **parameters)
+
+        bqm = BinaryQuadraticModel.from_ising(h, J, offset=offset)
+        return self.sample(bqm, scalar=scalar,
+                           bias_range=bias_range,
+                           quadratic_range=quadratic_range,
+                           ignored_variables=ignored_variables,
+                           ignored_interactions=ignored_interactions,
+                           ignore_offset=ignore_offset, **parameters)
 
 
 def _scale_back_response(bqm, response, scalar, ignored_interactions,
@@ -296,40 +300,5 @@ def _scaled_bqm(bqm, scalar, bias_range, quadratic_range,
     return bqm_copy
 
 
-def _scaled_hubo(h, j, offset, scalar, bias_range,
-                 quadratic_range,
-                 ignored_variables,
-                 ignored_interactions,
-                 ignore_offset):
-    """Helper function of sample_ising for scaling"""
-
-    if scalar is None:
-        scalar = _calc_norm_coeff(h, j, bias_range, quadratic_range,
-                                  ignored_variables, ignored_interactions)
-    h_sc = dict(h)
-    j_sc = dict(j)
-    offset_sc = offset
-    if not isinstance(scalar, Number):
-        raise TypeError("expected scalar to be a Number")
-
-    if scalar != 1:
-        if ignored_variables is None or ignored_interactions is None:
-            raise ValueError('ignored interactions or variables cannot be None')
-        j_sc = {}
-        for u, v in j.items():
-            if u in ignored_interactions:
-                j_sc[u] = v
-            else:
-                j_sc[u] = v * scalar
-
-        if not ignore_offset:
-            offset_sc = offset * scalar
-
-        h_sc = {}
-        for k, v in h.items():
-            if k in ignored_variables:
-                h_sc[k] = v
-            else:
-                h_sc[k] = v * scalar
-
-    return h_sc, j_sc, offset_sc
+def check_isin(key, key_list):
+    return sum(set(key) == set(key_tmp) for key_tmp in key_list)
diff --git a/docs/reference/bqm_functions.rst b/docs/reference/bqm_functions.rst
index 29f93774..5a8260d3 100644
--- a/docs/reference/bqm_functions.rst
+++ b/docs/reference/bqm_functions.rst
@@ -8,29 +8,6 @@ Functional interface to BQM methods and assorted utilities.
 
 .. currentmodule:: dimod
 
-Embedding
-=========
-
-See :ref:`embedding`.
-
-.. autosummary::
-   :toctree: generated/
-
-   embed_bqm
-
-
-Higher-order Polynomials
-========================
-
-.. autosummary::
-   :toctree: generated/
-
-   make_quadratic
-   poly_energy
-   poly_energies
-   create_poly
-
-
 Roof-duality
 ============
 
diff --git a/docs/reference/higherorder.rst b/docs/reference/higherorder.rst
new file mode 100644
index 00000000..4bba9d6d
--- /dev/null
+++ b/docs/reference/higherorder.rst
@@ -0,0 +1,104 @@
+.. _higherorder:
+
+=============================
+Solving Higher-order Problems
+=============================
+
+Sometimes it is nice to work with problems that are not restricted to quadratic
+interactions.
+
+Binary Polynomials
+==================
+
+.. automodule:: dimod.higherorder.polynomial
+
+.. autoclass:: BinaryPolynomial
+
+Methods
+~~~~~~~
+
+.. autosummary::
+   :toctree: generated/
+
+   BinaryPolynomial.copy
+   BinaryPolynomial.energies
+   BinaryPolynomial.energy
+   BinaryPolynomial.from_hising
+   BinaryPolynomial.from_hubo
+   BinaryPolynomial.normalize
+   BinaryPolynomial.relabel_variables
+   BinaryPolynomial.scale
+   BinaryPolynomial.to_binary
+   BinaryPolynomial.to_hising
+   BinaryPolynomial.to_hubo
+   BinaryPolynomial.to_spin
+
+Reducing to a Binary Quadratic Model
+====================================
+
+.. automodule:: dimod.higherorder.utils
+
+.. autosummary::
+   :toctree: generated/
+
+   make_quadratic
+
+
+Higher-order Composites
+=======================
+
+.. automodule:: dimod.reference.composites.higherordercomposites
+
+.. currentmodule:: dimod.reference.composites.higherordercomposites
+
+HigherOrderComposite
+--------------------
+
+.. autoclass:: HigherOrderComposite
+
+Properties
+~~~~~~~~~~
+
+.. autosummary::
+   :toctree: generated/
+
+   HigherOrderComposite.child
+   HigherOrderComposite.children
+   HigherOrderComposite.parameters
+   HigherOrderComposite.properties
+
+Methods
+~~~~~~~
+
+.. autosummary::
+   :toctree: generated/
+
+   HigherOrderComposite.sample_poly
+   HigherOrderComposite.sample_hising
+   HigherOrderComposite.sample_hubo
+
+PolyScaleComposite
+------------------
+
+.. autoclass:: PolyScaleComposite
+
+Properties
+~~~~~~~~~~
+
+.. autosummary::
+   :toctree: generated/
+
+   PolyScaleComposite.child
+   PolyScaleComposite.children
+   PolyScaleComposite.parameters
+   PolyScaleComposite.properties
+
+Methods
+~~~~~~~
+
+.. autosummary::
+   :toctree: generated/
+
+   PolyScaleComposite.sample_poly
+   PolyScaleComposite.sample_hising
+   PolyScaleComposite.sample_hubo
diff --git a/docs/reference/index.rst b/docs/reference/index.rst
index 9b3d2626..3b1fbd5c 100644
--- a/docs/reference/index.rst
+++ b/docs/reference/index.rst
@@ -16,6 +16,7 @@ Reference Documentation
    response
    vartypes
    bqm_functions
+   higherorder
    exceptions
    utilities
    testing
diff --git a/docs/reference/ref_samplers_composites.rst b/docs/reference/ref_samplers_composites.rst
index c61c3a19..fd6aeb91 100644
--- a/docs/reference/ref_samplers_composites.rst
+++ b/docs/reference/ref_samplers_composites.rst
@@ -114,25 +114,6 @@ 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
 ------------------------
 
diff --git a/docs/reference/samplers.rst b/docs/reference/samplers.rst
index a8a0e797..0f89f055 100644
--- a/docs/reference/samplers.rst
+++ b/docs/reference/samplers.rst
@@ -45,6 +45,16 @@ ABCs.
           :meth:`~.Sampler.sample`, :meth:`~.Sampler.sample_ising`, :meth:`~.Sampler.sample_qubo`
         - :meth:`~.Sampler.sample`, :meth:`~.Sampler.sample_ising`, :meth:`~.Sampler.sample_qubo`,
           :attr:`~.Composite.child`
+    *   - :class:`.PolySampler`
+        -
+        - :attr:`~.PolySampler.parameters`, :attr:`~.PolySampler.properties`
+        - :meth:`~.PolySampler.sample_poly`
+        - :meth:`~.PolySampler.sample_hising`, :meth:`~.PolySampler.sample_hubo`
+    *   - :class:`.ComposedPolySampler`
+        - :class:`.PolySampler`, :class:`.Composite`
+        - :attr:`~.PolySampler.parameters`, :attr:`~.PolySampler.properties`, :attr:`~.Composite.children`
+        - :meth:`~.Sampler.sample_poly`
+        - :meth:`~.Sampler.sample_hising`, :meth:`~.Sampler.sample_hubo`, :attr:`~.Composite.child`
 
 The table shows, for example, that the :class:`.Sampler` class requires that you implement
 the :attr:`~.Sampler.parameters` and :attr:`~.Sampler.properties` properties and at least
@@ -141,3 +151,43 @@ Mixin Properties
 
    Structured.adjacency
    Structured.structure
+
+
+Creating a Binary Polynomial Sampler
+====================================
+
+.. automodule:: dimod.core.polysampler
+
+.. currentmodule:: dimod
+.. autoclass:: PolySampler
+
+Abstract Properties
+-------------------
+
+.. autosummary::
+   :toctree: generated/
+
+   PolySampler.parameters
+   PolySampler.properties
+
+Abstract Methods
+----------------
+
+.. autosummary::
+   :toctree: generated/
+
+   PolySampler.sample_poly
+
+Mixin Methods
+-------------
+
+.. autosummary::
+   :toctree: generated/
+
+   PolySampler.sample_hising
+   PolySampler.sample_hubo
+
+Creating a Composed Binary Polynomial Sampler
+=============================================
+
+.. autoclass:: ComposedPolySampler
diff --git a/setup.py b/setup.py
index 2ac94a4b..6a2efd2f 100644
--- a/setup.py
+++ b/setup.py
@@ -46,6 +46,7 @@ extras_require = {'all': ['networkx>=2.0,<3.0',
 packages = ['dimod',
             'dimod.core',
             'dimod.generators',
+            'dimod.higherorder',
             'dimod.reference',
             'dimod.reference.composites',
             'dimod.reference.samplers',


Test patch:
diff --git a/tests/test_higherorder.py b/tests/test_higherorder.py
index 1dbcb381..c4cd80d1 100644
--- a/tests/test_higherorder.py
+++ b/tests/test_higherorder.py
@@ -24,7 +24,7 @@ class TestMakeQuadratic(unittest.TestCase):
 
     def test__spin_prod(self):
 
-        bqm = dimod.higherorder._spin_product(['a', 'b', 'p', 'aux'])
+        bqm = dimod.higherorder.utils._spin_product(['a', 'b', 'p', 'aux'])
 
         for v in ['a', 'b', 'p', 'aux']:
             self.assertIn(v, bqm)
@@ -41,7 +41,7 @@ class TestMakeQuadratic(unittest.TestCase):
 
         variables = ['a', 'b', 'p']
 
-        bqm = dimod.higherorder._binary_product(variables)
+        bqm = dimod.higherorder.utils._binary_product(variables)
 
         for v in variables:
             self.assertIn(v, bqm)
@@ -225,11 +225,9 @@ class TestMakeQuadratic(unittest.TestCase):
             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)
+        poly = dimod.BinaryPolynomial.from_hising(linear, j)
         samples = [[1, 1, -1], [1, -1, 1], [1, 1, 1], [-1, 1, -1]]
 
         en = poly_energies(samples, poly)
diff --git a/tests/test_higherorder_polynomial.py b/tests/test_higherorder_polynomial.py
new file mode 100644
index 00000000..504cf5f4
--- /dev/null
+++ b/tests/test_higherorder_polynomial.py
@@ -0,0 +1,268 @@
+# 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 unittest
+
+import numpy as np
+
+import dimod
+from dimod.higherorder import BinaryPolynomial
+
+
+class TestConstruction(unittest.TestCase):
+    # just that things don't fall down, we'll test correctness when
+    # testing other attributes
+    def test_from_dict(self):
+        BinaryPolynomial({'a': -1, tuple(): 1.3, 'bc': -1, ('a', 'b'): 1}, 'SPIN')
+
+    def test_from_iterator(self):
+        BinaryPolynomial(((term, -1) for term in itertools.combinations(range(100), 2)), 'SPIN')
+
+    def test_aggregation(self):
+        poly = BinaryPolynomial({'ab': 1, 'ba': 1, ('a', 'b'): 1, ('b', 'a'): 1}, 'SPIN')
+        self.assertEqual(poly, BinaryPolynomial({'ab': 4}, 'SPIN'))
+
+    def test_squares_binary(self):
+        poly = BinaryPolynomial({'aa': -1}, dimod.BINARY)
+        self.assertEqual(poly['a'], -1)
+
+    def test_squares_spin(self):
+        poly = BinaryPolynomial({'aa': -1}, dimod.SPIN)
+        self.assertEqual(poly[()], -1)
+
+    def test_cubes_spin(self):
+        poly = BinaryPolynomial({'aaa': -1}, dimod.SPIN)
+        self.assertEqual(poly['a'], -1)
+
+
+class Test__contains__(unittest.TestCase):
+    def test_single_term(self):
+        poly = BinaryPolynomial({('a', 'b'): 1}, 'SPIN')
+        self.assertIn('ab', poly)
+        self.assertIn('ba', poly)
+        self.assertIn(('a', 'b'), poly)
+        self.assertIn(('b', 'a'), poly)
+
+
+class Test__del__(unittest.TestCase):
+    def test_typical(self):
+        poly = BinaryPolynomial({'ab': -1, 'a': 1, 'b': 1}, 'BINARY')
+        self.assertEqual(len(poly), 3)
+        del poly[('a', 'b')]
+        self.assertEqual(len(poly), 2)
+
+
+class Test__eq__(unittest.TestCase):
+    def test_unlike_types(self):
+        polydict = {'ab': -1, 'a': 1, 'b': 1}
+        self.assertEqual(BinaryPolynomial(polydict, 'SPIN'), polydict)
+        self.assertNotEqual(BinaryPolynomial(polydict, 'SPIN'), 1)
+
+
+class Test__len__(unittest.TestCase):
+    def test_single_term(self):
+        poly = BinaryPolynomial({('a', 'b'): 1}, 'SPIN')
+        self.assertEqual(len(poly), 1)
+
+    def test_repeated_term(self):
+        poly = BinaryPolynomial({('a', 'b'): 1, 'ba': 1}, 'BINARY')
+        self.assertEqual(len(poly), 1)
+
+
+class Test__getitems__(unittest.TestCase):
+    def test_repeated_term(self):
+        poly = BinaryPolynomial({'ab': 1, 'ba': 1, ('a', 'b'): 1, ('b', 'a'): 1}, 'BINARY')
+        self.assertEqual(poly['ab'], 4)
+
+
+class Test_energies(unittest.TestCase):
+    def test_single_variable(self):
+        poly = BinaryPolynomial({'a': -1}, 'SPIN')
+
+        energies = poly.energies(([[-1], [1]], ['a']))
+        np.testing.assert_array_equal(energies, [1, -1])
+
+
+class TestDegree(unittest.TestCase):
+    def test_empty(self):
+        poly = BinaryPolynomial([], 'BINARY')
+        self.assertEqual(poly.degree, 0)
+
+    def test_degree0(self):
+        poly = BinaryPolynomial.from_hising({}, {}, 0)
+        self.assertEqual(poly.degree, 0)
+
+    def test_degree3(self):
+        poly = BinaryPolynomial.from_hubo({'abc': -1}, 0)
+        self.assertEqual(poly.degree, 3)
+
+
+class TestScale(unittest.TestCase):
+    def test_single_variable(self):
+        poly = BinaryPolynomial({'a': -1}, 'SPIN')
+        poly.scale(.5)
+        self.assertEqual(poly['a'], -.5)
+
+    def test_typical(self):
+        poly = BinaryPolynomial({'a': 1, 'ab': 1, '': 1}, 'BINARY')
+        poly.scale(2)
+        self.assertEqual(poly['a'], 2)
+        self.assertEqual(poly['ba'], 2)
+        self.assertEqual(poly[tuple()], 2)
+
+    def test_ignore_terms(self):
+        poly = BinaryPolynomial({'a': 1, 'ab': 1, '': 1}, 'BINARY')
+        poly.scale(2, ignored_terms=['', 'ba'])
+        self.assertEqual(poly['a'], 2)
+        self.assertEqual(poly['ba'], 1)
+        self.assertEqual(poly[tuple()], 1)
+
+    def test_scale_by_float(self):
+        poly = BinaryPolynomial({'a': 4}, 'SPIN')
+        poly.scale(.25)
+        self.assertEqual(poly['a'], 1)
+
+
+class TestNormalize(unittest.TestCase):
+    def test_empty(self):
+        poly = BinaryPolynomial({}, 'SPIN')
+        poly.normalize()
+
+    def test_typical(self):
+        poly = BinaryPolynomial({'a': 1, 'ab': 1, '': 1}, 'BINARY')
+        poly.normalize(.5)
+        self.assertEqual(poly['a'], .5)
+        self.assertEqual(poly['ba'], .5)
+        self.assertEqual(poly[tuple()], .5)
+
+    def test_int_division(self):
+        # we use division which previous caused issues in python2
+        poly = BinaryPolynomial({'a': 4}, 'SPIN')
+        poly.normalize(bias_range=1, poly_range=None, ignored_terms=[])
+        self.assertEqual(poly['a'], 1)
+
+    def test_ignored(self):
+        poly = BinaryPolynomial({'a': 4, 'b': 2}, 'SPIN')
+        poly.normalize(bias_range=1, poly_range=None, ignored_terms=['a'])
+        self.assertEqual(poly['a'], 4)  # should not be scaled or used in calcs
+        self.assertEqual(poly['b'], 1)  # should be normalized correctly
+
+
+class TestRelabel(unittest.TestCase):
+    def test_swap_partial(self):
+        poly = BinaryPolynomial({'abc': 1, 'bc': 1, 'ab': -1}, 'BINARY')
+        poly.relabel_variables({'b': 'c', 'c': 'b'})
+        self.assertEqual(poly, {'abc': 1, 'bc': 1, 'ac': -1})
+
+    def test_copy(self):
+        poly = BinaryPolynomial({'abc': 1, 'bc': 1, 'ab': -1}, 'BINARY')
+        new = poly.relabel_variables({'a': 'z'}, inplace=False)
+        self.assertEqual(new, {'zbc': 1, 'bc': 1, 'zb': -1})
+        self.assertEqual(poly, {'abc': 1, 'bc': 1, 'ab': -1})
+
+
+class Test_to_binary(unittest.TestCase):
+    def test_energy_equivalence_only_offset(self):
+        spipoly = BinaryPolynomial({'': 5}, 'SPIN')
+        binpoly = spipoly.to_binary()
+
+        variables = list(binpoly.variables)
+        for config in itertools.product((0, 1), repeat=len(variables)):
+            binary_sample = dict(zip(variables, config))
+            spin_sample = {v: 2*x - 1 for v, x in binary_sample.items()}
+            self.assertAlmostEqual(spipoly.energy(spin_sample), binpoly.energy(binary_sample))
+
+    def test_energy_equivalence_only_linear(self):
+        spipoly = BinaryPolynomial({'a': 5, 'b': -3}, 'SPIN')
+        binpoly = spipoly.to_binary()
+
+        variables = list(binpoly.variables)
+        for config in itertools.product((0, 1), repeat=len(variables)):
+            binary_sample = dict(zip(variables, config))
+            spin_sample = {v: 2*x - 1 for v, x in binary_sample.items()}
+            self.assertAlmostEqual(spipoly.energy(spin_sample), binpoly.energy(binary_sample))
+
+    def test_energy_equivalence(self):
+        spipoly = BinaryPolynomial({'abc': 1, 'bc': 1, 'ab': -1, '': 0}, 'SPIN')
+        binpoly = spipoly.to_binary()
+
+        variables = list(binpoly.variables)
+        for config in itertools.product((0, 1), repeat=len(variables)):
+            binary_sample = dict(zip(variables, config))
+            spin_sample = {v: 2*x - 1 for v, x in binary_sample.items()}
+            self.assertAlmostEqual(spipoly.energy(spin_sample), binpoly.energy(binary_sample))
+
+
+class Test_to_spin(unittest.TestCase):
+    def test_energy_equivalence_only_offset(self):
+        binpoly = BinaryPolynomial({'': 5}, 'BINARY')
+        spipoly = binpoly.to_spin()
+
+        variables = list(binpoly.variables)
+        for config in itertools.product((0, 1), repeat=len(variables)):
+            binary_sample = dict(zip(variables, config))
+            spin_sample = {v: 2*x - 1 for v, x in binary_sample.items()}
+            self.assertAlmostEqual(spipoly.energy(spin_sample), binpoly.energy(binary_sample))
+
+    def test_energy_equivalence_only_linear(self):
+        binpoly = BinaryPolynomial({'a': 5, 'b': -3}, 'BINARY')
+        spipoly = binpoly.to_spin()
+
+        variables = list(binpoly.variables)
+        for config in itertools.product((0, 1), repeat=len(variables)):
+            binary_sample = dict(zip(variables, config))
+            spin_sample = {v: 2*x - 1 for v, x in binary_sample.items()}
+            self.assertAlmostEqual(spipoly.energy(spin_sample), binpoly.energy(binary_sample))
+
+    def test_energy_equivalence(self):
+        binpoly = BinaryPolynomial({'abc': 1, 'bc': 1, 'ab': -1, '': 0}, 'BINARY')
+        spipoly = binpoly.to_spin()
+
+        variables = list(binpoly.variables)
+        for config in itertools.product((0, 1), repeat=len(variables)):
+            binary_sample = dict(zip(variables, config))
+            spin_sample = {v: 2*x - 1 for v, x in binary_sample.items()}
+            self.assertAlmostEqual(spipoly.energy(spin_sample), binpoly.energy(binary_sample))
+
+
+class TestHUBO(unittest.TestCase):
+    # .to_hubo and .from_hubo
+    def test_binary(self):
+        poly = BinaryPolynomial({'abc': 1, 'bc': 1, 'ab': -1, '': 0}, 'BINARY')
+        H, off = poly.to_hubo()
+        new = BinaryPolynomial.from_hubo(H, off)
+        self.assertEqual(poly, new)
+
+    def test_spin(self):
+        poly = BinaryPolynomial({'abc': 1, 'bc': 1, 'ab': -1, '': 0}, 'SPIN')
+        H, off = poly.to_hubo()
+        new = BinaryPolynomial.from_hubo(H, off)
+        self.assertEqual(poly, new.to_spin())
+
+
+class TestHising(unittest.TestCase):
+    # .to_hising and .from_hising
+    def test_binary(self):
+        poly = BinaryPolynomial({'abc': 1, 'bc': 1, 'ab': -1}, 'BINARY')
+        h, J, off = poly.to_hising()
+        new = BinaryPolynomial.from_hising(h, J, off)
+        self.assertEqual(poly, new.to_binary())
+
+    def test_spin(self):
+        poly = BinaryPolynomial({'abc': 1, 'bc': 1, 'ab': -1}, 'SPIN')
+        h, J, off = poly.to_hising()
+        new = BinaryPolynomial.from_hising(h, J, off)
+        self.assertEqual(poly, new)
diff --git a/tests/test_higherordercomposite.py b/tests/test_higherordercomposite.py
index 1b94c998..3712809a 100644
--- a/tests/test_higherordercomposite.py
+++ b/tests/test_higherordercomposite.py
@@ -24,11 +24,6 @@ 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}
@@ -63,7 +58,7 @@ class TestFixedVariableComposite(unittest.TestCase):
                                         keep_penalty_variables=True,
                                         discard_unsatisfied=True)
 
-        self.assertEqual(len(response.first.sample),5)
+        self.assertEqual(len(response.first.sample), 5)
         self.assertAlmostEqual(response.first.energy, -3.3)
         self.assertTrue(response.first.penalty_satisfaction)
 
diff --git a/tests/test_polyscalecomposite.py b/tests/test_polyscalecomposite.py
new file mode 100644
index 00000000..53213990
--- /dev/null
+++ b/tests/test_polyscalecomposite.py
@@ -0,0 +1,114 @@
+# 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 unittest
+
+import numpy as np
+
+import dimod
+
+from dimod import PolyScaleComposite, HigherOrderComposite, ExactSolver
+
+
+class RangeLimitedSampler(dimod.PolySampler):
+    parameters = None
+    properties = None
+
+    def sample_poly(self, poly, num_reads=1):
+
+        if any(bias > 1 for bias in poly.values()):
+            raise RuntimeError
+        if any(bias < -1 for bias in poly.values()):
+            raise RuntimeError
+
+        samples = np.ones((num_reads, len(poly.variables))), list(poly.variables)
+        sampleset = dimod.SampleSet.from_samples(samples, vartype=poly.vartype,
+                                                 energy=poly.energies(samples))
+        return sampleset
+
+
+class TestConstruction(unittest.TestCase):
+    def test_typical(self):
+        sampler = PolyScaleComposite(HigherOrderComposite(ExactSolver()))
+
+        self.assertTrue(hasattr(sampler, 'sample_poly'))
+        self.assertTrue(hasattr(sampler, 'sample_hising'))
+        self.assertTrue(hasattr(sampler, 'sample_hubo'))
+
+    def test_wrap_bqm(self):
+        with self.assertRaises(TypeError):
+            PolyScaleComposite(ExactSolver())
+
+
+class TestSampleHising(unittest.TestCase):
+    def test_all_zero(self):
+        sampler = PolyScaleComposite(RangeLimitedSampler())
+        sampler.sample_hising({'a': 0}, {'ab': 0, 'bc': 0, 'abc': 0})
+
+    def test_empty(self):
+        sampler = PolyScaleComposite(RangeLimitedSampler())
+        samples = sampler.sample_hising({}, {})
+        self.assertEqual(len(samples.variables), 0)
+
+    def test_normalizing(self):
+        sampler = PolyScaleComposite(RangeLimitedSampler())
+        samples = sampler.sample_hising({'a': 4}, {})
+        self.assertEqual(samples.first.energy, 4)
+
+    def test_scale(self):
+        sampler = PolyScaleComposite(RangeLimitedSampler())
+        samples = sampler.sample_hising({'a': 4}, {}, scalar=.25)
+        self.assertEqual(samples.first.energy, 4)
+
+    def test_fail_scale(self):
+        sampler = PolyScaleComposite(RangeLimitedSampler())
+        with self.assertRaises(RuntimeError):
+            sampler.sample_hising({'a': 4}, {}, scalar=1)
+
+    def test_all_energies(self):
+        sampler = PolyScaleComposite(HigherOrderComposite(ExactSolver()))
+
+        h = {'a': -1, 'b': 4}
+        J = {'abc': -1, 'ab': 1, 'aaa': .5}
+
+        sampleset = sampler.sample_hising(h, J, discard_unsatisfied=True)
+
+        for sample, energy in sampleset.data(['sample', 'energy']):
+            en = 0
+            for v, bias in h.items():
+                en += sample[v] * bias
+            for term, bias in J.items():
+                val = bias
+                for v in term:
+                    val *= sample[v]
+                en += val
+
+            self.assertAlmostEqual(energy, en)
+
+
+class TestSampleHubo(unittest.TestCase):
+    def test_empty(self):
+        sampler = PolyScaleComposite(RangeLimitedSampler())
+        samples = sampler.sample_hubo({})
+        self.assertEqual(len(samples.variables), 0)
+
+        sampler = PolyScaleComposite(RangeLimitedSampler())
+        samples = sampler.sample_hubo({'a': 4}, scalar=.25)
+        self.assertEqual(samples.first.energy, 4)
+
+    def test_fail_scale(self):
+        sampler = PolyScaleComposite(RangeLimitedSampler())
+        with self.assertRaises(RuntimeError):
+            sampler.sample_hubo({'a': 4}, scalar=1)
diff --git a/tests/test_scalecomposite.py b/tests/test_scalecomposite.py
index 754a6ae2..6a880e93 100644
--- a/tests/test_scalecomposite.py
+++ b/tests/test_scalecomposite.py
@@ -18,12 +18,53 @@ import unittest
 
 import dimod.testing as dtest
 from dimod import ExactSolver, ScaleComposite, HigherOrderComposite, \
-    BinaryQuadraticModel, Sampler
-from dimod.reference.composites.scalecomposite import _scaled_hubo, \
-    _check_params, _scaled_bqm
+    BinaryQuadraticModel, Sampler, PolySampler
+from dimod.reference.composites.scalecomposite import _check_params, _scaled_bqm, _calc_norm_coeff
+
+
+from numbers import Number
+
+
+def _scaled_hubo(h, j, offset, scalar, bias_range,
+                 quadratic_range,
+                 ignored_variables,
+                 ignored_interactions,
+                 ignore_offset):
+    """Helper function of sample_ising for scaling"""
+
+    if scalar is None:
+        scalar = _calc_norm_coeff(h, j, bias_range, quadratic_range,
+                                  ignored_variables, ignored_interactions)
+    h_sc = dict(h)
+    j_sc = dict(j)
+    offset_sc = offset
+    if not isinstance(scalar, Number):
+        raise TypeError("expected scalar to be a Number")
+
+    if scalar != 1:
+        if ignored_variables is None or ignored_interactions is None:
+            raise ValueError('ignored interactions or variables cannot be None')
+        j_sc = {}
+        for u, v in j.items():
+            if u in ignored_interactions:
+                j_sc[u] = v
+            else:
+                j_sc[u] = v * scalar
+
+        if not ignore_offset:
+            offset_sc = offset * scalar
+
+        h_sc = {}
+        for k, v in h.items():
+            if k in ignored_variables:
+                h_sc[k] = v
+            else:
+                h_sc[k] = v * scalar
 
+    return h_sc, j_sc, offset_sc
 
-class ScalingChecker(Sampler):
+
+class ScalingChecker(Sampler, PolySampler):
     def __init__(self, child_sampler, bqm=None, h=None, J=None, offset=0,
                  scalar=None, bias_range=1, quadratic_range=None,
                  ignored_variables=None, ignored_interactions=None,
@@ -61,6 +102,15 @@ class ScalingChecker(Sampler):
         assert self.offset == offset
         return self.child.sample_ising(h, J, offset=offset, **parameters)
 
+    def sample_poly(self, poly, **parameters):
+        h, J, offset = poly.to_hising()
+        assert self.h == h
+        selfJ = {frozenset(term): bias for term, bias in self.J.items()}
+        J = {frozenset(term): bias for term, bias in J.items()}
+        assert selfJ == J
+        assert self.offset == offset
+        return self.child.sample_poly(poly, **parameters)
+
     def parameters(self):
         return self.child.parameters()