search

odlgroup / odl

Operator Discretization Library https://odlgroup.github.io/odl/
Mozilla Public License 2.0
375 stars 105 forks source link

BroadcastOperator doesn't work with scipy > 1.8.1 #1626

Open Priscilla- opened 1 year ago

Priscilla- commented 1 year ago

Hello,

I'm reproducing the reconstruction example in https://odlgroup.github.io/odl/guide/pdhg_guide.html, and I found out that the BroadcastOperator raises an Error related to the format and sparsity of the data when computing the L operator. I'm on python3.10.1, scipy1.10 and numpy1.23.5 on an OSX. I've tried downgrading scipy to the v1.8.1, which solves the issue with the Broadcasting operator, but other problems with the data types to construct the phantom appear instead. I've tried hacking scipy commenting the lines 127-129 in "/scipy/sparse/_sputils.py":

    newdtype = np.dtype(dtype)
    # Commented the section below to fix problem with ODL's broadcasting operator
    #if newdtype == np.object_:
       # raise ValueError(
        #    "object dtype is not supported by sparse matrices"
        #)

And now the Broadcasting operator works.

Notice that scipy1.10.0 is installed with ODL.

sbanert commented 9 months ago

This seems to be related to this change in scipy: https://github.com/scipy/scipy/pull/15828 Apparently, ODL represents operators between product spaces with scipy's sparse matrix data structure (a sparse matrix of operators, that is). However, scipy only intends the usage for numeric data, which our operators not really are.

mehrhardt commented 9 months ago

What is the best resolution for ODL then? Change the operator representation to lists?

ozanoktem commented 9 months ago

This is a highly prioritised issue that needs to be resolved. The original solution (which works with scipy v1.6, the problem most likely appeared when scipy went from v1.8.1 to v1.9) generated a warning since the ODL way of using sparse matrices was not the usage scipy intended. Most likely we will need to refactor ODL as to not rely on scipy at all, e.g., consider operator representation as lists.

paulhausner commented 9 months ago

I took a quick look at this. Maybe the easiest way to resolve this is using a custom data strucutre for the product space operators which has the same signature as scipy's coo matrix. See the opened pull request for a first implementation