diff --git a/dimod/binary/binary_quadratic_model.py b/dimod/binary/binary_quadratic_model.py
index 7c2270411..628ba674b 100644
--- a/dimod/binary/binary_quadratic_model.py
+++ b/dimod/binary/binary_quadratic_model.py
@@ -12,6 +12,19 @@
 #    See the License for the specific language governing permissions and
 #    limitations under the License.
 
+r"""Binary quadratic models (BQMs) are problems of the form:
+
+.. math::
+
+    E(\bf{v})
+    = \sum_{i} a_i v_i
+    + \sum_{i<j} b_{i,j} v_i v_j
+    + c
+    \qquad\qquad v_i \in\{-1,+1\} \text{  or } \{0,1\}
+
+where :math:`a_{i}, b_{ij}, c` are real values.
+"""
+
 from __future__ import annotations
 
 import collections.abc as abc
@@ -70,18 +83,6 @@
 class BinaryQuadraticModel(QuadraticViewsMixin):
     r"""Binary quadratic model.
 
-    Binary quadratic models (BQMs) are problems of the form:
-
-    .. math::
-
-        E(\bf{v})
-        = \sum_{i=1} a_i v_i
-        + \sum_{i<j} b_{i,j} v_i v_j
-        + c
-        \qquad\qquad v_i \in\{-1,+1\} \text{  or } \{0,1\}
-
-    where :math:`a_{i}, b_{ij}, c` are real values.
-
     This class encodes Ising and quadratic unconstrained binary optimization
     (QUBO) models used by samplers such as the D-Wave system.
 
diff --git a/dimod/constrained/constrained.py b/dimod/constrained/constrained.py
index a1cf191ce..cf9e91d72 100644
--- a/dimod/constrained/constrained.py
+++ b/dimod/constrained/constrained.py
@@ -12,8 +12,34 @@
 #    See the License for the specific language governing permissions and
 #    limitations under the License.
 
-"""
-Constrained Quadratic Model class.
+r"""Constrained quadratic models are problems of the form:
+
+.. math::
+
+    \begin{align}
+        \text{Minimize an objective:} & \\
+        & \sum_{i} a_i x_i + \sum_{i \le j} b_{ij} x_i x_j + c, \\
+        \text{Subject to constraints:} & \\
+        & \sum_i a_i^{(m)} x_i + \sum_{i \le j} b_{ij}^{(m)} x_i x_j+ c^{(m)} \circ 0,
+        \quad m=1, \dots, M,
+    \end{align}
+
+where :math:`\{ x_i\}_{i=1, \dots, N}` can be binary\ [#]_, integer, or continuous
+variables, :math:`a_{i}, b_{ij}, c` are real values,
+:math:`\circ \in \{ \ge, \le, = \}` and  :math:`M` is the total number of constraints.
+
+.. [#]
+    For binary variables, the range of the quadratic-term summation is
+    :math:`i < j` because :math:`x^2 = x` for binary values :math:`\{0, 1\}`
+    and :math:`s^2 = 1` for spin values :math:`\{-1, 1\}`.
+
+Constraints can be categorized as either "hard" or "soft". Any hard constraint
+must be satisfied for a solution of the model to qualify as feasible. Soft
+constraints may be violated to achieve an overall good solution. By setting
+appropriate weights to soft constraints in comparison to the objective
+and to other soft constraints, you can express the relative importance of such
+constraints.
+
 """
 
 from __future__ import annotations
@@ -137,100 +163,16 @@ def __len__(self):
 class ConstrainedQuadraticModel(cyConstrainedQuadraticModel):
     r"""A constrained quadratic model.
 
-    Constrained quadratic models are problems of the form:
-
-    .. math::
-
-        \begin{align}
-            \text{Minimize an objective:} & \\
-            & \sum_{i} a_i x_i + \sum_{i \le j} b_{ij} x_i x_j + c, \\
-            \text{Subject to constraints:} & \\
-            & \sum_i a_i^{(m)} x_i + \sum_{i \le j} b_{ij}^{(m)} x_i x_j+ c^{(m)} \circ 0,
-            \quad m=1, \dots, M,
-        \end{align}
-
-    where :math:`\{ x_i\}_{i=1, \dots, N}` can be binary\ [#]_ or integer
-    variables, :math:`a_{i}, b_{ij}, c` are real values,
-    :math:`\circ \in \{ \ge, \le, = \}` and  :math:`M` is the total number of constraints.
-
-    .. [#]
-        For binary variables, the range of the quadratic-term summation is
-        :math:`i < j` because :math:`x^2 = x` for binary values :math:`\{0, 1\}`
-        and :math:`s^2 = 1` for spin values :math:`\{-1, 1\}`.
-
-    Constraints are often categorized as either "hard" or "soft". Any hard constraint
-    must be satisfied for a solution of the model to qualify as feasible. Soft
-    constraints may be violated to achieve an overall good solution. By setting
-    appropriate weights to soft constraints in comparison to the objective
-    and to other soft constraints, you can express the relative importance of such
-    constraints.
-
     The objective and constraints are encoded as either :class:`.QuadraticModel`
     or :class:`.BinaryQuadraticModel` depending on the variable types used.
 
     Example:
 
-        This example solves the simple `bin packing problem <https://w.wiki/3jz4>`_
-        of packing a set of items of different weights into the smallest
-        possible number of bins.
-
-        `dimod` provides a general :func:`~dimod.generators.random_bin_packing`
-        function to generate bin packing problems, and this example follows the
-        same naming conventions.
-
-        Consider four objects with weights between 0 and 1, and assume that each
-        bin has a capacity to hold up to a total weight of 1.
-
-        >>> weights = [.9, .7, .2, .1]
-        >>> capacity = 1
-
-        Variable :math:`y_j` indicates that bin :math:`j` is used. Clearly, no
-        more than four bins are needed.
-
-        >>> y = [dimod.Binary(f'y_{j}') for j in range(len(weights))]
-
-        Variable :math:`x_{i,j}` indicates that item :math:`i` is put in bin
-        :math:`j`.
-
-        >>> x = [[dimod.Binary(f'x_{i}_{j}') for j in range(len(weights))]
-        ...      for i in range(len(weights))]
-
         Create an empty constrained quadratic model ("empty" meaning that no
-        objective or constraints have set).
+        objective or constraints are set).
 
         >>> cqm = dimod.ConstrainedQuadraticModel()
 
-        The problem is to minimize the number of bins used. Therefore the objective
-        is to minimize the value of :math:`\sum_j y_j`.
-
-        >>> cqm.set_objective(sum(y))
-
-        Any feasible solution must meet the constraint that each item can only go
-        in one bin. You can express this constraint, for a given item :math:`i`,
-        with :math:`\sum_j x_{i, j} == 1`. Note that the label of each
-        constraint is returned so that you can access them in the future if
-        desired.
-
-        >>> for i in range(len(weights)):
-        ...     cqm.add_constraint(sum(x[i]) == 1, label=f'item_placing_{i}')
-        'item_placing_0'
-        'item_placing_1'
-        'item_placing_2'
-        'item_placing_3'
-
-        Finally, enforce the limits on each bin. You can express this constraint,
-        for a given bin :math:`j`, with :math:`\sum_i x_{i, j} * w_i <= c` where
-        :math:`w_i` is the weight of item :math:`i` and :math:`c` is the capacity.
-
-        >>> for j in range(len(weights)):
-        ...     cqm.add_constraint(
-        ...         sum(weights[i] * x[i][j] for i in range(len(weights))) - y[j] * capacity <= 0,
-        ...         label=f'capacity_bin_{j}')
-        'capacity_bin_0'
-        'capacity_bin_1'
-        'capacity_bin_2'
-        'capacity_bin_3'
-
     """
     def __init__(self):
         super().__init__()
diff --git a/dimod/include/dimod/abc.h b/dimod/include/dimod/abc.h
index c6e0c4f38..88f6f7560 100644
--- a/dimod/include/dimod/abc.h
+++ b/dimod/include/dimod/abc.h
@@ -140,13 +140,13 @@ class ConstQuadraticIterator {
 template <class Bias, class Index>
 class QuadraticModelBase {
  public:
-    /// The first template parameter (Bias).
+    /// The first template parameter (`Bias`).
     using bias_type = Bias;
 
-    /// The second template parameter (Index).
+    /// The second template parameter (`Index`).
     using index_type = Index;
 
-    /// Unsigned integral type that can represent non-negative values.
+    /// Unsigned integer type that can represent non-negative values.
     using size_type = std::size_t;
 
     /// @private  <- don't doc
@@ -176,8 +176,10 @@ class QuadraticModelBase {
     /// Add offset.
     void add_offset(bias_type bias);
 
+    /// Add interaction between variables `v` and `u`.
     void add_quadratic(index_type u, index_type v, bias_type bias);
 
+    /// Add interactions between row `row` and column `col`.
     void add_quadratic(std::initializer_list<index_type> row, std::initializer_list<index_type> col,
                        std::initializer_list<bias_type> biases);
 
@@ -186,7 +188,7 @@ class QuadraticModelBase {
                        index_type length);
 
     /**
-     * Add quadratic bias for the given variables at the end of eachother's neighborhoods.
+     * Add quadratic bias for the given variables at the end of each other's neighborhoods.
      *
      * # Parameters
      * - `u` - a variable.
@@ -216,14 +218,16 @@ class QuadraticModelBase {
     template <class T>
     void add_quadratic_from_dense(const T dense[], index_type num_variables);
 
+    /// Return an iterator to the beginning of the neighborhood of `v`.
     const_neighborhood_iterator cbegin_neighborhood(index_type v) const;
 
+    /// Return an iterator to the end of the neighborhood of `v`.
     const_neighborhood_iterator cend_neighborhood(index_type v) const;
 
-    /// todo
+    /// Return an iterator to the beginning of the quadratic interactions.
     const_quadratic_iterator cbegin_quadratic() const;
 
-    /// todo
+    /// Return an iterator to the end of the quadratic interactions.
     const_quadratic_iterator cend_quadratic() const;
 
     /// Remove the offset and all variables and interactions from the model.
@@ -232,7 +236,7 @@ class QuadraticModelBase {
     /**
      * Return the energy of the given sample.
      *
-     * The `sample_start` must be random access iterator pointing to the
+     * The `sample_start` must be a random access iterator pointing to the
      * beginning of the sample.
      *
      * The behavior of this function is undefined when the sample is not
@@ -250,14 +254,14 @@ class QuadraticModelBase {
     template <class T>
     void fix_variable(index_type v, T assignment);
 
-    /// Check whether u and v have an interaction
+    /// Check whether `u` and `v` have an interaction.
     bool has_interaction(index_type u, index_type v) const;
 
     /// Test whether two quadratic models are equal.
     template <class B, class I>
     bool is_equal(const QuadraticModelBase<B, I>& other) const;
 
-    /// Return True if the model has no quadratic biases.
+    /// Test whether the model has no quadratic biases.
     bool is_linear() const;
 
     /// The linear bias of variable `v`.
@@ -292,7 +296,7 @@ class QuadraticModelBase {
     /// Return the number of interactions in the quadratic model.
     size_type num_interactions() const;
 
-    /// The number of other variables `v` interacts with.
+    /// Return the number of other variables that `v` interacts with.
     size_type num_interactions(index_type v) const;
 
     /// Return the number of variables in the quadratic model.
@@ -302,24 +306,24 @@ class QuadraticModelBase {
     bias_type offset() const;
 
     /**
-     * Return the quadratic bias associated with `u`, `v`.
+     * Return the quadratic bias associated with `u` and `v`.
      *
      * If `u` and `v` do not have a quadratic bias, returns 0.
      *
-     * Note that this function does not return a reference, this is because
+     * Note that this function does not return a reference because
      * each quadratic bias is stored twice.
      *
      */
     bias_type quadratic(index_type u, index_type v) const;
 
     /**
-     * Return the quadratic bias associated with `u`, `v`.
+     * Return the quadratic bias associated with `u` and `v`.
      *
-     * Note that this function does not return a reference, this is because
+     * Note that this function does not return a reference because
      * each quadratic bias is stored twice.
      *
-     * Raises an `out_of_range` error if either `u` or `v` are not variables or
-     * if they do not have an interaction then the function throws an exception.
+     * Raises an `out_of_range` error if either `u` or `v` are not variables;
+     * if they do not have an interaction, the function throws an exception.
      */
     bias_type quadratic_at(index_type u, index_type v) const;
 
@@ -334,7 +338,7 @@ class QuadraticModelBase {
      */
     virtual void remove_variable(index_type v);
 
-    /// Multiply all biases 'scalar'
+    /// Multiply all biases by the value of `scalar`.
     void scale(bias_type scalar);
 
     /// Set the linear bias of variable `v`.
@@ -346,7 +350,7 @@ class QuadraticModelBase {
     /// Set the offset.
     void set_offset(bias_type offset);
 
-    /// Set the quadratic bias for the given variables.
+    /// Set the quadratic bias between variables `u` and `v`.
     void set_quadratic(index_type u, index_type v, bias_type bias);
 
     void substitute_variable(index_type v, bias_type multiplier, bias_type offset);
@@ -371,7 +375,7 @@ class QuadraticModelBase {
     /// Increase the size of the model by `n`. Returns the index of the first variable added.
     index_type add_variables(index_type n);
 
-    /// Return an empty neighborhood - useful when we don't have an adj
+    /// Return an empty neighborhood; useful when a variable does not have an adjacency.
     static const std::vector<OneVarTerm<bias_type, index_type>>& empty_neighborhood() {
         static std::vector<OneVarTerm<bias_type, index_type>> empty;
         return empty;
@@ -381,7 +385,7 @@ class QuadraticModelBase {
     void resize(index_type n);
 
     /// Protected version of vartype() that allows subclasses to distinguish
-    /// between the vartype_ called by mixin functions and the public API one.
+    /// between the `vartype_` called by mixin functions and the public API one.
     /// By default they are the same.
     virtual Vartype vartype_(index_type v) const { return vartype(v); }
 
diff --git a/dimod/include/dimod/binary_quadratic_model.h b/dimod/include/dimod/binary_quadratic_model.h
index 3d96f7436..6b412769e 100644
--- a/dimod/include/dimod/binary_quadratic_model.h
+++ b/dimod/include/dimod/binary_quadratic_model.h
@@ -19,23 +19,23 @@
 
 namespace dimod {
 
-/// A Binary Quadratic Model is a quadratic polynomial over binary variables.
+/// A binary quadratic model (BQM) is a quadratic polynomial over binary-valued variables.
 template <class Bias, class Index = int>
 class BinaryQuadraticModel : public abc::QuadraticModelBase<Bias, Index> {
  public:
-    /// The type of the base class.
+    /// Type of the base class.
     using base_type = abc::QuadraticModelBase<Bias, Index>;
 
-    /// The first template parameter (Bias).
+    /// First template parameter (`Bias`).
     using bias_type = Bias;
 
-    /// The second template parameter (Index).
+    /// Second template parameter (`Index`).
     using index_type = Index;
 
-    /// Unsigned integral that can represent non-negative values.
+    /// Unsigned integer that can represent non-negative values.
     using size_type = typename base_type::size_type;
 
-    /// Empty constructor. The vartype defaults to `Vartype::BINARY`.
+    /// Empty constructor. The variable type defaults to `Vartype::BINARY`.
     BinaryQuadraticModel();
 
     /// Create a BQM of the given `vartype`.
@@ -51,9 +51,9 @@ class BinaryQuadraticModel : public abc::QuadraticModelBase<Bias, Index> {
      *
      * Values on the diagonal are treated differently depending on the variable
      * type.
-     * If the BQM is SPIN-valued, then the values on the diagonal are
+     * If the BQM is SPIN-valued, values on the diagonal are
      * added to the offset.
-     * If the BQM is BINARY-valued, then the values on the diagonal are added
+     * If the BQM is BINARY-valued, values on the diagonal are added
      * as linear biases.
      *
      */
@@ -62,7 +62,7 @@ class BinaryQuadraticModel : public abc::QuadraticModelBase<Bias, Index> {
 
     using base_type::add_quadratic;
 
-    /*
+    /**
      * Construct a BQM from COO-formated iterators.
      *
      * A sparse BQM encoded in [COOrdinate] format is specified by three
@@ -70,7 +70,7 @@ class BinaryQuadraticModel : public abc::QuadraticModelBase<Bias, Index> {
      *
      * [COOrdinate]: https://w.wiki/n$L
      *
-     * `row_iterator` must be a random access iterator  pointing to the
+     * `row_iterator` must be a random access iterator pointing to the
      * beginning of the row data. `col_iterator` must be a random access
      * iterator pointing to the beginning of the column data. `bias_iterator`
      * must be a random access iterator pointing to the beginning of the bias
@@ -83,7 +83,7 @@ class BinaryQuadraticModel : public abc::QuadraticModelBase<Bias, Index> {
     /// Add one (disconnected) variable to the BQM and return its index.
     index_type add_variable();
 
-    /// Change the vartype of the binary quadratic model
+    /// Change the variable type of the BQM.
     void change_vartype(Vartype vartype);
 
     bias_type lower_bound() const;
diff --git a/dimod/include/dimod/constrained_quadratic_model.h b/dimod/include/dimod/constrained_quadratic_model.h
index f4120aa1b..8fdb15d9c 100644
--- a/dimod/include/dimod/constrained_quadratic_model.h
+++ b/dimod/include/dimod/constrained_quadratic_model.h
@@ -26,16 +26,17 @@
 
 namespace dimod {
 
+/// A constrained quadratic model (CQM) can include an objective and constraints as polynomials with one or two variables per term.
 template <class Bias, class Index = int>
 class ConstrainedQuadraticModel {
  public:
-    /// The first template parameter (Bias).
+    /// First template parameter (`Bias`).
     using bias_type = Bias;
 
-    /// The second template parameter (Index).
+    /// Second template parameter (`Index`).
     using index_type = Index;
 
-    /// Unsigned integral type that can represent non-negative values.
+    /// Unsigned integer type that can represent non-negative values.
     using size_type = std::size_t;
 
     ConstrainedQuadraticModel();
@@ -60,20 +61,27 @@ class ConstrainedQuadraticModel {
     /// If an exception is thrown, there are no changes to the model.
     index_type add_constraint(Constraint<bias_type, index_type> constraint);
 
+    /// Add `n` constraints.
     index_type add_constraints(index_type n);
 
+    /// Add a constraint with only linear coefficients.
     index_type add_linear_constraint(std::initializer_list<index_type> variables,
                                      std::initializer_list<bias_type> biases, Sense sense,
                                      bias_type rhs);
 
+    /// Add variable of type `vartype` with lower bound `lb` and upper bound `ub`.
     index_type add_variable(Vartype vartype, bias_type lb, bias_type ub);
 
+    /// Add variable of type `vartype`.
     index_type add_variable(Vartype vartype);
 
+    /// Add `n` variables of type `vartype`.
     index_type add_variables(Vartype vartype, index_type n);
 
+    /// Add `n` variables of type `vartype` with lower bound `lb` and upper bound `ub`.
     index_type add_variables(Vartype vartype, index_type n, bias_type lb, bias_type ub);
 
+    /// Change the variable type of variable `v` to `vartype`, updating the biases appropriately.
     void change_vartype(Vartype vartype, index_type v);
 
     void clear();
@@ -84,6 +92,7 @@ class ConstrainedQuadraticModel {
     std::weak_ptr<Constraint<bias_type, index_type>> constraint_weak_ptr(index_type c);
     std::weak_ptr<const Constraint<bias_type, index_type>> constraint_weak_ptr(index_type c) const;
 
+    /// Fix variable `v` in the model to value `assignment`.
     template <class T>
     void fix_variable(index_type v, T assignment);
 
@@ -102,17 +111,22 @@ class ConstrainedQuadraticModel {
     /// Remove a constraint from the model.
     void remove_constraint(index_type c);
 
+    /// Remove variable `v` from the model.
     void remove_variable(index_type v);
 
+    /// Set a lower bound of `lb` on variable `v`.
     void set_lower_bound(index_type v, bias_type lb);
 
+    /// Set an objective for the model.
     template <class B, class I>
     void set_objective(const abc::QuadraticModelBase<B, I>& objective);
 
+    /// Set an objective for the model using a mapping.
     template <class B, class I, class T>
     void set_objective(const abc::QuadraticModelBase<B, I>& objective,
                        const std::vector<T>& mapping);
 
+    /// Set an upper bound of `ub` on variable `v`.
     void set_upper_bound(index_type v, bias_type ub);
     void set_vartype(index_type v, Vartype vartype);
 
diff --git a/dimod/include/dimod/constraint.h b/dimod/include/dimod/constraint.h
index d32c64e9f..319b6e851 100644
--- a/dimod/include/dimod/constraint.h
+++ b/dimod/include/dimod/constraint.h
@@ -29,16 +29,17 @@ enum Sense { LE, GE, EQ };
 
 enum Penalty { LINEAR, QUADRATIC, CONSTANT };
 
+/// A constrained quadratic model (CQM) can include hard or soft constraints as polynomials with one or two variables per term.
 template <class Bias, class Index>
 class Constraint : public Expression<Bias, Index> {
  public:
-    /// The type of the base class.
+    /// Type of the base class.
     using base_type = Expression<Bias, Index>;
 
-    /// The first template parameter (Bias).
+    /// First template parameter (`Bias`).
     using bias_type = Bias;
 
-    /// The second template parameter (Index).
+    /// Second template parameter (`Index`).
     using index_type = Index;
 
     using size_type = typename base_type::size_type;
@@ -48,21 +49,44 @@ class Constraint : public Expression<Bias, Index> {
     Constraint();
     explicit Constraint(const parent_type* parent);
 
+    /// Return true for a one-hot constraint of discrete variables.
     bool is_onehot() const;
+
+    /// Return true for a soft constraint with a finite weight that can be violated.
     bool is_soft() const;
+
+    /// Mark the constraint as encoding a discrete variable.
     void mark_discrete(bool mark = true);
+
+    /// Return true if the constraint encodes a discrete variable.
     bool marked_discrete() const;
+
+    /// Return the penalty set for a soft constraint.
     Penalty penalty() const;
+
+    /// Return a constraint's right-hand side.
     bias_type rhs() const;
 
     // note: flips sign when negative
+    /// Scale by multiplying by `scalar`.
     void scale(bias_type scalar);
 
+    /// Sense (greater or equal, less or equal, equal) of a constraint.
     Sense sense() const;
+
+    /// Set the penalty for a soft constraint.
     void set_penalty(Penalty penalty);
+
+    /// Set a constraint's right-hand side.
     void set_rhs(bias_type rhs);
+
+    /// Set the sense (greater or equal, less or equal, equal) of a constraint.
     void set_sense(Sense sense);
+
+    /// Set the weight for a soft constraint.
     void set_weight(bias_type weight);
+
+    /// Return a soft constraint's weight.
     bias_type weight() const;
 
  private:
diff --git a/dimod/include/dimod/expression.h b/dimod/include/dimod/expression.h
index bdc021a5b..a4b71ec3e 100644
--- a/dimod/include/dimod/expression.h
+++ b/dimod/include/dimod/expression.h
@@ -40,7 +40,7 @@ class Expression : public abc::QuadraticModelBase<Bias, Index> {
     /// The second template parameter (Index).
     using index_type = Index;
 
-    /// Unsigned integral type that can represent non-negative values.
+    /// Unsigned integer type that can represent non-negative values.
     using size_type = std::size_t;
 
     using parent_type = ConstrainedQuadraticModel<bias_type, index_type>;
diff --git a/dimod/include/dimod/quadratic_model.h b/dimod/include/dimod/quadratic_model.h
index b53c73f8b..df80b14a5 100644
--- a/dimod/include/dimod/quadratic_model.h
+++ b/dimod/include/dimod/quadratic_model.h
@@ -24,19 +24,20 @@
 
 namespace dimod {
 
+/// A quadratic model (QM) is a polynomial with one or two variables per term.
 template <class Bias, class Index = int>
 class QuadraticModel : public abc::QuadraticModelBase<Bias, Index> {
  public:
     /// The type of the base class.
     using base_type = abc::QuadraticModelBase<Bias, Index>;
 
-    /// The first template parameter (Bias).
+    /// The first template parameter (`Bias`).
     using bias_type = Bias;
 
-    /// The second template parameter (Index).
+    /// The second template parameter (`Index`).
     using index_type = Index;
 
-    /// Unsigned integral that can represent non-negative values.
+    /// Unsigned integer that can represent non-negative values.
     using size_type = typename base_type::size_type;
 
     QuadraticModel();
@@ -46,12 +47,16 @@ class QuadraticModel : public abc::QuadraticModelBase<Bias, Index> {
     template <class B, class I>
     explicit QuadraticModel(const BinaryQuadraticModel<B, I>& bqm);
 
+    /// Add variable of type `vartype`.
     index_type add_variable(Vartype vartype);
 
+    /// Add `n` variables of type `vartype` with lower bound `lb` and upper bound `ub`.
     index_type add_variable(Vartype vartype, bias_type lb, bias_type ub);
 
+    /// Add `n` variables of type `vartype`.
     index_type add_variables(Vartype vartype, index_type n);
 
+    /// Add `n` variables of type `vartype` with lower bound `lb` and upper bound `ub`.
     index_type add_variables(Vartype vartype, index_type n, bias_type lb, bias_type ub);
 
     void clear();
@@ -63,13 +68,14 @@ class QuadraticModel : public abc::QuadraticModelBase<Bias, Index> {
     bias_type lower_bound(index_type v) const;
 
     /**
-     * Total bytes consumed by the biases, vartype info, bounds, and indices.
+     * Total bytes consumed by the biases, `vartype` info, bounds, and indices.
      *
      * If `capacity` is true, use the capacity of the underlying vectors rather
      * than the size.
      */
     size_type nbytes(bool capacity = false) const;
 
+    /// Remove variable `v`.
     void remove_variable(index_type v);
 
     // Resize the model to contain `n` variables.
@@ -78,23 +84,25 @@ class QuadraticModel : public abc::QuadraticModelBase<Bias, Index> {
     /**
      * Resize the model to contain `n` variables.
      *
-     * The `vartype` is used to any new variables added.
-     *
-     * The `vartype` must be `Vartype::BINARY` or `Vartype::SPIN`.
+     * Any added variables are of type `vartype` (value must be `Vartype::BINARY` or `Vartype::SPIN`)
      */
     void resize(index_type n, Vartype vartype);
 
     /**
      * Resize the model to contain `n` variables.
      *
-     * The `vartype` is used to any new variables added.
+     * Any added variables are of type `vartype` (value must be `Vartype::BINARY` or `Vartype::SPIN`)
+     * and have lower bound `lb` and upper bound `ub`.
      */
     void resize(index_type n, Vartype vartype, bias_type lb, bias_type ub);
 
+    /// Set a lower bound of `lb` on variable `v`.
     void set_lower_bound(index_type v, bias_type lb);
 
+    /// Set an upper bound of `ub` on variable `v`.
     void set_upper_bound(index_type v, bias_type ub);
 
+    /// Set the variable type of variable `v`.
     void set_vartype(index_type v, Vartype vartype);
 
     // todo: substitute_variable with vartype/bounds support
diff --git a/dimod/include/dimod/vartypes.h b/dimod/include/dimod/vartypes.h
index 792767ba9..95660a3d5 100644
--- a/dimod/include/dimod/vartypes.h
+++ b/dimod/include/dimod/vartypes.h
@@ -73,6 +73,7 @@ class vartype_limits<Bias, Vartype::REAL> {
 template <class Bias>
 class vartype_info {
  public:
+    /// Default upper bound
     static Bias default_max(Vartype vartype) {
         if (vartype == Vartype::BINARY) {
             return vartype_limits<Bias, Vartype::BINARY>::default_max();
@@ -86,6 +87,7 @@ class vartype_info {
             throw std::logic_error("unknown vartype");
         }
     }
+    /// Default lower bound
     static Bias default_min(Vartype vartype) {
         if (vartype == Vartype::BINARY) {
             return vartype_limits<Bias, Vartype::BINARY>::default_min();
@@ -99,6 +101,7 @@ class vartype_info {
             throw std::logic_error("unknown vartype");
         }
     }
+    /// Maximum supported value
     static Bias max(Vartype vartype) {
         if (vartype == Vartype::BINARY) {
             return vartype_limits<Bias, Vartype::BINARY>::max();
@@ -112,6 +115,7 @@ class vartype_info {
             throw std::logic_error("unknown vartype");
         }
     }
+    /// Minimum supported value
     static Bias min(Vartype vartype) {
         if (vartype == Vartype::BINARY) {
             return vartype_limits<Bias, Vartype::BINARY>::min();
diff --git a/dimod/quadratic/quadratic_model.py b/dimod/quadratic/quadratic_model.py
index 99a9d73e0..7e15ea192 100644
--- a/dimod/quadratic/quadratic_model.py
+++ b/dimod/quadratic/quadratic_model.py
@@ -12,6 +12,21 @@
 #    See the License for the specific language governing permissions and
 #    limitations under the License.
 
+r"""Quadratic models are problems of the form:
+
+.. math::
+
+    E(x) = \sum_i a_i x_i + \sum_{i \le j} b_{i, j} x_i x_j + c
+
+where :math:`\{ x_i\}_{i=1, \dots, N}` can be binary\ [#]_ or integer
+variables and :math:`a_{i}, b_{ij}, c` are real values.
+
+.. [#]
+    For binary variables, the range of the quadratic-term summation is
+    :math:`i < j` because :math:`x^2 = x` for binary values :math:`\{0, 1\}`
+    and :math:`s^2 = 1` for spin values :math:`\{-1, 1\}`.
+"""
+
 from __future__ import annotations
 
 import struct
@@ -121,23 +136,8 @@ def loads_data(self, data):
 
 
 class QuadraticModel(QuadraticViewsMixin):
-    r"""A quadratic model.
-
-    Quadratic models are problems of the form:
+    r"A quadratic model."
 
-    .. math::
-
-        E(x) = \sum_i a_i x_i + \sum_{i \le j} b_{i, j} x_i x_j + c
-
-    where :math:`\{ x_i\}_{i=1, \dots, N}` can be binary\ [#]_ or integer
-    variables and :math:`a_{i}, b_{ij}, c` are real values.
-
-    .. [#]
-        For binary variables, the range of the quadratic-term summation is
-        :math:`i < j` because :math:`x^2 = x` for binary values :math:`\{0, 1\}`
-        and :math:`s^2 = 1` for spin values :math:`\{-1, 1\}`.
-
-    """
     _DATA_CLASSES = {
         np.dtype(np.float32): cyQM_float32,
         np.dtype(np.float64): cyQM_float64,
diff --git a/docs/intro/intro_samplers.rst b/docs/intro/intro_samplers.rst
index 897ca0dbd..e762490e0 100644
--- a/docs/intro/intro_samplers.rst
+++ b/docs/intro/intro_samplers.rst
@@ -21,7 +21,7 @@ Example: Using a Reference Sampler
 
 To find solutions to the small four-node
 `maximum cut <https://en.wikipedia.org/wiki/Maximum_cut>`_
-BQM generated in the :ref:`intro_qm` section, shown again in the figure below,
+BQM generated in the :ref:`intro_models` section, shown again in the figure below,
 you can use one of dimod's reference samplers: its
 :class:`~dimod.reference.samplers.ExactSolver` test sampler, for example,
 calculates the energy of all possible samples.
diff --git a/docs/reference/cpp.rst b/docs/reference/cpp.rst
index c547b38ff..c73de9f77 100644
--- a/docs/reference/cpp.rst
+++ b/docs/reference/cpp.rst
@@ -4,23 +4,48 @@
 C++ API
 =======
 
-dimod
-=====
+This page describes the `dimod` package's C++ API.
 
-BinaryQuadraticModel
---------------------
+.. contents::
+    :local:
+    :depth: 3
+
+Models
+======
+
+For the Python API and descriptions of the various models, see :ref:`dimod_models`.
+
+Binary Quadratic Model (BQM)
+----------------------------
 
 .. doxygenclass:: dimod::BinaryQuadraticModel
     :members:
     :project: dimod
 
-QuadraticModel
---------------
+Constrained Quadratic Model (CQM)
+---------------------------------
+
+.. doxygenclass:: dimod::ConstrainedQuadraticModel
+    :members:
+    :project: dimod
+
+Constraints
+~~~~~~~~~~~
+
+.. doxygenclass:: dimod::Constraint
+    :members:
+    :project: dimod
+
+Quadratic Model (QM)
+--------------------
 
 .. doxygenclass:: dimod::QuadraticModel
     :members:
     :project: dimod
 
+Variable Type (Vartype)
+=======================
+
 Vartype
 -------
 
@@ -35,10 +60,19 @@ vartype_info
     :undoc-members:
     :project: dimod
 
+vartype_limits
+--------------
+
 .. Todo: vartype_limits. Getting it to look nice is possible but fiddly
 
-dimod::abc
-==========
+.. doxygenclass:: dimod::vartype_limits
+    :members:
+    :undoc-members:
+    :project: dimod
+
+
+dimod Abstract Base Class (`dimod::abc`)
+========================================
 
 .. doxygenclass:: dimod::abc::QuadraticModelBase
     :members:
@@ -47,8 +81,8 @@ dimod::abc
 
 .. Todo: dimod lp
 
-dimod::utils
-============
+dimod Utilities (`dimod::utils`)
+================================
 
 .. doxygenfunction:: zip_sort
    :project: dimod
diff --git a/docs/reference/models.rst b/docs/reference/models.rst
index 2b250901d..3297d9b3b 100644
--- a/docs/reference/models.rst
+++ b/docs/reference/models.rst
@@ -5,7 +5,7 @@ Models: BQM, CQM, QM, Others
 ============================
 
 This page describes the `dimod` package's quadratic models: classes,
-attributes, and methods. For an introduction and the data structure, see :ref:`intro_qm`.
+attributes, and methods. For an introduction and the data structure, see :ref:`intro_models`.
 
 .. contents::
     :local:
@@ -16,8 +16,11 @@ attributes, and methods. For an introduction and the data structure, see :ref:`i
 Binary Quadratic Models
 =======================
 
-For an introduction to BQMs, see
-:std:doc:`Concepts: Binary Quadratic Models <oceandocs:concepts/bqm>`.
+.. automodule:: dimod.binary.binary_quadratic_model
+
+For examples, see
+`Ocean's Getting Started examples <https://docs.ocean.dwavesys.com/en/stable/getting_started.html#examples>`_
+and the BQM examples in `D-Wave's collection of examples <https://github.com/dwave-examples>`_.
 
 BQM Class
 ---------
@@ -126,8 +129,11 @@ Generic BQM Constructor
 Constrained Quadratic Model
 ===========================
 
-For an introduction to BQMs, see
-:std:doc:`Concepts: Constrained Quadratic Models <oceandocs:concepts/cqm>`.
+.. automodule:: dimod.constrained.constrained
+
+For examples, see
+`Ocean's Getting Started examples <https://docs.ocean.dwavesys.com/en/stable/getting_started.html#examples>`_
+and the CQM examples in `D-Wave's collection of examples <https://github.com/dwave-examples>`_.
 
 CQM Class
 ---------
@@ -207,8 +213,12 @@ Sense Class
 Quadratic Models
 ================
 
-For an introduction to QMs, see
-:std:doc:`Concepts: Quadratic Models <oceandocs:concepts/qm>`.
+.. automodule:: dimod.quadratic.quadratic_model
+
+For examples, see
+`Ocean's Getting Started examples <https://docs.ocean.dwavesys.com/en/stable/getting_started.html#examples>`_
+and the examples in `D-Wave's collection of examples <https://github.com/dwave-examples>`_.
+
 
 QM Class
 --------