Unverified Commit e5021923 authored by Stephan Rave's avatar Stephan Rave Committed by GitHub
Browse files

Merge pull request #1454 from pymor/output_functional_not_none

Ensure that output_functional is never None
parents 8f8be37d d66cc39e
Pipeline #106216 failed with stages
in 41 minutes and 42 seconds
......@@ -3,7 +3,7 @@
# License: BSD 2-Clause License (https://opensource.org/licenses/BSD-2-Clause)
from pymor.algorithms.preassemble import preassemble
from pymor.algorithms.rules import RuleTable, match_class
from pymor.algorithms.rules import RuleTable, RuleNotMatchingError, match_class, match_always
from pymor.models.interface import Model
from pymor.operators.constructions import (AdjointOperator, AffineOperator, ConcatenationOperator,
FixedParameterOperator, LincombOperator,
......@@ -14,7 +14,7 @@ from pymor.operators.numpy import NumpyMatrixOperator
from pymor.vectorarrays.list import NumpyListVectorSpace
def convert_to_numpy_list_vector_array(obj):
def convert_to_numpy_list_vector_array(obj, space=None):
"""Use NumpyListVectorArrayMatrixOperator instead of NumpyMatrixOperator.
This simple function recursively converts |NumpyMatrixOperators| to corresponding
......@@ -25,19 +25,25 @@ def convert_to_numpy_list_vector_array(obj):
obj
Either an |Operator|, e.g. |NumpyMatrixOperator| or a |LincombOperator| of
|NumpyMatrixOperators|, or an entire |Model| that is to be converted.
space
The |VectorSpace| that is to be converted.
Returns
-------
The converted |Operator| or |Model|.
"""
if isinstance(obj, Model) and space is None:
space = obj.solution_space
assert space is not None
obj = preassemble(obj)
return ConvertToNumpyListVectorArrayRules().apply(obj)
return ConvertToNumpyListVectorArrayRules(space).apply(obj)
class ConvertToNumpyListVectorArrayRules(RuleTable):
def __init__(self):
def __init__(self, space):
super().__init__(use_caching=True)
self.space = space
@match_class(AdjointOperator, AffineOperator, ConcatenationOperator,
FixedParameterOperator, LincombOperator, SelectionOperator,
......@@ -45,23 +51,20 @@ class ConvertToNumpyListVectorArrayRules(RuleTable):
def action_recurse(self, op):
return self.replace_children(op)
@match_always
def action_only_range(self, op):
if not (op.range == self.space and op.source != self.space):
raise RuleNotMatchingError
range = NumpyListVectorSpace(op.range.dim, op.range.id)
return VectorArrayOperator(range.from_numpy(op.as_range_array().to_numpy()), adjoint=False, name=op.name)
@match_always
def action_only_source(self, op):
if not (op.range != self.space and op.source == self.space):
raise RuleNotMatchingError
source = NumpyListVectorSpace(op.source.dim, op.source.id)
return VectorArrayOperator(source.from_numpy(op.as_source_array().to_numpy()), adjoint=True, name=op.name)
@match_class(NumpyMatrixOperator)
def action_NumpyMatrixOperator(self, op):
vector = op.source.dim == 1
functional = op.range.dim == 1
if vector and functional:
raise NotImplementedError
if vector:
space = NumpyListVectorSpace(op.range.dim, op.range.id)
return VectorOperator(space.from_numpy(op.matrix.reshape((1, -1))), op.name)
elif functional:
space = NumpyListVectorSpace(op.source.dim, op.source.id)
return VectorFunctional(space.from_numpy(op.matrix.ravel()), op.name)
else:
return op.with_(new_type=NumpyListVectorArrayMatrixOperator)
@match_class(VectorArrayOperator)
def action_VectorArrayOperator(self, op):
space = NumpyListVectorSpace(op.array.dim, op.array.space.id)
return op.with_(new_type=VectorArrayOperator,
array=space.from_numpy(op.array.to_numpy()))
return op.with_(new_type=NumpyListVectorArrayMatrixOperator)
......@@ -62,17 +62,17 @@ class StationaryModel(Model):
if isinstance(rhs, VectorArray):
assert rhs in operator.range
rhs = VectorOperator(rhs, name='rhs')
output_functional = output_functional or ZeroOperator(NumpyVectorSpace(0), operator.source)
assert rhs.range == operator.range and rhs.source.is_scalar and rhs.linear
assert output_functional is None or output_functional.source == operator.source
assert output_functional.source == operator.source
super().__init__(products=products, error_estimator=error_estimator, visualizer=visualizer, name=name)
self.__auto_init(locals())
self.solution_space = operator.source
self.linear = operator.linear and (output_functional is None or output_functional.linear)
if output_functional is not None:
self.dim_output = output_functional.range.dim
self.linear = operator.linear and output_functional.linear
self.dim_output = output_functional.range.dim
def __str__(self):
return (
......@@ -121,7 +121,6 @@ class StationaryModel(Model):
if not use_adjoint:
return super()._compute_output_d_mu(solution, mu, return_array)
else:
assert self.output_functional is not None
assert self.operator.linear
assert self.output_functional.linear
dual_solutions = self.operator.range.empty()
......@@ -204,12 +203,12 @@ class StationaryModel(Model):
+ ConstantOperator(affine_shift, self.solution_space))
new_rhs = self.rhs
if self.output_functional is not None:
if not isinstance(self.output_functional, ZeroOperator):
new_output_functional = \
self.output_functional @ (IdentityOperator(self.solution_space)
+ ConstantOperator(affine_shift, self.solution_space))
else:
new_output_functional = None
new_output_functional = self.output_functional
return self.with_(operator=new_operator, rhs=new_rhs, output_functional=new_output_functional,
error_estimator=None)
......@@ -284,13 +283,14 @@ class InstationaryModel(Model):
initial_data = VectorOperator(initial_data, name='initial_data')
mass = mass or IdentityOperator(operator.source)
rhs = rhs or ZeroOperator(operator.source, NumpyVectorSpace(1))
output_functional = output_functional or ZeroOperator(NumpyVectorSpace(0), operator.source)
assert isinstance(time_stepper, TimeStepper)
assert initial_data.source.is_scalar
assert operator.source == initial_data.range
assert rhs.linear and rhs.range == operator.range and rhs.source.is_scalar
assert mass.linear and mass.source == mass.range == operator.source
assert output_functional is None or output_functional.source == operator.source
assert output_functional.source == operator.source
super().__init__(products=products, error_estimator=error_estimator, visualizer=visualizer, name=name)
......@@ -298,8 +298,7 @@ class InstationaryModel(Model):
self.__auto_init(locals())
self.solution_space = operator.source
self.linear = operator.linear and (output_functional is None or output_functional.linear)
if output_functional is not None:
self.dim_output = output_functional.range.dim
self.dim_output = output_functional.range.dim
def __str__(self):
return (
......@@ -336,8 +335,6 @@ class InstationaryModel(Model):
None -> D
self.mass -> E
"""
if self.output_functional is None:
raise ValueError('No output defined.')
A = - self.operator
B = self.rhs
C = self.output_functional
......
......@@ -103,10 +103,7 @@ class Model(CacheableObject, ParametricObject):
|NumPy array| with the computed output or a dict which at least
must contain the key `'output'`.
"""
if not getattr(self, 'output_functional', None):
return np.zeros((len(solution), 0))
else:
return self.output_functional.apply(solution, mu=mu).to_numpy()
return self.output_functional.apply(solution, mu=mu).to_numpy()
def _compute_solution_d_mu_single_direction(self, parameter, index, solution, mu=None, **kwargs):
"""Compute the partial derivative of the solution w.r.t. a parameter index
......
......@@ -21,6 +21,7 @@ if config.HAVE_TORCH:
from pymor.core.base import BasicObject
from pymor.models.interface import Model
from pymor.operators.constructions import ZeroOperator
from pymor.vectorarrays.numpy import NumpyVectorSpace
class NeuralNetworkModel(Model):
......@@ -71,8 +72,9 @@ if config.HAVE_TORCH:
self.__auto_init(locals())
self.solution_space = NumpyVectorSpace(neural_network.output_dimension)
if output_functional is not None:
self.dim_output = output_functional.range.dim
output_functional = output_functional or ZeroOperator(NumpyVectorSpace(0), self.solution_space)
assert output_functional.source == self.solution_space
self.dim_output = output_functional.range.dim
def _compute_solution(self, mu=None, **kwargs):
......
......@@ -178,7 +178,7 @@ class StationaryRBReductor(ProjectionBasedReductor):
'operator': project(fom.operator, RB, RB),
'rhs': project(fom.rhs, RB, None),
'products': {k: project(v, RB, RB) for k, v in fom.products.items()},
'output_functional': project(fom.output_functional, None, RB) if fom.output_functional else None
'output_functional': project(fom.output_functional, None, RB)
}
return projected_operators
......@@ -189,8 +189,7 @@ class StationaryRBReductor(ProjectionBasedReductor):
'operator': project_to_subbasis(rom.operator, dim, dim),
'rhs': project_to_subbasis(rom.rhs, dim, None),
'products': {k: project_to_subbasis(v, dim, dim) for k, v in rom.products.items()},
'output_functional': (project_to_subbasis(rom.output_functional, None, dim)
if rom.output_functional else None)
'output_functional': project_to_subbasis(rom.output_functional, None, dim)
}
return projected_operators
......@@ -256,7 +255,7 @@ class InstationaryRBReductor(ProjectionBasedReductor):
'rhs': project(fom.rhs, RB, None),
'initial_data': projected_initial_data,
'products': {k: project(v, RB, RB) for k, v in fom.products.items()},
'output_functional': project(fom.output_functional, None, RB) if fom.output_functional else None
'output_functional': project(fom.output_functional, None, RB)
}
return projected_operators
......@@ -283,8 +282,7 @@ class InstationaryRBReductor(ProjectionBasedReductor):
'rhs': project_to_subbasis(rom.rhs, dim, None),
'initial_data': projected_initial_data,
'products': {k: project_to_subbasis(v, dim, dim) for k, v in rom.products.items()},
'output_functional': (project_to_subbasis(rom.output_functional, None, dim)
if rom.output_functional else None)
'output_functional': project_to_subbasis(rom.output_functional, None, dim)
}
return projected_operators
......
......@@ -248,8 +248,7 @@ if config.HAVE_TORCH:
def _build_rom(self):
"""Construct the reduced order model."""
with self.logger.block('Building ROM ...'):
projected_output_functional = (project(self.fom.output_functional, None, self.reduced_basis)
if self.fom.output_functional else None)
projected_output_functional = project(self.fom.output_functional, None, self.reduced_basis)
rom = NeuralNetworkModel(self.neural_network, parameters=self.fom.parameters,
output_functional=projected_output_functional,
name=f'{self.fom.name}_reduced')
......@@ -439,8 +438,7 @@ if config.HAVE_TORCH:
def _build_rom(self):
"""Construct the reduced order model."""
with self.logger.block('Building ROM ...'):
projected_output_functional = (project(self.fom.output_functional, None, self.reduced_basis)
if self.fom.output_functional else None)
projected_output_functional = project(self.fom.output_functional, None, self.reduced_basis)
rom = NeuralNetworkInstationaryModel(self.fom.T, self.nt, self.neural_network,
parameters=self.fom.parameters,
output_functional=projected_output_functional,
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment