Problems computing degree of expressions with numpy data

[whart222]

The following script returns an error that indicates that the expression system things that the Objective constraint is a constant, when it clearly isn't. Digging around, it looks like there is some confusion generated by the NumPY data types.

import numpy as np
from pyomo.environ import *
opt = SolverFactory("glpk")
from pyomo.opt import SolverFactory

model = AbstractModel()

nsample = 500
nvariables = 20
X0 = np.ones([nsample,1])
model.X = np.random.uniform(0,10,([nsample,nvariables]))
X = np.concatenate([X0,model.X],axis = 1)

model.I = RangeSet(1,nsample) 
model.J = RangeSet(1,nvariables) 


error = np.random.normal(0,1,(nsample,1))
beta = np.random.randint(-5,5,size = ([nvariables+1,1]))
model.Y = np.dot(X,beta) + error

model.beta = Var(model.J)
model.beta0 = Var()

def obj_fun(model):
    
    return sum(abs(model.Y[i-1]-(model.beta0 + sum(model.X[i-1,j-1]*model.beta[j] for j in model.J) )) for i in model.I)

model.OBJ = Objective(rule = obj_fun, sense = minimize)



opt = SolverFactory('glpk')
instance = model.create_instance()
results = opt.solve(instance)
results.write()

This code generates the following error:

TypeError: Implicit conversion of Pyomo NumericValue type `<class 'pyomo.core.base.expr_coopr3._SumExpression'>' to a float is
disabled. This error is often the result of using Pyomo components as
arguments to one of the Python built-in math module functions when
defining expressions. Avoid this error by using Pyomo-provided math
functions.

And it generates the warning

WARNING: Constant objective detected, replacing with a placeholder to prevent solver failure.

In the cpxlp writer, the call

degree = canonical_degree(canonical_repn)

returns a value of zero for the degree!

[ghackebeil]

The problem is that the objective rule is returning a 1x1 array

 [<pyomo.core.base.expr_coopr3._SumExpression object at 0x10a895208>]

If you tweak the numpy data to use arrays of shape (n,) for vectors rather than (n,1), it fixes the problem.

error = np.random.normal(0,1, size=nsample)
beta = np.random.randint(-5,5, size=nvariables+1)
model.Y = np.dot(X,beta) + error
print(model.Y.shape)    # -> (500,)

In the original code, the shape of model.Y that is reported is (500,1).

[whart222]

How does this impact the expression generation? Why is Pyomo sensitive to this difference?

…

On Dec 18, 2016 4:16 PM, "Gabriel Hackebeil" ***@***.***> wrote: The problem is that the objective rule is returning a 1x1 array [<pyomo.core.base.expr_coopr3._SumExpression object at 0x10a895208>] If you tweak the numpy data to use arrays of shape (n,) for vectors rather than (n,1), it fixes the problem. error = np.random.normal(0,1, size=nsample) beta = np.random.randint(-5,5, size=nvariables+1) model.Y = np.dot(X,beta) + errorprint(model.Y.shape) # -> (500,) In the original code, the shape of model.Y that is reported is (500,1). — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#87 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAsb-MvHRM78dFVu-ZOcgYq_zf8KDaWoks5rJaKigaJpZM4LQOaR> .

[ghackebeil]

I don't think it impacted the expression generation, per se. The operator overloading still generates the expected expression, just inside an array. If you dot a zero-based indexed variable with a flat numpy array, you will get a sum of product expressions.

The fact that the original result here is a 1x1 array rather than a scalar just has to do with numpy returning an object of the correct shape based on the shapes of the arrays used in the expression.

The rest of the behavior simply has to do with us not expecting a list or array to be returned from a rule, so whatever happens after that point probably is not good.

[jsiirola]

@ghackebeil is right: the problem is not with expression generation, but rather with the error checking within the Objective.construct(). Actually, the problem is all the way down in as_numeric(). In there, we try to automatically detect and register "native" numeric types - that is, types that behave like POD numeric data - so that Pyomo can in theory work with other custom numeric types. The test we use is:

if obj.__class__ is (obj + 0).__class__:

The idea is that if adding 0 gives identity, then it must be a type that behaves like a numeric type. The problem is that the np.ndarry satisfies this test. Pyomo then adds the ndarray type to the list of known "native_numeric_types" ... and in turn, things like Constraint and Objective treat the ndarray as if it were constant (regardless of what it holds).

Now, what to do?

The one thing we can do is re-think native (POD) data detection occurs in as_numeric(). This seems like a pretty severe problem, as I never considered folks using generic containers that behaved like POD data, but could contain Pyomo expressions. I am not entirely sure how to detect the difference between something like a numpy.float64 and a numpy.ndarray. The former is something that we should consider to be a native_numeric_type, but the latter is NOT. Unfortunately, for both data types:

obj.__class__ is (obj + 0).__class__  # True
obj is not obj + 0  # True

For the record, this problem is related to both #31 and #68.