Checking updated PR...

No PEP8 issues.

Comment last updated at 2024-02-19 07:48:06 UTC

Has somebody investigated yet whether switching to scipy's dok_matrix instead would also work?

After some thinking (and talking with @JevgenijaAksjonova) I would go further and say it's probably most reasonable to not use a dedicated sparse-matrix class at all, but simply a dictionary associating index-tuples to operators. This makes __getitem__ trivial. The other needed operations are also easy, though they will be somewhat less efficient than with something based on lists of indices and values.

Why I hesitate about making a dedicated class: this won't be a proper sparse-matrix implementation unless a lot more work is put into it. One would in principle need to think about questions like how duplicate and out-of-order indices are to be handled. Getting this wrong could lead to strange and subtle bugs. With a plain dict, these issues don't arise.

Performance / memory concerns should not really matter here, assuming that ProductSpace is typically used with just a few spaces, accordingly just a few operators in the product operator, but these operators are so heavy as to amortize the overhead of the dict. That's a very different situation from scipy.sparse's use cases.

Maybe I'm wrong and some users need products over lots of spaces, though that would sound more like some tensor product space.

One setting where one might end up with a constructing operators that acts on products over lots of spaces is in Kaczmarz method. The idea is to apply the forward operator on subsets of the data, which in tomography is typically data acquired when the source is in a specific position. Thus, one easily ends up with products over hundreds (if not thousands) of spaces. ODL has an implementation of Kaczmarz, see here for the documentation and here for an example usage.

ProductSpaceOperator is not used in the example with Kaczmarz method.

@ozanoktem I get the point: the Kaczmarz method conceptually always inverts a large BroadcastOperator, which is under the hood a dense (!), single-column ProductSpaceOperator. Though, as Jevgenija said, the ODL Kaczmarz method does not currently express it this way, it makes sense to provide for that option, if only for convenience (as in the example, which uses the BroadcastOperator for a single forward pass).

In that case using a dict-based implementation does have a more significant overhead, compared to the cheap RayTransforms. I still tend to say that in cases where this matters something more fundamentally different should be done, and at any rate a COO sparse matrix is also not optimal for the dense BroadcastOperator, but it is certainly a better compromise.

There is still an alternative to be considered to having a custom COO class here: since the cases BroadcastOperator and ReductionOperator are always dense matrices, they could either just not use ProductSpaceOperator underneath, or ProductSpaceOperator could special-case n=1 and m=1 and use simple array storage for these.

All of that requires extra code though, so I tend to agree to going with the custom-class approach, then.

Some tests are failing on this branch:

space/pspace_test.py::test_element_setitem_fancy PASSED                                                                                                                                [ 54%]
space/pspace_test.py::test_element_setitem_broadcast PASSED                                                                                                                            [ 55%]
space/pspace_test.py::test_unary_ops FAILED                                                                                                                                            [ 56%]
space/pspace_test.py::test_operators[ op = '+' ] FAILED                                                                                                                                [ 57%]
space/pspace_test.py::test_operators[ op = '/' ] FAILED                                                                                                                                [ 58%]
space/pspace_test.py::test_operators[ op = '*' ] FAILED                                                                                                                                [ 60%]
space/pspace_test.py::test_operators[ op = '-' ] FAILED                                                                                                                                [ 61%]
space/pspace_test.py::test_operators[ op = '+=' ] FAILED                                                                                                                               [ 62%]
space/pspace_test.py::test_operators[ op = '/=' ] FAILED                                                                                                                               [ 63%]
space/pspace_test.py::test_operators[ op = '*=' ] FAILED                                                                                                                               [ 64%]
space/pspace_test.py::test_operators[ op = '-=' ] FAILED                                                                                                                               [ 65%]
space/pspace_test.py::test_ufuncs PASSED                                                                                                                                               [ 67%]
space/pspace_test.py::test_reductions PASSED                                                                                                                                           [ 68%]

I'm looking into why this happens, it's probably fairly trivial. Incidentally, space/pspace_test.py doesn't go through even in master (though the test cases above succeed there); there seem to be multiple issues at work.

Ok, the test discrepancies are purely down to the patch

--- a/odl/util/testutils.py
+++ b/odl/util/testutils.py
@@ -323,7 +323,7 @@ def noise_array(space):
     """
     from odl.space import ProductSpace
     if isinstance(space, ProductSpace):
-        return np.array([noise_array(si) for si in space])
+        return np.array([noise_array(si) for si in space], dtype=object)
     else:
         if space.dtype == bool:
             arr = np.random.randint(0, 2, size=space.shape, dtype=bool)

I think it's more the test suite that's at fault here, than the code being tested.

@leftaroundabout Overall, I agree with your assessment. I think the duplicate and out-of-order elements issues also arise with the scipy implementation so the suggested class-based approach should not break any existing code at least. Checks for duplicates and ordering could be implemented pretty cheaply here.

It might be misleading to call the class COOMatrix since as you pointed out it does not implement any proper matrix functionality (even though this might be of interest in the future) but maybe it should be renamed for now to avoid confusion.

@paulhausner

It might be misleading to call the class COOMatrix

I think the name is ok, as long as there's some documentation as to the fact that this is merely a bare-bones container data structure with only one use case.

Then I would support merging this; the pspace_tests need to be fixed separately. I haven't found anything else that would get broken by the change.

Besides making sure unit tests pass, one should also ensure that the examples listed in the links below (many of these rely on constructing product space operators) actually work:

• https://github.com/odlgroup/odl/tree/master/examples/solvers
• https://odlgroup.github.io/odl/getting_started/first_steps.html
• https://github.com/odlgroup/odlworkshop/blob/master/code/part4_spectral_ct.ipynb

I checked that the examples involving product space operators (pdhg_denoising_tgv.py and pdhg_tomography_tgv.py) run. That is, they make it through 10 iterations until crashing with an unrelated AttributeError: 'Colorbar' object has no attribute 'draw_all'; presumably that's #1636.

Can't vouch for correctness of computations. For now, so much in ODL just doesn't work at all that I'm afraid we first need to squash the issues that prevent getting results at all, then when all the examples run and plotting works again we have to scrutinise the results again.

Thanks to the high-level Operator class I'm actually fairly optimistic that if those examples run some iterations without error then they probably are also correct.

I think this should be merged. Although I'm not sure the custom class is really the best from a project-maintenance perspective, it is certainly close to what SciPy used to do, and it just works for the tests and examples we currently have.

There are some things in here that don't really belong to the sparse matrix issue, those are addressed in #1640 and #1638. It's probably easiest to first merge those PRs and then rebase this one on top of them.

Could you please rebase this PR?

Rebased version at https://github.com/leftaroundabout/odl/tree/feature/custom-coo-sparse

I have confirmed that this still succeeds the unit tests and relevant examples.

I can still not run the test example in https://github.com/odlgroup/odl/blob/master/examples/solvers/pdhg_tomography_tgv.py, get error message ValueError: object dtype is not supported by sparse matrices when running lines 74-75 in that example.

The code below for TV reconstruction using Douglas-Rachford that used to work earlier does not work. The problem is the line defining g_funcs, running that yields the error ValueError: shape of inp not equal to space shape: (360, 500) != (300, 300)

import numpy as np
import matplotlib.pyplot as plt
import odl

# --- Set up the forward operator (ray transform) --- #
      
# Reconstruction space: discretized functions on the rectangle
# [-20, 20]x[-20,20] with 300 samples per dimension.
reco_space = odl.uniform_discr(
    min_pt=[-20, -20], max_pt=[20, 20], shape=[300, 300], dtype='float32')

# Create a discrete Shepp-Logan phantom (modified version)
phantom = odl.phantom.shepp_logan(reco_space, modified=True)
# Show phantom
phantom.show(title='Phantom');

# Data space: Make a fan- geometry with flat detector
# Source positions: 360 positions uniformly spaced in 0, 2pi
source_pos = odl.uniform_partition(0, 2 * np.pi, 360)
#  Detector: 500 detectorelements uniformly distributed on [-68,68]
detector_pixels = odl.uniform_partition(-68, 68, 500)
# Aquisition geometry
geometry = odl.tomo.geometry.conebeam.FanBeamGeometry(source_pos, 
                                                      detector_pixels, 
                                                      src_radius=40, 
                                                      det_radius=40)

# Forward operator: Ray transform 
fwd_op = odl.tomo.RayTransform(reco_space, geometry)

# --- Generate artificial data --- #

# Create projection data by calling the ray transform on the phantom
ideal_data = fwd_op(phantom)
data = ideal_data + 0.5 * odl.phantom.white_noise(fwd_op.range) * np.mean(ideal_data)
    
# Show data
data.show(title='Projection Data (Sinogram)');

# --- Reconstruction: Total variation (with Douglas-Rachford) --- #

# Assemble all operators into a list.
alpha = 0.5
grad = odl.Gradient(fwd_op.domain)
lin_ops = [fwd_op, grad]

# Create functionals for the l2 distance and l1 norm.
g_funcs = [odl.solvers.L2NormSquared(fwd_op.domain).translated(data),
           alpha * odl.solvers.L1Norm(grad.range)]

# Functional of the bound constraint 0 <= f <= 1
f = odl.solvers.IndicatorBox(fwd_op.domain, 0, 1)

# Find scaling constants so that the solver converges.
# See the douglas_rachford_pd documentation for more information.
opnorm_fwdOp = odl.power_method_opnorm(fwd_op, xstart=data)
opnorm_grad = odl.power_method_opnorm(grad, xstart=data)
sigma = [1 / opnorm_fwdOp**2, 1 / opnorm_grad**2]
tau = 1.0

# Solve using the Douglas-Rachford Primal-Dual method
x = fwd_op.domain.zero()
odl.solvers.douglas_rachford_pd(x, f, g_funcs, lin_ops,
                                tau=tau, sigma=sigma, niter=100)

# --- Show results --- #
x.show('TV (with Douglas-Rachford)', force_show=True);

@ozanoktem

I can still not run the test example in https://github.com/odlgroup/odl/blob/master/examples/solvers/pdhg_tomography_tgv.py,

Were you using the rebased version? The PDHG examply works fine for me on that branch.

The branch here does not include the other fixes yet.

Pull request #1642 has been merged.