rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
follicle_rect_experiments.m
Go to the documentation of this file.
1 function res = follicle_rect_experiments(step,model,model_data,detailed_data)
2 %function follicle_rect_experiments(step[,model,model_data,detailed_data])
3 % experiments with the model of the human follicle growth
4 %
5 % model range is n=150, noutput = 50, nt = 500, time = 19.02 years
6 % implicit Euler time-discretization, accept relatively large error at
7 % initial time, for large overall time interval and few timesteps
8 % and "interactivity"
9 %
10 % For the results of the paper, simply run step 3.4, 5, 6, 7, 8
11 % in this order.
12 %
13 % (step = 1 simulation and visualization of detailed dynamical
14 % system, using filecaching => do not call directly)
15 % step = 1.1, n=150, initial data distributed
16 % over larger domain init_range = 20
17 % and different output functional 'expectation_M'
18 % n smaller, e.g. 120: decay of output!!!!
19 % (step = 2 reduced simulation and visualization of approximation
20 % => do not call directly)
21 % (step = 3 detailed and reduced simulation and visualization of
22 % error and estimator
23 % => do not call directly)
24 % step = 3.1, n=150, n_output = 50, lower states as reduced basis
25 % expected large error, but test of functionality.
26 % step = 3.2, n=150, n_output = 50, generate POD basis of
27 % single trajectory
28 % basis is saved for use in further steps
29 % and plotted
30 % step = 3.3, n=150, n_output , generate BT basis of
31 % single trajectory
32 % basis is saved for use in further steps
33 % and plotted
34 % step = 3.4, n=150, generate PODGreedy basis of
35 % 81 trajectories
36 % basis is saved for use in further steps
37 % and plotted
38 % step = 3.5 experiments with a part of POD-basis from 3.4
39 % step = 3.8 generate system matrices and export to file for
40 % remote BT generation. System A = 20301x20301
41 % => ss = sys(A,B,C,D) => 3 GB system ss!!!!!
42 % balreal takes about ...
43 % step 4.0: demo_rb_gui of POD basis from 3.4
44 % step 5.0: Optimization with RB model with POD basis from 3.4
45 % generation of optimization plots for initial submission
46 % of paper
47 % step 5.1: Optimization with RB model with POD basis from 3.4
48 % generation of optimization plots for revised paper
49 % step 6.0: Generation of trajectory plots for paper
50 % step 7.0: Generation of basis vector plots for paper
51 % step 8.0: Generation of POD-Greedy plots for paper
52 %
53 % some steps explicitly require a model and model_data
54 %
55 % Todo:
56 % generate BT basis
57 % compare quality of BT to POD on single trajectory
58 % compare quality of BT to POD on parameter corners
59 
60 % Bernard Haasdonk 21.4.2010
61 
62 res = [];
63 
64 if (nargin<1) || (isempty(step))
65  step = 1.1 % simulation and visualization of detailed dynamical system
66  %step = 2 % reduced simulation and visualization of approximation
67  %step =3 % detailed and reduced simulation and visualization of
68  % error and estimator
69 end;
70 
71 if (nargin < 2)
72  model = [];
73 end;
74 
75 if (nargin < 3)
76  model_data = [];
77 % model = follicle_rect_model([]);
78 % model_data = gen_model_data(model);
79 % model = follicle_rect_model([]);
80 % model_data = gen_model_data(model);
81 % model_data.RB = eye(model.dim_x,231); % i.e. first 231 states
82 % model.RB_generation_mode = 'from_detailed_data';
83 end;
84 
85 if (nargin < 4)
86  detailed_data = [];
87 end;
88 
89 switch step
90  case 1 % simulation and visualization
91  tic;
92  if isempty(model_data)
93  model_data = gen_model_data(model);
94  end;
95  sim_data = filecache_function(@detailed_simulation,model,model_data);
96 % sim_data = detailed_simulation(model,model_data);
97  t = toc;
98  mu = model.get_mu(model);
99  model.plot_title = ['state probabilities for mu=(',...
100  num2str(mu(1)),',',num2str(mu(2)),')'];
101  p = plot_sim_data(model,model_data,sim_data,model);
102 % model.plot_sim_data_state(model,model_data,sim_data,model);
103  %compute sum for plausibility check of time range/dimension
104 % l1norm = sum(sim_data.X);
105 % figure, plot(l1norm); xlabel('time'); title('l1norm of X over time');
106  fprintf('\n')
107  disp(['time for detailed_simulation :',num2str(t)]);
108  if model.debug
109  disp(['finished, please interact with data']);
110  keyboard;
111  end;
112  res = sim_data;
113 
114  case 1.1
115  tic;
116  params.range = 150;
117 % params.output_range = 50;
118  % params.output_plot_indices = [1,2];
119  params.init_range = 20;
120 % params.output_type = 'separatrix';
121  params.output_type = 'expectation_N';
122 % params.reg_epsilon = 0.00009;
123  params.reg_epsilon = 0.00000;
124  model = follicle_rect_model(params);
125  model.clim = [0,0.00005];
126  % model_data = gen_model_data(model);
127  % mu = (u1, u2), mu_ranges : {[0.005,0.02],[0.005,0.02]};
128 % mu = [0.005,0.005]; % => blob to 5e-5, yout to 3.7
129 % mu = [0.01,0.01]; % => blob to 5e-5, yout to 3.8
130 % mu = [0.02,0.02]; % => blob to 0, yout to almost 3.7,
131 % mu = [0.005,0.02]; % => blob to 0, yout to almost 9.2
132  mu = [0.02,0.005]; % => blob non existing, yout 0.3
133 
134  model = model.set_mu(model,mu);
135  % model.save_time_indices = 0:(model.nt/100):model.nt;
136  %model.save_time_indices = 1:100:1001;
137  %model.nt = 1000;
138  %model.T = 1000*model.nt;
139 % sim_data = follicle_rect_experiments(1,model,model_data);
140  sim_data = follicle_rect_experiments(1,model);
141  save(fullfile(rbmatlabtemp,'follicle_n200_example.mat'),...
142  'sim_data','model','model_data');
143  t = toc;
144  disp(['computation time t = ',num2str(t),' sec.']);
145  keyboard;
146 
147  case 2 % reduced simulation
148 
149  model.debug = 1;
150 
151  % produce high-dimensional data, e.g. reduced basis, etc.
152  % mu is not known
153  detailed_data = gen_detailed_data(model,model_data);
154  disp('gen_detailed_data successful');
155 
156  % produce low-dimensional data, e.g. gram matrices between
157  % reduced basis vectors, subgrid extraction, etc.
158  % mu is not known (should comprise previous offline and online steps)
159 
160  reduced_data = gen_reduced_data(model,detailed_data);
161  disp('gen_reduced_data successful');
162  model.N = reduced_data.N;
163 
164  % rb_sim_data = rb_simulation(reduced_data,model)
165  % perform reduced simulation, i.e. mu known
166  % produces some structure with suitable fields, e.g DOF vector
167 
168  disp('entering reduced simulation');
169  rb_sim_data = rb_simulation(model,reduced_data);
170  disp('reduced simulation successful');
171 
172  rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
173 
174  model.plot_title = 'reduced trajectory';
175  %model.plot_indices = 1:3; => default
176  p = plot_sim_data(model,model_data,rb_sim_data,[]);
177 
178  disp(['finished, please interact with data']);
179  keyboard;
180 
181  case 2.1
182 % par.range = 150;
183 % par.output_range = 20;
184  model = follicle_rect_model;
185  model_data = gen_model_data(model);
186  model_data.RB = eye(model.dim_x,231); % i.e. first 231 states
187  model.RB_generation_mode = 'from_detailed_data';
188  follicle_rect_experiments(2,model,model_data);
189 
190  case 3 % detailed and reduced simulation
191 
192  %model.x0_shift=0;
193  %model.u_function_ptr = @(model) 0; % define u constant 0
194  % => leads to error 0 :-) if x0_shift = 0
195 
196  %model.u_function_ptr = @(model) 0.1; % define u constant
197  %=> error linearly increasing
198  %model.u_function_ptr = @(model) 0.3*sin(2*model.t)+0.002*model.t;
199 
200  tic;
201 
202  % remove some fields such that caching works
203  tmp_model_data = model_data;
204  tmp_model = model;
205  if isfield(tmp_model_data,'RB')
206  tmp_model_data = rmfield(tmp_model_data,'RB');
207  end;
208  if isfield(tmp_model_data,'RB_generation_mode')
209  tmp_model = rmfield(tmp_model,'RB_generation_mode');
210  end;
211  sim_data = filecache_function(@detailed_simulation,tmp_model,tmp_model_data);
212  t1 = toc;
213  fprintf('\n');
214  disp(['time for detailed simulation: ',num2str(t1)]);
215 
216 % model.plot_title = 'full simulation';
217 % model.plot_sim_data_state(model,model_data,sim_data,model);
218 
219  if isempty(detailed_data)
220  detailed_data = gen_detailed_data(model,model_data);
221  end;
222 
223  reduced_data = gen_reduced_data(model,detailed_data);
224  model.N = reduced_data.N;
225  tic;
226  rb_sim_data = rb_simulation(model,reduced_data);
227  t2 = toc;
228  disp(['time for reduced simulation: ',num2str(t2)]);
229 
230  rb_sim_data = rb_reconstruction(model,detailed_data, ...
231  rb_sim_data);
232 
233 % model.plot_title = 'reduced simulation';
234 % model.plot_sim_data_state(model,model_data,rb_sim_data,model);
235 
236 % figure, plot3(sim_data.X(1,:),...
237 % sim_data.X(2,:),...
238 % sim_data.X(3,:), ...
239 % 'color','b' ...
240 % );
241 % hold on
242 % plot3(rb_sim_data.X(1,:),...
243 % rb_sim_data.X(2,:),...
244 % rb_sim_data.X(3,:),...
245 % 'color','r');
246 % axis equal;
247 % title('exact and reduced trajectories');
248 % legend({'exact','reduced'});
249 
250  er = sim_data;
251  er.X = sim_data.X-rb_sim_data.X;
252  er.Y = sim_data.Y-rb_sim_data.Y;
253 
254  model.plot_title = 'error';
255  model.plot_sim_data_state(model,model_data,er,model);
256 
257  err = sqrt(sum((detailed_data.G * er.X).*er.X,1));
258  max_t_err = err;
259 
260  % somehow replace by matrix version!!!
261  for i = 2:length(max_t_err)
262  if max_t_err(i)<max_t_err(i-1)
263  max_t_err(i) = max_t_err(i-1);
264  end;
265  end;
266 
267  % figure, plot(sim_data.time,[rb_sim_data.DeltaX(:),max_t_err(:),err(:)]),
268  if model.enable_error_estimator
269  figure, plot([rb_sim_data.DeltaX(:),max_t_err(:),err(:)]),
270  Ylim = get(gca,'Ylim');
271  set(gca,'Ylim',[0,max(Ylim)]);
272  title('rb state error and estimator');
273  legend({'estimator \Delta_x(t)','||e||_{L\infty([0,t],L2)}',...
274  '||e(t)||_{L2}'});
275  else
276  figure, plot([max_t_err(:),err(:)]),
277  Ylim = get(gca,'Ylim');
278  set(gca,'Ylim',[0,max(Ylim)]);
279  title('rb state error');
280  legend({'||e||_{L\infty([0,t],L2)}',...
281  '||e(t)||_{L2}'});
282  end;
283 
284  if model.dim_y > 1
285  er.Ynorm = sum(er.Y.^2)';
286  else
287  er.Ynorm = abs(er.Y(:));
288  end;
289 
290  if model.enable_error_estimator
291  figure, plot([rb_sim_data.DeltaY(:),er.Ynorm]);
292  title('rb output error and estimator');
293  Ylim = get(gca,'Ylim');
294  set(gca,'Ylim',[0,max(Ylim)]);
295  legend({'estimator','error'});
296  figure, plot(rb_sim_data.time_indices,rb_sim_data.Rnorms);
297  title('Residual norms')
298  else
299  figure, plot([er.Ynorm]);
300  title('rb output error');
301  Ylim = get(gca,'Ylim');
302  set(gca,'Ylim',[0,max(Ylim)]);
303  legend({'error'});
304  end;
305 
306 % plot_sim_data(model,model_data,sim_data,model);
307  if model.debug
308  disp(['finished, please interact with data']);
309  keyboard;
310  end;
311 
312  case 3.1 % single trajectory, lower states as reduced basis
313 
314  par = [];
315 % par.output_type = 'expectation_N';
316  model = follicle_rect_model(par);
317  model = model.set_mu(model,[0.01,0.01]);
318  model_data = gen_model_data(model);
319  % length i = 136 (range <=15)
320  i = find((model_data.Ns+model_data.Ms)<=15); % i.e. lower states
321 
322  E = speye(model.dim_x);
323  model_data.RB = E(:,i);
324  model.RB_generation_mode = 'from_detailed_data';
325  follicle_rect_experiments(3,model,model_data);
326  keyboard;
327 
328  case 3.2
329  % generate POD basis of single trajectory
330  params = [];
331  %params.range = 200;
332  %params.output_range = 50;
333  % params.output_plot_indices = [1,2];
334  %params.init_range = 20;
335 % params.output_type = 'separatrix';
336  %params.output_type = 'expectation_N';
337 % params.reg_epsilon = 0.00009;
338  %params.reg_epsilon = 0.00000;
339  model = follicle_rect_model(params);
340  % model.RB_generation_mode = 'PCA_trajectory';
341  % model.RB_stop_Nmax = 10;
342  model = model.set_mu(model,[0.01,0.01]);
343  model_data = gen_model_data(model);
344  model.RB_generation_mode = 'from_detailed_data';
345  tmpmodel = model;
346  tmpmodel.save_time_indices = 0:model.nt;
347  disp('computing trajectory');
348  sim_data = filecache_function(@detailed_simulation,tmpmodel,model_data);
349  disp('computing POD');
350  RBinit = orthonormalize_gram_schmidt(sim_data.X(:,1),model_data.G,1e-20);
351  Xorth = sim_data.X(:,2:end)- ...
352  RBinit * ((RBinit' * model_data.G) * sim_data.X(:,2:end));
353  % keyboard;
354  clear('sim_data');
355  [U,S,V,flag] = svds(Xorth,20); % single svds(..20) takes about
356  % 5 min!!
357  i = find(diag(S)>=eps);
358  U = U(:,i);
359  %[U,S] = svds_batch(sim_data.X,5,100);
360  % warning: basis not fully orthogonal!!!
361  % U = orth(U);
362  disp('computing gram-schmidt orthogonalization');
363  model_data.RB = orthonormalize_gram_schmidt([RBinit,U],...
364  model_data.G,eps);
365  % keyboard;
366  detailed_data = gen_detailed_data(model,model_data);
367  save(fullfile(rbmatlabtemp,'follicle_POD_detailed_data.mat'),...
368  'model','model_data','detailed_data');
369 
370  model.clim = [-0.00005,0.00005];
371  model.plot_detailed_data(model,detailed_data,[]);
372 % follicle_rect_experiments(3,model,model_data);
373 
374  case 3.3 % generate BT basis
375  % todo
376 
377  case 3.4 % generate POD-Greedy basis with projection error 9x9 grid
378 
379  % fn = fullfile(rbmatlabtemp,'follicle_rect_PODGreedy_detailed_data_proj_err_9x9.mat');
380  fn = fullfile(rbmatlabtemp,'follicle_rect_PODGreedy_detailed_data_proj_err_9x9_large_sorted.mat');
381  %fn = fullfile(rbmatlabtemp,'temp_detailed_data.mat');
382  if ~exist(fn)
383  % POD-Greedy basis with true projection error as indicator
384  tic;
385  par = [];
386  % par.range = 200;
387  % par.output_range = 50;
388  % par.init_range = 20;
389 
390  model = follicle_rect_model(par);
391  model.RB_generation_mode = 'greedy_log_uniform_fixed';
392  model.RB_detailed_train_savepath = ...
393  'follicle_rect_PODGreedy_train_data_proj_err_9x9';
394  %model.RB_error_indicator = 'projection-error';
395  model.RB_stop_epsilon = 1e-16;
396  model.RB_numintervals = [8,8];
397  model.RB_stop_timeout = 12*60*60; % 12 hours
398  model_data = gen_model_data(model);
399  detailed_data = gen_detailed_data(model,model_data);
400  save(fn,'model','model_data','detailed_data');
401  t = toc;
402  disp('computation time: ')
403  t
404 
405  end;
406  load(fn);
407  % length i = 136 (range <=15)
408  %i = find((model_data.Ns+model_data.Ms)<=15); % i.e. lower states
409 
410  model_data.RB = detailed_data.V;
411  model.RB_generation_mode = 'from_detailed_data';
412  follicle_rect_experiments(3,model,model_data);
413 
414  model.clim = [-0.00005,0.00005];
415  model.plot_detailed_data(model,detailed_data,[]);
416 
417  %%%%%%%%%
418  % old ...
419 
420  case 3.5
421 
422  % test POD basis, only part of basis vectors
423 % fn = fullfile(rbmatlabtemp,'follicle_multi_POD_detailed_data.mat');
424  fn = fullfile(rbmatlabtemp, ...
425  'follicle_PODgreedy_detailed_data_non_orth.mat');
426  if ~exist(fn)
427  error('please run step 3.4 first!!');
428  end;
429  load(fn);
430  model.clim = [-0.00005,0.00005];
431  plot_detailed_data(model,detailed_data,[]);
432  return;
433  % model_data.RB = model_data.RB(:,1:2); => Acceleration 1000 !!
434 % RB = orthonormalize_qr(model_data.RB(:,1:20));
435 % RB = orth(model_data.RB(:,1:2));
436  RB = model_data.RB(:,1:4);
437  model_data.RB = RB;
438  model.debug = 1;
439 % model.enable_error_estimator = 1;
440  model.enable_error_estimator = 1;
441  follicle_rect_experiments(3,model,model_data);
442  keyboard;
443 
444  case 3.7
445 
446  % test POD-Greedy basis
447 % fn = fullfile(rbmatlabtemp,'follicle_PODGreedy_detailed_data.mat');
448  fn = fullfile(rbmatlabtemp,'follicle_rect_PODGreedy_detailed_data_proj_err_9x9_large.mat');
449  if ~exist(fn)
450  error('please run step 3* (or any similar) first!!');
451  end;
452  load(fn);
453 % model_data.RB = model_data.RB(:,1:2); => Acceleration 1000 !!
454 % model_data.RB = model_data.RB(:,1:20);
455 % model.debug = 1;
456  model.enable_error_estimator = 0;
457  follicle_rect_experiments(3,model,model_data,detailed_data);
458  plot_detailed_data(model,detailed_data,[]);
459  keyboard;
460 
461  case 3.8 % generate BT export matrices
462 
463  params.range = 200;
464 % params.range = 70;
465  params.output_range = 50;
466  % params.output_plot_indices = [1,2];
467  model = follicle_rect_model(params);
468  model.clim = [0,0.00005];
469  model_data = gen_model_data(model);
470 
471  fn = 'follicle_rect_model_system_matrices.mat';
472  fn2 = 'follicle_rect_model_BT.m';
473  fn3 = 'follicle_rect_model_BT_basis.mat';
474  fullfn = fullfile(rbmatlabtemp, ...
475  fn);
476  fullfn2 = fullfile(rbmatlabtemp, ...
477  fn2);
478  if exist(fullfn)
479  delete(fullfn);
480 % error(['file ',fullfn,' existing, please remove.']);
481  end;
482  if exist(fullfn2)
483  delete(fullfn2);
484 % error(['file ',fullfn2,' existing, please remove.']);
485  end;
486 
487  % matrix export for remote BT basis computation
488  model.decomp_mode = 0;
489  A = filecache_function(model.A_function_ptr,model,model_data);
490  B = model.B_function_ptr(model,model_data);
491  C = model.C_function_ptr(model,model_data);
492  D = model.D_function_ptr(model,model_data);
493  save(fullfn,'model','A','B','C','D');
494  disp(['system matrices exported to file ',fullfn]);
495 
496  diary(fullfile(rbmatlabtemp,fn2));
497  disp( [' load(''',fn,'''); '...
498  ' sys = ss(A,B,C,D); ',...
499  ' tic; ', ...
500  ' [sysb,G,T,Tinv] = balreal(sys); ', ...
501  ' time = toc; ', ...
502  ' W = transpose(T(1:500,:)); ',...
503  ' V = Tinv(:,1:500); ',...
504  ' H = G(1:500); ',...
505  ' save([''~/',fn3,'''],''V'',''W'',''H'',''time'');']);
506  diary off;
507  disp(['BT script exported to file ',fullfn2]);
508 
509  disp(['please copy files to remote host into your matlab startup' ...
510  ' directory and execute'])
511  disp([' matlab -nodesktop -nosplash -r ',fn2]);
512 
513  disp(['then copy remote file ~/',fn3,' to local ',rbmatlabtemp,' and press enter'])
514 
515  keyboard;
516 
517  tmp = load(fullfile(rbmatlabtemp,fn3));
518  if isfield(tmp,'V') && isfield(tmp,'W')
519  disp('BT basis OK')
520  else
521  error('error in file generation');
522  end;
523 
524  case 3.9
525  % test BT basis
526  % todo
527 
528  case 4.0 % reduced basis in demo_rb_gui;
529 
530  % test POD basis only part of basis vectors
531 % fn = fullfile(rbmatlabtemp,'follicle_POD_detailed_data.mat');
532 % fn = fullfile(rbmatlabtemp,'follicle_PODgreedy_detailed_data_eps0.0186714.mat');
533 %fn = fullfile(rbmatlabtemp, ...
534 % 'follicle_PODGreedy_detailed_data_eps0.0186714.mat');
535 %fn = fullfile(rbmatlabtemp,'follicle_PODgreedy_detailed_data_eps0.mat');
536 %fn = fullfile(rbmatlabtemp,'follicle_PODgreedy_detailed_data_proj_err.mat');
537 fn = fullfile(rbmatlabtemp,'follicle_rect_PODgreedy_detailed_data_proj_err_9x9_large_sorted.mat');
538 % fn = fullfile(rbmatlabtemp,'follicle_multi_POD_detailed_data.mat');
539 % fn = fullfile(rbmatlabtemp,'follicle_PODGreedy_detailed_data_non_orth.mat');
540  if ~exist(fn)
541  error('please run step 3.4 (or similar) first!!');
542  end;
543  load(fn);
544  model.N = size(detailed_data.V,2);
545 % model.clim = [0.00000,0.00005];
546  model.clim = [0.00000,0.00000005];
547  demo_rb_gui(model,detailed_data);
548  keyboard;
549 
550  case 5.0
551  % solve optimization problem: target output curve for mu = (0.01,0.01)
552 
553  % load target output curve
554  params = [];
555  %params.range = 200;
556  %params.output_range = 50;
557  % params.output_plot_indices = [1,2];
558  %params.init_range = 20;
559 % params.output_type='separatrix';
560 % params.output_range = 50;
561  params.output_type='expectation_N';
562 % params.output_type='expectation_M';
563  model = follicle_rect_model(params);
564  model.clim = [0.00000,0.00000005];
565 % model.clim = [0,0.00005];
566  model_data = gen_model_data(model);
567 % mu = [0.013,0.003];
568  mu = [0.01,0.01];
569 % format long
570 % r = rand(1,2)*0.0005
571 % 1.0e-003 *
572 % 0.063493408146753 0.456687928069510
573 % mu0 = [0.019,0.0055] ;
574 % mu = mu0 + r
575 % mu = [0.019063493408147 0.005956687928070]
576 
577  model = model.set_mu(model,mu);
578 % mu = model.get_mu(model);
579  sim_data = filecache_function(@detailed_simulation,model,model_data);
580  Ytarget = sim_data.Y;
581 
582  % load reduced model
583 % fn = fullfile(rbmatlabtemp,'follicle_multi_POD_detailed_data.mat');
584 % fn = fullfile(rbmatlabtemp,'follicle_POD_detailed_data.mat');
585 % fn = fullfile(rbmatlabtemp,'follicle_PODGreedy_detailed_data_hand_orth.mat');
586 fn = fullfile(rbmatlabtemp,'follicle_rect_PODgreedy_detailed_data_proj_err_9x9_large.mat');
587 %fn = fullfile(rbmatlabtemp,'follicle_rect_PODgreedy_detailed_data_proj_err_9x9.mat');
588 % fn = fullfile(rbmatlabtemp,'follicle_PODgreedy_detailed_data_eps0.mat');
589 % fn = fullfile(rbmatlabtemp,'follicle_PODGreedy_detailed_data.mat');
590  if ~exist(fn)
591  error('please run basis generation step 3.* first!!');
592  end;
593 % keyboard;
594 % the following is required, if another output functional is desired
595  model_old = model;
596  model_data_old = model_data;
597  load(fn); % separatrix, 8x8
598  model = model_old; % expectation_N, 20x20
599  model_data.C_components = model_data_old.C_components;
600  detailed_data.C_components = model_data_old.C_components;
601 
602  sim_data_check = filecache_function(@detailed_simulation,model,model_data);
603  Ytarget_check = sim_data_check.Y;
604 
605  if ~isequal(Ytarget,Ytarget_check)
606  error('model and stored data inconsistent!');
607  end;
608  % model.C_function_ptr = model_old.C_function_ptr;
609  % keyboard;
610 % model.N = floor(size(detailed_data.V,2)/2);
611 % model.N = 12;
612  model.N = size(detailed_data.V,2);
613  reduced_data = gen_reduced_data(model,detailed_data);
614 
615  % set up optimization function
616  fun = @(mu) my_opt_target(mu,Ytarget,model,reduced_data)
617 
618  if 1
619  disp('generating surface');
620  % plot target surface
621  logr = log([0.005,0.02]);
622  logrange = logr(1):(logr(2)-logr(1))/20:logr(2);
623  range = exp(logrange)
624  [X,Y] = meshgrid(range,range);
625  Z = zeros(size(X));
626  for i = 1:length(X(:))
627  Z(i) = fun([X(i),Y(i)]);
628 % if (Z(i))>1e10
629 % keyboard;
630 % end;
631  end;
632  surf(X,Y,Z);
633  set(gca,'Zscale','log')
634  xlabel('u1');
635  ylabel('u2');
636  zlabel('J(u1,u2)','Fontsize',15);
637  title('optimization target J');
638  set(gca,'Fontsize',15);
639  yl = get(gca,'Ylabel');set(yl,'Fontsize',20);
640  xl = get(gca,'Xlabel');set(xl,'Fontsize',20);
641  ti = get(gca,'Title');set(ti,'Fontsize',20);
642  axis tight;
643  zl = get(gca,'Zlim'); % slightly stretch z-axes
644  zl = 1.1*(zl-sum(zl))+sum(zl);
645  set(gca,'Zlim',zl);
646  saveas(gcf,fullfile(rbmatlabhome,'optimization_target.eps'),'epsc');
647  saveas(gcf,fullfile(rbmatlabhome,'optimization_target.png'),'png');
648  close(gcf);
649  end;
650 
651 % mu0 = [0.005,0.005];
652 % mu0 = [0.02,0.005]; % good parameter for random target parameter
653 % mu0 = [0.02,0.02]; % => 6% rel error in mu, finds better local optimum!
654 % mu0 = [0.01,0.01]; % => rel error 3*10-5 = 0.003%
655 % mu0 = [0.005,0.02]; % => rel error 7.5%
656 % mu0 = [ 0.008435225652526 0.018314710503781]; %=> 6% rel error in mu
657 % mu0 = [ 0.017736939588032 0.019009898716363 ]; % => 11.5% rel error
658 % mu0 = [0.015181027322867 0.016366101958675]; % => 14% rel error
659  mu0 = [0.016146987021874 0.010883405293013]; % => 4.2 rel error
660  lb = [0.005,0.005];
661  ub = [0.02,0.02];
662 % options = [];
663  options = optimset();
664  options.MaxFunEvals = 1000;
665  options.TolFun = 1e-14;
666  options.TolX = 1e-14;
667  [mumin,fval,exitflag,output,lambda,grad] = fmincon(fun,mu0,[],[], ...
668  [],[],lb,ub,[],options);
669 % [mumin,fval,exitflag,output,lambda,grad] = fmincon(fun,mu0,[],[], ...
670 % [],[],lb,ub);
671  disp('the obtained optimal mu value is: ')
672  mumin
673 
674  disp('And target value J(mumin)')
675  fun(mumin)
676 
677  disp('The original mu value is')
678  mu
679 
680  disp('relative error in mu value:')
681  disp(norm(mu-mumin)/norm(mu));
682 
683  disp('And target value J(mu)')
684  fun(mu)
685 
686 % disp('optimization finished, please inspect workspace');
687  keyboard;
688 
689  %%%%%%%%%%%%%%%%%%%%%%%%%% single parameter optimization
690 
691  disp('starting 1d-optimization')
692  % set up optimization function
693  fun = @(mu1) my_opt_target_single_param(mu1,Ytarget,model,reduced_data)
694 
695  if 1
696  disp('generating surface');
697  % plot target surface
698  logr = log([0.005,0.02]);
699  logrange = logr(1):(logr(2)-logr(1))/100:logr(2);
700  X = exp(logrange)
701  Z = zeros(size(X),1);
702  for i = 1:length(X(:))
703  Z(i) = fun([X(i)]);
704 % if (Z(i))>1e10
705 % keyboard;
706 % end;
707  end;
708  plot(X,Z,'Linewidth',2);
709  set(gca,'Yscale','log')
710  xlabel('u1');
711  ylabel('J(u1,0.01)');
712  title('optimization target');
713  saveas(gcf,fullfile(rbmatlabhome,'optimization_target_single_param.eps'),'epsc');
714  saveas(gcf,fullfile(rbmatlabhome,'optimization_target_single_param.png'),'png');
715  keyboard;
716  close(gcf);
717  end;
718 
719 % mu0 = [0.005,0.005];
720 % mu0 = [0.02,0.005]; % good parameter for random target parameter
721 % mu0 = [0.02,0.02]; % => 6% rel error in mu, finds better local optimum!
722 % mu0 = [0.01,0.01]; % => rel error 3*10-5 = 0.003%
723 % mu0 = [0.005,0.02]; % => rel error 7.5%
724 % mu0 = [ 0.008435225652526 0.018314710503781]; %=> 6% rel error in mu
725 % mu0 = [ 0.017736939588032 0.019009898716363 ]; % => 11.5% rel error
726 % mu0 = [0.015181027322867 0.016366101958675]; % => 14% rel error
727  mu0 = 0.005;
728  lb = [0.005];
729  ub = [0.02];
730 % options = [];
731  options = optimset();
732  options.MaxFunEvals = 1000;
733  options.TolFun = 1e-14;
734  options.TolX = 1e-14;
735  [mumin,fval,exitflag,output,lambda,grad] = fmincon(fun,mu0,[],[], ...
736  [],[],lb,ub,[],options);
737 % [mumin,fval,exitflag,output,lambda,grad] = fmincon(fun,mu0,[],[], ...
738 % [],[],lb,ub);
739  disp('the obtained optimal mu1 value is: ')
740  mumin
741 
742  disp('And target value J(mumin,0.01)')
743  fun(mumin)
744 
745  disp('The original mu1 value is')
746  mu(1)
747 
748  disp('relative error in mu value:')
749  disp(norm(mu(1)-mumin)/norm(mu(1)));
750 
751  disp('And target value J(mu)')
752  fun(mu(1))
753 
754  % results of optimization:
755 
756  %mumin =
757  % 0.010525505381718 0.010279932143424
758  %And target value J(mumin)
759  %J(0.010526 0.01028).00011691
760  %ans =
761  % 1.169063479811699e-004
762  %The original mu value is
763  %mu =
764  % 0.010000000000000 0.010000000000000
765  %relative error in mu value:
766  % 0.042102132436309
767  %And target value J(mu)
768  %J(0.01 0.01).00019323
769  %ans =
770  % 1.932325185788862e-004
771  % iterations: 60
772  % funcCount: 383
773  % lssteplength: 1
774  % stepsize: 1.571820776306015e-014
775  % algorithm: 'medium-scale: SQP, Quasi-Newton, line-search'
776  % firstorderopt: 2.703745493153992e+002
777  % constrviolation: -0.005279932143424
778  % message: [1x773 char]
779 
780 
781  %the obtained optimal mu1 value is:
782  %mumin =
783  % 0.0100
784  %And target value J(mumin)
785  %J(0.0099999 0.01).00016254
786  % 1.6254e-004
787  %The original mu1 value is
788  % 0.0100
789  %relative error in mu value:
790  % 1.1433e-005
791  %And target value J(mu)
792  %J(0.01 0.01).00019323
793  % 1.9323e-004
794 
795  case 5.1 % new experiments for revised paper
796 
797  % solve optimization problem: target output curve for mu = (0.01,0.01)
798  % add noise!!
799 
800  % load target output curve
801  params = [];
802  %params.range = 200;
803  %params.output_range = 50;
804  % params.output_plot_indices = [1,2];
805  %params.init_range = 20;
806 % params.output_type='separatrix';
807 % params.output_range = 50;
808  params.output_type='expectation_N';
809 % params.output_type='expectation_M';
810  model = follicle_rect_model(params);
811  model.clim = [0.00000,0.00000005];
812 % model.clim = [0,0.00005];
813  model_data = gen_model_data(model);
814 % mu = [0.013,0.003];
815  mu = [0.01,0.01];
816 % format long
817 % r = rand(1,2)*0.0005
818 % 1.0e-003 *
819 % 0.063493408146753 0.456687928069510
820 % mu0 = [0.019,0.0055] ;
821 % mu = mu0 + r
822 % mu = [0.019063493408147 0.005956687928070]
823 
824  model = model.set_mu(model,mu);
825 % mu = model.get_mu(model);
826  sim_data = filecache_function(@detailed_simulation,model,model_data);
827  Ytarget_clean = sim_data.Y;
828 
829  % set random generator to deterministic value, noise 1% relative
830  % to Y
831  RandStream.setDefaultStream(RandStream('mt19937ar','seed',12345));
832  noise_level = 0.05;
833  noise = randn(size(Ytarget_clean)).*noise_level.*Ytarget_clean;
834  Ytarget = Ytarget_clean + noise;
835  p1 = plot(sim_data.time,Ytarget_clean,'k-','Linewidth',2);
836  hold on;
837  p2 = plot(sim_data.time,Ytarget,'b-');
838  title('clean and noisy output');
839  set(gca,'Fontsize',20);
840  title('clean and noisy reference output');
841  xlabel('time','Fontsize',20);
842  ylabel('output','Fontsize',20);
843  legend({'clean, y(t,\theta_{ref})','noisy, y_{meas}(t)'},'Location','SouthEast');
844  saveas(gcf,fullfile(rbmatlabhome,'clean_noisy_output.eps'),'epsc');
845  saveas(gcf,fullfile(rbmatlabhome,'clean_noisy_output.png'),'png');
846  disp('press enter to continue');
847  pause;
848  close(gcf);
849 
850  % load reduced model
851 % fn = fullfile(rbmatlabtemp,'follicle_multi_POD_detailed_data.mat');
852 % fn = fullfile(rbmatlabtemp,'follicle_POD_detailed_data.mat');
853 % fn = fullfile(rbmatlabtemp,'follicle_PODGreedy_detailed_data_hand_orth.mat');
854 fn = fullfile(rbmatlabtemp,'follicle_rect_PODgreedy_detailed_data_proj_err_9x9_large_sorted.mat');
855 %fn = fullfile(rbmatlabtemp,'follicle_rect_PODgreedy_detailed_data_proj_err_9x9.mat');
856 % fn = fullfile(rbmatlabtemp,'follicle_PODgreedy_detailed_data_eps0.mat');
857 % fn = fullfile(rbmatlabtemp,'follicle_PODGreedy_detailed_data.mat');
858  if ~exist(fn)
859  error('please run basis generation step 3.* first!!');
860  end;
861 % keyboard;
862 % the following is required, if another output functional is desired
863  model_old = model;
864  model_data_old = model_data;
865  load(fn); % separatrix, 8x8
866  model = model_old; % expectation_N, 20x20
867  model_data.C_components = model_data_old.C_components;
868  detailed_data.C_components = model_data_old.C_components;
869 
870  sim_data_check = filecache_function(@detailed_simulation,model,model_data);
871  Ytarget_check = sim_data_check.Y;
872 
873  if ~isequal(Ytarget_clean,Ytarget_check)
874  error('model and stored data inconsistent!');
875  end;
876  % model.C_function_ptr = model_old.C_function_ptr;
877  % keyboard;
878 % model.N = floor(size(detailed_data.V,2)/2);
879 % model.N = 12;
880  model.N = size(detailed_data.V,2);
881  reduced_data = gen_reduced_data(model,detailed_data);
882 
883  % set up optimization function
884  fun = @(mu) my_opt_target(mu,Ytarget,model,reduced_data)
885 
886  if 1
887  disp('generating surface');
888  % plot target surface
889  logr = log([0.005,0.02]);
890  logrange = logr(1):(logr(2)-logr(1))/20:logr(2);
891  range = exp(logrange)
892  [X,Y] = meshgrid(range,range);
893  Z = zeros(size(X));
894  for i = 1:length(X(:))
895  Z(i) = fun([X(i),Y(i)]);
896 % if (Z(i))>1e10
897 % keyboard;
898 % end;
899  end;
900  surf(X,Y,Z);
901  set(gca,'Zscale','log')
902  xlabel('u1');
903  ylabel('u2');
904  zlabel('J(u1,u2)','Fontsize',15);
905  title('optimization target J');
906  set(gca,'Fontsize',15);
907  yl = get(gca,'Ylabel');set(yl,'Fontsize',20);
908  xl = get(gca,'Xlabel');set(xl,'Fontsize',20);
909  ti = get(gca,'Title');set(ti,'Fontsize',20);
910  axis tight;
911  zl = get(gca,'Zlim'); % slightly stretch z-axes
912  zl = 1.1*(zl-sum(zl))+sum(zl);
913  set(gca,'Zlim',zl);
914  set(gca,'Clim',[-20000,200000])
915  saveas(gcf,fullfile(rbmatlabhome,'optimization_target_new.eps'),'epsc');
916  saveas(gcf,fullfile(rbmatlabhome,'optimization_target_new.png'),'png');
917  disp('press key to continue')
918  keyboard;
919  pause;
920  close(gcf);
921  end;
922 
923  % mu0 = [0.016146987021874 0.010883405293013]; % => 29% rel error
924  mu0 = [0.005,0.005]; % => 58% relative error
925  mu0 = [ 0.008435225652526 0.018314710503781]; %=> 46% rel error in mu
926 % mu0 = [0.02,0.005]; % good parameter for random target parameter
927 % mu0 = [0.02,0.02]; % => 6% rel error in mu, finds better local optimum!
928 % mu0 = [0.01,0.01]; % => rel error 3*10-5 = 0.003%
929 % mu0 = [0.005,0.02]; % => rel error 7.5%
930 % mu0 = [ 0.017736939588032 0.019009898716363 ]; % => 11.5% rel error
931 % mu0 = [0.015181027322867 0.016366101958675]; % => 14% rel error
932  lb = [0.005,0.005];
933  ub = [0.02,0.02];
934 % options = [];
935  options = optimset();
936  options.MaxFunEvals = 1000;
937  options.TolFun = 1e-14;
938  options.TolX = 1e-14;
939  [mumin,fval,exitflag,output,lambda,grad] = fmincon(fun,mu0,[],[], ...
940  [],[],lb,ub,[],options);
941 % [mumin,fval,exitflag,output,lambda,grad] = fmincon(fun,mu0,[],[], ...
942 % [],[],lb,ub);
943  disp('the obtained optimal mu value is: ')
944  mumin
945 
946  disp('And target value J(mumin)')
947  fun(mumin)
948 
949  disp('The original mu value is')
950  mu
951 
952  disp('relative error in mu value:')
953  disp(norm(mu-mumin)/norm(mu));
954 
955  disp('And target value J(mu)')
956  fun(mu)
957 
958 % disp('optimization finished, please inspect workspace');
959  keyboard;
960 
961  %%%%%%%%%%%%%%%%%%%%%%%%%% single parameter optimization
962 
963  disp('starting 1d-optimization')
964  % set up optimization function
965  fun = @(mu1) my_opt_target_single_param(mu1,Ytarget,model,reduced_data)
966 
967  if 1
968  disp('generating surface');
969  % plot target surface
970  logr = log([0.005,0.02]);
971  logrange = logr(1):(logr(2)-logr(1))/100:logr(2);
972  X = exp(logrange)
973  Z = zeros(size(X),1);
974  for i = 1:length(X(:))
975  Z(i) = fun([X(i)]);
976 % if (Z(i))>1e10
977 % keyboard;
978 % end;
979  end;
980  plot(X,Z,'Linewidth',2);
981  set(gca,'Yscale','log')
982  xlabel('u1');
983  ylabel('J(u1,0.01)');
984  title('optimization target');
985  saveas(gcf,fullfile(rbmatlabhome,'optimization_target_single_param_new.eps'),'epsc');
986  saveas(gcf,fullfile(rbmatlabhome,'optimization_target_single_param_new.png'),'png');
987  disp('press key to continue')
988  pause;
989  close(gcf);
990  end;
991 
992 % mu0 = [0.005,0.005];
993 % mu0 = [0.02,0.005]; % good parameter for random target parameter
994 % mu0 = [0.02,0.02]; % => 6% rel error in mu, finds better local optimum!
995 % mu0 = [0.01,0.01]; % => rel error 3*10-5 = 0.003%
996 % mu0 = [0.005,0.02]; % => rel error 7.5%
997 % mu0 = [ 0.008435225652526 0.018314710503781]; %=> 6% rel error in mu
998 % mu0 = [ 0.017736939588032 0.019009898716363 ]; % => 11.5% rel error
999 % mu0 = [0.015181027322867 0.016366101958675]; % => 14% rel error
1000  mu0 = 0.005;
1001  lb = [0.005];
1002  ub = [0.02];
1003 % options = [];
1004  options = optimset();
1005  options.MaxFunEvals = 1000;
1006  options.TolFun = 1e-14;
1007  options.TolX = 1e-14;
1008  [mumin,fval,exitflag,output,lambda,grad] = fmincon(fun,mu0,[],[], ...
1009  [],[],lb,ub,[],options);
1010 % [mumin,fval,exitflag,output,lambda,grad] = fmincon(fun,mu0,[],[], ...
1011 % [],[],lb,ub);
1012  disp('the obtained optimal mu1 value is: ')
1013  mumin
1014 
1015  disp('And target value J(mumin,0.01)')
1016  fun(mumin)
1017 
1018  disp('The original mu1 value is')
1019  mu(1)
1020 
1021  disp('relative error in mu value:')
1022  disp(norm(mu(1)-mumin)/norm(mu(1)));
1023 
1024  disp('And target value J(mu)')
1025  fun(mu(1))
1026 
1027 
1028 % model = set_mu(model,[mumin,mu(2)]);
1029 % rb_sim_data = rb_simulation(model,reduced_data);
1030 % Ytarget_recovered = rb_sim_data.Y;
1031 
1032 % plot(rb_sim_data.time,[Ytarget_clean; Ytarget; Ytarget_recovered]');
1033 % set(gca,'Fontsize',20);
1034 % title('Clean, noisy and recovered output');
1035 % xlabel('time','Fontsize',20);
1036 % ylabel('output','Fontsize',20);
1037  keyboard;
1038 
1039 
1040  case 6 % generate trajectory images for paper
1041 
1042  fn = fullfile(rbmatlabtemp,'follicle_rect_PODgreedy_detailed_data_proj_err_9x9_large_sorted.mat');
1043  load(fn);
1044 
1045  % model settings:
1046  % end time: 10 mio sec. = approx 19 years
1047  % dimension: 22801 = (151)^2
1048  % number time steps: 500
1049  % output: expectation_N
1050 
1051  % plot of some trajectories of model:
1052  % 3 snapshots over time for 2 parameters
1053  % and one output plot
1054 
1055  mus = {[0.02, 0.005],[0.005,0.02]};
1056  tis = [0,1,250,499];
1057 
1058  plot_params = [];
1059  plot_params.show_colorbar = 0;
1060  plot_params.axes_tight = 1;
1061  plot_params.no_lines = 1;
1062  plot_params.clim = [0.00000,0.00000005];
1063 % my_model = follicle_rect_model([]);
1064  plot_params.colormap = model.transformed_colormap;
1065 
1066  for i = 1:length(mus)
1067  mu = mus{i};
1068  model = model.set_mu(model,mu);
1069  tic, sim_data = filecache_function(@detailed_simulation,model, ...
1070  model_data);
1071  t = toc;
1072  disp(['detailed simulation: ',num2str(t),'seconds']);
1073  for ti = 1:length(tis)
1074  plot_params.timestep = tis(ti);
1075  if (ti==length(tis))
1076  plot_params.show_colorbar = 1;
1077  else
1078  plot_params.show_colorbar = 0;
1079  end;
1080  if (ti==1)
1081  plot_yticks = 1;
1082  else
1083  plot_yticks = 0;
1084  end;
1085  if (i==2)
1086  plot_xticks = 1;
1087  else
1088  plot_xticks = 0;
1089  end;
1090  model.plot_slice(model,model_data,sim_data,plot_params);
1091 % i
1092 % ti
1093  axis equal;
1094  axis tight;
1095  title('');
1096  if plot_params.show_colorbar
1097  set(gca,'Fontsize',20)
1098  end;
1099  if ~plot_xticks
1100  set(gca,'Xtick',[])
1101  xlabel('');
1102  else % enlarge font
1103  set(gca,'Fontsize',20)
1104  xl = get(gca,'Xlabel');set(xl,'Fontsize',25);
1105  end;
1106  if ~plot_yticks
1107  set(gca,'Ytick',[])
1108  ylabel('');
1109  else % enlarge font
1110  set(gca,'Fontsize',20)
1111  yl = get(gca,'Ylabel');set(yl,'Fontsize',25);
1112  end;
1113 
1114  saveas(gcf,fullfile(rbmatlabhome,...
1115  ['snapshot_mu',num2str(mu(1)),'_',...
1116  num2str(mu(2)),'_t',num2str(tis(ti)),'.eps']),'epsc');
1117  saveas(gcf,fullfile(rbmatlabhome,...
1118  ['snapshot_mu',num2str(mu(1)),'_',...
1119  num2str(mu(2)),'_t',num2str(tis(ti)),'.png']),'png');
1120  close(gcf);
1121  end;
1122  end;
1123 
1124  case 7 % generate further images for paper: basis vector plots
1125 
1126  % generate plots of basis vectors
1127 
1128  fn = fullfile(rbmatlabtemp,'follicle_rect_PODgreedy_detailed_data_proj_err_9x9_large_sorted.mat');
1129  load(fn);
1130 % load(fullfile(rbmatlabtemp,'follicle_temp_detailed_data'));
1131  model_data = gen_model_data(model);
1132  plot_params = [];
1133  plot_params.no_lines = 1;
1134  plot_params.timestep = 0;
1135  plot_params.show_colorbar = 0;
1136  plot_params.clim = [-0.0001,0.0001];
1137  plot_params.colormap = jet(256);
1138  sim_data =[];
1139  for i = 1:size(detailed_data.V,2);
1140  sim_data.X = detailed_data.V(:,i);
1141  model.plot_slice(model,model_data,sim_data,plot_params);
1142  colormap(jet(256));
1143  axis off
1144  axis equal
1145  title('');
1146  saveas(gcf,fullfile(rbmatlabhome,['basis_vector_',num2str(i),'.eps']),'epsc');
1147  saveas(gcf,fullfile(rbmatlabhome,['basis_vector_',num2str(i),'.png']),'png');
1148  close(gcf);
1149  end;
1150 
1151  case 8 %generate further plots: error decay curve and parameter selection
1152 
1153  fn = fullfile(rbmatlabtemp,'follicle_rect_PODgreedy_detailed_data_proj_err_9x9_large_sorted.mat');
1154  load(fn);
1155 
1156  % generate error decay curve
1157  p = plot(detailed_data.RB_info.max_err_sequence);
1158  set(gca,'Yscale','log');
1159  set(p,'Linewidth',2);
1160  title('POD-Greedy error decay');
1161  xlabel('number of basis vectors');
1162  ylabel('maximum error');
1163  set(gca,'Fontsize',15);
1164  yl = get(gca,'Ylabel');set(yl,'Fontsize',20);
1165  xl = get(gca,'Xlabel');set(xl,'Fontsize',20);
1166  ti = get(gca,'Title');set(ti,'Fontsize',20);
1167  saveas(gcf,fullfile(rbmatlabhome,'PODGreedy_error_decay.eps'),'epsc');
1168  saveas(gcf,fullfile(rbmatlabhome,'PODGreedy_error_decay.png'),'png');
1169  close(gcf);
1170 
1171  % generate plot of selected mu parameters
1172  p = plot(detailed_data.RB_info.mu_sequence(1,2:end),...
1173  detailed_data.RB_info.mu_sequence(2,2:end),'.', ...
1174  'Markersize',30);
1175  axis equal;
1176  set(gca,'Xlim',[0.004,0.021]);
1177  set(gca,'Ylim',[0.004,0.021]);
1178  title('PODGreedy selected parameters');
1179  xlabel('u1');
1180  ylabel('u2');
1181  set(gca,'Fontsize',15);
1182  yl = get(gca,'Ylabel');set(yl,'Fontsize',20);
1183  xl = get(gca,'Xlabel');set(xl,'Fontsize',20);
1184  ti = get(gca,'Title');set(ti,'Fontsize',20);
1185  saveas(gcf,fullfile(rbmatlabhome,'PODGreedy_selected_samples.eps'),'epsc');
1186  saveas(gcf,fullfile(rbmatlabhome,'PODGreedy_selected_samples.png'),'png');
1187  close(gcf);
1188 
1189  otherwise
1190  error('step unknown');
1191 end;
1192 
1193 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1194 
1195 function J = my_opt_target(mu,Ytarget,model,reduced_data);
1196 % function J = my_opt_target(mu,Ytarget,model,reduced_data);
1197 % my_opt_target
1198 
1199 tmpmodel = model.set_mu(model,mu);
1200 rb_sim_data = rb_simulation(tmpmodel,reduced_data);
1201 J = sum((rb_sim_data.Y-Ytarget).^2);
1202 disp(['J(',num2str(mu(:)'),')=',num2str(J)]);
1203 
1204 function J = my_opt_target_single_param(mu1,Ytarget,model,reduced_data);
1205 mu = model.get_mu(model);
1206 mu(1) = mu1;
1207 tmpmodel = model.set_mu(model,mu);
1208 rb_sim_data = rb_simulation(tmpmodel,reduced_data);
1209 J = sum((rb_sim_data.Y-Ytarget).^2);
1210 disp(['J(',num2str(mu(:)'),')=',num2str(J)]);
1211 
follicle_rect_model
function model = follicle_rect_model(params)
model of the human follicle growth
Definition: follicle_rect_model.m:17
filecache_function
function varargout = filecache_function(funcptr, varargin)
function used for file-caching other function calls.
Definition: filecache_function.m:17
demo_rb_gui
function demo_rb_gui(varargin)
reduced basis demo with sliders
Definition: demo_rb_gui.m:17
plot_sim_data
function p = plot_sim_data(model, model_data, sim_data, plot_params)
function performing the plot of the simulation results as specified in model.
Definition: plot_sim_data.m:17
follicle_rect_experiments
function res = follicle_rect_experiments(step, model, model_data, detailed_data)
experiments with the model of the human follicle growth
Definition: follicle_rect_experiments.m:17
Greedy
Customizable implementation of an abstract greedy algorithm.
Definition: DuneRBLeafNode.m:1
Generated on Fri Sep 30 2016 10:51:25 for rbmatlab by   doxygen 1.8.8