Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

iscomplexobj(BlockArray) after #259 #281

Closed
tbalke opened this issue Apr 23, 2022 · 0 comments · Fixed by #282
Closed

iscomplexobj(BlockArray) after #259 #281

tbalke opened this issue Apr 23, 2022 · 0 comments · Fixed by #282
Assignees
Labels
bug Something isn't working

Comments

@tbalke
Copy link
Contributor

tbalke commented Apr 23, 2022

I ran into an error when running examples/scripts/ct_astra_weighted_tv_admm.py.
It seems to be the case that in

if snp.iscomplexobj(v):
out = snp.exp(1j * snp.angle(v)) * tmp
else:
out = snp.sign(v) * tmp

the if-clause gets executed (because v is a BlockArray, see below) and the arrays get converted to complex64, causing the following output:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~/pythonModules/scico/examples/scripts/ct_astra_weighted_tv_admm.py in <module>
    123 )
    124 print(f"Solving on {device_info()}\n")
--> 125 admm_unweighted.solve()
    126 x_unweighted = postprocess(admm_unweighted.x)
    127 

~/pythonModules/scico/scico/optimize/admm.py in solve(self, callback)
    679         for self.itnum in range(self.itnum, self.itnum + self.maxiter):
    680             self.step()
--> 681             self.itstat_object.insert(self.itstat_insert_func(self))
    682             if callback:
    683                 self.timer.stop()

<string> in itstat_func(obj)

~/pythonModules/scico/scico/optimize/admm.py in norm_dual_residual(self)
    579         out = 0.0
    580         for zi, ziold, Ci in zip(self.z_list, self.z_list_old, self.C_list):
--> 581             out += norm(Ci.adj(zi - ziold)) ** 2
    582         return snp.sqrt(out)
    583 

~/pythonModules/scico/scico/_generic_operators.py in adj(self, y)
    600             return ComposedLinearOperator(self.H, y)
    601         if self.output_dtype != y.dtype:
--> 602             raise ValueError(f"dtype error: expected {self.output_dtype}, got {y.dtype}")
    603         if self.output_shape != y.shape:
    604             raise ValueError(

ValueError: dtype error: expected float32, got complex64

This is an issue that only arises after #259.

  • The first difference to earlier commits is that in
    zi = gi.prox(Cix + ui, 1 / rhoi, v0=zi)

    zi is a BlockArray, whereas in previous commits it is not. Why?
  • Secondly, the handling of the snp.iscomplexobj is not the same as shown below. Is this desired?

git checkout 4ac440a (before Simplify BlockArray implementation (#259))

In [22]: a = BlockArray.array([snp.zeros(2), snp.ones(2)])

In [23]: a
Out[23]: 
scico.blockarray.BlockArray: 
DeviceArray([0., 0., 1., 1.], dtype=float32)

In [24]: snp.iscomplexobj(a)
Out[24]: False

git checkout 4381be5 (after Simplify BlockArray implementation (#259))

In [2]: a = snp.blockarray([snp.zeros(2), snp.ones(2)])

In [3]: a
Out[3]: [DeviceArray([0., 0.], dtype=float32), DeviceArray([1., 1.], dtype=float32)]

In [4]: snp.iscomplexobj(a)
Out[4]: [DeviceArray(False, dtype=bool), DeviceArray(False, dtype=bool)]

In [5]: bool(snp.iscomplexobj(a))
Out[5]: True

@bwohlberg

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants