ENH: Add a tomography example for the adupdates solver.#1416
Conversation
|
Checking updated PR... No PEP8 issues. Comment last updated on August 29, 2018 at 11:28 Hours UTC |
|
I don’t understand the Travis failure. :( |
|
Seems to be this issue. Can you try older versions of pytest? Looks like we need to blacklist a few. |
|
I can confirm that Edit: On my workstation, |
|
Sounds very much like this issue: astra-toolbox/astra-toolbox#109 Try updating astra. |
|
Regarding the ASTRA issue: to me it looks like #684, which should have been fixed in #1376. But for the fix in #1376 to work, one need the ASTRA update that corresponds to the bug report astra-toolbox/astra-toolbox#109 that @adler-j mentioned. However, I think that this is still not in the latest release of ASTRA. |
|
You can check the ASTRA thing with a dev build from @wjp's conda channel: |
|
I'll try this later, since the pytest upgrade breaks my tests anyway. ;-) |
|
Regarding the pytest recursion issue: It does not occur in version 3.6.4 (I tested in on my workstation), which is the newest version available via conda previous to 3.7.1, and according to the link @kohr-h posted, it should be fixed in version 3.7.2. Edit: Corrected the version where the bug should be fixed. |
|
Thanks for investigating, @sbanert!
I guess you mean 3.7.2 (the upcoming release). In the meanwhile, we should blacklist the affected versions 3.7.0 and 3.7.1. Not that it's a terribly important core dependency, but the deps that are pulled in by default should work. There are two places where we have control over the version. Here: Line 147 in dd01b85 and here: Line 1 in dd01b85 For the syntax to blacklist versions you can check the However, I suggest that we deal with it in a separate issue, since I think we might want to bail out or at least warn in |
This prevents the pytest recursion bug described in pytest-dev/pytest#3771. For a discussion, see also pull request odlgroup#1416.
This prevents the pytest recursion bug described in pytest-dev/pytest#3771. For a discussion, see also pull request odlgroup#1416.
This prevents the pytest recursion bug described in pytest-dev/pytest#3771. For a discussion, see also pull request odlgroup#1416.
|
Regarding the ASTRA issue: After applying @kohr-h’s upgrade suggestion, the tests do no longer abort. However, I get 12 failures like this: __________________ test_projector[ geom='cone2d' - impl='astra_cpu' - angles='random' - in-place ] __________________
projector = RayTransform: uniform_discr([-20., -20.], [ 20., 20.], (100, 100), dtype='float32') -> DiscreteLp( FunctionSpace(I...
[-29.7, -29.1, -28.5, ..., 28.5, 29.1, 29.7]
),
rn((100, 100), dtype='float32', weighting=0.03651))
in_place = True
def test_projector(projector, in_place):
"""Test Ray transform forward projection."""
# TODO: this needs to be improved
# Accept 10% errors
places = 1
# Create Shepp-Logan phantom
vol = projector.domain.one()
# Calculate projection
if in_place:
proj = projector.range.zero()
projector(vol, out=proj)
else:
proj = projector(vol)
# We expect maximum value to be along diagonal
expected_max = projector.domain.partition.extent[0] * np.sqrt(2)
> assert almost_equal(proj.ufuncs.max(), expected_max, places=places)
E AssertionError: assert False
E + where False = almost_equal(33.670719, 56.568542494923804, places=1)
E + where 33.670719 = <bound method TensorSpaceUfuncs.max of <odl.util.ufuncs.TensorSpaceUfuncs object at 0x7ff89fbf6518>>()
E + where <bound method TensorSpaceUfuncs.max of <odl.util.ufuncs.TensorSpaceUfuncs object at 0x7ff89fbf6518>> = <odl.util.ufuncs.TensorSpaceUfuncs object at 0x7ff89fbf6518>.max
E + where <odl.util.ufuncs.TensorSpaceUfuncs object at 0x7ff89fbf6518> = DiscreteLp( FunctionSpace(IntervalProd([ 0.1072, -30. ], [ 6.1924, 30. ]), out_dtype='float32'),\n nonun... 11.42340565],\n [ 20.49684906, 24.00613022, 24.0072403 , ..., 14.09786224,\n 12.93267536, 11.78796482]]\n).ufuncs
odl/test/tomo/operators/ray_trafo_test.py:217: AssertionErrorMight these errors be due to #1415? |
Not exactly that, because we only look at forward projection here and check that we get close to the maximum ray intersection length. The change only affects the backprojection. |
This prevents the pytest recursion bug described in pytest-dev/pytest#3771. For a discussion, see also pull request odlgroup#1416.
d3ae519 to
b7dbd3a
Compare
kohr-h
left a comment
There was a problem hiding this comment.
Generally: Very useful, well-written and well-documented code.
But as an example, it misses its purpose by exposing overly many choices without a clear need, see the inline comments.
My suggestion is that you make this example more similar to the others by reducing the number of turning knobs, or by flagging this example as advanced and (perhaps) putting it into an advanced directory as per #1209
| # Some constants we shall use in the sequel to avoid magic numbers. | ||
| DIMENSIONS = 40.0 # Dimensions of reconstruction. | ||
| RESOLUTION = 1024 # Discretization level of reconstruction space. | ||
| SPLIT_METHOD = 'interlaced' # How to split the data? |
| EPSILON = 1e-6 # make sure that we do not divide by zero by adding this | ||
| STEPSIZE = 1. # Stepsize of the outer iterations | ||
| NITER = 5 # Number of iterations | ||
| RANDOM = True # Choose the oder of the inner iterations at random? |
There was a problem hiding this comment.
While I appreciate the flexibility that comes with this parametrization, I think that this level of flexibility tends to scare off people. There are just too many choices. The script misses its purpose as an example (aimed at simplicity and comprehensibility) and rather comes as a command line utility in disguise (aimed at flexibility and full-featuredness).
On top of that, not all of the parameters are really needed. Some are just "for implementer's convenience" and could be handled in a better way. Examples:
DIMENSIONS: Not needed. You can simply hardcode the extent of the reconstruction space and retrieve it later when you need it. That still gives you exactly one place to change it. (Apart from that, the name is slightly wrong:HALF_DIMENSIONSwould be more accurate.)RESOLUTION: Apart from the fact that most people would not think of resolution in this way, you don't actually need it as a module constant, just likeDIMENSIONS.
... and so on.
Ultimately, I think the only interesting knobs that I imagine somebody playing with this script would like to turn are related to the problem splitting, i.e., SPLIT_NUMBER AND SPLIT_METHOD, perhaps RANDOM. For the others: just do it as in the other examples and inline them.
|
|
||
| # Create the artificial data. | ||
| data_spaces = [op.range for op in ray_trafos] | ||
| data = [op(phantom) + NOISE_LEVEL * odl.phantom.white_noise(ds) |
There was a problem hiding this comment.
Here, for instance, there is no advantage of having NOISE_LEVEL as a constant above. The value is still entirely magical, and to change it, people will have to navigate to the top of the file.
| sampling_points1 = np.array(sampling_points1).T.tolist() | ||
| sampling_points2 = [list(p) for p in allpoints if p[j] % 2 == 1] | ||
| sampling_points2 = np.array(sampling_points2).T.tolist() | ||
| op1 = 2 * DIMENSIONS / RESOLUTION * odl.SamplingOperator( |
There was a problem hiding this comment.
2 * DIMENSIONS / RESOLUTION looks very much like reco_space.cell_sides[0] to me
| # data_fit_stepsizes = [1.0 / op.norm(estimate=True) ** 2 for op in ray_trafos] | ||
|
|
||
| # Now we build up the ingredients of our algorithm: | ||
| # Start at a black image, ... |
There was a problem hiding this comment.
Color only makes sense with a color map. Arrays are "zero", not "black".
|
Can I ignore the Python 2 failure? Otherwise, I think that I addressed all remarks. |
Yeah that's just some annoying conda issue, ignore it. I'll have a look at the code. |
kohr-h
left a comment
There was a problem hiding this comment.
A couple of nitpicks, but this looks good now.
| # --- Create simulated data (phantom) --- | ||
|
|
||
| # Reconstruction space: Set of two-dimensional quadratic images with | ||
| # given resolution. |
There was a problem hiding this comment.
This seems to be a comment based on the old version
| # Create the artificial data. | ||
| data_spaces = [op.range for op in ray_trafos] | ||
| data = [op(phantom) + 5.0 * odl.phantom.white_noise(ds) | ||
| for (op, ds) in zip(ray_trafos, data_spaces)] |
There was a problem hiding this comment.
It's fine like this, but if you want, you can use a fixed relative amount of noise (for some notion of "amount"), e.g.
noisefree_data = [op(phantom) for op in ray_trafos]
data = [proj + 0.05 * np.ptp(proj) * odl.phantom.white_noise(proj.space) for proj in noisefree_data](np.ptp is the one-function version of max - min that I didn't know until recently)
| # ... and for each dimension of the reconstruction space ... | ||
| reco_shape = reconstruction_space.shape | ||
| reco_dim = len(reco_shape) | ||
| for j in range(reco_dim): |
There was a problem hiding this comment.
j doesn't convey much, I'd prefer this:
for dim in range(reconstruction_space.ndim):
# Do stuff with `dim` instead of `j`We also use axis often in lieu of dim here.
| # respectively, in that dimension. | ||
| partial_der = odl.PartialDerivative( | ||
| reconstruction_space, j, pad_mode='order0') | ||
| allpoints = list(np.ndindex(reco_shape)) |
There was a problem hiding this comment.
Small thing, all_points reads nicer
| op1 = reconstruction_space.cell_sides[j] * odl.SamplingOperator( | ||
| reconstruction_space, sampling_points1) * partial_der | ||
| op2 = reconstruction_space.cell_sides[j] * odl.SamplingOperator( | ||
| reconstruction_space, sampling_points2) * partial_der |
There was a problem hiding this comment.
IMO the reconstruction_space name takes up a bit too much space, why not use the shorter reco_space as in the other examples?
| sampling_points1 = [list(p) for p in allpoints if p[j] % 2 == 0] | ||
| sampling_points1 = np.array(sampling_points1).T.tolist() | ||
| sampling_points2 = [list(p) for p in allpoints if p[j] % 2 == 1] | ||
| sampling_points2 = np.array(sampling_points2).T.tolist() |
There was a problem hiding this comment.
"1" and "2" don't convey meaning, use "even" and "odd"
Also, since the values are consumed immediately, shorter names are okay, e.g. odd_pts and even_pts.
| + \| \partial_{i, 2} x \|_1. | ||
|
|
||
| Here, :math:`\partial_{i, 1}` and :math:`\partial_{i, 2}` contain the even and | ||
| odd components, respectively, of the discretized :math:`i`th partial |
There was a problem hiding this comment.
Not sure, but
:math:`i`th
may not be rendered correctly -- the backtick must be followed by a whitespace or punctuation, as far as I know.
:math:`i`-th
should work, though.
| # In the stepsizes, we avoid the possible division by zero by adding a small | ||
| # positive value. The matrix corresponding to the operator `op` has only | ||
| # nonnegative entries, which makes this ensures that the final results are | ||
| # positve. |
There was a problem hiding this comment.
I guess you wanna delete "makes this"
| # Start at a zero image, ... | ||
| x = reconstruction_space.zero() | ||
|
|
||
| # collect all the functionals, ... |
There was a problem hiding this comment.
Above you also used ... at the front, which made the continuation more obvious.
|
Just a general comment w.r.t the git flow here. The commit messages here are typically of this type: "MAINT: Remove algorithm parameters and use keyword arguments." While this is locally informative (e.g. me following this PR, its obvious what you mean), being locally informative is not really what we're after down the line. I'd personally prefer if there was at least a hint as to where the change was done, e.g. "MAINT: Remove algorithm parameters and use kwargs in adupdates example." That way, one year from now when you're confused as to why things are as they are, it's a lot easier to scroll through the commits. Another option (actually the best option if you had infinite time available) is that you use these "local" commit messages during development, and then before the branch is merged, you rebase the branch and merge commits appropriately, also updating the commit messages. For example, this whole PR could quite frankly be squashed into a single commit with all the intermediate commit messages ("MAINT: Add relative instead of absolute noise.") removed. |
|
I can do, as you suggested, an interactive rebase before we merge. My reason is basically that it is easier to see what I changed after the reviews and that it’s easier to revert specific changes if someone disagrees with them. |
|
That's a very reasonable general workflow @sbanert. So please go ahead and take any final steps left to do, ping me, and I'll merge it. |
6b2c84c to
54da540
Compare
This pull request adds an example how to use the alternating dual updates method (#1243) and vector-valued stepsizes (#1166) in TV regularization for tomography.