This is much faster than fitting the features, and can be used
by trapping to replace some constraint index calculations.

Also:
Add docstring/type annotations
Rename a variable so it doesn't shadow newly-imported name

There are several major differences from the version in example 8:
* Coupled with PolynomialLibrary logic using PolynomialLibrary.powers_.
Previously, this coupling was implicit and actually tied to the specific order
in Example 8, not the order users would expect.
* Instead of handling pre-flattened 2D array, constraints are built with
target and feature axes separately, letting the caller handle formatting
the constraints as "target" or "feature".  This gives us flexibility to
later clean up the constraint logic outside this function without needing to
return here and removes the need for striding calculations like:
n_tgts * (n_tgts + k - 1) + i + n_tgts * int(j * (2 * n_tgts - j - 3) / 2.0)
Since constraints are built directly from the terms and exponents in the
library, its much easier to read and debug the code.
* Renamed "r" as "n_tgts", etc.  Single-letter index variables are
removed.  When they occur in homogenous iterators, they are concatenated with
what they're iterating (e.g. tgt_i, tgt_j)
* Broke out the different kinds of constraints into helper functions to make
it easier to test.

This allows follow on work to auto-attach constraints to a problem.

(other than degree, which has to be 2)