3 % experiments with the model of the human follicle growth
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
10 % For the results of the paper, simply run step 3.4, 5, 6, 7, 8
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
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
28 % basis is saved
for use in further steps
30 % step = 3.3, n=150, n_output , generate BT basis of
32 % basis is saved
for use in further steps
34 % step = 3.4, n=150, generate PODGreedy basis of
36 % basis is saved
for use in further steps
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 ...
44 % step 5.0: Optimization with RB model with POD basis from 3.4
45 % generation of optimization plots
for initial submission
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
53 % some steps explicitly require a model and model_data
57 % compare quality of BT to POD on single trajectory
58 % compare quality of BT to POD on parameter corners
60 % Bernard Haasdonk 21.4.2010
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
78 % model_data = gen_model_data(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';
90 case 1 % simulation and visualization
92 if isempty(model_data)
93 model_data = gen_model_data(model);
96 % sim_data = detailed_simulation(model,model_data);
98 mu = model.get_mu(model);
99 model.plot_title = ['state probabilities for mu=(',...
100 num2str(mu(1)),',',num2str(mu(2)),')'];
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');
107 disp(['time for detailed_simulation :',num2str(t)]);
109 disp(['finished, please interact with data']);
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;
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
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;
138 %model.T = 1000*model.nt;
141 save(fullfile(rbmatlabtemp,
'follicle_n200_example.mat'),...
142 'sim_data',
'model',
'model_data');
144 disp([
'computation time t = ',num2str(t),
' sec.']);
147 case 2 % reduced simulation
151 % produce high-dimensional data, e.g. reduced basis, etc.
153 detailed_data = gen_detailed_data(model,model_data);
154 disp(
'gen_detailed_data successful');
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)
160 reduced_data = gen_reduced_data(model,detailed_data);
161 disp(
'gen_reduced_data successful');
162 model.N = reduced_data.N;
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
168 disp(
'entering reduced simulation');
169 rb_sim_data = rb_simulation(model,reduced_data);
170 disp(
'reduced simulation successful');
172 rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
174 model.plot_title =
'reduced trajectory';
175 %model.plot_indices = 1:3; =>
default
178 disp([
'finished, please interact with data']);
183 % par.output_range = 20;
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';
190 case 3 % detailed and reduced simulation
193 %model.u_function_ptr = @(model) 0; % define u constant 0
194 % => leads to error 0 :-)
if x0_shift = 0
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;
202 %
remove some fields such that caching works
203 tmp_model_data = model_data;
205 if isfield(tmp_model_data,
'RB')
206 tmp_model_data = rmfield(tmp_model_data,'RB');
208 if isfield(tmp_model_data,'RB_generation_mode')
209 tmp_model = rmfield(tmp_model,'RB_generation_mode');
214 disp(['time for detailed simulation: ',num2str(t1)]);
216 % model.plot_title = 'full simulation';
217 % model.plot_sim_data_state(model,model_data,sim_data,model);
219 if isempty(detailed_data)
220 detailed_data = gen_detailed_data(model,model_data);
223 reduced_data = gen_reduced_data(model,detailed_data);
224 model.N = reduced_data.N;
226 rb_sim_data = rb_simulation(model,reduced_data);
228 disp(['time for reduced simulation: ',num2str(t2)]);
230 rb_sim_data = rb_reconstruction(model,detailed_data, ...
233 % model.plot_title = 'reduced simulation';
234 % model.plot_sim_data_state(model,model_data,rb_sim_data,model);
236 % figure, plot3(sim_data.X(1,:),...
237 % sim_data.X(2,:),...
238 % sim_data.X(3,:), ...
242 % plot3(rb_sim_data.X(1,:),...
243 % rb_sim_data.X(2,:),...
244 % rb_sim_data.X(3,:),...
247 % title('exact and reduced trajectories');
248 % legend({
'exact',
'reduced'});
251 er.X = sim_data.X-rb_sim_data.X;
252 er.Y = sim_data.Y-rb_sim_data.Y;
254 model.plot_title =
'error';
255 model.plot_sim_data_state(model,model_data,er,model);
257 err = sqrt(sum((detailed_data.G * er.X).*er.X,1));
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);
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)}',...
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)}',...
285 er.Ynorm = sum(er.Y.^2)
';
287 er.Ynorm = abs(er.Y(:));
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
')
299 figure, plot([er.Ynorm]);
300 title('rb output error
');
301 Ylim = get(gca,'Ylim
');
302 set(gca,'Ylim
',[0,max(Ylim)]);
306 % plot_sim_data(model,model_data,sim_data,model);
308 disp(['finished, please interact with data
']);
312 case 3.1 % single trajectory, lower states as reduced basis
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
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);
329 % generate POD basis of single trajectory
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
';
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));
355 [U,S,V,flag] = svds(Xorth,20); % single svds(..20) takes about
357 i = find(diag(S)>=eps);
359 %[U,S] = svds_batch(sim_data.X,5,100);
360 % warning: basis not fully orthogonal!!!
362 disp('computing gram-schmidt orthogonalization');
363 model_data.RB = orthonormalize_gram_schmidt([RBinit,U],...
366 detailed_data = gen_detailed_data(model,model_data);
367 save(fullfile(rbmatlabtemp,'follicle_POD_detailed_data.mat'),...
368 'model','model_data','detailed_data');
370 model.clim = [-0.00005,0.00005];
371 model.plot_detailed_data(model,detailed_data,[]);
374 case 3.3 % generate BT basis
377 case 3.4 % generate POD-
Greedy basis with projection error 9x9 grid
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');
383 % POD-
Greedy basis with true projection error as indicator
387 % par.output_range = 50;
388 % par.init_range = 20;
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');
402 disp('computation time: ')
407 % length i = 136 (range <=15)
408 %i = find((model_data.Ns+model_data.Ms)<=15); % i.e. lower states
410 model_data.RB = detailed_data.V;
411 model.RB_generation_mode = 'from_detailed_data';
414 model.clim = [-0.00005,0.00005];
415 model.plot_detailed_data(model,detailed_data,[]);
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');
427 error('please run step 3.4 first!!');
430 model.clim = [-0.00005,0.00005];
431 plot_detailed_data(model,detailed_data,[]);
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);
439 % model.enable_error_estimator = 1;
440 model.enable_error_estimator = 1;
447 % fn = fullfile(rbmatlabtemp,'follicle_PODGreedy_detailed_data.mat');
448 fn = fullfile(rbmatlabtemp,'follicle_rect_PODGreedy_detailed_data_proj_err_9x9_large.mat');
450 error('please run step 3* (or any similar) first!!');
453 % model_data.RB = model_data.RB(:,1:2); => Acceleration 1000 !!
454 % model_data.RB = model_data.RB(:,1:20);
456 model.enable_error_estimator = 0;
458 plot_detailed_data(model,detailed_data,[]);
461 case 3.8 % generate BT export matrices
465 params.output_range = 50;
466 % params.output_plot_indices = [1,2];
468 model.clim = [0,0.00005];
469 model_data = gen_model_data(model);
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, ...
476 fullfn2 = fullfile(rbmatlabtemp, ...
480 % error(['file ',fullfn,' existing, please remove.']);
484 % error(['file ',fullfn2,' existing, please remove.']);
487 % matrix export for remote BT basis computation
488 model.decomp_mode = 0;
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]);
496 diary(fullfile(rbmatlabtemp,fn2));
497 disp( [' load(''',fn,'''); '...
498 ' sys = ss(A,B,C,D); ',...
500 ' [sysb,G,T,Tinv] = balreal(sys); ', ...
502 ' W = transpose(T(1:500,:)); ',...
503 ' V = Tinv(:,1:500); ',...
504 ' H = G(1:500); ',...
505 ' save([''~/',fn3,'''],''V'',''W'',''H'',''time'');']);
507 disp(['BT script exported to file ',fullfn2]);
509 disp(['please copy files to remote host into your matlab startup' ...
510 ' directory and execute'])
511 disp([' matlab -nodesktop -nosplash -r ',fn2]);
513 disp(['then copy remote file ~/',fn3,' to local ',rbmatlabtemp,' and press enter'])
517 tmp = load(fullfile(rbmatlabtemp,fn3));
518 if isfield(tmp,'V') && isfield(tmp,'W')
521 error('error in file generation');
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');
541 error('please run step 3.4 (or similar) first!!');
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);
551 % solve optimization problem: target output curve for mu = (0.01,0.01)
553 % load target output curve
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';
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];
570 % r = rand(1,2)*0.0005
572 % 0.063493408146753 0.456687928069510
573 % mu0 = [0.019,0.0055] ;
575 % mu = [0.019063493408147 0.005956687928070]
577 model = model.set_mu(model,mu);
578 % mu = model.get_mu(model);
580 Ytarget = sim_data.Y;
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');
591 error('please run basis generation step 3.* first!!');
594 % the following is required, if another output functional is desired
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;
603 Ytarget_check = sim_data_check.Y;
605 if ~isequal(Ytarget,Ytarget_check)
606 error('model and stored data inconsistent!');
608 % model.C_function_ptr = model_old.C_function_ptr;
610 % model.N = floor(size(detailed_data.V,2)/2);
612 model.N = size(detailed_data.V,2);
613 reduced_data = gen_reduced_data(model,detailed_data);
615 % set up optimization function
616 fun = @(mu) my_opt_target(mu,Ytarget,model,reduced_data)
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);
626 for i = 1:length(X(:))
627 Z(i) = fun([X(i),Y(i)]);
633 set(gca,'Zscale','log')
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);
643 zl = get(gca,'Zlim'); % slightly stretch z-axes
644 zl = 1.1*(zl-sum(zl))+sum(zl);
646 saveas(gcf,fullfile(rbmatlabhome,'optimization_target.eps'),'epsc');
647 saveas(gcf,fullfile(rbmatlabhome,'optimization_target.png'),'png');
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
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,[],[], ...
671 disp('the obtained optimal mu value is: ')
674 disp('And target value J(mumin)')
677 disp('The original mu value is')
680 disp('relative error in mu value:')
681 disp(norm(mu-mumin)/norm(mu));
683 disp('And target value J(mu)')
686 % disp('optimization finished, please inspect workspace');
689 %%%%%%%%%%%%%%%%%%%%%%%%%% single parameter optimization
691 disp('starting 1d-optimization')
692 % set up optimization function
693 fun = @(mu1) my_opt_target_single_param(mu1,Ytarget,model,reduced_data)
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);
701 Z = zeros(size(X),1);
702 for i = 1:length(X(:))
708 plot(X,Z,'Linewidth',2);
709 set(gca,'Yscale','log')
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');
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
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,[],[], ...
739 disp('the obtained optimal mu1 value is: ')
742 disp('And target value J(mumin,0.01)')
745 disp('The original mu1 value is')
748 disp('relative error in mu value:')
749 disp(norm(mu(1)-mumin)/norm(mu(1)));
751 disp('And target value J(mu)')
754 % results of optimization:
757 % 0.010525505381718 0.010279932143424
758 %And target value J(mumin)
759 %J(0.010526 0.01028).00011691
761 % 1.169063479811699e-004
762 %The original mu value is
764 % 0.010000000000000 0.010000000000000
765 %relative error in mu value:
767 %And target value J(mu)
768 %J(0.01 0.01).00019323
770 % 1.932325185788862e-004
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]
781 %the obtained optimal mu1 value is:
784 %And target value J(mumin)
785 %J(0.0099999 0.01).00016254
787 %The original mu1 value is
789 %relative error in mu value:
791 %And target value J(mu)
792 %J(0.01 0.01).00019323
795 case 5.1 % new experiments for revised paper
797 % solve optimization problem: target output curve for mu = (0.01,0.01)
800 % load target output curve
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';
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];
817 % r = rand(1,2)*0.0005
819 % 0.063493408146753 0.456687928069510
820 % mu0 = [0.019,0.0055] ;
822 % mu = [0.019063493408147 0.005956687928070]
824 model = model.set_mu(model,mu);
825 % mu = model.get_mu(model);
827 Ytarget_clean = sim_data.Y;
829 % set random generator to deterministic value, noise 1% relative
831 RandStream.setDefaultStream(RandStream('mt19937ar','seed',12345));
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);
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');
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');
859 error('please run basis generation step 3.* first!!');
862 % the following is required, if another output functional is desired
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;
871 Ytarget_check = sim_data_check.Y;
873 if ~isequal(Ytarget_clean,Ytarget_check)
874 error('model and stored data inconsistent!');
876 % model.C_function_ptr = model_old.C_function_ptr;
878 % model.N = floor(size(detailed_data.V,2)/2);
880 model.N = size(detailed_data.V,2);
881 reduced_data = gen_reduced_data(model,detailed_data);
883 % set up optimization function
884 fun = @(mu) my_opt_target(mu,Ytarget,model,reduced_data)
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);
894 for i = 1:length(X(:))
895 Z(i) = fun([X(i),Y(i)]);
901 set(gca,'Zscale','log')
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);
911 zl = get(gca,'Zlim'); % slightly stretch z-axes
912 zl = 1.1*(zl-sum(zl))+sum(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')
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
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,[],[], ...
943 disp('the obtained optimal mu value is: ')
946 disp('And target value J(mumin)')
949 disp('The original mu value is')
952 disp('relative error in mu value:')
953 disp(norm(mu-mumin)/norm(mu));
955 disp('And target value J(mu)')
958 % disp('optimization finished, please inspect workspace');
961 %%%%%%%%%%%%%%%%%%%%%%%%%% single parameter optimization
963 disp('starting 1d-optimization')
964 % set up optimization function
965 fun = @(mu1) my_opt_target_single_param(mu1,Ytarget,model,reduced_data)
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);
973 Z = zeros(size(X),1);
974 for i = 1:length(X(:))
980 plot(X,Z,'Linewidth',2);
981 set(gca,'Yscale','log')
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')
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
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,[],[], ...
1012 disp('the obtained optimal mu1 value is: ')
1015 disp('And target value J(mumin,0.01)')
1018 disp('The original mu1 value is')
1021 disp('relative error in mu value:')
1022 disp(norm(mu(1)-mumin)/norm(mu(1)));
1024 disp('And target value J(mu)')
1028 % model = set_mu(model,[mumin,mu(2)]);
1029 % rb_sim_data = rb_simulation(model,reduced_data);
1030 % Ytarget_recovered = rb_sim_data.Y;
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);
1040 case 6 % generate trajectory images for paper
1042 fn = fullfile(rbmatlabtemp,'follicle_rect_PODgreedy_detailed_data_proj_err_9x9_large_sorted.mat');
1046 % end time: 10 mio sec. = approx 19 years
1047 % dimension: 22801 = (151)^2
1048 % number time steps: 500
1049 % output: expectation_N
1051 % plot of some trajectories of model:
1052 % 3 snapshots over time for 2 parameters
1053 % and one output plot
1055 mus = {[0.02, 0.005],[0.005,0.02]};
1056 tis = [0,1,250,499];
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];
1064 plot_params.colormap = model.transformed_colormap;
1066 for i = 1:length(mus)
1068 model = model.set_mu(model,mu);
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;
1078 plot_params.show_colorbar = 0;
1090 model.plot_slice(model,model_data,sim_data,plot_params);
1096 if plot_params.show_colorbar
1097 set(gca,'Fontsize',20)
1103 set(gca,'Fontsize',20)
1104 xl = get(gca,'Xlabel');set(xl,'Fontsize',25);
1110 set(gca,'Fontsize',20)
1111 yl = get(gca,'Ylabel');set(yl,'Fontsize',25);
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');
1124 case 7 % generate further images for paper: basis vector plots
1126 % generate plots of basis vectors
1128 fn = fullfile(rbmatlabtemp,'follicle_rect_PODgreedy_detailed_data_proj_err_9x9_large_sorted.mat');
1130 % load(fullfile(rbmatlabtemp,'follicle_temp_detailed_data'));
1131 model_data = gen_model_data(model);
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);
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);
1146 saveas(gcf,fullfile(rbmatlabhome,['basis_vector_',num2str(i),'.eps']),'epsc');
1147 saveas(gcf,fullfile(rbmatlabhome,['basis_vector_',num2str(i),'.png']),'png');
1151 case 8 %generate further plots: error decay curve and parameter selection
1153 fn = fullfile(rbmatlabtemp,'follicle_rect_PODgreedy_detailed_data_proj_err_9x9_large_sorted.mat');
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');
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),'.', ...
1176 set(gca,'Xlim',[0.004,0.021]);
1177 set(gca,'Ylim',[0.004,0.021]);
1178 title('PODGreedy selected parameters');
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');
1190 error('step unknown');
1193 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1195 function J = my_opt_target(mu,Ytarget,model,reduced_data);
1196 % function J = my_opt_target(mu,Ytarget,model,reduced_data);
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)]);
1204 function J = my_opt_target_single_param(mu1,Ytarget,model,reduced_data);
1205 mu = model.get_mu(model);
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)]);
function model = follicle_rect_model(params)
model of the human follicle growth
function varargout = filecache_function(funcptr, varargin)
function used for file-caching other function calls.
function demo_rb_gui(varargin)
reduced basis demo with sliders
function p = plot_sim_data(model, model_data, sim_data, plot_params)
function performing the plot of the simulation results as specified in model.
function res = follicle_rect_experiments(step, model, model_data, detailed_data)
experiments with the model of the human follicle growth
Customizable implementation of an abstract greedy algorithm.