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

examples: invoke tti example with --constant argument #1914

Merged
merged 16 commits into from
Aug 4, 2022

Conversation

ofmla
Copy link
Contributor

@ofmla ofmla commented Apr 29, 2022

I just realized that python ~/devito/examples/seimic/tti/tti_example.py --constant does not use the constant-tti model. This PR adds a very simple solution which will work for the constant-tti and layers-tti models

Copy link
Contributor

@FabioLuporini FabioLuporini left a comment

Choose a reason for hiding this comment

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

Thanks a lot! I just left a few comments

# With a new m as Constant
v0 = Constant(name="v", value=2.0, dtype=np.float32)
solver.forward(save=save, vp=v0)
solver.model.update('vp', v0.value)
Copy link
Contributor

Choose a reason for hiding this comment

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

this shouldn't be required strictly speaking?

Copy link
Contributor

Choose a reason for hiding this comment

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

I prefer the runtime argument vp=v0 to avoid keeping track of what was changed in the model

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Well, v0 is a constant and passing it directly to solver.forward won't work. It gives the following error:

Traceback (most recent call last):
  File "devito/examples/seismic/tti/tti_example.py", line 103, in <module>
    run(shape=shape, spacing=spacing, nbl=args.nbl, tn=tn,
  File "devito/examples/seismic/tti/tti_example.py", line 48, in run
    solver.forward(vp=2., save=save)
  File "/home/oscar/devito/examples/seismic/tti/wavesolver.py", line 155, in forward
    summary = self.op_fwd(save).apply(src=src, rec=rec, u=u, v=v,
  File "/home/oscar/devito/devito/operator/operator.py", line 708, in apply
    args = self.arguments(**kwargs)
  File "/home/oscar/devito/devito/operator/operator.py", line 557, in arguments
    args = self._prepare_arguments(**kwargs)
  File "/home/oscar/devito/devito/operator/operator.py", line 449, in _prepare_arguments
    args.update(p._arg_values(**kwargs))
  File "/home/oscar/devito/devito/types/dense.py", line 843, in _arg_values
    for i, s in zip(self.dimensions, new.shape):
AttributeError: 'float' object has no attribute 'shape'

and for that reason I used model.update

# With a new vp as a scalar value
solver.forward(save=save, vp=2.0)
solver.model.update('vp', 2.0)
Copy link
Contributor

Choose a reason for hiding this comment

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

perhaps the update should be performed internally to forward?

if you had multiple calls to forward and each time you had to update one or more parameters "by hand", it would be painful

Copy link
Contributor

Choose a reason for hiding this comment

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

Should probably yes be internal so the dt directly computed from runtime vo

Copy link
Contributor Author

Choose a reason for hiding this comment

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

So, the same applies here: passing an scalar to solver.forward() would lead to an error

Traceback (most recent call last):
  File "devito/examples/seismic/tti/tti_example.py", line 103, in <module>
    run(shape=shape, spacing=spacing, nbl=args.nbl, tn=tn,
  File "devito/examples/seismic/tti/tti_example.py", line 48, in run
    solver.forward(vp=2., save=save)
  File "/home/oscar/devito/examples/seismic/tti/wavesolver.py", line 155, in forward
    summary = self.op_fwd(save).apply(src=src, rec=rec, u=u, v=v,
  File "/home/oscar/devito/devito/operator/operator.py", line 708, in apply
    args = self.arguments(**kwargs)
  File "/home/oscar/devito/devito/operator/operator.py", line 557, in arguments
    args = self._prepare_arguments(**kwargs)
  File "/home/oscar/devito/devito/operator/operator.py", line 449, in _prepare_arguments
    args.update(p._arg_values(**kwargs))
  File "/home/oscar/devito/devito/types/dense.py", line 843, in _arg_values
    for i, s in zip(self.dimensions, new.shape):
AttributeError: 'float' object has no attribute 'shape'

If you planned for this to be possible, then it is likely that there is some error somewhere else. If this is not the case, you need to use model.update before

Copy link
Contributor

Choose a reason for hiding this comment

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

This looks like it ends up in the wrong spot (dense) instead of constant.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

So, are you saying you planned that the user could set a physical parameter i.e. velocity as a Constant or scalar and pass it to the forward method of the solver, and it should work smoothly?

if args.constant:
preset = 'constant-tti'
elif args.azi:
preset = 'layers-tti-noazimuth'
Copy link
Contributor

Choose a reason for hiding this comment

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

else:
    preset = 'layers-tii'

Copy link
Contributor

Choose a reason for hiding this comment

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

Don't particularly like these additional presets. There is also no reason constant would not be compatible with no azimuth

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I wrote it in that way because there is not a'constant-noazimuth'preset model. In any case you would have to play with conditionals I think

@@ -375,9 +375,9 @@ def update(self, name, value):
param = getattr(self, name)
except AttributeError:
# No physical parameter with tha name, create it
setattr(self, name, self._gen_phys_param(name, value, self.space_order))
setattr(self, name, self._gen_phys_param(value, name, self.space_order))
Copy link
Contributor

Choose a reason for hiding this comment

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

We need to put this under CI somehow

Copy link
Contributor Author

Choose a reason for hiding this comment

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

All that comes into my mind is

from examples.seismic import demo_model
import numpy as np


def test_model_updt(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0),
                    space_order=4, nbl=10):

    # Layered model (tti)
    tti_model = demo_model('layers-tti', shape=shape, spacing=spacing,
                           space_order=space_order, nbl=nbl)
    tpl1_set = set(tti_model.physical_parameters)

    # Layered model (isotropic)
    iso_model = demo_model('layers-isotropic', shape=shape, spacing=spacing,
                           space_order=space_order, nbl=nbl)
    tpl2_set = set(iso_model.physical_parameters)

    # Physical parameters in either set but not in the intersection.
    diff_phys_par = tuple(tpl1_set ^ tpl2_set)

    # Convert iso model in tti model
    for i in diff_phys_par:
        iso_model.update(i, getattr(tti_model, i).data[nbl:-nbl, nbl:-nbl, nbl:-nbl])
        assert np.array_equal(getattr(iso_model, i).data, getattr(tti_model, i).data)

I can add this script to the tests folder if you agree

Copy link
Contributor

@mloubout mloubout left a comment

Choose a reason for hiding this comment

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

I overall don't have an issue with this as it adds test/example cases but I think it should be done as a proper constant tti as currently it just slightly modifies the current VP constant only

@ofmla
Copy link
Contributor Author

ofmla commented May 2, 2022

I overall don't have an issue with this as it adds test/example cases but I think it should be done as a proper constant tti as currently it just slightly modifies the current VP constant only

I tried to fix a mistake (constant case was not running) without changing the test design. The test only update the vp parameter, so I left the things as they were.

@FabioLuporini FabioLuporini added the examples examples label May 24, 2022
Copy link
Contributor

@FabioLuporini FabioLuporini left a comment

Choose a reason for hiding this comment

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

What's still missing:

  • CI test for the constant case
  • (maybe) moving the "model.update" inside solver.forward

preset = 'layers-tti-noazimuth' if args.azi else 'layers-tti'
# if user defines constant and noazimuth arguments at the same time raise exception
if args.constant and args.azi:
parser.error('argument --noazimuth: not allowed with argument --constant')
Copy link
Contributor

Choose a reason for hiding this comment

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

Does this result in an actual Python exception?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

That is what the documentation says:

ArgumentParser.error(message)
This method prints a usage message including the message to the standard error and terminates the program with a status code of 2.

@@ -81,8 +83,15 @@ def test_tti_stability(shape, kernel):
help="Choice of finite-difference kernel")
args = parser.parse_args()

# Switch to TTI kernel if input is acoustic kernel
preset = 'layers-tti-noazimuth' if args.azi else 'layers-tti'
# if user defines constant and noazimuth arguments at the same time raise exception
Copy link
Contributor

Choose a reason for hiding this comment

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

If user (always start comments with capital letter)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Is it worth spending a commit here? 😅

@ofmla
Copy link
Contributor Author

ofmla commented Jun 14, 2022

What's still missing:

  • CI test for the constant case
  • (maybe) moving the "model.update" inside solver.forward

Well, the constant case is already tested in tests/test_tti.py, I mean the forward modeling for TTI case is executed using a constant model. It seems to me that there is not any wave-equation_example.py (i.e. tti_example.py, acoustic_example.py) in CI and that's OK. The point of this PR is that the tti_example.py script uses a layered model by default and although the user passes the argument --constant the exemple won't run on a constant model. The PR changes this.

On the other hand the solver.forward gets the physical parameters from model or kwargs. If I pass vp as argument to solver.forward it is passed to self.op_fwd(save).apply(...,**kwargs) even updating model via model.update within solver.forward. As result I get the same errors I mentioned in the coments above. That is why I use model.update before. If vp is a numpy array there is no problem passing it to solver.forward but if it is a scalar or Constant breakage will occur.

@@ -81,8 +74,15 @@ def test_tti_stability(shape, kernel):
help="Choice of finite-difference kernel")
args = parser.parse_args()

# Switch to TTI kernel if input is acoustic kernel
preset = 'layers-tti-noazimuth' if args.azi else 'layers-tti'
preset = 'layers-tti'
Copy link
Contributor

Choose a reason for hiding this comment

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

I would maybe rather move this inside the nested else

if ...
    ...
else:
    if args.azi:
        ....
    else:
        preset = 'layers-tti'

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Ok

@@ -0,0 +1,26 @@
from examples.seismic import demo_model
Copy link
Contributor

Choose a reason for hiding this comment

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

We shouldn't test examples inside tests because in theory the examples folder shouldn't (couldn't) even be accessible.

So I've been thinking about this...

How about we add it somewhere here? https://github.com/devitocodes/devito/blob/master/examples/seismic/tti/tti_example.py#L72

It should then be captured by CI through this line?
https://github.com/devitocodes/devito/blob/master/.github/workflows/examples.yml#L49
(@mloubout )

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I moved the model update test to tti_exemple.py as suggested

@@ -37,13 +37,6 @@ def run(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=250.0,
# Define receiver geometry (spread across `x, y`` just below surface)
rec, u, v, summary = solver.forward(save=save, autotune=autotune)

if preset == 'constant':
Copy link
Contributor

Choose a reason for hiding this comment

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

sorry I'm confused now... why is this being entirely dropped?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

As I see it, those lines are not longer necessary

Copy link
Contributor

@FabioLuporini FabioLuporini left a comment

Choose a reason for hiding this comment

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

LGTM. Thanks!

Copy link
Contributor

@mloubout mloubout left a comment

Choose a reason for hiding this comment

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

Thanks for the patience and updates. Left a minor comment on the test.

I still have one main point of discussion. So the constant TTI case is not Constant right, it's Function with constant values.
The main idea behind the constant isotropic was to have the Constant case to see the perf difference when there is no variation. I don't think there is much interest in a constant model if using Function personally (@FabioLuporini feel free to shit me down) since it's just a normal one with different values. So I'd rather make this an actual true constant model or remove it if there is no need for it.

Would still keep the update to Model and test OC.

@@ -71,6 +64,29 @@ def test_tti_stability(shape, kernel):
assert np.isfinite(norm(rec))


@pytest.mark.parametrize('shape', [(51, 51), (16, 16, 16)])
def test_model_updt(shape):
Copy link
Contributor

Choose a reason for hiding this comment

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

Some nitpicking:

  • Can this be moved in model.py
  • Can we test all cases, i.e 1) replacing the value of existing param by an array, 2) replacing existing param with a new Function 3) new param array, 4) new param Function

Copy link
Contributor Author

Choose a reason for hiding this comment

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

done

@@ -81,8 +97,16 @@ def test_tti_stability(shape, kernel):
help="Choice of finite-difference kernel")
args = parser.parse_args()

# Switch to TTI kernel if input is acoustic kernel
preset = 'layers-tti-noazimuth' if args.azi else 'layers-tti'
if args.constant:
Copy link
Contributor

Choose a reason for hiding this comment

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

Can probably make it a dict, but fine like that too.

@ofmla
Copy link
Contributor Author

ofmla commented Jun 22, 2022

Thanks for the patience and updates. Left a minor comment on the test.

I still have one main point of discussion. So the constant TTI case is not Constant right, it's Function with constant values. The main idea behind the constant isotropic was to have the Constant case to see the perf difference when there is no variation. I don't think there is much interest in a constant model if using Function personally (@FabioLuporini feel free to shit me down) since it's just a normal one with different values. So I'd rather make this an actual true constant model or remove it if there is no need for it.

Would still keep the update to Model and test OC.

You are right. I reversed the changes made and kept the model updating, thx

@mloubout mloubout force-pushed the tti_example_constant branch from a19f6eb to b2ea438 Compare July 17, 2022 01:22
@codecov
Copy link

codecov bot commented Jul 17, 2022

Codecov Report

Merging #1914 (9f7803c) into master (c860a49) will decrease coverage by 0.07%.
The diff coverage is 22.64%.

@@            Coverage Diff             @@
##           master    #1914      +/-   ##
==========================================
- Coverage   87.89%   87.81%   -0.08%     
==========================================
  Files         214      214              
  Lines       36251    36286      +35     
  Branches     5468     5477       +9     
==========================================
+ Hits        31864    31866       +2     
- Misses       3871     3903      +32     
- Partials      516      517       +1     
Impacted Files Coverage Δ
examples/seismic/acoustic/acoustic_example.py 50.00% <0.00%> (ø)
examples/seismic/tti/tti_example.py 33.87% <5.00%> (-9.88%) ⬇️
examples/seismic/model.py 83.33% <12.50%> (-8.13%) ⬇️
examples/seismic/preset_models.py 64.11% <100.00%> (-0.42%) ⬇️
tests/test_tti.py 100.00% <100.00%> (ø)
devito/ir/support/basic.py 92.40% <0.00%> (+0.20%) ⬆️

Help us with your feedback. Take ten seconds to tell us how you rate us.

@mloubout
Copy link
Contributor

Looks like this broke some notebooks and flake8 not passing.

@review-notebook-app
Copy link

Check out this pull request on  ReviewNB

See visual diffs & provide feedback on Jupyter Notebooks.


Powered by ReviewNB

@ofmla
Copy link
Contributor Author

ofmla commented Jul 17, 2022

Looks like this broke some notebooks and flake8 not passing.

I made some changes. I hope that it will work.

@mloubout mloubout force-pushed the tti_example_constant branch from c8108bb to 87e7d99 Compare July 18, 2022 00:14
@@ -39,7 +39,7 @@ def run(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=1000.0,
# Define receiver geometry (spread across x, just below surface)
rec, u, summary = solver.forward(save=save, autotune=autotune)

if preset == 'constant':
if preset == 'constant-isotropic':
# With a new m as Constant
Copy link
Contributor

Choose a reason for hiding this comment

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

double gap if you could fix

Copy link
Contributor

Choose a reason for hiding this comment

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

I will merge it as is because it's been around for too long, but we should check if there's a pep8 rule for flake8 to intercept it

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I don't see any problem here

Copy link
Contributor

Choose a reason for hiding this comment

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

He means there is an extra space between With and a

Copy link
Contributor Author

Choose a reason for hiding this comment

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

ok

@@ -81,8 +92,16 @@ def test_tti_stability(shape, kernel):
help="Choice of finite-difference kernel")
args = parser.parse_args()

# Switch to TTI kernel if input is acoustic kernel
preset = 'layers-tti-noazimuth' if args.azi else 'layers-tti'
if args.constant:
Copy link
Contributor

Choose a reason for hiding this comment

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

redundant if-else?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Well, it is never easy to satisfy all parties. I had written without else but one of you wanted an if-else. That is a matter of taste.

solver.forward(save=save, vp=v0)
# With a new vp as a scalar value
solver.forward(save=save, vp=2.0)
if preset in ['constant-tti', 'constant-tti-noazimuth']:
Copy link
Contributor

Choose a reason for hiding this comment

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

@ofmla wait... so this way it will run twice if preset=constant-tti for example ? right before this if and then inside?

Copy link
Contributor Author

@ofmla ofmla Jul 18, 2022

Choose a reason for hiding this comment

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

Yes, the intention is to run the forward modeling again passing scalars as Constant physical parameters

d = {'vp': vp, 'epsilon': epsilon, 'delta': delta, 'theta': theta, 'phi': phi}
for k, v in d.items():
v = Constant(name=k, value=v, dtype=np.float32)
solver.forward(save=save, vp=vp, epsilon=epsilon, delta=delta,
Copy link
Contributor

Choose a reason for hiding this comment

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

I don't think this does the Constant the k,v loop has no effect on this call

Copy link
Contributor

Choose a reason for hiding this comment

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

It should be something like

k, v in d.items():
            d[k] = Constant(name=k, value=v, dtype=np.float32)
solver.forward(save=save, **d)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

ok

@@ -39,7 +39,7 @@ def run(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=1000.0,
# Define receiver geometry (spread across x, just below surface)
rec, u, summary = solver.forward(save=save, autotune=autotune)

if preset == 'constant':
if preset == 'constant-isotropic':
# With a new m as Constant
Copy link
Contributor

Choose a reason for hiding this comment

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

He means there is an extra space between With and a

@mloubout
Copy link
Contributor

mloubout commented Aug 1, 2022

Needs flake8, rebase and GTG

@ofmla
Copy link
Contributor Author

ofmla commented Aug 1, 2022

Needs flake8, rebase and GTG

The flake8 errors are in files that I didn't touch. It is weird

@mloubout
Copy link
Contributor

mloubout commented Aug 1, 2022

The flake8 errors are in files that I didn't touch. It is weird

Must be the new standards, is fixed in #1979 so I'd say wait for it to be merged then rebase and should be all good

# With new physical parameters as Constants
d = {'vp': vp, 'epsilon': epsilon, 'delta': delta, 'theta': theta, 'phi': phi}
for k, v in d.items():
d[k] = Constant(name=k, value=v, dtype=np.float32)
Copy link
Contributor

Choose a reason for hiding this comment

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

shuold be a dict comprehension

@ofmla ofmla force-pushed the tti_example_constant branch from 2204804 to bfe49d0 Compare August 3, 2022 02:34
@mloubout mloubout force-pushed the tti_example_constant branch from bfe49d0 to 9f7803c Compare August 4, 2022 03:12
@mloubout mloubout merged commit db8202d into devitocodes:master Aug 4, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
examples examples
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants