Inconsistency on A = trans(B) mapping

[ptillet]

Hey!

I'm getting more and more familiar with PyViennaCL, and I am doing the necessary changes to make PyViennaCL compatible with the latest ViennaCL development version (which makes row_major a dynamic parameter inside matrix_base).

Anyway, I've noticed inside dense_matrix.h (inside the matrix_base wrapper) the following:

.add_property("trans", pyvcl_do_1ary_op<vcl::matrix<TYPE, F>, vcl::matrix_base<TYPE,F>&, op_trans, 0>)

This is slightly inconsistent, if I'm not mistaken, with the original ViennaCL operator:

  self_type & matrix_base::operator=(const matrix_expression< const self_type, const self_type, op_trans> & proxy)

Which seems to allow ranges and slices to be assigned to transpositions. I see two fixes here, we could either have a trans() method in each matrix{ /_range/_slice} wrapper, or have an external p.trans(A) function. I think that the former option is better, since it doesn't induce any interface change

[tsmithe]

You're quite right, and I agree that the former option is better -- but of course, it doesn't really matter, since there are a couple of odd functions that are provided as 'external' p._viennacl.do_x(...) functions, rather than as class methods, and these are always wrapped by the thin python layer on top of _viennacl. It doesn't really matter if _viennacl changes a bit, because its only consumer is pyviennacl itself :)

Anyhow, thanks for the work -- I'm a bit swamped by other things at the moment, but will try and get back to your other report soon!

[tsmithe]

Also note that assignations in Python such as

A = B

are not assignations, but reference re-assignments. So the equals operator doesn't work like it does in C++, and you'd have to do something like

A[:] = B

to get what you mean by A = B in C++.

[tsmithe]

(This also means that to allow range and slice transposition assignments a bit of code might be necessary in the Python layer of PyViennaCL, but I would need to investigate that a bit more..)

[ptillet]

Actually, the problem seems to be deeper,. While matrix_base can be constructed from a matrix_expression<>, there is no such mechanism for matrix_range and matrix_slice. An additional constructor

matrix_range::matrix_range(matrix_expression<> expr) : base_type(expr){ }

Could be easily added, but I'm not sure whether it's safe or not. Maybe @karlrupp could help on this :P

[karlrupp]

Adding such a constructor would break the current semantics of a matrix_range and matrix_slice: Right now they represent a range/slice of a physically existing memory buffer, without owning the buffer. For an expression, this cannot be fulfilled.

If we want to change the semantics, then this would have quite a large impact on how expressions are mapped to kernels. I suggest to do this only after the scheduler has been fully integrated.

[ptillet]

Okay, so we should probably leave it as it is now, and let ranges and slices be un-usable as LHSes in python. This is not a major problem anyway.

[karlrupp]

Hmm, this doesn't sound like a good solution either. Philippe, could you please explain the issue again? Your introductory code snippet seems to be incomplete, some template arguments are probably interpreted as HTML.

[tsmithe]

Currently, it should be possible to use ranges and slices as LHSes -- I'm sure I have tested this previously, so I'm not entirely sure where the fault is. As I mentioned before, you have to do assignment in the slightly roundabout way

A[:] = B

I just tested this again and there seems to be a pair of (different) bugs anyway: one in using trans; and another in computing the sizes for an assignment as above:

>>> b[:] = c
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python2.7/dist-packages/pyviennacl/pycore.py", line 851, in __setitem__
    (item.shape, value.shape))
TypeError: Cannot assign across different shapes! ((3, 1) and (3, 2))

I'll get to these in the next day or so!

[ptillet]

Yes, it seems like ranges can be used as LHSes, I just didn't have investigated the code enough. My initial problem was with pyvcl_op::do_op() being unable to use a range as ReturnT (since it constructs an instance of ReturnT from a template matrix/vector_expression), while range = matrix_expression<>, seems to be allowed in C++. I couldn't find how range = matrix_expression<> is performed in python, but it seems like it is implemented indeed. Is there any additional copy going on?

2014-05-09 11:44 GMT+02:00 Toby Smithe notifications@github.com:

Currently, it should be possible to use ranges and slices as LHSes -- I'm sure I have tested this previously, so I'm not entirely sure where the fault is. As I mentioned before, you have to do assignment in the slightly roundabout way

A[:] = B

but I just test this and there seems to be a pair of (different) bugs anyway: one in using trans; and another in computing the sizes for an assignment as above:

b[:] = c Traceback (most recent call last): File "", line 1, in File "/usr/lib/python2.7/dist-packages/pyviennacl/pycore.py", line 851, in setitem (item.shape, value.shape)) TypeError: Cannot assign across different shapes! ((3, 1) and (3, 2))

I'll get to these in the next day or so!

— Reply to this email directly or view it on GitHubhttps://github.com/viennacl/pyviennacl-dev/issues/16#issuecomment-42649475 .

[tsmithe]

Ah, I wouldn't worry too much about the pyvcl_op calls -- they're only used for ViennaCL functions that don't have a scheduler call, such as solvers, and various miscellaneous functions (see extra_functions.cpp). For things like assignment, PyViennaCL uses the scheduler, constructing an 'Assign' class, which then takes any expression on the RHS and constructs an expression statement object. So the ReturnT and operator= don't matter in this case :-)

[ptillet]

Okay ! Perfect then, I can close the issue.

2014-05-09 12:24 GMT+02:00 Toby Smithe notifications@github.com:

Ah, I wouldn't worry too much about the pyvcl_op calls -- they're only used for ViennaCL functions that don't have a scheduler call, such as solvers, and various miscellaneous functions (see extra_functions.cpp). For things like assignment, PyViennaCL uses the scheduler, constructing an 'Assign' class, which then takes any expression on the RHS and constructs an expression statement object. So the ReturnT and operator= don't matter in this case :-)

— Reply to this email directly or view it on GitHubhttps://github.com/viennacl/pyviennacl-dev/issues/16#issuecomment-42652334 .

[tsmithe]

Excellent. It's been quite a productive discussion, too!

[tsmithe]

FYI, now I've fixed issue #18, you can again use a (dense) range / slice object as an LHS in an assignment and see the changes in the original object:

>>> import pyviennacl as p
>>> a = p.Matrix(5,5,0.2)
>>> b = a[0:3, 2:4]
>>> c = p.Matrix(3,2, 66.6)
>>> b[:] = c
>>> a
array([[  0.2,   0.2,  66.6,  66.6,   0.2],
       [  0.2,   0.2,  66.6,  66.6,   0.2],
       [  0.2,   0.2,  66.6,  66.6,   0.2],
       [  0.2,   0.2,   0.2,   0.2,   0.2],
       [  0.2,   0.2,   0.2,   0.2,   0.2]])

Once we've got the scheduler trans issue fixed and lanczos building on Windows, I'll put together a new point release. "Release early, release often," they say.. Hopefully this time I'll manage to get everything coordinated sufficiently to make a proper announcement!

At the beginning of my GSoC work, I'll probably branch current master into a 1.0 branch, as per usual practice, but I don't envisage making any big code changes till then.

[ptillet]

Great!

Tell me when you will branch the current master inside pyviennacl, I have already done the proper compatibility fixes on my personnal pyviennacl fork ;-)

2014-05-09 23:00 GMT+02:00 Toby Smithe notifications@github.com:

FYI, now I've fixed issue 18, you can now use a range / slice object as an LHS in an assignment and see the changes in the original object:

import pyviennacl as p>>> a = p.Matrix(5,5,0.2)>>> b = a[0:3, 2:4]>>> c = p.Matrix(3,2, 66.6)>>> b[:] = c>>> aarray([[ 0.2, 0.2, 66.6, 66.6, 0.2], [ 0.2, 0.2, 66.6, 66.6, 0.2], [ 0.2, 0.2, 66.6, 66.6, 0.2], [ 0.2, 0.2, 0.2, 0.2, 0.2], [ 0.2, 0.2, 0.2, 0.2, 0.2]])

Once we've got the scheduler trans issue fixed and lanczos building on Windows, I'll put together a new point release. "Release early, release often," they say.. Hopefully this time I'll manage to get everything coordinated sufficiently to make a proper announcement!

At the beginning of my GSoC work, I'll probably branch current master into a 1.0 branch, as per usual practice, but I don't envisage making any big code changes till then.

— Reply to this email directly or view it on GitHubhttps://github.com/viennacl/pyviennacl-dev/issues/16#issuecomment-42713342 .

[tsmithe]

I will. I want to keep the 1.0 branch following the 1.5 branch of ViennaCL -- though I might backport ViennaCL fixes if there aren't any more actual releases. @karlrupp -- do you think a ViennaCL 1.5.2 release would in order, including the build-Lanczos-on-MSVC and trans scheduler fixes? Or should I just backport those fixes manually?

[karlrupp]

@tsmithe Yes, a bugfix release is appropriate to have PyViennaCL built on top of a release version (i.e. 'defined state'). I'll take care of this - just allow for a couple of days as I'm about to leave for a conference.

[tsmithe]

Great -- thanks!