source: coopr.pyomo/trunk/coopr/pyomo/io/ampl/cAmpl/handlers/expression/prodexp.c @ 4439

Revision 4439, 15.3 KB checked in by tekl, 3 years ago (diff)

Fix bug in C NL writer for handling certain products

Line 
1#include "prodexp.h"
2#include "../../cAmpl.h"
3
4int _handle_prodexp(PyObject * context, PyObject * exp, PyObject ** ampl_repn) {
5    // tmp
6    PyObject * _result = NULL;
7    PyObject * _result2 = NULL;
8    Py_ssize_t i;
9
10    // denom=1.0
11    PyObject * denom = Py_BuildValue("f", 1.0);
12   
13    // for e in exp._denominator:
14    PyObject * exp__denominator = PyObject_GetAttrString(exp, "_denominator");
15    Py_ssize_t len_exp__denominator = PySequence_Length(exp__denominator);
16    for(i = 0; i < len_exp__denominator; i++) {
17        PyObject * e = PySequence_GetItem(exp__denominator, i);
18
19        // if e.fixed_value():
20        PyObject * e_fixed_value = PyObject_CallMethod(e, "fixed_value", NULL);
21        if(PyObject_IsTrue(e_fixed_value) == TRUE) {
22            // denom *= e.value
23            PyObject * e_value = PyObject_GetAttrString(e, "value");
24            _result = PyNumber_InPlaceMultiply(denom, e_value);
25            Py_DECREF(denom);
26            denom = _result;
27            _result = NULL;
28            Py_DECREF(e_value);
29
30        // elif e.is_constant():
31        } else {
32            PyObject * e_is_constant = PyObject_CallMethod(e, "is_constant", "()");
33            if(PyObject_IsTrue(e_is_constant) == TRUE) {
34                // denom *= e()
35                PyObject * e_ = PyObject_CallObject(e, NULL);
36                _result = PyNumber_Multiply(denom, e_);
37                Py_DECREF(denom);
38                denom = _result;
39                _result = NULL;
40                Py_DECREF(e_);
41
42            // else:
43            } else {
44                // ampl_repn._nonlinear_expr = exp
45                PyObject_SetAttrString(*ampl_repn, "_nonlinear_expr", exp);
46
47                // break
48                Py_DECREF(e_is_constant);
49                Py_DECREF(e_fixed_value);
50                Py_DECREF(e);
51                break;
52            }
53            Py_DECREF(e_is_constant);
54        }
55        Py_DECREF(e_fixed_value);
56
57        // if denom == 0.0:
58        PyObject * zero = Py_BuildValue("f", 0.0);
59        if(PyObject_RichCompareBool(denom, zero, Py_EQ) == TRUE) {
60            // print "Divide-by-zero error - offending sub-expression:"
61            PySys_WriteStdout("Divide-by-zero error - offending sub-expression:\n");
62
63            // e.pprint()
64            PyObject_CallMethod(e, "pprint", "()");
65
66            // raise ZeroDivisionError
67            PyErr_SetString(PyExc_ZeroDivisionError, "");
68            Py_DECREF(zero);
69            Py_DECREF(e);
70            return ERROR;
71        }
72
73        Py_DECREF(zero);
74        Py_DECREF(e);
75    }
76
77    // if not ampl_repn._nonlinear_expr is None:
78    PyObject * ampl_repn__nle = PyObject_GetAttrString(*ampl_repn, "_nonlinear_expr");
79    if(ampl_repn__nle != Py_None) {
80        // opt
81        PyObject * _update = PyString_FromString("update");
82
83        // we have a nonlinear expression ... build up all the vars
84        // for e in exp._denominator:
85        PyObject * exp__denominator = PyObject_GetAttrString(exp, "_denominator");
86        Py_ssize_t len_exp__denominator = PySequence_Length(exp__denominator);
87        for(i = 0; i < len_exp__denominator; i++) {
88            PyObject * e = PySequence_GetItem(exp__denominator, i);
89
90            // arg_repn = generate_ampl_repn(e)
91            PyObject * arg_repn = recursive_generate_ampl_repn(context, e);
92            Py_DECREF(e);
93            if(arg_repn == NULL) return ERROR;
94
95            // opt
96            PyObject * ampl_repn__nlv = PyObject_GetAttrString(*ampl_repn, "_nonlinear_vars");
97
98            // ampl_repn._nonlinear_vars.update(arg_repn._linear_terms_var)
99            PyObject * arg_repn__ltv = PyObject_GetAttrString(arg_repn, "_linear_terms_var");
100            PyObject_CallMethodObjArgs(ampl_repn__nlv, _update, arg_repn__ltv, NULL);
101            Py_DECREF(arg_repn__ltv);
102
103            // ampl_repn._nonlinear_vars.update(arg_repn._nonlinear_vars)
104            PyObject * arg_repn__nlv = PyObject_GetAttrString(arg_repn, "_nonlinear_vars");
105            PyObject_CallMethodObjArgs(ampl_repn__nlv, _update, arg_repn__nlv, NULL);
106            Py_DECREF(arg_repn__nlv);
107
108            Py_DECREF(ampl_repn__nlv);
109            Py_DECREF(arg_repn);
110        }
111        Py_DECREF(exp__denominator);
112
113        // for e in exp._numerator;
114        PyObject * exp__numerator = PyObject_GetAttrString(exp, "_numerator");
115        Py_ssize_t len_exp__numerator = PySequence_Length(exp__numerator);
116        Py_ssize_t i;
117        for(i = 0; i < len_exp__numerator; i++) {
118            PyObject * e = PySequence_GetItem(exp__numerator, i);
119
120            // arg_repn = generate_ampl_repn(e)
121            PyObject * arg_repn = recursive_generate_ampl_repn(context, e);
122            if(arg_repn == NULL) return ERROR;
123            Py_DECREF(e);
124
125            // opt
126            PyObject * ampl_repn__nlv = PyObject_GetAttrString(*ampl_repn, "_nonlinear_vars");
127
128            // ampl_repn._nonlinear_vars.update(arg_repn._linear_terms_var)
129            PyObject * arg_repn__ltv = PyObject_GetAttrString(arg_repn, "_linear_terms_var");
130            PyObject_CallMethodObjArgs(ampl_repn__nlv, _update, arg_repn__ltv, NULL);
131            Py_DECREF(arg_repn__ltv);
132
133            // ampl_repn._nonlinear_vars.update(arg_repn._nonlinear_vars)
134            PyObject * arg_repn__nlv = PyObject_GetAttrString(arg_repn, "_nonlinear_vars");
135            PyObject_CallMethodObjArgs(ampl_repn__nlv, _update, arg_repn__nlv, NULL);
136            Py_DECREF(arg_repn__nlv);
137
138            Py_DECREF(ampl_repn__nlv);
139            Py_DECREF(arg_repn);
140        }
141        Py_DECREF(exp__numerator);
142
143        Py_DECREF(ampl_repn__nle);
144        Py_DECREF(_update);
145        return TRUE;
146    }
147    Py_DECREF(ampl_repn__nle);
148
149    // OK, the denominator is a constant
150    // build up the ampl_repns for the numerator
151    // C NOTE: we break from Python objects a bit here for the counters
152    // n_linear_args = 0
153    int n_linear_args = 0;
154
155    // n_nonlinear_args = 0
156    int n_nonlinear_args = 0;
157
158    // arg_repns = list();
159    PyObject * arg_repns = Py_BuildValue("[]");
160
161    // for i in xrange(len(exp._numerator)):
162    // C NOTE: don't actually care about the xrange
163    PyObject * exp__numerator = PyObject_GetAttrString(exp, "_numerator");
164    Py_ssize_t len_exp__numerator = PySequence_Length(exp__numerator);
165    for(i = 0; i < len_exp__numerator; i++) {
166        // e = exp._numerator[i]
167        PyObject * e = PySequence_GetItem(exp__numerator, i);
168
169        // e_repn = generate_ampl_repn(e)
170        PyObject * e_repn = recursive_generate_ampl_repn(context, e);
171        Py_DECREF(e);
172        if(e_repn == NULL) return ERROR;
173
174        // arg_repns.append(e_repn)
175        PyList_Append(arg_repns, e_repn);
176
177        // check if the expression is not nonlinear else it is nonlinear
178        // if not e_repn._nonlinear_expr is None:
179        PyObject * e_repn__nle = PyObject_GetAttrString(e_repn, "_nonlinear_expr");
180        if(e_repn__nle != Py_None) {
181            // n_nonlinear_args += 1
182            n_nonlinear_args += 1;
183
184        // check whether the expression is constant or else it is linear
185        // elif not ((len(e_repn._linear_terms_var) == 0) and (e_repn._nonlinear_expr is None)):
186        } else {
187            PyObject * e_repn__ltv = PyObject_GetAttrString(e_repn, "_linear_terms_var");
188            Py_ssize_t len_e_repn__ltv = PyDict_Size(e_repn__ltv);
189            if(!(len_e_repn__ltv == 0 && e_repn__nle == Py_None)) {
190                // n_linear_args += 1
191                n_linear_args += 1;
192            }
193            Py_DECREF(e_repn__ltv);
194        }
195        Py_DECREF(e_repn__nle);
196    }
197    Py_DECREF(exp__numerator);
198
199    // is_nonlinear = False;
200    // C NOTE: more trickery to avoid extraneous PyObjects
201    int is_nonlinear = FALSE;
202
203    // if n_linear_args > 1 or n_nonlinear_args > 0:
204    if(n_linear_args > 1 || n_nonlinear_args > 0) {
205        is_nonlinear = TRUE;
206    }
207
208    // if is_nonlinear is True:
209    if(is_nonlinear == TRUE) {
210        // do like AMPL and simply return the expression
211        // without extracting the potentially linear part
212        // ampl_repn = ampl_representation()
213        *ampl_repn = new_ampl_representation();
214
215        // ampl_repn._nonlinear_expr = exp
216        PyObject_SetAttrString(*ampl_repn, "_nonlinear_expr", exp);
217
218        // for repn in arg_repns:
219        Py_ssize_t len_arg_repns = PySequence_Length(arg_repns);
220        for(i = 0; i < len_arg_repns; i++) {
221            PyObject * repn = PySequence_GetItem(arg_repns, i);
222
223            // opt
224            PyObject * ampl_repn__nlv = PyObject_GetAttrString(*ampl_repn, "_nonlinear_vars");
225
226            // ampl_repn._nonlinear_vars.update(repn._linear_terms_var)
227            PyObject * repn__ltv = PyObject_GetAttrString(repn, "_linear_terms_var");
228            PyObject_CallMethod(ampl_repn__nlv, "update", "(O)", repn__ltv);
229            Py_DECREF(repn__ltv);
230
231            // ampl_repn._nonlinear_vars.update(repn._nonlinear_vars)
232            PyObject * repn__nlv = PyObject_GetAttrString(repn, "_nonlinear_vars");
233            PyObject_CallMethod(ampl_repn__nlv, "update", "(O)", repn__nlv);
234            Py_DECREF(repn__nlv);
235        }
236
237        // return ampl_repn
238        Py_DECREF(arg_repns);
239        return TRUE;
240
241    // else:
242    } else {
243        // is linear or constant
244        // ampl_repn = current_repn = arg_repns[0]
245        PyObject * current_repn = PySequence_GetItem(arg_repns, 0);
246        *ampl_repn = current_repn;
247
248        // for i in xrange(1, len(arg_repns)):
249        // C NOTE: don't care about the xrange
250        Py_ssize_t len_arg_repns = PySequence_Length(arg_repns);
251        for(i = 1; i < len_arg_repns; i++) {
252            // e_repn = arg_repns[i]
253            PyObject * e_repn = PySequence_GetItem(arg_repns, i);
254
255            // TODO verify this works in different import scenarios
256            // ampl_repn = ampl_representation()
257            *ampl_repn = new_ampl_representation();
258
259            // const_c * const_e
260            // ampl_repn._constant = current_repn._constant * e_repn._constant
261            PyObject * current_repn__constant = PyObject_GetAttrString(current_repn, "_constant");
262            PyObject * e_repn__constant = PyObject_GetAttrString(e_repn, "_constant");
263            _result = PyNumber_Multiply(current_repn__constant, e_repn__constant);
264            PyObject_SetAttrString(*ampl_repn, "_constant", _result);
265            Py_DECREF(_result);
266
267            // const_e * L_c
268            // if e_repn._constant != 0.0:
269            PyObject * zero = Py_BuildValue("f", 0.0);
270            if(PyObject_RichCompareBool(e_repn__constant, zero, Py_NE) == TRUE) {
271                // opt
272                PyObject * current_repn__ltv = PyObject_GetAttrString(current_repn, "_linear_terms_var");
273                PyObject * current_repn__ltc = PyObject_GetAttrString(current_repn, "_linear_terms_coef");
274                PyObject * ampl_repn__ltv = PyObject_GetAttrString(*ampl_repn, "_linear_terms_var");
275                PyObject * ampl_repn__ltc = PyObject_GetAttrString(*ampl_repn, "_linear_terms_coef");
276
277                // for(var_name, var) in current_repn._linear_terms_var.iteritems():
278                PyObject * var_name, * var;
279                i = 0;
280                while(PyDict_Next(current_repn__ltv, &i, &var_name, &var)) {
281                    // ampl_repn._linear_terms_coef[var_name] = current_repn._linear_terms_coef[var_name] * e_repn._constant
282                    PyObject * cr__ltc_vn = PyDict_GetItem(current_repn__ltc, var_name);
283                    _result = PyNumber_Multiply(cr__ltc_vn, e_repn__constant);
284                    Py_DECREF(cr__ltc_vn);
285                    PyDict_SetItem(ampl_repn__ltc, var_name, _result);
286                    Py_DECREF(_result); _result = NULL;
287
288                    // ampl_repn._linear_terms_var[var_name] = var
289                    PyDict_SetItem(ampl_repn__ltv, var_name, var);
290                }
291
292                Py_DECREF(current_repn__ltv);
293                Py_DECREF(current_repn__ltc);
294                Py_DECREF(ampl_repn__ltv);
295                Py_DECREF(ampl_repn__ltc);
296            }
297
298            // const_c * L_e
299            // if current_repn._constant != 0.0:
300            if(PyObject_RichCompareBool(current_repn__constant, zero, Py_NE) == TRUE) {
301                // opt
302                PyObject * e_repn__ltv = PyObject_GetAttrString(e_repn, "_linear_terms_var");
303                PyObject * e_repn__ltc = PyObject_GetAttrString(e_repn, "_linear_terms_coef");
304                PyObject * ampl_repn__ltv = PyObject_GetAttrString(*ampl_repn, "_linear_terms_var");
305                PyObject * ampl_repn__ltc = PyObject_GetAttrString(*ampl_repn, "_linear_terms_coef");
306
307                // for (e_var_name,e_var) in e_repn._linear_terms_var.iteritems():
308                PyObject * e_var_name, * e_var;
309                i = 0;
310                while(PyDict_Next(e_repn__ltv, &i, &e_var_name, &e_var)) {
311                    // opt
312                    PyObject * current_repn__constant = PyObject_GetAttrString(current_repn, "_constant");
313                    PyObject * er__ltc_evn = PyDict_GetItem(e_repn__ltc, e_var_name);
314                    _result = PyNumber_Multiply(current_repn__constant, er__ltc_evn);
315                    Py_DECREF(current_repn__constant);
316                    Py_DECREF(er__ltc_evn);
317
318                    // if e_var_name in ampl_repn._linear_terms_var:
319                    if(PyMapping_HasKey(ampl_repn__ltv, e_var_name) == TRUE) {
320                        // ampl_repn._linear_terms_coef[e_var_name] += current_repn._constant * e_repn._linear_terms_coef[e_var_name]
321                        PyObject * _tmp = PyDict_GetItem(ampl_repn__ltc, e_var_name);
322                        _result2 = PyNumber_Add(_tmp, _result);
323                        PyDict_SetItem(ampl_repn__ltc, e_var_name, _result2);
324                        Py_DECREF(_tmp);
325                        Py_DECREF(_result2); _result2 = NULL;
326
327                    // else:
328                    } else {
329                        // ampl_repn._linear_terms_coef[e_var_name] = current_repn._constant * e_repn._linear_terms_coef[e_var_name]
330                        PyDict_SetItem(ampl_repn__ltc, e_var_name, _result);
331                    }
332                    Py_DECREF(_result); _result = NULL;
333                }
334            }
335
336            Py_DECREF(e_repn__constant);
337            Py_DECREF(current_repn__constant);
338            Py_DECREF(zero);
339        }
340
341        // current_repn = ampl_repn
342        current_repn = *ampl_repn;
343    }
344
345    // now deal with the product expression's coefficient that needs
346    // to be applied to all parts of the ampl_repn
347    // opt
348    PyObject * exp_coef = PyObject_GetAttrString(exp, "coef");
349    PyObject * quotient = PyNumber_Divide(exp_coef, denom);
350    Py_DECREF(exp_coef);
351
352    // ampl_repn._constant *= exp.coef/denom
353    PyObject * ampl_repn__constant = PyObject_GetAttrString(*ampl_repn, "_constant");
354    _result = PyNumber_Multiply(ampl_repn__constant, quotient);
355    Py_DECREF(ampl_repn__constant);
356    PyObject_SetAttrString(*ampl_repn, "_constant", _result);
357    Py_DECREF(_result); _result = NULL;
358
359    // for var_name in ampl_repn._linear_terms_coef:
360    PyObject * ampl_repn__ltc = PyObject_GetAttrString(*ampl_repn, "_linear_terms_coef");
361    PyObject * var_name, * _var;
362    i = 0;
363    // C NOTE: TODO this is scary. is this mutation allowed?
364    while(PyDict_Next(ampl_repn__ltc, &i, &var_name, &_var)) {
365        _result = PyNumber_Multiply(_var, quotient);
366        PyDict_SetItem(ampl_repn__ltc, var_name, _result);
367        Py_DECREF(_result); _result = NULL;
368    }
369
370    // return ampl_repn
371    return TRUE;
372}
Note: See TracBrowser for help on using the repository browser.