GitHub

romtools.workflows.greedy.run_greedy

The greedy procedure iteratively constructs a reduced basis ROM until it reaches a desired tolerance.

The algorithm is as follows:

We generate a parameter training set, $\mathcal{D}_{\mathrm{train}}, |\mathcal{D}_{\mathrm{train}} | = N_{\mathrm{train}}$

We select an initial sample, $\mu_1 \in \mathcal{D}_{\mathrm{train}} \text{ and set } \mathcal{D}_{\mathrm{train}}=\mathcal{D}_{\mathrm{train}} - {\mu_1}$

We then solve the FOM to obtain the solution, $\mathbf{u}(\mu_1)$.

We select a second sample, $\mu_2 \in \mathcal{D}_{\mathrm{train}} \text{ and set } \mathcal{D}_{\mathrm{train}}=\mathcal{D}_{\mathrm{train}} - {\mu_2}$

We then solve the FOM to obtain the solution, $\mathbf{u}(\mu_2)$.

We employ the first two solutions to compute the trial space, e.g., $$\boldsymbol \Phi = \mathrm{orthogonalize}(\mathbf{u}(\mu_2)) - \mathbf{u}_{\mathrm{shift}}, \; \mathbf{u}_{\mathrm{shift}} = \mathbf{u}(\mu_1)$$

We then solve the resulting ROM for the remaining parameter samples $\mathcal{D}_{\mathrm{train}}$ to generate approximate solutions $\mathbf{u}(\mu), \mu \in \mathcal{D}_{\mathrm{train}}$

For each ROM solution $\mathbf{u}(\mu), \mu \in \mathcal{D}_{\mathrm{train}}$ we compute an error estimate, $ e \left(\mu \right) $

If the maximum error estimate is less than some tolerance, we exit the algorithm. If not, we:

Set $ \mu^* = \underset{ \mu \in \mathcal{D}_{\mathrm{train}} }{ \mathrm{arg\; max} } \; e \left(\mu \right) $

Remove $\mu^{*}$ from the training set, $\mathcal{D}_{\mathrm{train}}= \mathcal{D}_{\mathrm{train}} - {\mu^{*}}$

Solve the FOM for $\mathbf{u}(\mu^{*})$

Set the reduced basis to be $\boldsymbol \Phi = \mathrm{orthogonalize}(\boldsymbol \Phi, \mathbf{u}(\mu^{*})) - \mathbf{u}_{\mathrm{shift}}$

Go back to step 6 and continue until convergence.

This function implements the basic greedy workflow. In addition, we enable an adaptive error estimate based on a QoI. This is based on the fact that, throughout the greedy algorithm, we have a history of error indicators as well as a set of FOM and ROM runs for the same parameter instance. We leverage this data to improve the error estimate. Specifically, for the max error estimates over the first $j$ iterations, $e_1,\ldots,e_j$, we additionally can access QoI errors $e^q_1,\ldots,e^q_j$ from our FOM-ROM runs, and we define a scaling based on $$C = \frac{\sum_{i=1}^j | e^q_j| }{ \sum_{i=1}^j | e_j | }.$$ The error at the $j+1$ iteration can be approximated by $C e(\mu)$. This will adaptively scale the error estimate to match the QoI error as closely as possible, which can be helpful for defining exit criterion.

Alternatively, if calibrated_error is set to False, then the scaling factor $C$ will always be $1$.

View Source

  1#
  2# ************************************************************************
  3#
  4#                         ROM Tools and Workflows
  5# Copyright 2019 National Technology & Engineering Solutions of Sandia,LLC
  6#                              (NTESS)
  7#
  8# Under the terms of Contract DE-NA0003525 with NTESS, the
  9# U.S. Government retains certain rights in this software.
 10#
 11# ROM Tools and Workflows is licensed under BSD-3-Clause terms of use:
 12#
 13# Redistribution and use in source and binary forms, with or without
 14# modification, are permitted provided that the following conditions
 15# are met:
 16#
 17# 1. Redistributions of source code must retain the above copyright
 18# notice, this list of conditions and the following disclaimer.
 19#
 20# 2. Redistributions in binary form must reproduce the above copyright
 21# notice, this list of conditions and the following disclaimer in the
 22# documentation and/or other materials provided with the distribution.
 23#
 24# 3. Neither the name of the copyright holder nor the names of its
 25# contributors may be used to endorse or promote products derived
 26# from this software without specific prior written permission.
 27#
 28# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 29# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 30# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 31# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
 32# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 33# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 34# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 35# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 36# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
 37# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
 38# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 39# POSSIBILITY OF SUCH DAMAGE.
 40#
 41# Questions? Contact Eric Parish (ejparis@sandia.gov)
 42#
 43# ************************************************************************
 44#
 45
 46'''
 47The greedy procedure iteratively constructs a reduced basis ROM until it reaches a desired tolerance.
 48The algorithm is as follows:
 49 1. We generate a parameter training set, $\\mathcal{D}_{\\mathrm{train}}, |\\mathcal{D}_{\\mathrm{train}}
 50    | = N_{\\mathrm{train}}$
 51 2. We select an initial sample, $\\mu_1 \\in \\mathcal{D}_{\\mathrm{train}} \\text{ and set }
 52    \\mathcal{D}_{\\mathrm{train}}=\\mathcal{D}_{\\mathrm{train}} - \\{\\mu_1\\}$
 53 3. We then solve the FOM to obtain the solution, $\\mathbf{u}(\\mu_1)$.
 54 4. We select a second sample, $\\mu_2 \\in \\mathcal{D}_{\\mathrm{train}} \\text{ and set }
 55    \\mathcal{D}_{\\mathrm{train}}=\\mathcal{D}_{\\mathrm{train}} - \\{\\mu_2\\}$
 56 4. We then solve the FOM to obtain the solution, $\\mathbf{u}(\\mu_2)$.
 57 5. We employ the first two solutions to compute the trial space, e.g.,
 58    $$\\boldsymbol \\Phi = \\mathrm{orthogonalize}(\\mathbf{u}(\\mu_2)) -  \\mathbf{u}_{\\mathrm{shift}},
 59    \\; \\mathbf{u}_{\\mathrm{shift}} = \\mathbf{u}(\\mu_1)$$
 60 6. We then solve the resulting ROM for the remaining parameter samples $\\mathcal{D}_{\\mathrm{train}}$ to
 61    generate approximate solutions $\\mathbf{u}(\\mu), \\mu \\in \\mathcal{D}_{\\mathrm{train}}$
 62 7. For each ROM solution $\\mathbf{u}(\\mu), \\mu \\in \\mathcal{D}_{\\mathrm{train}}$ we compute an error
 63    estimate, $ e \\left(\\mu \\right) $
 64 8. If the maximum error estimate is less than some tolerance, we exit the algorithm. If not, we:
 65   - Set $ \\mu^* = \\underset{ \\mu \\in \\mathcal{D}_{\\mathrm{train}} }{ \\mathrm{arg\\; max} } \\; e
 66        \\left(\\mu \\right) $
 67   - Remove $\\mu^{\\*}$ from the training set, $\\mathcal{D}_{\\mathrm{train}}=
 68        \\mathcal{D}_{\\mathrm{train}} - \\{\\mu^{\\*}\\}$
 69   - Solve the FOM for $\\mathbf{u}(\\mu^{\\*})$
 70   - Set the reduced basis to be $\\boldsymbol \\Phi = \\mathrm{orthogonalize}(\\boldsymbol \\Phi,
 71      \\mathbf{u}(\\mu^{\\*})) -  \\mathbf{u}_{\\mathrm{shift}}$
 72   - Go back to step 6 and continue until convergence.
 73
 74This function implements the basic greedy workflow. In addition, we enable an adaptive
 75error estimate based on a QoI. This is based on the fact that, throughout the greedy algorithm, we have a
 76history of error indicators as well as a set of FOM and ROM runs for the same parameter instance. We leverage this data
 77to improve the error estimate. Specifically, for the max error estimates over the first $j$ iterations,
 78$e_1,\\ldots,e_j$, we additionally can access QoI errors $e^q_1,\\ldots,e^q_j$ from our FOM-ROM runs, and we define a
 79scaling based on $$C = \\frac{\\sum_{i=1}^j | e^q_j| }{ \\sum_{i=1}^j | e_j | }.$$
 80The error at the $j+1$ iteration can be approximated by $C e(\\mu)$.
 81This will adaptively scale the error estimate to match the QoI error as
 82closely as possible, which can be helpful for defining exit criterion.
 83
 84Alternatively, if `calibrated_error` is set to `False`, then the scaling factor $C$ will always be $1$.
 85'''
 86
 87import time
 88import numpy as np
 89
 90import romtools.linalg.linalg as la
 91from romtools.workflows.models import QoiModel
 92from romtools.workflows.parameter_spaces import ParameterSpace
 93from romtools.workflows.workflow_utils import create_empty_dir
 94from romtools.workflows.model_builders import QoiModelWithErrorEstimateBuilder
 95
 96
 97def _create_parameter_dict(parameter_names, parameter_values):
 98    return dict(zip(parameter_names, parameter_values))
 99
100def run_greedy(fom_model: QoiModel,
101               rom_model_builder: QoiModelWithErrorEstimateBuilder,
102               parameter_space: ParameterSpace,
103               absolute_greedy_work_directory: str,
104               tolerance: float,
105               testing_sample_size: int = 10,
106               random_seed: int = 1,
107               calibrated_error: bool=True):
108    '''
109    Main implementation of the greedy algorithm.
110    '''
111    greedy_directory = absolute_greedy_work_directory
112    create_empty_dir(greedy_directory)
113    offline_directory_prefix = 'offline_data'
114
115    run_directory_prefix = "run_"
116    greedy_file = open(f"{greedy_directory}/greedy_status.log", "w", encoding="utf-8") # pylint: disable=consider-using-with
117    greedy_file.write("Greedy reduced basis status \n")
118    fom_time = 0.
119    rom_time = 0.
120    basis_time = 0.
121
122    np.random.seed(random_seed)
123
124    # create parameter samples
125    parameter_samples = parameter_space.generate_samples(testing_sample_size)
126    parameter_names = parameter_space.get_names()
127
128    # Make FOM/ROM directories
129    for sample_index, sample in enumerate(parameter_samples):
130        fom_run_directory =  f'{greedy_directory}/fom/{run_directory_prefix}{sample_index}'
131
132        create_empty_dir(fom_run_directory)
133        parameter_dict = _create_parameter_dict(parameter_names, sample)
134        fom_model.populate_run_directory(fom_run_directory, parameter_dict)
135
136    training_samples = np.array([0, 1], dtype='int')
137    samples_left = np.arange(2, testing_sample_size)
138
139    # Run FOM training cases
140    t0 = time.time()
141    training_dirs = []
142
143    for i in training_samples:
144        greedy_file.write(f"Running FOM sample {i} \n")
145        sample_index = i
146        parameter_dict = _create_parameter_dict(parameter_names,parameter_samples[sample_index])
147        fom_run_directory = f'{greedy_directory}/fom/{run_directory_prefix}{sample_index}'
148        create_empty_dir(fom_run_directory)
149        fom_model.populate_run_directory(fom_run_directory,parameter_dict)
150        fom_model.run_model(fom_run_directory,parameter_dict)
151        fom_qoi = fom_model.compute_qoi(fom_run_directory,parameter_dict)
152        if i == 0:
153            fom_qois = fom_qoi[None]
154        else:
155            fom_qois = np.append(fom_qois,fom_qoi[None],axis=0)
156        training_dirs.append(fom_run_directory)
157    fom_time += time.time() - t0
158
159    # Create ROM bases
160    t0 = time.time()
161    greedy_file.write("Creating ROM bases \n")
162
163    outer_loop_counter = 0
164    updated_offline_data_dir = f'{greedy_directory}/rom_iteration_{outer_loop_counter}/{offline_directory_prefix}/'
165
166    create_empty_dir(updated_offline_data_dir)
167    rom_model = rom_model_builder.build_from_training_dirs(updated_offline_data_dir,training_dirs)
168
169    basis_time += time.time() - t0
170
171    converged = False
172    max_error_indicators = np.zeros(0)
173    reg = QoIvsErrorIndicatorRegressor(calibrated_error)
174    qoi_errors = np.zeros(0)
175    predicted_qoi_errors = np.zeros(0)
176    greedy_file.flush()
177
178    while converged is False:
179        print(f'Greedy iteration # {outer_loop_counter}')
180        error_indicators = np.zeros(samples_left.size)
181
182        t0 = time.time()
183        greedy_file.write(f"Greedy iteration # {outer_loop_counter}\n" +
184                          f"Parameter samples: \n {parameter_samples}\n")
185        greedy_file.flush()
186        for counter, sample_index in enumerate(samples_left):
187            greedy_file.write(f"    Running ROM sample {sample_index}\n")
188            greedy_file.flush()
189            parameter_dict = _create_parameter_dict(parameter_names,parameter_samples[sample_index])
190            rom_run_directory = f'{greedy_directory}/rom_iteration_{outer_loop_counter}/{run_directory_prefix}{sample_index}'
191            create_empty_dir(rom_run_directory)
192            rom_model.populate_run_directory(rom_run_directory,parameter_dict)
193            rom_model.run_model(rom_run_directory,parameter_dict)
194            error_indicators[counter] = rom_model.compute_error_estimate(rom_run_directory,parameter_dict)
195
196        rom_time += time.time() - t0
197
198        sample_with_highest_error_indicator = samples_left[la.argmax(error_indicators)]
199        max_error_indicators = np.append(max_error_indicators,
200                                         np.amax(error_indicators))
201        greedy_file.write(f"Sample {sample_with_highest_error_indicator}"
202                          " had the highest error indicator of"
203                          f" {max_error_indicators[-1]} \n")
204
205        outer_loop_counter += 1
206        if outer_loop_counter > 1:
207            predicted_max_qoi_error = reg.predict(error_indicators[la.argmax(error_indicators)])
208            greedy_file.write("Our MLEM error estimate is "
209                              f"{predicted_max_qoi_error}\n")
210            greedy_file.flush()
211            predicted_qoi_errors = np.append(predicted_qoi_errors,
212                                             predicted_max_qoi_error)
213            if np.amax(predicted_max_qoi_error < tolerance):
214                converged = True
215                print('Run converged, max approximate qoi error'
216                      f' = {predicted_max_qoi_error}')
217                break
218
219        t0 = time.time()
220        greedy_file.write("Running FOM sample"
221                          f" {sample_with_highest_error_indicator}\n")
222        greedy_file.flush()
223
224
225        ## Identify the sample with the highest error and run FOM
226        fom_run_directory = f'{greedy_directory}/fom/{run_directory_prefix}{sample_with_highest_error_indicator}'
227        parameter_dict = _create_parameter_dict(parameter_names,parameter_samples[sample_with_highest_error_indicator])
228        fom_model.run_model(fom_run_directory,parameter_dict)
229        fom_qoi = fom_model.compute_qoi(fom_run_directory,parameter_dict)
230
231        fom_time += time.time() - t0
232
233        # Get ROM QoI to calibrate our error estimator
234        outer_loop_counter_m_1 = outer_loop_counter - 1
235        rom_run_directory = (
236            f'{greedy_directory}/rom_iteration_{outer_loop_counter_m_1}/'
237            f'{run_directory_prefix}{sample_with_highest_error_indicator}'
238        )
239        rom_qoi = rom_model.compute_qoi(rom_run_directory,parameter_dict)
240        qoi_error = np.linalg.norm(rom_qoi - fom_qoi) / np.linalg.norm(fom_qoi)
241        greedy_file.write(f"Sample {sample_with_highest_error_indicator} had "
242                          f"an error of {qoi_error}\n")
243        qoi_errors = np.append(qoi_errors, qoi_error)
244        reg.fit(max_error_indicators, qoi_errors)
245
246        ## Update our samples
247        training_samples = np.append(training_samples,
248                                     sample_with_highest_error_indicator)
249        samples_left = np.delete(samples_left, la.argmax(error_indicators))
250
251        # Update ROM basis
252        t0 = time.time()
253        greedy_file.write("Computing ROM bases \n")
254        greedy_file.flush()
255        training_dirs = [f'{greedy_directory}/fom/run_{i}' for i in training_samples]
256
257        updated_offline_data_dir = f'{greedy_directory}/rom_iteration_{outer_loop_counter}/{offline_directory_prefix}/'
258        create_empty_dir(updated_offline_data_dir)
259        rom_model = rom_model_builder.build_from_training_dirs(updated_offline_data_dir,training_dirs)
260
261        basis_time += time.time() - t0
262
263        # Add a new sample
264        new_parameter_sample = parameter_space.generate_samples(1)
265        parameter_samples = np.append(parameter_samples,
266                                      new_parameter_sample, axis=0)
267        new_sample_number = testing_sample_size + outer_loop_counter - 1
268        rom_run_directory = f'{greedy_directory}/rom_iteration_{outer_loop_counter}/{run_directory_prefix}{new_sample_number}'
269        fom_run_directory = f'{greedy_directory}/fom/{run_directory_prefix}{new_sample_number}'
270        create_empty_dir(rom_run_directory)
271        create_empty_dir(fom_run_directory)
272        parameter_dict = _create_parameter_dict(parameter_names,new_parameter_sample[0])
273
274        rom_model.populate_run_directory(rom_run_directory,parameter_dict)
275        fom_model.populate_run_directory(fom_run_directory,parameter_dict)
276
277        samples_left = np.append(samples_left, new_sample_number)
278        greedy_file.flush()
279        np.savez(f'{greedy_directory}/greedy_stats',
280                 max_error_indicators=max_error_indicators,
281                 qoi_errors=qoi_errors,
282                 predicted_qoi_errors=predicted_qoi_errors,
283                 training_samples=training_samples,
284                 fom_time=fom_time,
285                 rom_time=rom_time,
286                 basis_time=basis_time)
287
288    greedy_file.close()
289
290
291class QoIvsErrorIndicatorRegressor:
292    '''
293    Regressor for modeling the relationship between QoI error and error
294    indicator.
295
296    This class provides a simple linear regression model for estimating the
297    scaling factor (c) between QoI error and the error indicator.
298    '''
299    def __init__(self, calibrated_error: bool=True):
300        '''
301        Initializes the regressor with a default scaling factor of 1.0.
302        '''
303        self.__c = 1.
304        self.__calibrated_error = calibrated_error
305
306    def fit(self, x, y):
307        '''
308        If calibrated_error was specified, this function fits the regression
309        model to the provided data points (x, y) to estimate the scaling factor.
310
311        Args:
312            x: Error indicator values.
313            y: Corresponding QoI error values.
314        '''
315        if self.__calibrated_error:
316            self.__c = np.mean(y) / np.mean(x)
317
318    def predict(self, x):
319        '''
320        Predicts QoI error based on the error indicator using the estimated
321        scaling factor.
322
323        Args:
324            x: Error indicator value(s) for prediction.
325        '''
326        return self.__c*x

def run_greedy( fom_model: romtools.workflows.models.QoiModel, rom_model_builder: romtools.workflows.model_builders.QoiModelWithErrorEstimateBuilder, parameter_space: romtools.workflows.parameter_spaces.ParameterSpace, absolute_greedy_work_directory: str, tolerance: float, testing_sample_size: int = 10, random_seed: int = 1, calibrated_error: bool = True): View Source

101def run_greedy(fom_model: QoiModel,
102               rom_model_builder: QoiModelWithErrorEstimateBuilder,
103               parameter_space: ParameterSpace,
104               absolute_greedy_work_directory: str,
105               tolerance: float,
106               testing_sample_size: int = 10,
107               random_seed: int = 1,
108               calibrated_error: bool=True):
109    '''
110    Main implementation of the greedy algorithm.
111    '''
112    greedy_directory = absolute_greedy_work_directory
113    create_empty_dir(greedy_directory)
114    offline_directory_prefix = 'offline_data'
115
116    run_directory_prefix = "run_"
117    greedy_file = open(f"{greedy_directory}/greedy_status.log", "w", encoding="utf-8") # pylint: disable=consider-using-with
118    greedy_file.write("Greedy reduced basis status \n")
119    fom_time = 0.
120    rom_time = 0.
121    basis_time = 0.
122
123    np.random.seed(random_seed)
124
125    # create parameter samples
126    parameter_samples = parameter_space.generate_samples(testing_sample_size)
127    parameter_names = parameter_space.get_names()
128
129    # Make FOM/ROM directories
130    for sample_index, sample in enumerate(parameter_samples):
131        fom_run_directory =  f'{greedy_directory}/fom/{run_directory_prefix}{sample_index}'
132
133        create_empty_dir(fom_run_directory)
134        parameter_dict = _create_parameter_dict(parameter_names, sample)
135        fom_model.populate_run_directory(fom_run_directory, parameter_dict)
136
137    training_samples = np.array([0, 1], dtype='int')
138    samples_left = np.arange(2, testing_sample_size)
139
140    # Run FOM training cases
141    t0 = time.time()
142    training_dirs = []
143
144    for i in training_samples:
145        greedy_file.write(f"Running FOM sample {i} \n")
146        sample_index = i
147        parameter_dict = _create_parameter_dict(parameter_names,parameter_samples[sample_index])
148        fom_run_directory = f'{greedy_directory}/fom/{run_directory_prefix}{sample_index}'
149        create_empty_dir(fom_run_directory)
150        fom_model.populate_run_directory(fom_run_directory,parameter_dict)
151        fom_model.run_model(fom_run_directory,parameter_dict)
152        fom_qoi = fom_model.compute_qoi(fom_run_directory,parameter_dict)
153        if i == 0:
154            fom_qois = fom_qoi[None]
155        else:
156            fom_qois = np.append(fom_qois,fom_qoi[None],axis=0)
157        training_dirs.append(fom_run_directory)
158    fom_time += time.time() - t0
159
160    # Create ROM bases
161    t0 = time.time()
162    greedy_file.write("Creating ROM bases \n")
163
164    outer_loop_counter = 0
165    updated_offline_data_dir = f'{greedy_directory}/rom_iteration_{outer_loop_counter}/{offline_directory_prefix}/'
166
167    create_empty_dir(updated_offline_data_dir)
168    rom_model = rom_model_builder.build_from_training_dirs(updated_offline_data_dir,training_dirs)
169
170    basis_time += time.time() - t0
171
172    converged = False
173    max_error_indicators = np.zeros(0)
174    reg = QoIvsErrorIndicatorRegressor(calibrated_error)
175    qoi_errors = np.zeros(0)
176    predicted_qoi_errors = np.zeros(0)
177    greedy_file.flush()
178
179    while converged is False:
180        print(f'Greedy iteration # {outer_loop_counter}')
181        error_indicators = np.zeros(samples_left.size)
182
183        t0 = time.time()
184        greedy_file.write(f"Greedy iteration # {outer_loop_counter}\n" +
185                          f"Parameter samples: \n {parameter_samples}\n")
186        greedy_file.flush()
187        for counter, sample_index in enumerate(samples_left):
188            greedy_file.write(f"    Running ROM sample {sample_index}\n")
189            greedy_file.flush()
190            parameter_dict = _create_parameter_dict(parameter_names,parameter_samples[sample_index])
191            rom_run_directory = f'{greedy_directory}/rom_iteration_{outer_loop_counter}/{run_directory_prefix}{sample_index}'
192            create_empty_dir(rom_run_directory)
193            rom_model.populate_run_directory(rom_run_directory,parameter_dict)
194            rom_model.run_model(rom_run_directory,parameter_dict)
195            error_indicators[counter] = rom_model.compute_error_estimate(rom_run_directory,parameter_dict)
196
197        rom_time += time.time() - t0
198
199        sample_with_highest_error_indicator = samples_left[la.argmax(error_indicators)]
200        max_error_indicators = np.append(max_error_indicators,
201                                         np.amax(error_indicators))
202        greedy_file.write(f"Sample {sample_with_highest_error_indicator}"
203                          " had the highest error indicator of"
204                          f" {max_error_indicators[-1]} \n")
205
206        outer_loop_counter += 1
207        if outer_loop_counter > 1:
208            predicted_max_qoi_error = reg.predict(error_indicators[la.argmax(error_indicators)])
209            greedy_file.write("Our MLEM error estimate is "
210                              f"{predicted_max_qoi_error}\n")
211            greedy_file.flush()
212            predicted_qoi_errors = np.append(predicted_qoi_errors,
213                                             predicted_max_qoi_error)
214            if np.amax(predicted_max_qoi_error < tolerance):
215                converged = True
216                print('Run converged, max approximate qoi error'
217                      f' = {predicted_max_qoi_error}')
218                break
219
220        t0 = time.time()
221        greedy_file.write("Running FOM sample"
222                          f" {sample_with_highest_error_indicator}\n")
223        greedy_file.flush()
224
225
226        ## Identify the sample with the highest error and run FOM
227        fom_run_directory = f'{greedy_directory}/fom/{run_directory_prefix}{sample_with_highest_error_indicator}'
228        parameter_dict = _create_parameter_dict(parameter_names,parameter_samples[sample_with_highest_error_indicator])
229        fom_model.run_model(fom_run_directory,parameter_dict)
230        fom_qoi = fom_model.compute_qoi(fom_run_directory,parameter_dict)
231
232        fom_time += time.time() - t0
233
234        # Get ROM QoI to calibrate our error estimator
235        outer_loop_counter_m_1 = outer_loop_counter - 1
236        rom_run_directory = (
237            f'{greedy_directory}/rom_iteration_{outer_loop_counter_m_1}/'
238            f'{run_directory_prefix}{sample_with_highest_error_indicator}'
239        )
240        rom_qoi = rom_model.compute_qoi(rom_run_directory,parameter_dict)
241        qoi_error = np.linalg.norm(rom_qoi - fom_qoi) / np.linalg.norm(fom_qoi)
242        greedy_file.write(f"Sample {sample_with_highest_error_indicator} had "
243                          f"an error of {qoi_error}\n")
244        qoi_errors = np.append(qoi_errors, qoi_error)
245        reg.fit(max_error_indicators, qoi_errors)
246
247        ## Update our samples
248        training_samples = np.append(training_samples,
249                                     sample_with_highest_error_indicator)
250        samples_left = np.delete(samples_left, la.argmax(error_indicators))
251
252        # Update ROM basis
253        t0 = time.time()
254        greedy_file.write("Computing ROM bases \n")
255        greedy_file.flush()
256        training_dirs = [f'{greedy_directory}/fom/run_{i}' for i in training_samples]
257
258        updated_offline_data_dir = f'{greedy_directory}/rom_iteration_{outer_loop_counter}/{offline_directory_prefix}/'
259        create_empty_dir(updated_offline_data_dir)
260        rom_model = rom_model_builder.build_from_training_dirs(updated_offline_data_dir,training_dirs)
261
262        basis_time += time.time() - t0
263
264        # Add a new sample
265        new_parameter_sample = parameter_space.generate_samples(1)
266        parameter_samples = np.append(parameter_samples,
267                                      new_parameter_sample, axis=0)
268        new_sample_number = testing_sample_size + outer_loop_counter - 1
269        rom_run_directory = f'{greedy_directory}/rom_iteration_{outer_loop_counter}/{run_directory_prefix}{new_sample_number}'
270        fom_run_directory = f'{greedy_directory}/fom/{run_directory_prefix}{new_sample_number}'
271        create_empty_dir(rom_run_directory)
272        create_empty_dir(fom_run_directory)
273        parameter_dict = _create_parameter_dict(parameter_names,new_parameter_sample[0])
274
275        rom_model.populate_run_directory(rom_run_directory,parameter_dict)
276        fom_model.populate_run_directory(fom_run_directory,parameter_dict)
277
278        samples_left = np.append(samples_left, new_sample_number)
279        greedy_file.flush()
280        np.savez(f'{greedy_directory}/greedy_stats',
281                 max_error_indicators=max_error_indicators,
282                 qoi_errors=qoi_errors,
283                 predicted_qoi_errors=predicted_qoi_errors,
284                 training_samples=training_samples,
285                 fom_time=fom_time,
286                 rom_time=rom_time,
287                 basis_time=basis_time)
288
289    greedy_file.close()

Main implementation of the greedy algorithm.

class QoIvsErrorIndicatorRegressor: View Source

292class QoIvsErrorIndicatorRegressor:
293    '''
294    Regressor for modeling the relationship between QoI error and error
295    indicator.
296
297    This class provides a simple linear regression model for estimating the
298    scaling factor (c) between QoI error and the error indicator.
299    '''
300    def __init__(self, calibrated_error: bool=True):
301        '''
302        Initializes the regressor with a default scaling factor of 1.0.
303        '''
304        self.__c = 1.
305        self.__calibrated_error = calibrated_error
306
307    def fit(self, x, y):
308        '''
309        If calibrated_error was specified, this function fits the regression
310        model to the provided data points (x, y) to estimate the scaling factor.
311
312        Args:
313            x: Error indicator values.
314            y: Corresponding QoI error values.
315        '''
316        if self.__calibrated_error:
317            self.__c = np.mean(y) / np.mean(x)
318
319    def predict(self, x):
320        '''
321        Predicts QoI error based on the error indicator using the estimated
322        scaling factor.
323
324        Args:
325            x: Error indicator value(s) for prediction.
326        '''
327        return self.__c*x

Regressor for modeling the relationship between QoI error and error indicator.

This class provides a simple linear regression model for estimating the scaling factor (c) between QoI error and the error indicator.

QoIvsErrorIndicatorRegressor(calibrated_error: bool = True) View Source

300    def __init__(self, calibrated_error: bool=True):
301        '''
302        Initializes the regressor with a default scaling factor of 1.0.
303        '''
304        self.__c = 1.
305        self.__calibrated_error = calibrated_error

Initializes the regressor with a default scaling factor of 1.0.

def fit(self, x, y): View Source

307    def fit(self, x, y):
308        '''
309        If calibrated_error was specified, this function fits the regression
310        model to the provided data points (x, y) to estimate the scaling factor.
311
312        Args:
313            x: Error indicator values.
314            y: Corresponding QoI error values.
315        '''
316        if self.__calibrated_error:
317            self.__c = np.mean(y) / np.mean(x)

If calibrated_error was specified, this function fits the regression model to the provided data points (x, y) to estimate the scaling factor.

Arguments:

x: Error indicator values.
y: Corresponding QoI error values.

def predict(self, x): View Source

319    def predict(self, x):
320        '''
321        Predicts QoI error based on the error indicator using the estimated
322        scaling factor.
323
324        Args:
325            x: Error indicator value(s) for prediction.
326        '''
327        return self.__c*x

Predicts QoI error based on the error indicator using the estimated scaling factor.

Arguments:

x: Error indicator value(s) for prediction.