From 69b9ff8a25627c1a12676d4770f45acbfa32c0d3 Mon Sep 17 00:00:00 2001
From: Andreas Buhr <andreas@andreasbuhr.de>
Date: Mon, 2 Mar 2015 16:32:08 +0100
Subject: [PATCH] [operators.constructions] add selection operator

---
 src/pymor/operators/constructions.py | 65 ++++++++++++++++++++++++++++
 src/pymortests/operators.py          | 35 +++++++++++++++
 2 files changed, 100 insertions(+)

diff --git a/src/pymor/operators/constructions.py b/src/pymor/operators/constructions.py
index d0eab0f5d0..8394f7ec9f 100644
--- a/src/pymor/operators/constructions.py
+++ b/src/pymor/operators/constructions.py
@@ -608,3 +608,68 @@ def invert_options(self):
 
     def apply_inverse(self, U, ind=None, mu=None, options=None):
         return self.operator.apply_inverse(U, ind=ind, mu=self.mu, options=options)
+
+
+class SelectionOperator(OperatorBase):
+    """An |Operator| selecting one out of a list of |Operators|.
+
+    operators[i] is used
+    if parameterfunctional(mu) is less or equal than boundaries[i]
+    and greater than boundaries[i-1]::
+
+        -infty ------- boundaries[i] ---------- boundaries[i+1] ------- infty
+                            |                        |
+        --- operators[i] ---|---- operators[i+1] ----|---- operators[i+2]
+                            |                        |
+
+    Parameters
+    ----------
+    operators
+        List of |Operators| from which one |Operator| is
+        selected based on a parameter.
+    parameterfunctional
+        A |ParameterFunctional| used for the selection of one |Operator|..
+    name
+        Name of the operator.
+
+    """
+    def __init__(self, operators, parameterfunctional, boundaries, name=None):
+        assert len(operators) > 0
+
+        assert len(boundaries) == len(operators) - 1
+        # check that boundaries are ascending:
+        for i in range(len(boundaries)-1):
+            assert boundaries[i] <= boundaries[i+1]
+
+        assert all(isinstance(op, OperatorInterface) for op in operators)
+
+        assert all(op.source == operators[0].source for op in operators[1:])
+        assert all(op.range == operators[0].range for op in operators[1:])
+        self.source = operators[0].source
+        self.range = operators[0].range
+        self.operators = operators
+        self.linear = all(op.linear for op in operators)
+
+        self.name = name
+        self.build_parameter_type(inherits=list(operators) + [parameterfunctional])
+        self._try_assemble = not self.parametric
+
+        self.boundaries = boundaries
+        self.parameterfunctional = parameterfunctional
+
+    def _get_operator_number(self, mu):
+        value = self.parameterfunctional.evaluate(mu)
+        for i in range(len(self.boundaries)):
+            if self.boundaries[i] >= value:
+                return i
+        return len(self.boundaries)
+
+    def apply(self, U, ind=None, mu=None):
+        mu = self.parse_parameter(mu)
+        operator_number = self._get_operator_number(mu)
+        return self.operators[operator_number].apply(U,ind=ind, mu=mu)
+
+    def as_vector(self, mu=None):
+        mu = self.parse_parameter(mu)
+        operator_number = self._get_operator_number(mu)
+        return self.operators[operator_number].as_vector(mu=mu)
diff --git a/src/pymortests/operators.py b/src/pymortests/operators.py
index ebdb20bd1c..5917af9a89 100644
--- a/src/pymortests/operators.py
+++ b/src/pymortests/operators.py
@@ -8,6 +8,9 @@
 import numpy as np
 import pytest
 
+from pymor.operators.constructions import SelectionOperator
+from pymor.parameters.base import ParameterType
+from pymor.parameters.functionals import GenericParameterFunctional
 from pymor.core.exceptions import InversionError
 from pymor.la.numpyvectorarray import NumpyVectorArray
 from pymor.tools.floatcmp import float_cmp_all
@@ -17,6 +20,37 @@
 from pymortests.pickle import assert_picklable, assert_picklable_without_dumps_function
 
 
+def test_selection_op():
+    p1 = MonomOperator(1)
+    select_rhs_functional = GenericParameterFunctional(
+        lambda x: round(float(x["nrrhs"])), 
+        ParameterType({"nrrhs" : tuple()})
+    )
+    s1 = SelectionOperator(
+        operators = [p1], 
+        boundaries = [], 
+        parameterfunctional = select_rhs_functional,
+        name = "foo"
+    )
+    x = np.linspace(-1., 1., num=3)
+    vx = NumpyVectorArray(x[:, np.newaxis])
+    assert np.allclose(p1.apply(vx,mu=0).data, s1.apply(vx,mu=0).data)
+
+    s2 = SelectionOperator(
+        operators = [p1,p1,p1,p1],
+        boundaries = [-3, 3, 7],
+        parameterfunctional = select_rhs_functional,
+        name = "Bar"
+    )
+
+    assert s2._get_operator_number({"nrrhs":-4}) == 0
+    assert s2._get_operator_number({"nrrhs":-3}) == 0
+    assert s2._get_operator_number({"nrrhs":-2}) == 1
+    assert s2._get_operator_number({"nrrhs":3}) == 1
+    assert s2._get_operator_number({"nrrhs":4}) == 2
+    assert s2._get_operator_number({"nrrhs":7}) == 2
+    assert s2._get_operator_number({"nrrhs":9}) == 3
+
 def test_lincomb_op():
     p1 = MonomOperator(1)
     p2 = MonomOperator(2)
@@ -250,3 +284,4 @@ def test_apply_inverse_wrong_ind(operator_with_arrays):
     for ind in invalid_inds(V):
         with pytest.raises(Exception):
             op.apply_inverse(V, mu=mu, ind=ind)
+