cvxgrp / qcml

A Python parser for generating Python/C/Matlab solver interfaces
Other
42 stars 9 forks source link

C codegen: use of undeclared identifier 'result' #44

Closed alexvoronov closed 10 years ago

alexvoronov commented 10 years ago

This is a continuation of a discussion in #42.

I have this issue when I try to compile C code generated for the following example:

#!/usr/bin/env python
from qcml import QCML

if __name__ == '__main__':

    s = ("""variable x(2)
        parameter D(2,2) positive
        #minimize sum(2*D*(-x) + (x))
        minimize sum(2*D*x)
        """)

    p = QCML(debug=False)
    p.parse(s)

    p.canonicalize()
    p.codegen("C")
    p.save("bug5a")

Generated code is as following:

#include <stdlib.h>
#include "bug5a.h"

/* ----------------------- BEGIN GENERATED CODE --------------------------- */
qc_socp * qc_bug5a2socp(const bug5a_params * params, const bug5a_dims * dims)
{
    /*
     * maps 'params' into the C socp data type
     * 'params' ought to contain:
     *   'D' has shape Matrix(2,2)
     */

    /* all local variables */
    long i;  /* loop index */
    long *q_ptr;
    long *A_row_ptr, *A_col_ptr;
    long *G_row_ptr, *G_col_ptr;
    double *A_data_ptr, *G_data_ptr;
    long nnzA, nnzG;
    qc_matrix *G_csc, *A_csc;  /* possibly un-used */
    qc_matrix G_coo, A_coo;    /* possibly un-used */

    /* allocate socp data structure */
    qc_socp * data = (qc_socp *) calloc(1, sizeof(qc_socp));
    if (!data) return qc_socp_free(data);

    /* allocate problem dimensions */
    data->p = 2;
    data->m = 0;
    data->n = 4;

    /* allocate c vector */
    data->c = (double *) calloc(data->n, sizeof(double));
    if (!data->c) return qc_socp_free(data);

    /* allocate h vector */
    data->h = NULL;

    /* allocate b vector */
    data->b = (double *) calloc(data->p, sizeof(double));
    if (!data->b) return qc_socp_free(data);

    /* allocate G matrix */
    nnzG = 0;
    data->Gx = NULL;
    data->Gp = NULL;
    data->Gi = NULL;
    G_data_ptr = data->Gx;
    G_row_ptr = data->Gi;
    G_col_ptr = data->Gp;

    /* allocate A matrix */
    nnzA = 1*result->nnz + 2;
    data->Ax = (double *) malloc(nnzA * sizeof(double));
    data->Ap = (long *) malloc(nnzA * sizeof(long));
    data->Ai = (long *) malloc(nnzA * sizeof(long));
    if ((!data->Ax) || (!data->Ap) || (!data->Ai)) return qc_socp_free(data);
    A_data_ptr = data->Ax;
    A_row_ptr = data->Ai;
    A_col_ptr = data->Ap;

    /* allocate the cone sizes */
    data->l = 0;
    data->nsoc = 0;
    data->q = NULL;

    /* stuffing the objective vector */
    /* for the constraint _t0 + -1*2.0*D*x == 0 */
    for(i = 0; i < 2; ++i) *A_row_ptr++ = i;
    for(i = 2; i < 4; ++i) *A_col_ptr++ = i;
    for(i = 0; i < 2; ++i) *A_data_ptr++ = 1;
    for(i = 0; i < result->nnz; ++i) *A_row_ptr++ = result->i[i];
    for(i = 0; i < result->nnz; ++i) *A_col_ptr++ = result->j[i];
    for(i = 0; i < result->nnz; ++i) *A_data_ptr++ = result->v[i];

    for(i = 0; i < 2; ++i) data->c[i + 2] = 1;

    /* convert G and A ptrs into a qc_matrix */
    A_coo.m = data->p; A_coo.n = data->n; A_coo.nnz = nnzA;
    A_coo.i = data->Ai;
    A_coo.j = data->Ap;
    A_coo.v = data->Ax;

    /* convert the matrices to column compressed form */
    A_csc = qc_compress(&A_coo);
    if (!A_csc) return qc_socp_free(data);
    /* reassign into data, pointer now owned by data */
    data->Ai = A_csc->i;
    data->Ap = A_csc->j;
    data->Ax = A_csc->v;

    return data;
}

void qc_socp2bug5a(double * x, bug5a_vars * vars, const bug5a_dims * dims)
{
    /*
     * recovers the problem variables from the solver variable 'x'
     * assumes the variables struct is externally allocated
     * the user must keep track of the variable length;
     */
    vars->x = x + 0;  /* length 2 */
}
/* ------------------------ END GENERATED CODE ---------------------------- */

In the generated code, variable result is used, but it is not declared anywhere.

alexvoronov commented 10 years ago

I tried a few more examples, with more parameters. It seems like the parameters that are multiplied by constants get result variable, and other parameters 'get' proper variables in the generated code.

Here is an example with two parameters, D being multiplied by 2, and C not:

#!/usr/bin/env python
from qcml import QCML

if __name__ == '__main__':

    s = ("""variable x(2)
        parameter D(2,2) positive
        parameter C(2,2) positive
        minimize sum(2*D*x + C*x)
        """)

    p = QCML(debug=False)
    p.parse(s)

    p.canonicalize()
    p.codegen("C")
    p.save("bug5b")

In the generated code, D is referred to as result, and C as params->C:

    /* allocate A matrix */
    nnzA = 1*result->nnz + 1*params->C->nnz + 4;
    data->Ax = (double *) malloc(nnzA * sizeof(double));
    data->Ap = (long *) malloc(nnzA * sizeof(long));
    data->Ai = (long *) malloc(nnzA * sizeof(long));
    if ((!data->Ax) || (!data->Ap) || (!data->Ai)) return qc_socp_free(data);
    A_data_ptr = data->Ax;
    A_row_ptr = data->Ai;
    A_col_ptr = data->Ap;

    /* allocate the cone sizes */
    data->l = 0;
    data->nsoc = 0;
    data->q = NULL;

    /* stuffing the objective vector */
    /* for the constraint _t0 + -1*2.0*D*x == 0 */
    for(i = 0; i < 2; ++i) *A_row_ptr++ = i;
    for(i = 2; i < 4; ++i) *A_col_ptr++ = i;
    for(i = 0; i < 2; ++i) *A_data_ptr++ = 1;
    for(i = 0; i < result->nnz; ++i) *A_row_ptr++ = result->i[i];
    for(i = 0; i < result->nnz; ++i) *A_col_ptr++ = result->j[i];
    for(i = 0; i < result->nnz; ++i) *A_data_ptr++ = result->v[i];

    /* for the constraint _t1 + -1*C*x == 0 */
    for(i = 2; i < 4; ++i) *A_row_ptr++ = i;
    for(i = 4; i < 6; ++i) *A_col_ptr++ = i;
    for(i = 0; i < 2; ++i) *A_data_ptr++ = 1;
    for(i = 0; i < params->C->nnz; ++i) *A_row_ptr++ = 2 + 1*(params->C->i[i]);
    for(i = 0; i < params->C->nnz; ++i) *A_col_ptr++ = params->C->j[i];
    for(i = 0; i < params->C->nnz; ++i) *A_data_ptr++ = -params->C->v[i];

Another example, with three parameters:

#!/usr/bin/env python
from qcml import QCML

if __name__ == '__main__':

    s = ("""variable x(2)
        parameter D(2,2) positive
        parameter C(2,2) positive
        parameter B(2,2) positive
    minimize sum(2*D*x + C*x + B*x)
        """)

    p = QCML(debug=False)
    p.parse(s)

    p.canonicalize()
    p.codegen("C")
    p.save("bug5c")

Generated code:

    /* allocate A matrix */
    nnzA = 1*params->B->nnz + 1*result->nnz + 1*params->C->nnz + 6;
    data->Ax = (double *) malloc(nnzA * sizeof(double));
    data->Ap = (long *) malloc(nnzA * sizeof(long));
    data->Ai = (long *) malloc(nnzA * sizeof(long));
    if ((!data->Ax) || (!data->Ap) || (!data->Ai)) return qc_socp_free(data);
    A_data_ptr = data->Ax;
    A_row_ptr = data->Ai;
    A_col_ptr = data->Ap;

    /* allocate the cone sizes */
    data->l = 0;
    data->nsoc = 0;
    data->q = NULL;

    /* stuffing the objective vector */
    /* for the constraint _t0 + -1*C*x == 0 */
    for(i = 0; i < 2; ++i) *A_row_ptr++ = i;
    for(i = 2; i < 4; ++i) *A_col_ptr++ = i;
    for(i = 0; i < 2; ++i) *A_data_ptr++ = 1;
    for(i = 0; i < params->C->nnz; ++i) *A_row_ptr++ = params->C->i[i];
    for(i = 0; i < params->C->nnz; ++i) *A_col_ptr++ = params->C->j[i];
    for(i = 0; i < params->C->nnz; ++i) *A_data_ptr++ = -params->C->v[i];

    /* for the constraint _t1 + -1*B*x == 0 */
    for(i = 2; i < 4; ++i) *A_row_ptr++ = i;
    for(i = 4; i < 6; ++i) *A_col_ptr++ = i;
    for(i = 0; i < 2; ++i) *A_data_ptr++ = 1;
    for(i = 0; i < params->B->nnz; ++i) *A_row_ptr++ = 2 + 1*(params->B->i[i]);
    for(i = 0; i < params->B->nnz; ++i) *A_col_ptr++ = params->B->j[i];
    for(i = 0; i < params->B->nnz; ++i) *A_data_ptr++ = -params->B->v[i];

    /* for the constraint _t2 + -1*2.0*D*x == 0 */
    for(i = 0; i < result->nnz; ++i) *A_row_ptr++ = 4 + 1*(result->i[i]);
    for(i = 0; i < result->nnz; ++i) *A_col_ptr++ = result->j[i];
    for(i = 0; i < result->nnz; ++i) *A_data_ptr++ = result->v[i];
    for(i = 4; i < 6; ++i) *A_row_ptr++ = i;
    for(i = 6; i < 8; ++i) *A_col_ptr++ = i;
    for(i = 0; i < 2; ++i) *A_data_ptr++ = 1;

In both cases, matrix D is being referred to as result.

The order in which I declare matrices B, C, D does not matter. What matters seems to be the multiplication by a constant. Here is an example with two matrices being multiplied by a constant:

#!/usr/bin/env python
from qcml import QCML

if __name__ == '__main__':

    s = ("""variable x(2)
        parameter C(2,2) positive
        parameter B(2,2) positive
        parameter D(2,2) positive
        minimize sum(2 * D * x + C * x + 3 * B * x)
        """)

    p = QCML(debug=False)
    p.parse(s)

    p.canonicalize()
    p.codegen("C")
    p.save("bug5e")

Generated code:

    /* allocate A matrix */
    nnzA = 2*result->nnz + 1*params->C->nnz + 6;
    data->Ax = (double *) malloc(nnzA * sizeof(double));
    data->Ap = (long *) malloc(nnzA * sizeof(long));
    data->Ai = (long *) malloc(nnzA * sizeof(long));
    if ((!data->Ax) || (!data->Ap) || (!data->Ai)) return qc_socp_free(data);
    A_data_ptr = data->Ax;
    A_row_ptr = data->Ai;
    A_col_ptr = data->Ap;

    /* allocate the cone sizes */
    data->l = 0;
    data->nsoc = 0;
    data->q = NULL;

    /* stuffing the objective vector */
    /* for the constraint _t0 + -1*C*x == 0 */
    for(i = 0; i < 2; ++i) *A_row_ptr++ = i;
    for(i = 2; i < 4; ++i) *A_col_ptr++ = i;
    for(i = 0; i < 2; ++i) *A_data_ptr++ = 1;
    for(i = 0; i < params->C->nnz; ++i) *A_row_ptr++ = params->C->i[i];
    for(i = 0; i < params->C->nnz; ++i) *A_col_ptr++ = params->C->j[i];
    for(i = 0; i < params->C->nnz; ++i) *A_data_ptr++ = -params->C->v[i];

    /* for the constraint _t1 + -1*3.0*B*x == 0 */
    for(i = 2; i < 4; ++i) *A_row_ptr++ = i;
    for(i = 4; i < 6; ++i) *A_col_ptr++ = i;
    for(i = 0; i < 2; ++i) *A_data_ptr++ = 1;
    for(i = 0; i < result->nnz; ++i) *A_row_ptr++ = 2 + 1*(result->i[i]);
    for(i = 0; i < result->nnz; ++i) *A_col_ptr++ = result->j[i];
    for(i = 0; i < result->nnz; ++i) *A_data_ptr++ = result->v[i];

    /* for the constraint _t2 + -1*2.0*D*x == 0 */
    for(i = 0; i < result->nnz; ++i) *A_row_ptr++ = 4 + 1*(result->i[i]);
    for(i = 0; i < result->nnz; ++i) *A_col_ptr++ = result->j[i];
    for(i = 0; i < result->nnz; ++i) *A_data_ptr++ = result->v[i];
    for(i = 4; i < 6; ++i) *A_row_ptr++ = i;
    for(i = 6; i < 8; ++i) *A_col_ptr++ = i;
    for(i = 0; i < 2; ++i) *A_data_ptr++ = 1;

Here, both B and D have result in the generated code.

echu commented 10 years ago

okay. give it a shot now. you were right, for the generated code, any expression of the form const*D*x would cause problems. in fact, it was anything of the form scalar*D*x (since i only introduce new variables when multiplying two matrix expressions).

this now works:

variable x(2)
parameter D(2,2) positive
parameter C(2,2) positive
parameter B(2,2) positive
parameter a positive
minimize sum(a*D*x + C*x + B*x)

the generated code:

1    qc_socp * qc_{name}2socp(const {name}_params * params, const {name}_dims * dims)
2    {
3        /*
4         * maps 'params' into the C socp data type
5         * 'params' ought to contain:
6         *   'a' has shape Scalar()
7         *   'C' has shape Matrix(2,2)
8         *   'B' has shape Matrix(2,2)
9         *   'D' has shape Matrix(2,2)
10        */
11   
12       /* all local variables */
13       long i;  /* loop index */
14       long *q_ptr;
15       long *A_row_ptr, *A_col_ptr;
16       long *G_row_ptr, *G_col_ptr;
17       double *A_data_ptr, *G_data_ptr;
18       long nnzA, nnzG;
19       qc_matrix *G_csc, *A_csc;  /* possibly un-used */
20       qc_matrix G_coo, A_coo;    /* possibly un-used */
21   
22       /* allocate socp data structure */
23       qc_socp * data = (qc_socp *) calloc(1, sizeof(qc_socp));
24       if (!data) return qc_socp_free(data);
25   
26       /* allocate problem dimensions */
27       data->p = 6;
28       data->m = 0;
29       data->n = 8;
30   
31       /* allocate c vector */
32       data->c = (double *) calloc(data->n, sizeof(double));
33       if (!data->c) return qc_socp_free(data);
34   
35       /* allocate h vector */
36       data->h = NULL;
37   
38       /* allocate b vector */
39       data->b = (double *) calloc(data->p, sizeof(double));
40       if (!data->b) return qc_socp_free(data);
41   
42       /* allocate G matrix */
43       nnzG = 0;
44       data->Gx = NULL;
45       data->Gp = NULL;
46       data->Gi = NULL;
47       G_data_ptr = data->Gx;
48       G_row_ptr = data->Gi;
49       G_col_ptr = data->Gp;
50   
51       /* allocate A matrix */
52       nnzA = 1*params->B->nnz + 1*params->D->nnz + 1*params->C->nnz + 6;
53       data->Ax = (double *) malloc(nnzA * sizeof(double));
54       data->Ap = (long *) malloc(nnzA * sizeof(long));
55       data->Ai = (long *) malloc(nnzA * sizeof(long));
56       if ((!data->Ax) || (!data->Ap) || (!data->Ai)) return qc_socp_free(data);
57       A_data_ptr = data->Ax;
58       A_row_ptr = data->Ai;
59       A_col_ptr = data->Ap;
60   
61       /* allocate the cone sizes */
62       data->l = 0;
63       data->nsoc = 0;
64       data->q = NULL;
65   
66       /* stuffing the objective vector */
67       /* for the constraint _t0 + -1*C*x == 0 */
68       for(i = 0; i < 2; ++i) *A_row_ptr++ = i;
69       for(i = 2; i < 4; ++i) *A_col_ptr++ = i;
70       for(i = 0; i < 2; ++i) *A_data_ptr++ = 1;
71       for(i = 0; i < params->C->nnz; ++i) *A_row_ptr++ = params->C->i[i];
72       for(i = 0; i < params->C->nnz; ++i) *A_col_ptr++ = params->C->j[i];
73       for(i = 0; i < params->C->nnz; ++i) *A_data_ptr++ = -params->C->v[i];
74   
75       /* for the constraint _t1 + -1*B*x == 0 */
76       for(i = 2; i < 4; ++i) *A_row_ptr++ = i;
77       for(i = 4; i < 6; ++i) *A_col_ptr++ = i;
78       for(i = 0; i < 2; ++i) *A_data_ptr++ = 1;
79       for(i = 0; i < params->B->nnz; ++i) *A_row_ptr++ = 2 + 1*(params->B->i[i]);
80       for(i = 0; i < params->B->nnz; ++i) *A_col_ptr++ = params->B->j[i];
81       for(i = 0; i < params->B->nnz; ++i) *A_data_ptr++ = -params->B->v[i];
82   
83       /* for the constraint _t2 + -1*D*x == 0 */
84       for(i = 0; i < params->D->nnz; ++i) *A_row_ptr++ = 4 + 1*(params->D->i[i]);
85       for(i = 0; i < params->D->nnz; ++i) *A_col_ptr++ = params->D->j[i];
86       for(i = 0; i < params->D->nnz; ++i) *A_data_ptr++ = -params->D->v[i];
87       for(i = 4; i < 6; ++i) *A_row_ptr++ = i;
88       for(i = 6; i < 8; ++i) *A_col_ptr++ = i;
89       for(i = 0; i < 2; ++i) *A_data_ptr++ = 1;
90   
91       for(i = 0; i < 2; ++i) data->c[i + 2] = 1;
92       for(i = 0; i < 2; ++i) data->c[i + 4] = 1;
93       for(i = 0; i < 2; ++i) data->c[i + 6] = params->a;
94   
95       /* convert G and A ptrs into a qc_matrix */
96       A_coo.m = data->p; A_coo.n = data->n; A_coo.nnz = nnzA;
97       A_coo.i = data->Ai;
98       A_coo.j = data->Ap;
99       A_coo.v = data->Ax;
100  
101      /* convert the matrices to column compressed form */
102      A_csc = qc_compress(&A_coo);
103      if (!A_csc) return qc_socp_free(data);
104      /* reassign into data, pointer now owned by data */
105      data->Ai = A_csc->i;
106      data->Ap = A_csc->j;
107      data->Ax = A_csc->v;
108  
109      return data;
110  }

this is mostly fixed, but i'm aware that it's still broken for something that looks like:

variable x(2)
parameter D(2,2) positive
parameter C(2,2) positive
parameter B(2,2) positive
parameter a positive
minimize sum(a*D*x + C*x + B*x)
2*C*x == 0

(The constraint causes some problems.) I'll be working on a fix over the next week.

echu commented 10 years ago

Guess I had some extra time. :)

i think i've fixed all these issues now. let me know if it still persists.

alexvoronov commented 10 years ago

Thanks a lot! Those multiplications work now!

But I found another one :)

you will not believe how much can be revealed by reimplementing an old exercise in a new system :)

So here is an example that will also get result in the generated code:

#!/usr/bin/env python
from qcml import QCML

if __name__ == '__main__':

    s = ("""variable x(2)
        parameter D(2,2) positive
        minimize sum(x)
        subject to
        x - D * x == 0
        """)

    p = QCML(debug=False)
    p.parse(s)

    p.canonicalize()
    p.codegen("C")
    p.save("gen_bug")

Part of generated code:

    /* allocate A matrix */
    nnzA = 1*result->nnz;
    data->Ax = (double *) malloc(nnzA * sizeof(double));
    data->Ap = (long *) malloc(nnzA * sizeof(long));
    data->Ai = (long *) malloc(nnzA * sizeof(long));
    if ((!data->Ax) || (!data->Ap) || (!data->Ai)) return qc_socp_free(data);
    A_data_ptr = data->Ax;
    A_row_ptr = data->Ai;
    A_col_ptr = data->Ap;

    /* allocate the cone sizes */
    data->l = 0;
    data->nsoc = 0;
    data->q = NULL;

    /* stuffing the objective vector */
    for(i = 0; i < 2; ++i) data->c[i + 0] = 1;

    /* for the constraint x + -1*D*x == 0 */
    for(i = 0; i < result->nnz; ++i) *A_row_ptr++ = result->i[i];
    for(i = 0; i < result->nnz; ++i) *A_col_ptr++ = result->j[i];
    for(i = 0; i < result->nnz; ++i) *A_data_ptr++ = result->v[i];
echu commented 10 years ago

I added another test case for this. Give it a try now. Let me know if you find any other issues! Thanks so much!

alexvoronov commented 10 years ago

Thanks, the generated code for my example compiles now!

I get runtime errors now (segmentation fault: EXC_BAD_ACCESS in qcml_utils.c:137), but I suspect they come from my side, probably I did something wrong with qc_matrix. I'll open a new issue if I'm stuck :)

Edit: yep, I had a bug in my memory allocation for qc_matrix.