rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
scm_infsup_greedy.m
Go to the documentation of this file.
1 function SCM_data=scm_infsup_greedy(model, params, operator, innerproduct_matrix)
2 %model=scm_greedy(model, model_data)
3 %
4 % SCM greedy function
5 % computes necessary values for the calculation of the infsup constant
6 % requires model and model data structure
7 %
8 %
9 % SCM Properties:
10 %
11 %
12 % SCM_data.Ma : number of parameters with an already
13 % computed constant near to mu
14 % SCM_data.Mp : number of parameters without computed
15 % constant near to mu
16 % SCM_data.NumberOfParameters : size of the test set used in the SCM
17 % Greedy procedure
18 % SCM_data.eps : error boundary of the Greedy procedure
19 % SCM_data.max_iterations : maximal number of Greedy iterations
20 % SCM_data.parallel : = 0 non parallel version
21 % = 1 parallel (pool stays active after
22 % computaton)
23 % = 2 parallel (closing pool afterwards)
24 % David Kreplin 05.01.2015
25 
26 if isempty(params)
27  warning('No SCM properties are set -> using standart values')
28  %SCM standart values
29  params.Ma = inf;
30  params.Mp = 0;
31  params.NumberOfParameters = 100;
32  params.eps=0.1;
33  params.max_iterations=inf;
34  params.parallel=0;
35 end
36 
37 if ~isfield(params, 'Ma')
38  params.Ma = inf;
39 end
40 
41 if ~isfield(params, 'Mp')
42  params.Mp = 0;
43 end
44 
45 if ~isfield(params, 'NumberOfParameters')
46  warning('SCM: No Number of Parameters is set -> setting NumberOfParameters = 100');
47  params.NumberOfParameters = 100;
48 end
49 
50 if ~isfield(params, 'eps')
51  params.eps = 0.1;
52 end
53 
54 if ~isfield(params, 'max_iterations')
55  params.max_iterations = inf;
56 end
57 
58 if ~isfield(params, 'parallel')
59  params.parallel = 0;
60 end
61 
62 if ~isfield(params, 'para_rand_seed')
63  params.para_rand_seed = 1234;
64 end
65 
66 
67 %% Parameter set
68 if ~isfield(params, 'ParameterSet')
69  rand('state',params.para_rand_seed);
70  params.ParameterSet = rand_uniform(params.NumberOfParameters,model.mu_ranges);
71 end
72 NumberOfParameters = length(params.ParameterSet(1,:));
73 
74 SCM_data.params = params;
75 
76 ind =@ (i,j,q) j + (i-1)*q - i*(i-1)/2;
77 
78 
79 
80 
81 %% Boundary BOX for YLB
82 model.decomp_mode=1;
83 [A,~] = operator(model);
84 A=remove_dirichlet(A,1,SCM_data.params.dirichlet_gids );
85 
86 Q = length(A);
87 SCM_data.params.Q=Q;
88 SCM_data.lower_bound_infsup = zeros(Q*(Q+1)/2,1);
89 SCM_data.upper_bound_infsup = zeros(Q*(Q+1)/2,1);
90 
91 disp('-- SCM: computing boundary box');
92 %remove dirchlet
93 %innerproduct_matrix = model.get_inner_product_matrix(model_data);
94 innerproduct_matrix(SCM_data.params.dirichlet_gids ,:)=[];
95 innerproduct_matrix(:,SCM_data.params.dirichlet_gids )=[];
96 innerproduct_matrix = 0.5*(innerproduct_matrix + innerproduct_matrix');
97 for i = 1:Q
98  for j = i:Q
99 
100  fun = @(x) A{1,i}'*(innerproduct_matrix\(A{1,j}*x));
101 
102  [lambda1, success] = safeEig(fun, size(A{i}, 1), innerproduct_matrix, 'lr');
103  if success
104  [lambda2, success] = safeEig(fun, size(A{i}, 1), innerproduct_matrix, 'sr');
105  end
106 
107  if ~success
108  warning('scm_infsup_greedy: eigs did not converge, trying eig');
109  spectrum = eig(full(A{1,i}'*(innerproduct_matrix'\A{1,j})), full(innerproduct_matrix));
110  if any(abs(imag(spectrum)) > 1e-9)
111  warning('scm_infsup_greedy: ignoring complex eigenvalues');
112  end
113  spectrum = real(spectrum);
114  lambda1 = max(spectrum);
115  lambda2 = min(spectrum);
116  end
117 
118  SCM_data.upper_bound_infsup(ind(i,j,Q)) = lambda1;
119  SCM_data.lower_bound_infsup(ind(i,j,Q)) = lambda2;
120  end;
121 end
122 
123 % SCM_data.lower_bound_infsup=[0;0;-7.3747e-16];
124 % SCM_data.upper_bound_infsup=[1;0;1.267];
125 
126 %keyboard;
127 %% Initalize C
128 SCM_data.C=[];
129 SCM_data.C=AddC(params.ParameterSet(:,1),model,operator,innerproduct_matrix,SCM_data);
130 SCM_data.params.ParameterSet(:,1)=[];
131 NumberOfParameters = NumberOfParameters - 1;
132 % SCM Greedy
133 n=0;
134 eta=params.eps+1;
135 
136 %% SCM Greedy - parallel version
137 if SCM_data.params.parallel >= 1
138  % s = matlabpool('size'); %OLDER MATLAB VERSION
139  % if s == 0
140  % matlabpool('open');
141  % end;
142  if isempty(gcp('nocreate'))
143  parpool;
144  end;
145  %reduce parallel data amount
146  ParameterSet = SCM_data.params.ParameterSet;
147  while SCM_data.params.eps <= eta
148  n=n+1;
149  eta_temp=zeros(NumberOfParameters,1);
150  parfor i = 1:NumberOfParameters
151  alpha_UB = scm_infsup_ub(ParameterSet(:,i),model,operator,SCM_data);
152  alpha_LB = scm_infsup_lb(ParameterSet(:,i),model,operator,SCM_data);
153  eta_temp(i) = (sqrt(alpha_UB)-sqrt(alpha_LB))/sqrt(alpha_UB);
154  end
155 
156  [eta,index]=max(eta_temp);
157 
158  SCM_data.C=AddC(ParameterSet(:,index),model,operator,innerproduct_matrix,SCM_data);
159  ParameterSet(:,index) = [];
160  NumberOfParameters = NumberOfParameters-1;
161 
162  disp(['-- SCM: step ' num2str(n) ' - error: ' num2str(eta)]);
163  if n >= SCM_data.params.max_iterations
164  break
165  end;
166  end
167  SCM_data.params.ParameterSet=ParameterSet;
168  if params.parallel == 2
169  % matlabpool('close'); %OLDER MATLAB VERSION
170  delete(gcp('nocreate'))
171  end
172 else
173 
174  %% non-parallel version
175  while SCM_data.params.eps <= eta
176  n=n+1;
177  eta_temp=zeros(NumberOfParameters,1);
178  for i = 1:NumberOfParameters
179  alpha_UB = scm_infsup_ub(SCM_data.params.ParameterSet(:,i),model,operator,SCM_data);
180  alpha_LB = scm_infsup_lb(SCM_data.params.ParameterSet(:,i),model,operator,SCM_data);
181  eta_temp(i) = abs(sqrt(alpha_UB)-sqrt(alpha_LB))/sqrt(alpha_UB);
182  end
183  [eta,index]=max(eta_temp);
184 
185  SCM_data.C=AddC(SCM_data.params.ParameterSet(:,index),model,operator,innerproduct_matrix,SCM_data);
186 
187  SCM_data.params.ParameterSet(:,index) = [];
188  NumberOfParameters = NumberOfParameters-1;
189 
190  disp(['-- SCM: step ' num2str(n) ' - error: ' num2str(eta)]);
191  if n >= SCM_data.params.max_iterations
192  break
193  end;
194  end
195 
196 end
197 
198 
199 
200 end
201 
202 %%
203 function [lambda, success] = safeEig(fun, sizeA, B, sigma)
204 
205 lambda = 0;
206 success = 0;
207 opts.tol = 1e-8;
208 
209 while ~success && opts.tol < 1e-2
210 
211  try
212  lambda = eigs(fun, sizeA, B, 1, sigma, opts);
213  success = 1;
214  catch
215  warning('safeEig: eigs exception, lowering the tolerance');
216  opts.tol = 100*opts.tol;
217  end
218 end
219 
220 if abs(imag(lambda)) > 1e-9
221  warning('safeEig: ignoring complex eigenvalues');
222 end
223 lambda = real(lambda);
224 
225 end
226 
227 
228 function C=AddC(mu,model,operator,innerproduct_matrix,SCM_data)
229 ind =@ (i,j,q) j + (i-1)*q - i*(i-1)/2;
230 % innerproduct_matrix = model.get_inner_product_matrix(model_data);
231 % innerproduct_matrix(SCM_data.params.dirichlet_gids ,:)=[];
232 % innerproduct_matrix(:,SCM_data.params.dirichlet_gids )=[];
233 
234 if isempty(SCM_data.C)
235  C.Set = mu;
236  model = set_mu(model,mu);
237  model.decomp_mode=0;
238  [A,~] = operator(model);
239  A=remove_dirichlet(A,0,SCM_data.params.dirichlet_gids);
240  %[V,D] = eig(full(A')*inv(full(innerproduct_matrix))'*full(A),full(innerproduct_matrix));
241  %[v,d] = fast_eigs(A'*(innerproduct_matrix'\A),innerproduct_matrix,'sm' );
242  A = 0.5*(A+A');
243  fun = @(X) A'\(innerproduct_matrix*(A\X));
244  [v,d] = eigs(fun, size(A, 1), innerproduct_matrix, 1, 'sm');
245  % a = min(max(D));
246  C.infsup =d; %infsup;
247 
248  %v=V(:,i);
249  model.decomp_mode=1;
250  [A_Matrix,~] = operator(model);
251  A_Matrix=remove_dirichlet(A_Matrix,1,SCM_data.params.dirichlet_gids );
252  y_UB=zeros(SCM_data.params.Q*(SCM_data.params.Q+1)/2,1);
253  norm_value = v'*innerproduct_matrix*v;
254  for i=1:SCM_data.params.Q
255  for j=i:SCM_data.params.Q
256  y_UB(ind(i,j,SCM_data.params.Q)) = v'*A_Matrix{1,i}'*(innerproduct_matrix'\(A_Matrix{1,j}*v))/norm_value;
257  end
258  end
259 
260 
261  model.decomp_mode=2;
262  [A,~] = operator(model);
263  ParameterVector=zeros(0,(SCM_data.params.Q)*(SCM_data.params.Q+1)/2);
264  for i=1:SCM_data.params.Q
265  for j=i:SCM_data.params.Q
266  ParameterVector(ind(i,j,SCM_data.params.Q)) = -(2-(i==j))*(A(i)*A(j));
267  end
268  end
269  C.ParameterVector = ParameterVector;
270 
271  %Parameter for lower bound computation
272 
273  C.UBVector = y_UB;
274 else
275  C.Set = [SCM_data.C.Set mu];
276  model = set_mu(model,mu);
277  model.decomp_mode=0;
278  [A,~] = operator(model);
279  A=remove_dirichlet(A,0,SCM_data.params.dirichlet_gids);
280  % [V,D] = eig(full(A')*inv(full(innerproduct_matrix))'*full(A),full(innerproduct_matrix));
281  %[v,d] = fast_eigs(A'*(innerproduct_matrix'\A),innerproduct_matrix,'sm' );
282  A = 0.5*(A+A');
283  fun = @(X) A'\(innerproduct_matrix*(A\X));
284  [v,d] = eigs(fun, size(A, 1), innerproduct_matrix, 1, 'sm');
285  % a = min(max(D));
286  %C.infsup = [C.infsup;infsup];
287  C.infsup = [SCM_data.C.infsup;d];
288 
289  %v=V(:,i);
290  model.decomp_mode=1;
291  [A_Matrix,~] = operator(model);
292  A_Matrix=remove_dirichlet(A_Matrix,1,SCM_data.params.dirichlet_gids);
293  y_UB=zeros(SCM_data.params.Q*(SCM_data.params.Q+1)/2,1);
294  norm_value = v'*innerproduct_matrix*v;
295  for i=1:SCM_data.params.Q
296  for j=i:SCM_data.params.Q
297  y_UB(ind(i,j,SCM_data.params.Q)) = v'*A_Matrix{1,i}'*(innerproduct_matrix'\(A_Matrix{1,j}*v))/norm_value;
298  end
299  end
300 
301  model.decomp_mode=2;
302  [A,~] = operator(model);
303  ParameterVector=zeros(0,(SCM_data.params.Q)*(SCM_data.params.Q+1)/2);
304  for i=1:SCM_data.params.Q
305  for j=i:SCM_data.params.Q
306  ParameterVector(ind(i,j,SCM_data.params.Q)) = -(2-(i==j))*(A(i)*A(j));
307  end
308  end
309  C.ParameterVector = [SCM_data.C.ParameterVector;ParameterVector];
310 
311  C.UBVector = [SCM_data.C.UBVector y_UB];
312 
313 end
314 
315 end
316 
317 %% Upperbound computation
318 
319 function alpha_UB = scm_infsup_ub(mu,model,operator,SCM_data)
320 ind =@ (i,j,q) j + (i-1)*q - i*(i-1)/2;
321 
322 model = set_mu(model,mu);
323 model.decomp_mode=2;
324 [A,~] = operator(model);
325 factor= zeros(1,(SCM_data.params.Q)*(SCM_data.params.Q+1)/2);
326 
327 for i=1:SCM_data.params.Q
328  for j=i:SCM_data.params.Q
329  factor(ind(i,j,SCM_data.params.Q)) = (2-(i==j))*(A(i)*A(j));
330  end
331 end
332 
333 alpha_UB=min(factor*SCM_data.C.UBVector);
334 
335 end
336 
337 
338 
339 
340 
341 function A=remove_dirichlet(A,mode,I)
342 if mode == 0
343  A(I,:) = [];
344  A(:,I)= [];
345 else
346  for i = 1:length(A)
347  A{1,i}(I,:) = [];
348  A{1,i}(:,I) = [];
349  end
350 end
351 end
Customizable implementation of an abstract greedy algorithm.
Definition: DuneRBLeafNode.m:1
function SCM_data = scm_infsup_greedy(model, params, operator, innerproduct_matrix)
model=scm_greedy(model, model_data)