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

IO Slicing of LTIModel #1955

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

IO Slicing of LTIModel #1955

wants to merge 1 commit into from

Conversation

artpelling
Copy link
Member

@artpelling artpelling commented Feb 13, 2023

As a follow up to #1767 and #1916, I added a basic method to slice LTIModels based on numpy arrays. I would like to do this for the CSE tutorial, because it is a bit clunky in the notebook. I was not sure how to do it for arbitrary Operators. Please let me know, if you have an idea. Maybe we can use a RuleTable?

@artpelling artpelling requested a review from pmli February 13, 2023 19:20
@artpelling artpelling added the pr:new-feature Introduces a new feature label Feb 13, 2023
@artpelling artpelling added this to the 2023.1 milestone Feb 13, 2023
@@ -624,6 +625,40 @@ def to_abcde_files(self, files_basename):
if E is not None:
_mmwrite(Path(files_basename + '.E'), E)

def slice(self, input_indices=None, output_indices=None):
Copy link
Member

Choose a reason for hiding this comment

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

It would be nice to use the special method __getitem__.

Copy link
Member

Choose a reason for hiding this comment

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

I'd agree, in particular because it allows you to pass proper slices idiomatically. However, we should discuss, whether this operation is common enough and whether the user would expect Models to be sliceable and what the expected behavior is.
As inputs and outputs are well-defined for each Model (at least currently ...), we should also think about a generic implementation which does not assume knowledge of the internal structure.

Copy link
Member

Choose a reason for hiding this comment

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

I can't comment on the commonality of this operation (I never had a use for it myself).
Regarding other Models, for StationaryModels it would make sense to slice by the number of right-hand sides and outputs. For the general, nonlinear InstationaryModel

$$ \begin{align*} M(\mu) \dot{x}(t) & = f(x(t), u(t), \mu) \\ y(t) & = g(x(t), u(t), \mu) \end{align*} $$

the result should be an InstationaryModel

$$ \begin{align*} M(\mu) \dot{x}(t) & = f(x(t), S_1 u(t), \mu) \\ y(t) & = S_2^T g(x(t), u(t), \mu) \end{align*} $$

where $S_1$ and $S_2$ are some columns of the identity matrix. I'm not sure how to represent $S_1$ without defining a new Operator class for $f$ with a restricted input parameter. At least, $S_2$ should be easy when output functionals are separated, or it can be the ComponentProjectionOperator.

Comment on lines +658 to +660
return self.with_(B=NumpyMatrixOperator(to_matrix(self.B)[:, input_indices].reshape(-1, m)),
C=NumpyMatrixOperator(to_matrix(self.C)[output_indices, :].reshape(p, -1)),
D=NumpyMatrixOperator(to_matrix(self.D)[output_indices, input_indices].reshape(p, m)))
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, to allow parametric B, C, and D, we would need to use ConcatenationOperators, e.g.,

self.B @ NumpyMatrixOperator(np.eye(self.dim_input)[:, input_indices])

Additionally, pymor.algorithms.simplify.contract could be called to do the multiplication in the case of non-parametric Operators.

Copy link
Member

Choose a reason for hiding this comment

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

I'd agree with @pmli that this should be the fallback. This not only affects the parametric models but also models coming from an external pde solver. (It would be slightly better to create a VectorArrayOperator instead of a NumpyMatrixOperator as contract will convert it anyway.) One could still have this as a special case when B, C are NumpyMatrixOperators. But maybe it's not worth performance-wise to have extra code.

@pmli
Copy link
Member

pmli commented Feb 13, 2023

I would like to do this for the CSE tutorial

Note that we usually don't backport new features to old releases. You could install pyMOR from a git branch, but I feel that it would be good to only use an official release in the tutorial.

@sdrave sdrave force-pushed the main branch 4 times, most recently from 18e66f2 to e80291b Compare March 7, 2023 20:03
@sdrave sdrave removed the autoupdate label May 12, 2023
@sdrave
Copy link
Member

sdrave commented Jun 7, 2023

@artpelling, June 22 is hard freeze for the next release. Do you want to finish this by then or should we remove it from the milestone?

@artpelling
Copy link
Member Author

@artpelling, June 22 is hard freeze for the next release. Do you want to finish this by then or should we remove it from the milestone?

I thought about this one. There are too many aspects that still need conceptual decisions to be made. So we should move it.

@artpelling artpelling modified the milestones: 2023.1, 2023.2 Jun 7, 2023
@pmli pmli modified the milestones: 2023.2, 2024.1 Nov 2, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
pr:new-feature Introduces a new feature
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants