-
Notifications
You must be signed in to change notification settings - Fork 1
/
script_Gauss_Genz.m
174 lines (152 loc) · 6.09 KB
/
script_Gauss_Genz.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
%% script_Gauss_Genz
%
% Description:
% Script to perform a stability and error analysis of the Gaussian RBF for Genz' test functions
%
% Author: Jan Glaubitz
% Date: June 22, 2021
clc, clear
%warning('off','all')
%% Free parameters
N = 20^2; % number of data points
CC = 100; % number of tests for Genz
noise_level = 0; % amount of uniform noise (0 means no noise, while a>0 mean 10^(-a))
points = 'equid'; % type of data points (equid, Halton, random)
kernel = 'G'; % kernel (G, MQ, IQ, Wendland, TPS, cubic)
order = 2; % order (for Wendland function)
d = 0; % polynomial degree
%% prepare script
% dimension and precision
a = 0; b = 1; dim = 2; % domain is [0,1]^2
precision = 32; % use usual double precision
% domain, data points, and RBb
X = generate_points( dim, a, b, N, points); % generate data points
DM = DistanceMatrix( X, X ); % matrix containing the distances between points
h = min( DM + 42*eye(N), [], 'all' ); % measure of distance between data points
rbf = initialize_RBF( kernel, dim, order ); % initialize RBF
% set up shape parameters
mid = ceil((1/h)^(1/10)); % exponent corresponding to h
exp_start = mid-6; exp_end = mid+1; n = 100;
Ep = logspace(exp_start,exp_end,n)'; % generates a vector of n logarithmically spaced points
% values of interest
cond_nr = []; % condition number
opt = []; % optimal values
s = []; % stability values
err = []; % errors
for i=1:n
%% Update shape parameter
ep = Ep(i); % vector of shape parameters
[i, n]
%% Compute the condition number
cond_number = Cond( a, b, rbf, ep, X, d );
%% Compute moments of the RBFs
m_RBF = RBF_moments( a, b, kernel, rbf, ep, X );
%% Compute RBF-CF weights (in a different manner)
w = compute_weights( a, b, rbf, ep, X, m_RBF, d, precision );
%% Compute the stability measure s_N and otimal value C_N[1]
if d < 0
opt_value = sum(w); % optimal value
else
opt_value = (b-a)^2; % optimal value
end
stab_measure = sum(abs(w)); % stability measure
%% Genz test functions
error = Genz( a, b, X, w, CC, noise_level ); % average errors
%% Store values
cond_nr = [cond_nr; cond_number]; % condition number
opt = [opt; opt_value]; % optimal values
s = [s; stab_measure]; % stability values
err = [err; error]; % error of CF
end
%format shortE
%[min_error,index] = min( err(:,1) );
%[p_count, d, 1, min_error, Ep(index), s(index)]
%[min_error,index] = min( err(:,4) );
%[p_count, d, 4, min_error, Ep(index), s(index)]
%% Illustrate results - Genz 1
figure(1)
p = plot( Ep,s,'b-.', Ep,opt,'k:', Ep,err(:,1),'r-');
set(p, 'LineWidth',3.5)
set(gca, 'FontSize', 26) % Increasing ticks fontsize
xlim([ Ep(1); Ep(end) ])
xlabel('$\varepsilon$','Interpreter','latex')
xticks(10.^[exp_start+1 0.5*(exp_start+1+exp_end) exp_end]);
set(gca, 'XScale', 'log')
set(gca, 'YScale', 'log')
if d >= 0
lgnd = legend('$\|C_N\|_{\infty}$','$\|I\|_{\infty}$','$| C_N[g_1] - I[g_1] |$','Interpreter','latex','Location','best');
else
lgnd = legend('$\|C_N\|_{\infty}$','$C_N[1]$','$| C_N[g_1] - I[g_1] |$','Interpreter','latex','Location','best');
end
set(lgnd, 'Interpreter','latex', 'FontSize',26, 'color','none')
grid on
%% Illustrate results - Genz 2
figure(2)
p = plot( Ep,s,'b-.', Ep,opt,'k:', Ep,err(:,2),'r-');
set(p, 'LineWidth',3.5)
set(gca, 'FontSize', 26) % Increasing ticks fontsize
xlim([ Ep(1); Ep(end) ])
xlabel('$\varepsilon$','Interpreter','latex')
xticks(10.^[exp_start+1 0.5*(exp_start+exp_end+1) exp_end]);
set(gca, 'XScale', 'log')
set(gca, 'YScale', 'log')
if d >= 0
lgnd = legend('$\|C_N\|_{\infty}$','$\|I\|_{\infty}$','$| C_N[g_2] - I[g_2] |$','Interpreter','latex','Location','best');
else
lgnd = legend('$\|C_N\|_{\infty}$','$C_N[1]$','$| C_N[g_2] - I[g_2] |$','Interpreter','latex','Location','best');
end
set(lgnd, 'Interpreter','latex', 'FontSize',26, 'color','none')
grid on
%% Illustrate results - Genz 3
figure(3)
p = plot( Ep,s,'b-.', Ep,opt,'k:', Ep,err(:,3),'r-');
set(p, 'LineWidth',3.5)
set(gca, 'FontSize', 26) % Increasing ticks fontsize
xlim([ Ep(1); Ep(end) ])
xlabel('$\varepsilon$','Interpreter','latex')
xticks(10.^[exp_start+1 0.5*(exp_start+exp_end+1) exp_end]);
set(gca, 'XScale', 'log')
set(gca, 'YScale', 'log')
if d >= 0
lgnd = legend('$\|C_N\|_{\infty}$','$\|I\|_{\infty}$','$| C_N[g_3] - I[g_3] |$','Interpreter','latex','Location','best');
else
lgnd = legend('$\|C_N\|_{\infty}$','$C_N[1]$','$| C_N[g_3] - I[g_3] |$','Interpreter','latex','Location','best');
end
set(lgnd, 'Interpreter','latex', 'FontSize',26, 'color','none')
grid on
%% Illustrate results - Genz 4
figure(4)
p = plot( Ep,s,'b-.', Ep,opt,'k:', Ep,err(:,4),'r-');
set(p, 'LineWidth',3.5)
set(gca, 'FontSize', 26) % Increasing ticks fontsize
xlim([ Ep(1); Ep(end) ])
xlabel('$\varepsilon$','Interpreter','latex')
xticks(10.^[exp_start+1 0.5*(exp_start+exp_end+1) exp_end]);
set(gca, 'XScale', 'log')
set(gca, 'YScale', 'log')
if d >= 0
lgnd = legend('$\|C_N\|_{\infty}$','$\|I\|_{\infty}$','$| C_N[g_4] - I[g_4] |$','Interpreter','latex','Location','best');
else
lgnd = legend('$\|C_N\|_{\infty}$','$C_N[1]$','$| C_N[g_4] - I[g_4] |$','Interpreter','latex','Location','best');
end
set(lgnd, 'Interpreter','latex', 'FontSize',26, 'color','none')
grid on
%% Illustrate results - Genz 1 & 4 combined
figure(5)
p = plot( Ep,s,'b-.', Ep,err(:,1), 'k--', Ep,err(:,4),'r-');
line2 = yline(opt(1),':','$\|I\|_{\infty}$','Interpreter','latex');
set(line2, 'LineWidth',3.5, 'FontSize',26 );
set(p, 'LineWidth',3.5)
set(gca, 'FontSize', 26) % Increasing ticks fontsize
xlim([ Ep(1); Ep(end) ])
xlabel('$\varepsilon$','Interpreter','latex')
xticks(10.^[exp_start+1 0.5*(exp_start+exp_end+1) exp_end]);
set(gca, 'XScale', 'log')
set(gca, 'YScale', 'log')
if d >= 0
lgnd = legend('$\|C_N\|_{\infty}$','$| C_N[g_1] - I[g_1] |$','$| C_N[g_4] - I[g_4] |$','Interpreter','latex','Location','southwest');
else
lgnd = legend('$\|C_N\|_{\infty}$','$| C_N[g_1] - I[g_1] |$','$| C_N[g_4] - I[g_4] |$','Interpreter','latex','Location','southwest');
end
set(lgnd, 'Interpreter','latex', 'FontSize',26, 'color','none')
grid on