Successive constraints method #1989
Conversation
|
Seems like I should complain more often ... Maybe we can discuss things offline next week? |
Yes of course, let's discuss this next week. |
HenKlei
force-pushed
the
successive-constraints-method
branch
from
55a155b
to
61d1f23
Compare
|
@sdrave Do you think you will have time to discuss and look into this before the release? |
|
|
HenKlei
force-pushed
the
successive-constraints-method
branch
from
61d1f23
to
1c0bcd7
Compare
HenKlei
force-pushed
the
successive-constraints-method
branch
from
15ac296
to
65ab174
Compare
|
It would be great if you could take a look at this @sdrave when you find time. |
|
The pipelines are failing due to a very old scipy version (1.5.4) where the 'highs' method (which is now the default) was not available as solver of a linear program in scipy. 'simplex' would be a method that is supported in probably all versions but is considered legacy in current versions. So maybe also not the best option to use that one as default. How about pinning scipy to >=1.7.0? |
Debian oldstable has 1.6.0. I think that should be enough as well? Any reasons to require 1.7? |
There was a problem hiding this comment.
Hi @HenKlei, I think the code generally looks nice.
My biggest concern is that you have to solve the same eigenproblems ('SM') for both the upper and lower bound functional (and the upper bound functional does not allow passing the solutions via __init__). I wonder whether it makes more sense to refactor the code into an algorithms.scm module, with a function which instantiates both functionals. That would be somewhat less elegant, but would avoid significant offline costs.
Further remarks:
- Have you considered also implementing the greedy algorithm from the paper? Should be fairly straightforward.
- I have never looked into this personally, but how does the situation look like for inf-sup stable problems? Can we easily support this as well?
- As already discussed, it would be good to extend the demo/add another demo to validate the performance of the method (only looking at the coercivity constants).
| Parameters | ||
| ---------- | ||
| operator | ||
| |LincombOperator| for which to provide a lower bound on the coercivity constant. |
There was a problem hiding this comment.
| |LincombOperator| for which to provide a lower bound on the coercivity constant. | |
| |LincombOperator| for which to provide a upper bound on the coercivity constant. |
| def __init__(self, operator, constraint_parameters, | ||
| method='highs', options={}, M=None, bounds=None, coercivity_constants=None): | ||
| from pymor.operators.constructions import LincombOperator | ||
| assert isinstance(operator, LincombOperator) |
There was a problem hiding this comment.
Also check that all of operator.operators are non-parametric and linear.
| for f in operator.coefficients) | ||
|
|
||
| if self.M is not None: | ||
| self.M = min(self.M, len(self.constraint_parameters)) |
There was a problem hiding this comment.
Maybe output a warning if M is clipped?
| |LincombOperator| for which to provide a lower bound on the coercivity constant. | ||
| constraint_parameters | ||
| List of |Parameters| used to construct the constraints. | ||
| method |
| If `None`, all parameters from `constraint_parameters` are used. | ||
| bounds | ||
| Either `None` or a list of tuples containing lower and upper bounds | ||
| for the design variables. |
There was a problem hiding this comment.
Explain what design variables means in this context.
| self.logger.info('Computing coercivity constants for parameters by solving eigenvalue problems ...') | ||
| self.coercivity_constants = [] | ||
| for mu in constraint_parameters: | ||
| fixed_parameter_op = FixedParameterOperator(operator, mu=mu) |
There was a problem hiding this comment.
Why not
| fixed_parameter_op = FixedParameterOperator(operator, mu=mu) | |
| fixed_parameter_op = operator.assemble(mu) |
| _, indices = self.kdtree.query(mu.to_numpy(), k=self.M) | ||
| selected_parameters = [self.constraint_parameters[i] for i in list(indices)] | ||
| else: | ||
| indices = list(range(len(self.constraint_parameters))) |
There was a problem hiding this comment.
| indices = list(range(len(self.constraint_parameters))) | |
| indices = np.arange(len(self.constraint_parameters)) |
?
|
|
||
| def __init__(self, operator, constraint_parameters): | ||
| from pymor.operators.constructions import FixedParameterOperator, LincombOperator | ||
| assert isinstance(operator, LincombOperator) |
| self.minimizers = [] | ||
| from pymor.algorithms.eigs import eigs | ||
| for mu in self.constraint_parameters: | ||
| fixed_parameter_op = FixedParameterOperator(operator, mu=mu) |
| # select reduction algorithm with error estimator | ||
| ################################################# | ||
| bounds = None | ||
| coercivity_constants = None |
You're right, 1.6.0 would be sufficient. |
Since @sdrave was complaining about the missing successive constraints method for a few times now, I started implementing something in that direction. However, this is far from being finished at the moment.