-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathscript_PHS_Genz.m
166 lines (146 loc) · 5.58 KB
/
script_PHS_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
%% script_PHS_Genz
%
% Description:
% Script to perform a stability and error analysis of PHS for Genz' test functions
%
% Author: Jan Glaubitz
% Date: June 22, 2021
clc, clear
%% Free parameters
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 = 'TPS'; % kernel (G, MQ, IQ, Wendland, TPS, cubic, quintic)
ep = 1; % shape parameter
order = 2; % order (for Wendland function)
d = 1; % polynomial degree
%% prepare script
% dimension and precision
a = 0; b = 1; dim = 2; % domain is [0,1]^2
precision = 32; % use usual double precision
NN = (2:10:62).^2;
% values of interest
cond_nr = []; % condition number
opt = []; % optimal values
s = []; % stability values
err = []; % errors
for i=1:length(NN)
%% Update data points
N = NN(i); % number of data points
[i, length(NN)]
% data points, and RBb
X = generate_points( dim, a, b, N, points); % generate data points
rbf = initialize_RBF( kernel, dim, order ); % initialize RBF
%% 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 );
%m_RBF2 = RBF_moments( a, b, 'numint', 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
%% Illustrate results - Genz 1
figure(1)
p = plot( NN,s,'b-.', NN,opt,'k:', NN,err(:,1),'r-');
set(p, 'LineWidth',3.5)
set(gca, 'FontSize', 26) % Increasing ticks fontsize
xlim([ NN(1); NN(end) ])
xlabel('$N$','Interpreter','latex')
xticks(10.^[1 2 3]);
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( NN,s,'b-.', NN,opt,'k:', NN,err(:,2),'r-');
set(p, 'LineWidth',3.5)
set(gca, 'FontSize', 26) % Increasing ticks fontsize
xlim([ NN(1); NN(end) ])
xlabel('$N$','Interpreter','latex')
xticks(10.^[1 2 3 4]);
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( NN,s,'b-.', NN,opt,'k:', NN,err(:,3),'r-');
set(p, 'LineWidth',3.5)
set(gca, 'FontSize', 26) % Increasing ticks fontsize
xlim([ NN(1); NN(end) ])
xlabel('$N$','Interpreter','latex')
xticks(10.^[1 2 3 4]);
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( NN,s,'b-.', NN,opt,'k:', NN,err(:,4),'r-');
set(p, 'LineWidth',3.5)
set(gca, 'FontSize', 26) % Increasing ticks fontsize
xlim([ NN(1); NN(end) ])
xlabel('$N$','Interpreter','latex')
xticks(10.^[1 2 3 4]);
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( NN,s,'b-.', NN,err(:,1), 'k--', NN,err(:,4),'r-');
set(p, 'LineWidth',3.5)
% yline for optimal stability value
aux_line = yline(opt(1),':','$\|I\|_{\infty}$','Interpreter','latex');
set(aux_line, 'LineWidth',3.5, 'FontSize',26 );
aux_line.LabelVerticalAlignment = 'bottom';
lgnd = legend('$\|C_N\|_{\infty}$','$| C_N[g_1] - I[g_1] |$','$| C_N[g_4] - I[g_4] |$','Interpreter','latex','Location','southwest');
% axes
set(gca, 'FontSize', 26) % Increasing ticks fontsize
xlim([ NN(1); NN(end) ])
xlabel('$N$','Interpreter','latex')
xticks(10.^[1 2 3 4]);
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