diff options
Diffstat (limited to 'inst/examples/02_design_of_experiments/stk_example_doe05.m')
| -rw-r--r-- | inst/examples/02_design_of_experiments/stk_example_doe05.m | 47 |
1 files changed, 33 insertions, 14 deletions
diff --git a/inst/examples/02_design_of_experiments/stk_example_doe05.m b/inst/examples/02_design_of_experiments/stk_example_doe05.m index 3ea5500..84be3e2 100644 --- a/inst/examples/02_design_of_experiments/stk_example_doe05.m +++ b/inst/examples/02_design_of_experiments/stk_example_doe05.m @@ -11,7 +11,7 @@ % Copyright Notice % -% Copyright (C) 2015-2017 CentraleSupelec +% Copyright (C) 2015-2017, 2020 CentraleSupelec % Copyright (C) 2016 EDF R&D % Copyright (C) 2013, 2014 SUPELEC % @@ -23,7 +23,7 @@ % This file is part of % % STK: a Small (Matlab/Octave) Toolbox for Kriging -% (http://sourceforge.net/projects/kriging) +% (https://github.com/stk-kriging/stk/) % % STK is free software: you can redistribute it and/or modify it under % the terms of the GNU General Public License as published by the Free @@ -38,8 +38,7 @@ % You should have received a copy of the GNU General Public License % along with STK. If not, see <http://www.gnu.org/licenses/>. -stk_disp_examplewelcome; -stk_figure ('stk_example_eqi_test'); +stk_disp_examplewelcome; stk_figure ('stk_example_doe05'); %% Problem definition @@ -47,7 +46,8 @@ stk_figure ('stk_example_eqi_test'); % The goal is to find the minimum of f on the domain BOX. % 1D test function -f = @(x)(x .* sin (x)); % Define a 1D test function +f = @(x)(x .* sin (x)); % Define a 1D test function... +% f = @(x)((x - 6) .^ 2) % ...or another one? DIM = 1; % Dimension of the factor space BOX = stk_hrect ([0; 12], {'x'}); % Factor space (hyper-rectangle object) @@ -58,7 +58,7 @@ NOISE_VARIANCE = 2 ^ 2; GRID_SIZE = 200; % Number of points in the grid xg = stk_sampling_regulargrid (GRID_SIZE, DIM, BOX); -% Give names explicit names to the points of the grid +% Give explicit names to the points of the grid xg.rownames = arrayfun ... (@(i)(sprintf ('xg(%03d)', i)), 1:GRID_SIZE, 'UniformOutput', false)'; @@ -74,7 +74,7 @@ BUDGET = 100; % Total evaluation budget REESTIM_PERIOD = 10; % How often should we re-estimate the cov parameters ? SAMPCRIT_NAME = 'AKG'; % Choose a particular sampling criterion -% Note: the two criteria proposed here compute an "expected improvement" of +% Note: the three criteria proposed here compute an "expected improvement" of % some kind. As such, they return positive values, and must be maximized. switch SAMPCRIT_NAME @@ -83,7 +83,14 @@ switch SAMPCRIT_NAME POINT_BATCH_SIZE = @(x, n) BUDGET - n; sampcrit = stk_sampcrit_eqi ([], QUANTILE_ORDER, POINT_BATCH_SIZE); case 'AKG' + % Use the "approximate KG" criterion, with Scott's original method + % for the construction of the reference grid (i.e., taking past + % observations points plus the candidate point as the reference grid). sampcrit = stk_sampcrit_akg (); + case 'KG' + % Use the "exact KG" criterion, with can be obtained by taking the + % reference grid equal to the entire input grid + sampcrit = stk_sampcrit_akg ([], xg); end @@ -92,7 +99,7 @@ end % Construction of the initial design x0 = stk_sampling_regulargrid (N0, DIM, BOX); -% Give names explicit names to the points in the initial design +% Give explicit names to the points in the initial design x0.rownames = arrayfun ... (@(i)(sprintf ('init%03d', i)), 1:N0, 'UniformOutput', false)'; @@ -103,7 +110,7 @@ z0 = z0 + sqrt (NOISE_VARIANCE) * randn (size (z0)); %% Specification of the model (Gaussian process prior) -model0 = stk_model ('stk_materncov52_iso'); +model0 = stk_model (@stk_materncov52_iso); % Assume that the variance of the observation noise is known model0.lognoisevariance = log (NOISE_VARIANCE); @@ -139,22 +146,28 @@ for iter = 1:(BUDGET - N0) % Instanciate sampling criterion object with model sampcrit.model = model; - % Compute the EQI criterion on the grid + % Compute the sampling criterion on the grid [crit_val, z_pred] = sampcrit (xg); if mod (iter, PLOT_PERIOD) == 1 + % Figure: upper panel - stk_subplot (2, 1, 1); hold off; % CG#12 + stk_subplot (2, 1, 1); cla; stk_plot1d ([],[], xg, zg, z_pred); hold on; plot (x, z, 'k.'); + title (sprintf ('n = %d + %d = %d.\n\n', N0, iter - 1, N0 + iter - 1)); + % Figure: lower panel stk_subplot (2, 1, 2); cla; plot (xg, crit_val); xlim (BOX); ylabel (SAMPCRIT_NAME); + end - if all (crit_val == 0), break, end + % Stop if the criterion becomes equal to zero everywhere + if all (crit_val == 0), break; end - % Pick the point where the EQI is maximum as our next evaluation point + % Pick the point where the sampling criterion is maximum + % as our next evaluation point [crit_max, i_max] = max (crit_val); x_new = xg(i_max, :); @@ -177,8 +190,14 @@ end % Display the final DoE data = stk_dataframe ([x z], {'x', 'z'}); disp (data); +% Premature stopping ? +if iter < (BUDGET - N0) + warning ('The algorithm stopped prematurely.'); + iter = iter - 1; +end + % Total number of evaluations ? -fprintf ('\nNumber of evaluations: %d + %d = %d.\n\n', N0, BUDGET - N0, BUDGET); +fprintf ('\nNumber of evaluations: %d + %d = %d.\n\n', N0, iter, N0 + iter); %#ok<*AGROW> |