error.py 16 KB
Newer Older
1
# This file is part of the pyMOR project (http://www.pymor.org).
2
# Copyright 2013-2020 pyMOR developers and contributors. All rights reserved.
3 4 5 6 7 8 9 10
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)

from numbers import Number
import time

import numpy as np

from pymor.core.logger import getLogger
Stephan Rave's avatar
Stephan Rave committed
11
from pymor.models.basic import StationaryModel
12 13 14
from pymor.parallel.dummy import dummy_pool


15 16
def reduction_error_analysis(rom, fom, reductor, test_mus,
                             basis_sizes=0,
Tim Keil's avatar
Tim Keil committed
17 18 19
                             error_estimator=True, condition=False, error_norms=(),
                             error_norm_names=None,
                             error_estimator_norm_index=0, custom=(),
20 21 22 23 24
                             plot=False, plot_custom_logarithmic=True,
                             pool=dummy_pool):
    """Analyze the model reduction error.

    The maximum model reduction error is estimated by solving the reduced
Stephan Rave's avatar
Stephan Rave committed
25
    |Model| for given random |Parameters|.
26 27 28

    Parameters
    ----------
Stephan Rave's avatar
Stephan Rave committed
29
    rom
Stephan Rave's avatar
Stephan Rave committed
30
        The reduced |Model|.
Stephan Rave's avatar
Stephan Rave committed
31
    fom
Stephan Rave's avatar
Stephan Rave committed
32
        The high-dimensional |Model|.
33
    reductor
Stephan Rave's avatar
Stephan Rave committed
34
        The reductor which has created `rom`.
35
    test_mus
36
        List of |Parameters| to compute the errors for.
37 38 39 40
    basis_sizes
        Either a list of reduced basis dimensions to consider, or
        the number of dimensions (which are then selected equidistantly,
        always including the maximum reduced space dimension).
41 42
        The dimensions are input for the `dim`-Parameter of
        `reductor.reduce()`.
Tim Keil's avatar
Tim Keil committed
43
    error_estimator
Stephan Rave's avatar
Stephan Rave committed
44
        If `True` evaluate the error estimator of `rom`
45 46 47
        on the test |Parameters|.
    condition
        If `True`, compute the condition of the reduced system matrix
Stephan Rave's avatar
Stephan Rave committed
48
        for the given test |Parameters| (can only be specified if
Stephan Rave's avatar
Stephan Rave committed
49 50
        `rom` is an instance of |StationaryModel|
        and `rom.operator` is linear).
51 52 53 54 55
    error_norms
        List of norms in which to compute the model reduction error.
    error_norm_names
        Names of the norms given by `error_norms`. If `None`, the
        `name` attributes of the given norms are used.
Tim Keil's avatar
Tim Keil committed
56 57
    error_estimator_norm_index
        When `error_estimator` is `True` and `error_norms` are specified,
58
        this is the index of the norm in `error_norms` w.r.t. which
Tim Keil's avatar
Tim Keil committed
59
        to compute the effectivity of the error estimator.
60
    custom
61 62 63
        List of custom functions which are evaluated for each test
        |parameter values| and basis size. The functions must have
        the signature ::
64

Stephan Rave's avatar
Stephan Rave committed
65
            def custom_value(rom, fom, reductor, mu, dim):
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
                pass

    plot
        If `True`, generate a plot of the computed quantities w.r.t.
        the basis size.
    plot_custom_logarithmic
        If `True`, use a logarithmic y-axis to plot the computed custom
        values.
    pool
        If not `None`, the |WorkerPool| to use for parallelization.

    Returns
    -------
    Dict with the following fields:

Tim Keil's avatar
Tim Keil committed
81
        :mus:                       The test |Parameters| which have been considered.
Stephan Rave's avatar
Stephan Rave committed
82

Tim Keil's avatar
Tim Keil committed
83
        :basis_sizes:               The reduced basis dimensions which have been considered.
Stephan Rave's avatar
Stephan Rave committed
84

Tim Keil's avatar
Tim Keil committed
85 86 87 88
        :norms:                     |Array| of the norms of the high-dimensional solutions
                                    w.r.t. all given test |Parameters|, reduced basis
                                    dimensions and norms in `error_norms`.
                                    (Only present when `error_norms` has been specified.)
Stephan Rave's avatar
Stephan Rave committed
89

Tim Keil's avatar
Tim Keil committed
90
        :max_norms:                 Maxima of `norms` over the given test |Parameters|.
Stephan Rave's avatar
Stephan Rave committed
91

Tim Keil's avatar
Tim Keil committed
92
        :max_norm_mus:              |Parameters| corresponding to `max_norms`.
Stephan Rave's avatar
Stephan Rave committed
93

Tim Keil's avatar
Tim Keil committed
94 95 96 97
        :errors:                    |Array| of the norms of the model reduction errors
                                    w.r.t. all given test |Parameters|, reduced basis
                                    dimensions and norms in `error_norms`.
                                    (Only present when `error_norms` has been specified.)
Stephan Rave's avatar
Stephan Rave committed
98

Tim Keil's avatar
Tim Keil committed
99
        :max_errors:                Maxima of `errors` over the given test |Parameters|.
Stephan Rave's avatar
Stephan Rave committed
100

Tim Keil's avatar
Tim Keil committed
101
        :max_error_mus:             |Parameters| corresponding to `max_errors`.
Stephan Rave's avatar
Stephan Rave committed
102

Tim Keil's avatar
Tim Keil committed
103 104
        :rel_errors:                `errors` divided by `norms`.
                                    (Only present when `error_norms` has been specified.)
Stephan Rave's avatar
Stephan Rave committed
105

Tim Keil's avatar
Tim Keil committed
106
        :max_rel_errors:            Maxima of `rel_errors` over the given test |Parameters|.
Stephan Rave's avatar
Stephan Rave committed
107

Tim Keil's avatar
Tim Keil committed
108
        :max_rel_error_mus:         |Parameters| corresponding to `max_rel_errors`.
Stephan Rave's avatar
Stephan Rave committed
109

Tim Keil's avatar
Tim Keil committed
110 111
        :error_norm_names:          Names of the given `error_norms`.
                                    (Only present when `error_norms` has been specified.)
Stephan Rave's avatar
Stephan Rave committed
112

Tim Keil's avatar
Tim Keil committed
113 114 115 116
        :error_estimates:           |Array| of the model reduction error estimates
                                    w.r.t. all given test |Parameters| and reduced basis
                                    dimensions.
                                    (Only present when `error_estimator` is `True`.)
Stephan Rave's avatar
Stephan Rave committed
117

Tim Keil's avatar
Tim Keil committed
118
        :max_error_estimate:        Maxima of `error_estimates` over the given test |Parameters|.
Stephan Rave's avatar
Stephan Rave committed
119

Tim Keil's avatar
Tim Keil committed
120
        :max_error_estimate_mus:    |Parameters| corresponding to `max_error_estimates`.
Stephan Rave's avatar
Stephan Rave committed
121

Tim Keil's avatar
Tim Keil committed
122 123 124
        :effectivities:             `errors` divided by `error_estimates`.
                                    (Only present when `error_estimator` is `True` and `error_norms`
                                    has been specified.)
Stephan Rave's avatar
Stephan Rave committed
125

Tim Keil's avatar
Tim Keil committed
126
        :min_effectivities:         Minima of `effectivities` over the given test |Parameters|.
Stephan Rave's avatar
Stephan Rave committed
127

Tim Keil's avatar
Tim Keil committed
128
        :min_effectivity_mus:       |Parameters| corresponding to `min_effectivities`.
Stephan Rave's avatar
Stephan Rave committed
129

Tim Keil's avatar
Tim Keil committed
130
        :max_effectivities:         Maxima of `effectivities` over the given test |Parameters|.
Stephan Rave's avatar
Stephan Rave committed
131

Tim Keil's avatar
Tim Keil committed
132
        :max_effectivity_mus:       |Parameters| corresponding to `max_effectivities`.
Stephan Rave's avatar
Stephan Rave committed
133

Tim Keil's avatar
Tim Keil committed
134 135 136 137
        :errors:                    |Array| of the reduced system matrix conditions
                                    w.r.t. all given test |Parameters| and reduced basis
                                    dimensions.
                                    (Only present when `conditions` is `True`.)
Stephan Rave's avatar
Stephan Rave committed
138

Tim Keil's avatar
Tim Keil committed
139
        :max_conditions:            Maxima of `conditions` over the given test |Parameters|.
Stephan Rave's avatar
Stephan Rave committed
140

Tim Keil's avatar
Tim Keil committed
141
        :max_condition_mus:         |Parameters| corresponding to `max_conditions`.
Stephan Rave's avatar
Stephan Rave committed
142

Tim Keil's avatar
Tim Keil committed
143 144 145 146
        :custom_values:             |Array| of custom function evaluations
                                    w.r.t. all given test |Parameters|, reduced basis
                                    dimensions and functions in `custom`.
                                    (Only present when `custom` has been specified.)
Stephan Rave's avatar
Stephan Rave committed
147

Tim Keil's avatar
Tim Keil committed
148
        :max_custom_values:         Maxima of `custom_values` over the given test |Parameters|.
Stephan Rave's avatar
Stephan Rave committed
149

Tim Keil's avatar
Tim Keil committed
150
        :max_custom_values_mus:     |Parameters| corresponding to `max_custom_values`.
Stephan Rave's avatar
Stephan Rave committed
151

Tim Keil's avatar
Tim Keil committed
152
        :time:                      Time (in seconds) needed for the error analysis.
Stephan Rave's avatar
Stephan Rave committed
153

Tim Keil's avatar
Tim Keil committed
154 155
        :summary:                   String containing a summary of all computed quantities for
                                    the largest (last) considered basis size.
Stephan Rave's avatar
Stephan Rave committed
156

Tim Keil's avatar
Tim Keil committed
157 158
        :figure:                    The figure containing the generated plots.
                                    (Only present when `plot` is `True`.)
159 160
    """

Stephan Rave's avatar
Stephan Rave committed
161
    assert not error_norms or (fom and reductor)
162 163
    assert error_norm_names is None or len(error_norm_names) == len(error_norms)
    assert not condition \
Stephan Rave's avatar
Stephan Rave committed
164
        or isinstance(rom, StationaryModel) and rom.operator.linear
165 166 167 168 169

    logger = getLogger('pymor.algorithms.error')
    if pool is None or pool is dummy_pool:
        pool = dummy_pool
    else:
170
        logger.info(f'Using pool of {len(pool)} workers for error analysis')
171 172 173 174 175

    tic = time.time()

    if isinstance(basis_sizes, Number):
        if basis_sizes == 1:
Stephan Rave's avatar
Stephan Rave committed
176
            basis_sizes = [rom.solution_space.dim]
177 178
        else:
            if basis_sizes == 0:
Stephan Rave's avatar
Stephan Rave committed
179 180 181
                basis_sizes = rom.solution_space.dim + 1
            basis_sizes = min(rom.solution_space.dim + 1, basis_sizes)
            basis_sizes = np.linspace(0, rom.solution_space.dim, basis_sizes).astype(int)
182 183 184
    if error_norm_names is None:
        error_norm_names = tuple(norm.name for norm in error_norms)

Tim Keil's avatar
Tim Keil committed
185 186
    norms, error_estimates, errors, conditions, custom_values = \
        list(zip(*pool.map(_compute_errors, test_mus, fom=fom, reductor=reductor, error_estimator=error_estimator,
Stephan Rave's avatar
Stephan Rave committed
187
                           error_norms=error_norms, condition=condition, custom=custom, basis_sizes=basis_sizes)))
188 189 190 191 192 193
    print()

    result = {}
    result['mus'] = test_mus = np.array(test_mus)
    result['basis_sizes'] = basis_sizes

194
    summary = [('number of samples', str(len(test_mus)))]
195 196 197 198 199 200 201 202 203 204 205 206 207 208

    if error_norms:
        result['norms'] = norms = np.array(norms)
        result['max_norms'] = max_norms = np.max(norms, axis=0)
        result['max_norm_mus'] = max_norm_mus = test_mus[np.argmax(norms, axis=0)]
        result['errors'] = errors = np.array(errors)
        result['max_errors'] = max_errors = np.max(errors, axis=0)
        result['max_error_mus'] = max_error_mus = test_mus[np.argmax(errors, axis=0)]
        result['rel_errors'] = rel_errors = errors / norms[:, :, np.newaxis]
        result['max_rel_errors'] = np.max(rel_errors, axis=0)
        result['max_rel_error_mus'] = test_mus[np.argmax(rel_errors, axis=0)]
        for name, norm, norm_mu, error, error_mu in zip(error_norm_names,
                                                        max_norms, max_norm_mus,
                                                        max_errors[:, -1], max_error_mus[:, -1]):
209 210 211 212
            summary.append((f'maximum {name}-norm',
                            f'{norm:.7e} (mu = {error_mu})'))
            summary.append((f'maximum {name}-error',
                            f'{error:.7e} (mu = {error_mu})'))
213 214
        result['error_norm_names'] = error_norm_names

Tim Keil's avatar
Tim Keil committed
215 216 217 218
    if error_estimator:
        result['error_estimates'] = error_estimates = np.array(error_estimates)
        result['max_error_estimates'] = max_error_estimates = np.max(error_estimates, axis=0)
        result['max_error_estimate_mus'] = max_error_estimate_mus = test_mus[np.argmax(error_estimates, axis=0)]
219
        summary.append(('maximum estimated error',
Tim Keil's avatar
Tim Keil committed
220
                        f'{max_error_estimates[-1]:.7e} (mu = {max_error_estimate_mus[-1]})'))
221

Tim Keil's avatar
Tim Keil committed
222 223
    if error_estimator and error_norms:
        result['effectivities'] = effectivities = errors[:, error_estimator_norm_index, :] / error_estimates
224 225 226 227
        result['max_effectivities'] = max_effectivities = np.max(effectivities, axis=0)
        result['max_effectivity_mus'] = max_effectivity_mus = test_mus[np.argmax(effectivities, axis=0)]
        result['min_effectivities'] = min_effectivities = np.min(effectivities, axis=0)
        result['min_effectivity_mus'] = min_effectivity_mus = test_mus[np.argmin(effectivities, axis=0)]
Tim Keil's avatar
Tim Keil committed
228
        summary.append(('minimum error estimator effectivity',
Stephan Rave's avatar
Stephan Rave committed
229
                        f'{min_effectivities[-1]:.7e} (mu = {min_effectivity_mus[-1]})'))
Tim Keil's avatar
Tim Keil committed
230
        summary.append(('maximum error estimator effectivity',
Stephan Rave's avatar
Stephan Rave committed
231
                        f'{max_effectivities[-1]:.7e} (mu = {max_effectivity_mus[-1]})'))
232 233 234 235 236 237

    if condition:
        result['conditions'] = conditions = np.array(conditions)
        result['max_conditions'] = max_conditions = np.max(conditions, axis=0)
        result['max_condition_mus'] = max_condition_mus = test_mus[np.argmax(conditions, axis=0)]
        summary.append(('maximum system matrix condition',
Stephan Rave's avatar
Stephan Rave committed
238
                        f'{max_conditions[-1]:.7e} (mu = {max_condition_mus[-1]})'))
239 240 241 242 243 244

    if custom:
        result['custom_values'] = custom_values = np.array(custom_values)
        result['max_custom_values'] = max_custom_values = np.max(custom_values, axis=0)
        result['max_custom_values_mus'] = max_custom_values_mus = test_mus[np.argmax(custom_values, axis=0)]
        for i, (value, mu) in enumerate(zip(max_custom_values[:, -1], max_custom_values_mus[:, -1])):
245 246
            summary.append((f'maximum custom value {i}',
                            f'{value:.7e} (mu = {mu})'))
247 248 249 250 251

    toc = time.time()
    result['time'] = toc - tic
    summary.append(('elapsed time', str(toc - tic)))

Stephan Rave's avatar
Stephan Rave committed
252 253
    summary_fields, summary_values = list(zip(*summary))
    summary_field_width = np.max(list(map(len, summary_fields))) + 2
Stephan Rave's avatar
Stephan Rave committed
254
    summary_lines = [f'    {field+":":{summary_field_width}} {value}'
255 256 257 258 259 260 261
                     for field, value in zip(summary_fields, summary_values)]
    summary = 'Stochastic error estimation:\n' + '\n'.join(summary_lines)
    result['summary'] = summary

    if plot:
        import matplotlib.pyplot as plt
        fig = plt.figure()
Tim Keil's avatar
Tim Keil committed
262
        num_plots = (int(bool(error_norms) or error_estimator) + int(bool(error_norms) and error_estimator)
263
                     + int(condition) + int(bool(custom)))
264 265
        current_plot = 1

Tim Keil's avatar
Tim Keil committed
266
        if bool(error_norms) or error_estimator:
267 268 269 270 271 272
            ax = fig.add_subplot(1, num_plots, current_plot)
            legend = []
            if error_norms:
                for name, errors in zip(error_norm_names, max_errors):
                    ax.semilogy(basis_sizes, errors)
                    legend.append(name)
Tim Keil's avatar
Tim Keil committed
273 274 275
            if error_estimator:
                ax.semilogy(basis_sizes, max_error_estimates)
                legend.append('error estimator')
276
            ax.legend(legend)
277
            ax.set_title('maximum errors')
278 279
            current_plot += 1

Tim Keil's avatar
Tim Keil committed
280
        if bool(error_norms) and error_estimator:
281
            ax = fig.add_subplot(1, num_plots, current_plot)
282 283
            ax.semilogy(basis_sizes, min_effectivities)
            ax.semilogy(basis_sizes, max_effectivities)
284
            ax.legend(('min', 'max'))
Tim Keil's avatar
Tim Keil committed
285
            ax.set_title('error estimator effectivities')
286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311
            current_plot += 1

        if condition:
            ax = fig.add_subplot(1, num_plots, current_plot)
            ax.semilogy(basis_sizes, max_conditions)
            ax.set_title('maximum condition')
            current_plot += 1

        if custom:
            ax = fig.add_subplot(1, num_plots, current_plot)
            legend = []
            for i, values in enumerate(custom_values):
                if plot_custom_logarithmic:
                    ax.semilogy(basis_sizes, values)
                else:
                    ax.plot(basis_sizes, values)
                legend.append('value ' + str(i))
            ax.legend(legend)
            ax.set_title('maximum custom values')
            current_plot += 1

        result['figure'] = fig

    return result


Tim Keil's avatar
Tim Keil committed
312
def _compute_errors(mu, fom, reductor, error_estimator, error_norms, condition, custom, basis_sizes):
313 314 315 316 317
    import sys

    print('.', end='')
    sys.stdout.flush()

Tim Keil's avatar
Tim Keil committed
318
    error_estimates = np.empty(len(basis_sizes)) if error_estimator else None
319 320 321 322 323
    norms = np.empty(len(error_norms))
    errors = np.empty((len(error_norms), len(basis_sizes)))
    conditions = np.empty(len(basis_sizes)) if condition else None
    custom_values = np.empty((len(custom), len(basis_sizes)))

Stephan Rave's avatar
Stephan Rave committed
324 325 326 327 328
    if fom:
        logging_disabled = fom.logging_disabled
        fom.disable_logging()
        U = fom.solve(mu)
        fom.disable_logging(logging_disabled)
329 330 331 332 333 334
        for i_norm, norm in enumerate(error_norms):
            n = norm(U)
            n = n[0] if hasattr(n, '__len__') else n
            norms[i_norm] = n

    for i_N, N in enumerate(basis_sizes):
335
        rom = reductor.reduce(dims={k: N for k in reductor.bases})
336 337
        result = rom.compute(solution=True, solution_error_estimate=error_estimator, mu=mu)
        u = result['solution']
Tim Keil's avatar
Tim Keil committed
338
        if error_estimator:
339
            e = result['solution_error_estimate']
340
            e = e[0] if hasattr(e, '__len__') else e
Tim Keil's avatar
Tim Keil committed
341
            error_estimates[i_N] = e
Stephan Rave's avatar
Stephan Rave committed
342
        if fom and reductor:
343
            URB = reductor.reconstruct(u)
344 345 346 347 348
            for i_norm, norm in enumerate(error_norms):
                e = norm(U - URB)
                e = e[0] if hasattr(e, '__len__') else e
                errors[i_norm, i_N] = e
        if condition:
Stephan Rave's avatar
Stephan Rave committed
349
            conditions[i_N] = np.linalg.cond(rom.operator.assemble(mu).matrix) if N > 0 else 0.
350
        for i_custom, cust in enumerate(custom):
Stephan Rave's avatar
Stephan Rave committed
351
            c = cust(rom=rom,
Stephan Rave's avatar
Stephan Rave committed
352
                     fom=fom,
353
                     reductor=reductor,
354 355 356 357 358
                     mu=mu,
                     dim=N)
            c = c[0] if hasattr(c, '__len__') else c
            custom_values[i_custom, i_N] = c

Tim Keil's avatar
Tim Keil committed
359
    return norms, error_estimates, errors, conditions, custom_values