Blatt 08, 03.07.2019¶

Aktivieren Sie wie gewohnt ihre Arbeitsumgebung und starten Sie den Jupyter Notebook server, siehe zB Blatt 1, Aufgabe 0.
Importieren Sie numpy und pymor.basic und machen Sie matplotlib für das Notebook nutzbar.

Wir betrachten das Diffusionsproblem von Blatt 04, Aufgabe 1 mit dem Ausgabe Funktional

$$ s(v) := \frac{1}{|\Omega_\text{out}|} \int_{\Omega_\text{out}} v\; \text{d}x \quad\quad\text{mit } \Omega_\text{out} := (0, 1) \times (0.9, 1) $$

und das folgende Minimierungsproblem: finde das Minimum

$$ \mu* := \underset{\mu \in \mathcal{P}}{\text{argmin}} J(\mu) $$

wobei $J(\mu) := s(u_\mu)$, und $u_\mu$ für $\mu \in \mathcal{P}$ die Lösung des Diffusionsproblems sei; $J$ wird oft Paramter-zu-Ausgabe Abbilding oder input-output map genannt.

Aufgabe 1: der schwache greedy Algorithmus¶

Erstellen Sie ein analytisches Problem unter Verwedung von problem_B4_A1_parametric und fügen Sie dem Problem das Ausgabefunktional mit der id output hinzu (siehe auch API Dokumentation von StationaryProblem).

Diskretisieren Sie das Problem mit Hilfe des cg discretizers und werten Sie $J$ für den Paramter $\mu = ${'R': 1, 'B': 1} aus.
- extrahieren Sie dazu das Ausgabefunktional s aus dem operators-dict der Diskretisierung
- definieren Sie $J$

Erstellen Sie einen CoerciveRBReductor mit entsprechenden product und coercivity_estimator

Erstellen Sie eine reduzierte Basis mit Hilfe eines (weak discrete) greedy Algorithmus.
- Machen Sie sich mit dem greedy Algorithmus in pyMOR vertraut. Dieser Nutzt den Residuenbasierten offline/online zerlegbaren Fehlerschätzer (weak) des CoerciveRBReductor, um den Modellreduktionsfehler online-effizient für eine große Menge an Trainingsparametern (discrete) zu schätzen. Es werden also für jedes Element der reduzierten Basis nur eine hochdimensionale Lösung der PDE benötigt.
- Verwenden Sie eine Trainingsmenge von 10000 gleichverteilten Paramtern.
- Stellen Sie sicher, dass der maximale geschätzte Modellreduktionsfehler (über die Trainingsmenge) 1e-3 nicht überschreitet.
- Stellen Sie sicher, dass die reduzierte Basis mit Hilfe eines Gram-Schmidt Verfahrens orthonormalisiert wird, setzen Sie hierzu extension_params={'method': 'gram_schmidt'}.

Stellen Sie den Fehlerabfall während der Basisgenerierung mit Hilfe der vom greedy Algorithmus zurück gegebenen Daten grafisch dar.

Aufgabe 2: PDE-constrained minimization¶

Definieren Sie eine Funktione
```
def parameter_to_output(R, B):
    ...
```
welche mit Hilfe der reduzierten Diskretisierung und des reduzierten Ausgabefunktionals eine RB-Approximation von $J$ berechnet.
- Extrahieren Sie dazu die reduzierte Diskretisierung aus den vom greedy Algorithmus übergebenen Daten.

def plot_3d(x, y, f, xlabel='X', ylabel='Y'):
    # compute_value_matrix
    f_of_x = np.zeros((len(x), len(y)))
    import time
    t = time.time()
    for ii in range(len(x)):
        for jj in range(len(y)):
            f_of_x[jj][ii] = f(x[ii], y[jj]) # ii and jj are switched on purpose!
    x, y = np.meshgrid(x, y)
    print('{} evaluations took {}s'.format(len(x)*len(y), time.time() - t))
    # plot
    from mpl_toolkits.mplot3d import Axes3D # required for 3d plots
    from matplotlib import cm # required for colormaps
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.plot_surface(x, y, f_of_x, cmap=cm.viridis)
    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel)
    ax.invert_xaxis()
    return ax

def plot_point(ax, point, color):
    x, y = point[0], point[1]
    ax.plot([x, x], [y, y], ax.get_zlim(), color, zorder=10)

Erstellen Sie einen 3d Plot der RB-Approximation von $J$ mit Hilfe von plot_3d.
- Geben Sie als x-Werte 100 gleichverteilte Werte im Definitionsbereich der Parameterkomponente R an, entsprechend als y-Werte B. Nutzen Sie dazu den Parameterraum der reduzierten Diskretisierung.

ax = plot_3d(#x = 
             #Y = 
             parameter_to_output,
             xlabel='R', ylabel='B')
xlocs, xlabels = plt.xticks()
ylocs, ylabels = plt.yticks()

Machen Sie sich mit der minimize Funktion in scipy vertraut und vervollständigen Sie die Definition von find_minimum, um das Minimierungsproblem zu lösen.
- wählen Sie als Minimierer den BFGS Algorithmus für "Bound-Constrained minimization"
- lassen Sie den Gradienten von $J$ mit Hilfe von Finiten Differenzen Approximieren
- nutzen Sie den Paramterraum der reduzierten Diskretisierung um die "Bound-Constraints" zu definieren
- übergeben Sie options={'ftol': 1e-15}

def find_minimum(initial_guess):
    from scipy.optimize import minimize
    return minimize(#...
                   )

def print_report(result):
    if result.success:
        print('succeded after {} iterations (required {} PDE solves)'.format(result.nit, result.nfev))
        print('J({}) = {}'.format(result.x, result.fun))
    else:
        print('failed!')

Lassen Sie die Minimierung für vier verschiedene Anfangswerte mit Hilfe des folgenden Codes druchlaufen. Was können Sie beobachten?

for mu, color in zip(rd.parameter_space.sample_uniformly(2), ('red', 'yellow', 'orange', 'purple')):
    initial_guess = [mu['R'], mu['B']]    
    print('minimizing from {}...'.format(initial_guess))
    result = find_minimum(initial_guess)
    print_report(result)
    plot_point(ax, result.x, color)
    print('')

Erstellen Sie einen 3d plot der RB-Approximation von $J$ wie oben, wobei die $x$ und $y$ Punkte nicht gleichverteilt, sonder mit Hilfe des Logarithmus zur Basis 10 verteilt sind.

ax = plot_3d(# X = 
             # Y = 
             # function to plot
             , xlabel='R', ylabel='B')
xlocs, _ = plt.xticks()
ylocs, _ = plt.yticks()
plt.xticks(np.linspace(min(xlocs), max(xlocs), len(xlabels)), xlabels)
plt.yticks(np.linspace(min(ylocs), max(ylocs), len(ylabels)), ylabels)

Wiederholen Sie auch die Minimierung und grafische Darstellung der gefundenen Minima. Passen Sie letzteren entsprechend an.

for mu, color in zip(rd.parameter_space.sample_uniformly(2), ('red', 'yellow', 'orange', 'purple')):
    initial_guess = [mu['R'], mu['B']]    
    print('minimizing from {}...'.format(initial_guess))
    result = find_minimum(initial_guess)
    print_report(result)
    plot_point(ax,
               # point
               , color)
    print('')