Unexpected behaviour for dimod.make_quadratic

Description dimod.make_quadratic has a keyword argument strength. According to the docs, this variable should take positive values. However, positive values seem to return incorrect results, while negative values do work.

To Reproduce

import numpy as np
import matplotlib.pyplot as pt
import dimod

def bqm_to_dict(bqm, vars):
    ''' convert BQM to dictionary '''
    d = {}
    for i,j,v in bqm.iter_quadratic():
        d[(vars.index(i),vars.index(j))] = v
    for i,v in bqm.iter_linear():
        d[(vars.index(i),)] = v
    d[()] = bqm.offset
    return d

def eval_energy(merged_dict, X):
    ''' evaluate energy of PUSO/BUSO '''
    energy = 0.
    for idx, val in merged_dict.items():
        energy -= val * np.prod([X[int(i)] for i in idx])
    return energy

# generate random PUSO problem with 3 variables
coeffs = np.random.rand(7)*2-1
dict_puso = {
 ('0',): coeffs[0],
 ('1',): coeffs[1],
 ('2',): coeffs[2],
 ('0', '1'): coeffs[3],
 ('0', '2'): coeffs[4],
 ('1', '2'): coeffs[5],
 ('0', '1', '2'): coeffs[6]}

# quadratization
bqm = dimod.make_quadratic(poly=dict_puso, strength=10, vartype=dimod.SPIN)

# reorder variables
vars_head = ['0','1','2']
vars_tail = [v for v in bqm.variables if not v in vars_head]
vars = vars_head + vars_tail

# get dictionary representation of BQM
dict_buso = bqm_to_dict(bqm, vars=vars)

# visualize energy landscape over possible spin configurations
xvals, yvals_puso, yvals_buso = [], [], []
for x in [-1,1]:
    for y in [-1,1]:
        for z in [-1,1]:
            E_puso = eval_energy(dict_puso, np.array([x,y,z]))

            E_buso = np.inf
            for aux1 in [-1,1]:
                for aux2 in [-1,1]:
                    temp = eval_energy(dict_buso, np.array([x,y,z,aux1,aux2]))
                    if temp < E_buso:
                        E_buso = temp

            xvals.append(f'({x},{y},{z})')
            yvals_puso.append(E_puso)
            yvals_buso.append(E_buso)

pt.plot(xvals, yvals_puso)
pt.title('PUSO')
pt.show()

pt.plot(xvals, yvals_buso)
pt.title('BUSO')
pt.show()

Expected behavior The 2 plots generated by the code above should be similar (as they represent the energy landscape before and after quadratization). This isn't the case for positive values of strength, while it is for negative values

Environment:

dimod: 0.12.14

dwavesystems / dwave-ocean-sdk

Unexpected behaviour for dimod.make_quadratic #305