Issue 49 matrix representation of linear operator

[aringh]

Added the function for matrix representation of a linear operator; see issue #49

[aringh]

There seems to be two errors in this build, both being more or less them same. However I'm not sure of what the problem is... where is the test supposed to be put? I thought that the structure of the tests was supposed to be the same as for the implementation, and thus I put it in test/discr/utility_test.py.

Can either @adler-j or @kohr-h have a look at this, whenever you have time?

[kohr-h]

It obviously doesn't like the two identical file names utility_test.py under discr/ and util/. Why don't you put the code in discr_mappings.py? It wasn't intended for that from the beginning, but I think it fits there. As there is no test file for that module (I think), you can just rename your test file.

[aringh]

I renamed it to operator_utilities. Looks ok?

[adler-j]

A short module header goes here. Basically a few rows about what this module contains.

[adler-j]

Almost done, still need discreate -> discrete in quite a few places, and this should be in odl/operator imo, not odl/discr.

[kohr-h]

Almost done, still need discreate -> discrete in quite a few places, and this should be in odl/operator imo, not odl/discr.

Maybe you're right. I suggested discr when I was starting to think that the operator alone would not suffice as information to calculate the matrix representation, but up to implementation errors (not consistent), this is possible. So it can go to operators. Just needs a better name. I'm not too creative now, cannot think of anything.

[aringh]

I moved it to operator (and moved the test accordingly), changed the typos, and changed to use the MatVecOperator instead (also in the operator_test file).

The code works for FnBase only, so for product operators one need to use the methods suggested in #49. The following code works also for ProducSpace (as long as it is not a ProductSpace where one component is another ProductSpace...) Last call: should we keep the simple one or switch for this more complicated code?

def my_matrix_representation(op):
    assert op.is_linear
    assert (isinstance(op.domain, odl.space.base_ntuples.FnBase) or
            isinstance(op.domain, odl.set.pspace.ProductSpace))
    assert (isinstance(op.range, odl.space.base_ntuples.FnBase) or
            isinstance(op.range, odl.set.pspace.ProductSpace))

    # Get the size of the range, and handle ProductSpace
    n = []
    op_ran = op.range
    if isinstance(op_ran, odl.set.pspace.ProductSpace):
        num_ran = op_ran.size
        for i in range(num_ran):
            n.append(op_ran[i].size)
    else:
        num_ran = 1
        n.appned(op_ran.size)

    # Get the size of the domain, and handle ProductSpace
    m = []
    op_dom = op.domain
    if isinstance(op_dom, odl.set.pspace.ProductSpace):
        num_dom = op_dom.size
        for i in range(num_dom):
            m.append(op_dom[i].size)
    else:
        num_dom = 1
        m.appned(op_dom.size)

    # Generate the matrix
    matrix = np.zeros([np.sum(n), np.sum(m)])
    tmp = op_ran.element()
    v = op_dom.element()
    index = 0
    for i in range(num_dom):
        for j in range(m[i]):
            v.set_zero()
            v[i][j] = 1.0
            tmp_two = op(v, out=tmp)
            tmp_three = []
            for k in range(num_ran):
                tmp_three = np.concatenate((tmp_three, tmp_two[k].asarray()))
            matrix[:, index] = tmp_three
            index += 1

    return matrix

[adler-j]

Quite sure that code fails for FnBase, what would v[0][0] be in that case?

Anyway I also see the purpose of having this for ProductSpace, and since we have put some effort in it, why not. Just make sure to add a test for it.

[aringh]

I have rebased this, in order to be able to use ProductSpaceOperator. This new implementation can take operators defined from ProductSpace (or FnBase) to ProductSpace (or FnBase), and there are a few test to try the code in different cases. More complicated code, but easier for the end-user.

[kohr-h]

This can be one instance check: isinstance(op.domain, (FnBase, ProductSpace))

[landscape-bot]

Repository health increased by 0.27% when pulling 872df47 on issue-49__matrix_representation_of_linear_operator into 5d8832f on master.

• 12 new problems were found (including 0 errors and 0 code smells).
• No problems were fixed.

[adler-j]

Its quite obnoxious, is it not? Maybe disabling the comment function (no idea why it was enabled suddenly).

[kohr-h]

Yes, if it gives one comment per push, that's a bit too much. In principle, it can be handy, so we can keep it, I guess. Or maybe codacy is an alternative.

[adler-j]

Fixed the most un-needed warnings. Now it's actually quite useful, if only it wouldn't write actual messages.

[kohr-h]

Almost good, just one blank missing - it will glue like "[...] nor ProductSpacewith only [...]"

[aringh]

Done!

[kohr-h]

Fixed the most un-needed warnings. Now it's actually quite useful, if only it wouldn't write actual messages.

Comment on an old diff? I don't see any more print statements.

[adler-j]

This is on the somewhat parallel landscape messages :)

[kohr-h]

Got it! ;-)

[aringh]

Ok to merge?

[kohr-h]

Looks good, let's do it!

[kohr-h]

Wait, let me fix two small things.

[kohr-h]

Got a nice speedup by 3.3 with these optimizations. Basically trashed the tmp_result and instead of setting the complete tmp_dom back to zero, I only reset the one value I know I changed in the last step.