dimod.higherorder.utils.make_quadratic gives inconsistent / incorrect results

**Description**
[D-Wave has indicated at least two approaches](https://docs.dwavesys.com/docs/latest/c_handbook_3.html) for reducing higher order terms to quadratic, namely Reduction by Substitution and Reduction by Minimum Selection. Both of these are shown to be mathematically correct approaches.

The issue is with the [make_quadratic](https://docs.ocean.dwavesys.com/projects/dimod/en/latest/reference/generated/dimod.higherorder.utils.make_quadratic.html#dimod.higherorder.utils.make_quadratic) function. The function has at least one clear bug, which is that it returns completely different results if vartype = 'BINARY' or dimod.BINARY, which according to the [vartype documentation](https://docs.ocean.dwavesys.com/projects/dimod/en/latest/reference/vartypes.html#dimod.Vartype) those should both be identical.

For the trivial example case higher order polynomial:
`
poly = {(0, 1, 2): 10}
`

If vartype = dimod.BINARY, we get the following output:
`
BinaryQuadraticModel({0: 0.0, 1: 0.0, '0*1': 3.0, 2: 0}, {(0, 1): 1.0, (0, '0*1'): -2.0, (1, '0*1'): -2.0, ('0*1', 2): 10}, 0.0, Vartype.BINARY)
`
This is exactly what you would expect if we were using the Reduction by Substitution approach. However, if vartype='BINARY', we get the following output:
`
BinaryQuadraticModel({0: -5.0, 1: -5.0, '0*1': -5.0, 'aux0,1': -8.0, 2: 0}, {(0, 1): 2.0, (0, '0*1'): 2.0, (0, 'aux0,1'): 4.0, (1, '0*1'): 2.0, (1, 'aux0,1'): 4.0, ('0*1', 'aux0,1'): 4.0, ('0*1', 2): 10}, 9.0, Vartype.BINARY)
`
**Steps To Reproduce**
```
import dimod

poly = {(0, 1, 2): 10}
bqm = dimod.make_quadratic(poly, 1.0, 'SPIN')
print(bqm)
print('\n')
bqm = dimod.make_quadratic(poly, 1.0, 'BINARY')
print(bqm)
print('\n')
bqm = dimod.make_quadratic(poly, 1.0, dimod.SPIN)
print(bqm)
print('\n')
bqm = dimod.make_quadratic(poly, 1.0, dimod.BINARY)
print(bqm)
print('\n')

assert dimod.make_quadratic(poly, 1.0, 'SPIN') == dimod.make_quadratic(poly, 1.0, dimod.SPIN)
assert dimod.make_quadratic(poly, 1.0, 'BINARY') == dimod.make_quadratic(poly, 1.0, dimod.BINARY)
```

**Expected Behavior**
I would expect 3 things that are currently missing in the code:

1. Consistency: the two assert statements in the code above should always be true.
2. Clarity: documentation should make clear either the option to choose the reduction method, or clarify which reduction method is being used explicitly.
3. Correctness: a unit test is required that validates the exact expected output for several different input QUBOs / Ising Hamiltonians.

A final note: it seems inconsistent that somehow make_quadratic requires only 1 ancilla qubit for BINARY and 2 ancilla qubits for SPIN, as converting between the two expressions require no ancilla qubits whatsoever.

**Environment**
- OS: Mac High Sierra 10.13.6
- Python version: Python 3.7.4

**Additional Context**
No additional context required.