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

Min-theta approach #861

Merged
merged 8 commits into from
Feb 17, 2020
Merged

Min-theta approach #861

merged 8 commits into from
Feb 17, 2020

Conversation

ftalbrecht
Copy link
Contributor

Adds support for a variant of Bernards {min,max}-theta approach. However, we need an energy product which is not simply fom.operator.assemble(mu=mu_star), since we also need to clear columns.

This is of interest to @TiKeil

Copy link
Member

@sdrave sdrave left a comment

Choose a reason for hiding this comment

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

I think this PR should be split into two PRs. One for the discretizer and one for the functionals. Also, tests should be added.

If not `None`, a |Parameter| `mu` for which to assemble the symmetric part of the |Operator| of the resulting
|Model| `fom`, such that `fom.operator.assemble(mu) == fom.products['energy']` (apart from the fact that the
latter includes Dirichlet unit rows and columns, while the former only includes Dirichlet rows), if
`fom.operator` is already symmetric. Only available if `preassemble == False`.
Copy link
Member

Choose a reason for hiding this comment

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

I guess you mean 'if preassemble is True'. You can use FixedParameterOperator to make this work also when preassemble is False.

if energy_product:
eLi += [DiffusionOperator(grid, boundary_info, diffusion_function=df, dirichlet_clear_diag=True,
dirichlet_clear_columns=True)
for i, df in enumerate(p.diffusion.functions)]
Copy link
Member

Choose a reason for hiding this comment

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

Since you only want to assemble for a fixed parameter, it does not make sense to iterate over the coefficients. Just instantiate DiffusionOperator with p.diffusion.

if energy_product:
eLi += [AdvectionOperator(grid, boundary_info, advection_function=af, dirichlet_clear_diag=True,
dirichlet_clear_columns=True)
for i, af in enumerate(p.advection.functions)]
Copy link
Member

Choose a reason for hiding this comment

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

I don't think it makes much sense to have the asymmetrical advection part in the energy product.

coefficients += list(p.advection.coefficients)
elif p.advection is not None:
Li += [AdvectionOperator(grid, boundary_info, advection_function=p.advection,
dirichlet_clear_diag=True, name='advection')]
if energy_product:
eLi += [AdvectionOperator(grid, boundary_info, advection_function=p.advection,
dirichlet_clear_diag=True, dirichlet_clear_columns=True)]
Copy link
Member

Choose a reason for hiding this comment

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

See above.

if energy_product:
eLi += [ReactionOperator(grid, boundary_info, coefficient_function=rf, dirichlet_clear_diag=True,
dirichlet_clear_columns=True)
for i, rf in enumerate(p.reaction.functions)]
Copy link
Member

Choose a reason for hiding this comment

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

As above, we do not need the affine decomposition here.

src/pymor/parameters/functionals.py Show resolved Hide resolved
thetas = [ConstantParameterFunctional(theta) if not isinstance(theta, ParameterFunctionalInterface) else theta
for theta in thetas]
assert all([isinstance(f, ParameterFunctionalInterface) for f in thetas])
self.build_parameter_type(*chain(thetas))
Copy link
Member

Choose a reason for hiding this comment

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

Why is chain used here? Since it only gets one argument, it should not do anything.

assert all([isinstance(theta, (Number, ParameterFunctionalInterface)) for theta in thetas])
thetas = [ConstantParameterFunctional(theta) if not isinstance(theta, ParameterFunctionalInterface) else theta
for theta in thetas]
assert all([isinstance(f, ParameterFunctionalInterface) for f in thetas])
Copy link
Member

Choose a reason for hiding this comment

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

After the line above this assertion does not make much sense. Probably you want to assert before that you only have Numbers and ParameterFunctionals.

thetas = [ConstantParameterFunctional(f) if not isinstance(f, ParameterFunctionalInterface) else f
for f in thetas]
assert all([isinstance(f, ParameterFunctionalInterface) for f in thetas])
self.build_parameter_type(*chain(thetas))
Copy link
Member

Choose a reason for hiding this comment

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

See above.

assert isinstance(thetas, (list, tuple))
thetas = [ConstantParameterFunctional(f) if not isinstance(f, ParameterFunctionalInterface) else f
for f in thetas]
assert all([isinstance(f, ParameterFunctionalInterface) for f in thetas])
Copy link
Member

Choose a reason for hiding this comment

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

See above.

@ftalbrecht ftalbrecht marked this pull request as ready for review February 12, 2020 13:50
@ftalbrecht
Copy link
Contributor Author

Split energy product to #875, added docs, added simple tests...

@TiKeil
Copy link

TiKeil commented Feb 12, 2020

I would like to have a base functional as a generalization for the max and min theta approaches such that I can choose theta_q in the numerator.
This is for example required if I want to have an estimator for the partial derivatives w.r.t to parameters.
Would that be possible? I can also do that later in a new pull request ofc.

@ftalbrecht
Copy link
Contributor Author

I would like to have a base functional as a generalization for the max and min theta approaches such that I can choose theta_q in the numerator.

Not sure what you mean by that. Note, however, that the max functional evaluates differently than the min functional, since the former also allows negative thetas...

@TiKeil
Copy link

TiKeil commented Feb 12, 2020

I would like to have a base functional as a generalization for the max and min theta approaches such that I can choose theta_q in the numerator.

Not sure what you mean by that. Note, however, that the max functional evaluates differently than the min functional, since the former also allows negative thetas...

I am not sure about the min-theta but at least the max theta approach can be generalized to any bilinear or linear functional which has the same non parametric parts (which you denoted by a_q or L_q ). Lets suppose I want to estimate a function which can be written as

b(u,v, \mu) = sum_{q} anotherTheta_q(\mu) a_q(u, v).

Then I can at least estimate the continuity constant of b w.r.t the mu_bar energy norm by

max_{mu} ( | anotherTheta_q ( mu) | / theta_q(mu_bar) ) * 1

Thus, it would be nice to have a class where I can have theta_q and anotherTheta_q .
Or am I wrong here?

@@ -3,6 +3,7 @@
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)

from numbers import Number
from itertools import chain
Copy link
Member

Choose a reason for hiding this comment

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

Unused?

"""|ParameterFunctional| implementing the min-theta approach from [Haa17]_ (Proposition 2.35).

Let V denote a Hilbert space and let a: V x V -> K denote a parametric coercive bilinear form with affine
decomposition
Copy link
Member

Choose a reason for hiding this comment

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

:: ?


a(u, v, mu) = sum_{q = 1}^Q theta_q(mu) a_q(u, v),

for Q positive coefficient |ParameterFunctional|s theta_1, ..., theta_Q and positive semi-definit component
Copy link
Member

Choose a reason for hiding this comment

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

definite


for Q positive coefficient |ParameterFunctional|s theta_1, ..., theta_Q and positive semi-definit component
bilinear forms a_1, ..., a_Q: V x V -> K. Let mu_bar be a parameter with respect to which the coercivity constant
of a(., ., mu_bar) is known, i.e. we known alpha_mu_bar > 0, s.t.
Copy link
Member

Choose a reason for hiding this comment

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

:: ?

alpha_mu_bar |u|_V^2 <= a(u, u, mu=mu_bar).

The min-theta approach from [Haa17]_ (Proposition 2.35) allows to obtain a computable bound for the coercivity
constant of a(., ., mu) for arbitrary parameters mu, since
Copy link
Member

Choose a reason for hiding this comment

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

:: ?


if all theta_q(mu_bar) != 0.

Given a list of the thetas, the parameter mu_bar and the constant gamma_mu_bar, this functional thus evaluates
Copy link
Member

Choose a reason for hiding this comment

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

|Parameter|

assert len(thetas) > 0
assert all([isinstance(theta, (Number, ParameterFunctional)) for theta in thetas])
thetas = [ConstantParameterFunctional(f) if not isinstance(f, ParameterFunctional) else f
for f in thetas]
Copy link
Member

Choose a reason for hiding this comment

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

tuple?

self.build_parameter_type(*[t for t in thetas if t.parametric])
mu_bar = self.parse_parameter(mu_bar)
thetas_mu_bar = np.array([theta(mu_bar) for theta in thetas])
assert np.all(np.logical_or(thetas_mu_bar < 0, thetas_mu_bar > 0))
Copy link
Member

Choose a reason for hiding this comment

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

thetas_mu_bar != 0 ?

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 used this due to float cmp problems, using float_cmp now instead.

assert isinstance(gamma_mu_bar, Number)
assert gamma_mu_bar > 0
self.__auto_init(locals())
self.thetas_mu_bar = thetas_mu_bar # why is this required after __auto_init?
Copy link
Member

Choose a reason for hiding this comment

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

See above.





Copy link
Member

Choose a reason for hiding this comment

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

Remove whitespace.

@codecov
Copy link

codecov bot commented Feb 12, 2020

Codecov Report

Merging #861 into master will increase coverage by 0.06%.
The diff coverage is 95.34%.

Impacted Files Coverage Δ
src/pymor/parameters/functionals.py 92.97% <95.34%> (+0.71%) ⬆️

@ftalbrecht
Copy link
Contributor Author

Tried to address your points...

@ftalbrecht ftalbrecht merged commit 5a21744 into master Feb 17, 2020
@ftalbrecht ftalbrecht deleted the min-theta-approach branch February 17, 2020 08:29
@sdrave sdrave added this to the 2020.1 milestone Jul 21, 2020
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.

4 participants