Institut für Numerische und Angewandte Mathematik
FB Mathematik und Informatik der Universität Münster
Prof. Dr. Mario Ohlberger und Dr. Stephan Rave

Übung zum Programmierpraktikum Numerische Mehrskalenmethoden und Modellreduktion

Sommersemester 2017

Musterlösung zu Blatt 3

In [1]:
%matplotlib notebook
from pymor.basic import *
set_log_levels({'pymor': 'WARN',
                'pymor.functions.bitmap': 'CRITICAL'})

Aufgabe 1 a)

Liegt die Datei problems.py:

In [2]:
%cat problems.py

from pymor.basic import *


def problem_1_1_a():
    return StationaryProblem(
        domain=RectDomain([[0., 0.], [1., 1.]]),

        diffusion=ConstantFunction(1, 2),

        rhs=ExpressionFunction('(np.sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 1.', 2, ()),

        name='Blatt1_Aufgabe_1a'
    )


def problem_1_1_b():
    return StationaryProblem(
        domain=RectDomain(bottom='neumann'),

        neumann_data=ConstantFunction(-1., 2),

        diffusion=ExpressionFunction('1. - (np.sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 0.999', 2, ()),

        name='Blatt1_Aufgabe_1b'
    )


def problem_1_1_c():
    return StationaryProblem(
        domain=RectDomain(bottom='neumann'),

        neumann_data=ConstantFunction(-1., 2),

        diffusion=ExpressionFunction(
            '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',
            2, (),
            values={'K': 10}
        ),

        name='Blatt1_Aufgabe_1c'
    )


def problem_1_1_d():
    return StationaryProblem(
        domain=RectDomain(bottom='neumann'),

        neumann_data=ConstantFunction(-1., 2),

        diffusion=BitmapFunction('RB.png', range=[0.001, 1]),

        name='Blatt1_Aufgabe_1d'
    )


def problem_1_2_a():
    return StationaryProblem(
        domain=RectDomain(bottom='neumann'),

        neumann_data=ExpressionFunction('-cos(pi*x[...,0])**2*neum', 2, (), parameter_type={'neum': ()}),

        diffusion=ExpressionFunction(
            '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',
            2, (),
            values={'K': 10}
        ),

        parameter_space=CubicParameterSpace({'neum': ()}, -10, 10),

        name='Blatt1_Aufgabe_2a'
    )


def problem_1_2_b():
    return StationaryProblem(
        domain=RectDomain(bottom='neumann'),

        neumann_data=ExpressionFunction('-cos(pi*x[...,0])**2*neum', 2, (), parameter_type={'neum': ()}),

        diffusion=ExpressionFunction(
            '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * (1 - diffu)',
            2, (),
            values={'K': 10},
            parameter_type={'diffu': ()}
        ),

        parameter_space=CubicParameterSpace(
            {'neum': (), 'diffu': ()},
            ranges={'diffu': [0.0001, 1.], 'neum': [-10, 10]}
        ),

        name='Blatt1_Aufgabe_2b'
    )


def problem_1_2_c():
    return StationaryProblem(
        domain=RectDomain(bottom='neumann'),

        neumann_data=ConstantFunction(-1., 2),

        diffusion=LincombFunction(
            [ConstantFunction(1., 2),
             BitmapFunction('R.png', range=[1, 0]),
             BitmapFunction('B.png', range=[1, 0])],
            [1.,
             ExpressionParameterFunctional('-(1 - R)', {'R': ()}),
             ExpressionParameterFunctional('-(1 - B)', {'B': ()})]
        ),

        parameter_space=CubicParameterSpace({'R': (), 'B': ()}, 0.0001, 1.),

        name='Blatt1_Aufgabe_2c'
    )

im aktuellen Arbeitsverzeichnis, so können wir diese als Python-Modul importieren:

In [3]:
import problems

Die darin enthaltenen Methoden können wir dann wie gewohnt verwenden:

In [4]:
for f in [problems.problem_1_1_a, problems.problem_1_1_b,
          problems.problem_1_1_c, problems.problem_1_1_d]:
    p = f()
    d, _ = discretize_stationary_cg(p, diameter=1/100)
    d.visualize(d.solve(), legend=d.name)

Den parametrischen Problemen haben wir in problems.py bereits einen ParameterSpace zugewiesen, der automatisch an die Diskretisierung vererbt wird:

In [5]:
for f in [problems.problem_1_2_a, problems.problem_1_2_b,
          problems.problem_1_2_c]:
    p = f()
    d, _ = discretize_stationary_cg(p, diameter=1/100)
    mu = d.parameter_space.sample_randomly(1)[0]
    d.visualize(d.solve(mu), legend=str(mu))

Aufgabe 1 b)

In algorithms.py ist, wie in Aufgabe 2 b) verlangt, bereits auch die Orthonormalisierung der Basis in greedy_approximation implementiert:

In [6]:
%cat algorithms.py

import numpy as np
from pymor.algorithms.gram_schmidt import gram_schmidt


def orthogonal_projection(U, basis, product=None,
                          basis_is_orthogonal=False,
                          return_projection_matrix=False):

    rhs = product.apply2(basis, U) if product else basis.dot(U)
    if basis_is_orthogonal:
        M = np.eye(len(basis))
        v = rhs
    else:
        M = product.apply2(basis, basis) if product else basis.gramian()
        try:
            v = np.linalg.solve(M, rhs)
        except np.linalg.LinAlgError:
            v = np.zeros(rhs.shape)

    U_proj = basis.lincomb(v.T)

    if return_projection_matrix:
        return U_proj, M
    else:
        return U_proj


def greedy_approximation(U, basis_size, rtol=0, product=None, orthogonalize=False):
    # initialization
    basis = U.space.empty()
    max_errors = []

    # compute norms of vectors in U for relative approximation error computation
    norms = np.sqrt(product.pairwise_apply2(U, U) if product else U.pairwise_dot(U))

    # main loop
    for n in range(basis_size):

        # compute projections
        # since the greedy basis is hierarchical, calling orthogonal_projection is
        # suboptimal since previous results could be reused
        U_proj = orthogonal_projection(U, basis, product=product, basis_is_orthogonal=orthogonalize)

        # compute projection errors
        U_err = U - U_proj
        errors = np.sqrt(product.pairwise_apply2(U_err, U_err) if product else U_err.pairwise_dot(U_err))

        if np.all(errors / norms < rtol):
            break  # relative approximation error threshold reached

        # select vector for basis extension
        max_errors.append(np.max(errors))
        basis.append(U[np.argmax(errors)])

        if orthogonalize:
            gram_schmidt(basis, offset=n, product=product, copy=False)
            if len(basis) == n:
                break  # Gram-Schmidt has removed newly added vector since it is numerically linear dependent

    return basis, max_errors

Natürlich kann für eigene Module auch from ... import ... verwendet werden:

In [7]:
from algorithms import orthogonal_projection, greedy_approximation

Aufgabe 2 a)

Im Folgenden verwenden wir problem_1_2_b, bei dem die Effekte der unterschiedlichen Basisgenerierungsverfahen am besten beobachtet werden können:

In [8]:
d, _ = discretize_stationary_cg(problems.problem_1_2_b(), diameter=1/100)

Zunächst berechnen wir wieder Snapshots für den Approximationsraum (U) und zur Berechnung des Approximationsfehlers (V):

In [9]:
U = d.solution_space.empty()
print('U: ', end='')
for mu in d.parameter_space.sample_randomly(30):
    U.append(d.solve(mu))
    print('.', end='', flush=True)
    
V = d.solution_space.empty()
print('\nV: ', end='')
for mu in d.parameter_space.sample_randomly(100, seed=99):
    V.append(d.solve(mu))
    print('.', end='', flush=True)

U: ..............................
V: ....................................................................................................

Den Projektionsfehler berechnen wir wie gewohnt, verwenden dabei aber unsere neue orthogonal_projection-Methode:

In [10]:
errors = []
for N in range(len(U)):
    V_proj = orthogonal_projection(V, U[:N], product=d.h1_0_semi_product)
    errors.append(np.max(d.h1_0_semi_norm(V - V_proj)))

Zur Berechnung einer Orthonormalbasis verwenden wir die pymor.algorithms.gram_schmidt.gram_schmidt-Methode. Wenn dabei nicht copy=False als Parameter angegeben wird, bleibt das Ursprungs-Array U unverändert:

In [11]:
basis = gram_schmidt(U, product=d.h1_0_semi_product)
errors_orth = []
for N in range(len(U)):
    V_proj = orthogonal_projection(V, basis[:N], product=d.h1_0_semi_product, basis_is_orthogonal=True)
    errors_orth.append(np.max(d.h1_0_semi_norm(V - V_proj)))

Schließlich plotten wir die Ergebnisse. Dabei zeigt sich klar die erhöhte numerische Stabilität (Kondition der Projektionsmatrix ist 1), wenn die Projektionsbasis zunächst orthonormalisiert wird:

In [12]:
from matplotlib import pyplot as plt
plt.figure()
plt.semilogy(errors_orth, marker='o', label='w. orth')
plt.semilogy(errors, label='wo. orth')
plt.legend()
plt.show()

Aufgabe 2 b)

Zusätzlich berechnen wir Approximationsräume mit dem Greedy-Algorithmus (mit und ohne Orthonormalisierung), den wir in greedy_approximation implementiert haben. Wieder zeigt sich eine höhere numerische Stabilität, wenn die Basis im Greedy-Algorithmus zusätzlich orthonormalisiert wird:

In [13]:
greedy_basis, _ = greedy_approximation(U, 30, product=d.h1_0_semi_product, orthogonalize=False)
greedy_orth_basis, _ = greedy_approximation(U, 30, product=d.h1_0_semi_product, orthogonalize=True)

errors_greedy = []
errors_greedy_orth = []
for N in range(len(U)):
    V_proj = orthogonal_projection(V, greedy_basis[:N], product=d.h1_0_semi_product)
    errors_greedy.append(np.max(d.h1_0_semi_norm(V - V_proj)))
    
    V_proj = orthogonal_projection(V, greedy_orth_basis[:N], product=d.h1_0_semi_product)
    errors_greedy_orth.append(np.max(d.h1_0_semi_norm(V - V_proj)))

plt.figure()
plt.semilogy(errors_orth, marker='o', label='w. orth')
plt.semilogy(errors, label='wo. orth')
plt.semilogy(errors_greedy_orth, marker='s', label='greedy w. orth')
plt.semilogy(errors_greedy, label='greedy wo. orth')
plt.legend()
plt.show()

Zusatzaufgabe 1

Zur Lösung dieser Aufgabe benötigen wir die Image-, ImageDraw- und ImageFont-Submodule von PIL:

In [14]:
from PIL import Image, ImageDraw, ImageFont

Mit PIL können wir programmatisch Bitmap-Grafiken generieren. Da BitmapFunction gegenwärtig jedoch immer den Namen einer Datei mit der entsprechenden Grafik erwartet, müssen wir die generierten Grafik in einer Datei zwischenspeichern. Hierzu können wir das tempfile-Modul der Python-Standardbibliothek verwenden, um leicht temporäre, automatisch wieder gelöschte Dateien anzulegen:

In [15]:
from tempfile import NamedTemporaryFile

Nun definieren wir die word_problem-Methode:

In [16]:
def word_problem(text):
    font = ImageFont.truetype('DejaVuSans-Bold.ttf', 64)  # load some font from file of given size
    size = font.getsize(text)                             # compute width and height of rendered text
    size = (size[0] + 20, size[1] + 20)                   # add a border of 10 pixels around the text
    
    def make_bitmap_function(char_num):                   # we need to genereate a BitmapFunction for each character
        img = Image.new('L', size)                        # create new Image object of given dimensions
        d = ImageDraw.Draw(img)                           # create ImageDraw object for the given Image
        
        # in order to position the character correctly, we first draw all characters from the first
        # up to the wanted character
        d.text((10, 10), text[:char_num + 1], font=font, fill=255)
        
        # next we erase all previous character by drawing a black rectangle
        d.rectangle(((0,0), (font.getsize(text[:char_num])[0] + 10, size[1])), fill=0, outline=0)
        
        # open a new temporary file
        with NamedTemporaryFile(suffix='.png') as f:  # after leaving this 'with' block, the temporary
                                                      # file is automatically deleted
            img.save(f, format='png')
            return BitmapFunction(f.name, bounding_box=[(0,0), size], range=[0., 1.])
    
    # create BitmapFunctions for each character
    dfs = [make_bitmap_function(n) for n in range(len(text))]
    
    # create an indicator function for the background
    background = ConstantFunction(1., 2) - LincombFunction(dfs, np.ones(len(dfs)))
    
    # form the linear combination
    dfs = [background,] + dfs
    coefficients = [1,] + [ProjectionParameterFunctional('chars', (len(text),), (i,)) for i in range(len(text))]
    diffusion = LincombFunction(dfs, coefficients)
    
    return StationaryProblem(
        domain=RectDomain(dfs[1].bounding_box, bottom='neumann'),
        neumann_data=ConstantFunction(-1., 2),
        diffusion=diffusion,
        parameter_space=CubicParameterSpace(diffusion.parameter_type, 0.01, 1.)
    )

Und testen diese mit dem String 'pyMOR':

In [17]:
d, _ = discretize_stationary_cg(word_problem('pyMOR'), diameter=1.)
d.visualize(d.solve([0.01, 0.1, 0.01, 0.1, 0.01]))