Skip to content

ENH: Add a tomography example for the adupdates solver.#1416

Merged
kohr-h merged 1 commit into
odlgroup:masterfrom
sbanert:adupdates_example
Aug 29, 2018
Merged

ENH: Add a tomography example for the adupdates solver.#1416
kohr-h merged 1 commit into
odlgroup:masterfrom
sbanert:adupdates_example

Conversation

@sbanert

@sbanert sbanert commented Aug 17, 2018

Copy link
Copy Markdown
Contributor

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.

@pep8speaks

pep8speaks commented Aug 17, 2018

Copy link
Copy Markdown

Checking updated PR...

No PEP8 issues.

Comment last updated on August 29, 2018 at 11:28 Hours UTC

@sbanert

sbanert commented Aug 17, 2018

Copy link
Copy Markdown
Contributor Author

I don’t understand the Travis failure. :(

@kohr-h

kohr-h commented Aug 17, 2018

Copy link
Copy Markdown
Member

Seems to be this issue. Can you try older versions of pytest? Looks like we need to blacklist a few.

@sbanert

sbanert commented Aug 20, 2018

Copy link
Copy Markdown
Contributor Author

I can confirm that pytest 3.6.1 works on my laptop until it hits odl/test/tomo/backends/astra_cuda_test.py .[1] 5393 abort (core dumped) pytest --pep8 (my computer does not have an Nvidia GPU), whereas pytest 3.7.1 breaks with the recursion problem.

Edit: On my workstation, pytest 3.6.2 works fine (and it has got an Nvidia GPU), whereas pytest 3.7.1 breaks again.

@adler-j

adler-j commented Aug 20, 2018

Copy link
Copy Markdown
Member

Sounds very much like this issue:

astra-toolbox/astra-toolbox#109

Try updating astra.

@aringh

aringh commented Aug 20, 2018

Copy link
Copy Markdown
Member

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.

@kohr-h

kohr-h commented Aug 20, 2018

Copy link
Copy Markdown
Member

You can check the ASTRA thing with a dev build from @wjp's conda channel:

conda install -c wjpalenstijn astra-toolbox=1.9.0dev

@sbanert

sbanert commented Aug 20, 2018

Copy link
Copy Markdown
Contributor Author

I'll try this later, since the pytest upgrade breaks my tests anyway. ;-)

@sbanert

sbanert commented Aug 20, 2018

Copy link
Copy Markdown
Contributor Author

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.

@kohr-h

kohr-h commented Aug 20, 2018

Copy link
Copy Markdown
Member

Thanks for investigating, @sbanert!

should be fixed in version 3.7.1

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:

odl/setup.py

Line 147 in dd01b85

install_requires=[requires],

and here:

pytest >= 3.0.3

For the syntax to blacklist versions you can check the numpy line in the requirements.txt file.

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 odl.test() as well.

sbanert pushed a commit to sbanert/odl that referenced this pull request Aug 21, 2018
This prevents the pytest recursion bug described in
pytest-dev/pytest#3771. For a discussion, see also
pull request odlgroup#1416.
sbanert pushed a commit to sbanert/odl that referenced this pull request Aug 21, 2018
This prevents the pytest recursion bug described in
pytest-dev/pytest#3771. For a discussion, see also
pull request odlgroup#1416.
sbanert pushed a commit to sbanert/odl that referenced this pull request Aug 21, 2018
This prevents the pytest recursion bug described in
pytest-dev/pytest#3771. For a discussion, see also
pull request odlgroup#1416.
@sbanert

sbanert commented Aug 21, 2018

Copy link
Copy Markdown
Contributor Author

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: AssertionError

Might these errors be due to #1415?

@kohr-h

kohr-h commented Aug 21, 2018

Copy link
Copy Markdown
Member

Might 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.
However, that's not to say that no other scaling factors changed 😉

sbanert pushed a commit to sbanert/odl that referenced this pull request Aug 21, 2018
This prevents the pytest recursion bug described in
pytest-dev/pytest#3771. For a discussion, see also
pull request odlgroup#1416.
@sbanert sbanert force-pushed the adupdates_example branch from d3ae519 to b7dbd3a Compare August 21, 2018 16:51

@kohr-h kohr-h left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Choices?

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?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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_DIMENSIONS would 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 like DIMENSIONS.

... 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)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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(

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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, ...

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Color only makes sense with a color map. Arrays are "zero", not "black".

@sbanert

sbanert commented Aug 28, 2018

Copy link
Copy Markdown
Contributor Author

Can I ignore the Python 2 failure? Otherwise, I think that I addressed all remarks.

@kohr-h

kohr-h commented Aug 28, 2018

Copy link
Copy Markdown
Member

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 kohr-h left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A couple of nitpicks, but this looks good now.

# --- Create simulated data (phantom) ---

# Reconstruction space: Set of two-dimensional quadratic images with
# given resolution.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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)]

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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):

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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))

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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()

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"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

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess you wanna delete "makes this"

# Start at a zero image, ...
x = reconstruction_space.zero()

# collect all the functionals, ...

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Above you also used ... at the front, which made the continuation more obvious.

@adler-j

adler-j commented Aug 29, 2018

Copy link
Copy Markdown
Member

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.

@sbanert

sbanert commented Aug 29, 2018

Copy link
Copy Markdown
Contributor Author

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.

@kohr-h

kohr-h commented Aug 29, 2018

Copy link
Copy Markdown
Member

That's a very reasonable general workflow @sbanert.
In this case, though, the PR was already approved, so the only thing left to do was "whatever you want done before merge". When that's done I'd just merge it after having a last glance.

So please go ahead and take any final steps left to do, ping me, and I'll merge it.

@sbanert sbanert force-pushed the adupdates_example branch from 6b2c84c to 54da540 Compare August 29, 2018 11:28
@kohr-h kohr-h merged commit 43fefdd into odlgroup:master Aug 29, 2018
@sbanert sbanert deleted the adupdates_example branch August 29, 2018 11:57
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

5 participants