follicle_experiments.m

Go to the documentation of this file.

function res = follicle_experiments(step,model,model_data,detailed_data)
%function follicle_experiments(step[,model,model_data,detailed_data])
% experiments with the model of the human follicle growth
%
% model range is n=200, noutput = 50, nt = 500, time = 19.02 years
% implicit Euler time-discretization, accept relatively large error at
% initial time, for large overall time interval and few timesteps
% and "interactivity"
%
% (step = 1  simulation and visualization of detailed dynamical
%            system, using filecaching => do not call directly)
% step = 1.1 simulation and visualization of detailed dynamical
%            system  => n=20,n_output=20 seems too small 
% step = 1.2 simulation and visualization of detailed dynamical
%            system  => n=50, n_output=20 seems better
% step = 1.3 simulation and visualization of detailed dynamical
%            system  => n=50, n_output=50 seems better, no big
%            difference to 1.2
% step = 1.4 simulation and visualization of detailed dynamical
%            system  => n=70, n_output=50 no big difference to
%            1.2, 1.3
% step = 1.5 simulation and visualization of detailed dynamical
%            system  => n=200, n_output=50 seems better
%            takes about 2 minutes. Results saved for step 1.6
%            parameter variation explicit by setting mu.
% step 1.6:  check influence of different nt on accuracy
%            => relative error to nt=1000 pointwise <= 4% so OK for me     
% (step = 1.7 
%            load precomputed example for range = 200
%            show second heap around 70x70 by suitable color
%            scaling. Generate Movie of data => Does not yet work!!!)
% NEW: step = 1.8, n=200, n_output = 50, initial data distributed
%            over larger domain init_range = 20 
%            and different output functional 'expectation_M'
% (step = 2   reduced simulation and visualization of approximation
%            => do not call directly)
% step = 2.1 n=20, n_output = 20, full reduced basis. should be
%            identical to 1.1
% (step = 3   detailed and reduced simulation and visualization of
%            error and estimator
%            => do not call directly)
% step = 3.1 n=20, n_output = 20, full reduced basis. error should
%            be almost 0, but therefore large overestimation (but
%            still small estimator)
% step = 3.2 n=20, n_output = 20, partial reduced basis (stupid basis)
%            (first 136 states). error should
%            be larger than 0, but overestimation reduced only
%            factor 3
% step = 3.3 n=20, n_output = 20, partial reduced basis
%            (lower 136 triangular states). error should
%            be larger than 0, but overestimation reduced only
%            factor 3. Error & estimator less than 3.2, so more
%            clever basis
% step = 3.4, n=200, n_output = 50, generate POD basis of
%            single trajectory
%            basis is saved for use in further steps
%            and plotted
% step = 3.41, n=200, n_output = 50, generate POD basis of
%            9 trajectories
%            basis is saved for use in further steps
%            and plotted
% step = 3.42, n=200, n_output = 50, generate POD basis of
%            25 trajectories
%            basis is saved for use in further steps
%            and plotted
% step = 3.5 experiments with a part of POD-basis from 3.4
% step = 3.6 n=200, n_output = 50, generate POD-Greedy basis with 
%            true error as indicator, 5x5 p-points logarithmically 
%            equidistant distributed
% NEW: case 3.62 % search for several mu largest required regularization parameter
% step = 3.7 test POD-greedy-basis from step 3.6

% step = 3.8 generate system matrices and export to file for
%             remote BT generation. System A = 20301x20301
%             => ss = sys(A,B,C,D) => 3 GB system ss!!!!!
%             balreal takes about ...
% step 4.0: demo_rb_gui of POD basis from 3.4
% step 5.0: Optimization with RB model with POD basis from 3.4
%
% some steps explicitly require a model and model_data 
%
% Done:
% check different dynamics in parameter range => 1.5 with differing mu
% check correct nt size: double and error determination
% visualization of basis
% generate multiple trajectory POD basis
% distributed initial data -> 1.8
% choice of output functions: 'average','Ex','Ey' -> select in model
%
% Done 7.9.2011:
%   (again) generate single trajectory POD basis => 3.4 (15 minutes)
%    generate POD-Greedy, true-error basis overnight, epsilon 0
%         => 3.6 running => after 2 hours early stop: 
%            did not find PCA extension!
%    generate POD-Greedy, projection-error basis overnight 9x9, epsilon 0
%         => 3.63 remote on am100, 21 Basis vectors, good basis!
%    opitmization with basis and "expectation N"
%          => 4.0. landscape: perfect: error 0.268 very localized
%             optimum value: instead of mu= (0.01,0.01)
%             J(0.01009    0.010123).16455
%             sehr schoen
% make demo_rb_gui work
% solve optimization problem with reduced model
% profiler matrix assembly, double assembly problem
%   => OK, solved, model initialization now without model_data.
% correct output functional constant for separatrix :
% 
% Todo:
% generate POD multiple trajectory basis => 3.41, now 25 trajectories
%         => hangup... => am100 started 8.9.2011, 17:10h
% test different bases
% check stability issue => generate case, that demonstrates problem
% compute output functional constant for expectations

% select best basis for optimization
% check error estimator in first step: is implicit computation of
%       error estimator OK? Why not 0? terrible overestimation
% check possible size for BT from control toolbox: adjust n small enough!!!!!
%  => N= 70 works, but strange results
% generate BT basis
% compare quality of BT to POD on single trajectory
% compare quality of BT to POD on parameter corners
%
% movie of default trajectory => moviegen does not work
% step = 3.9 compute BT-basis (compute remotely on am100)
% step = 3.9 test BT-basis (computed remotely on am100)
%  - POD-Greedy basis using error estimator
%  - POD-Greedy basis by adaptive grid

res = [];

% Bernard Haasdonk 21.4.2010

if (nargin<1) || (isempty(step))
  step = 1.1 % simulation and visualization of detailed dynamical system
           %step = 2 % reduced simulation and visualization of approximation
           %step =3 % detailed and reduced simulation and visualization of
           % error and estimator
end;

if (nargin < 2)  
  model = [];
end;

if (nargin < 3)  
  model_data = [];
%  model = follicle_model([]);
%  model_data = gen_model_data(model);
%  model = follicle_model([]);
%  model_data = gen_model_data(model);
%  model_data.RB = eye(model.dim_x,231); % i.e. first 231 states 
%  model.RB_generation_mode = 'from_detailed_data';
end;

if (nargin < 4)
  detailed_data = [];
end;

%switch basis
% case 1
%  model_data.RB = eye(model.dim_x,231); % i.e. first 231 states 
%  model.RB_generation_mode = 'from_detailed_data';
% case 2 
%  model_data.RB = eye(model.dim_x,200); % i.e. first 200 states 
%  model.RB_generation_mode = 'from_detailed_data';
% case 3
%  model_data.RB = eye(model.dim_x,20); % i.e. first 200 states 
%  model.RB_generation_mode = 'from_detailed_data';
% case 4 
% otherwise
%  error('basis not implemented!');
%end;
%model_data.RB = eye(model.dim_x,model.dim_x); % i.e. non-reduced 
%model_data.RB = eye(model.dim_x,200); % i.e. first 200 states 

switch step
 case 1 % simulation and visualization
  tic;  
  if isempty(model_data)
    model_data = gen_model_data(model);
  end;
  sim_data = filecache_function(@detailed_simulation,model,model_data);
%  sim_data = detailed_simulation(model,model_data);
  t = toc;
  mu = model.get_mu(model);
  model.plot_title = ['state probabilities for mu=(',...
                      num2str(mu(1)),',',num2str(mu(2)),')'];
  p = plot_sim_data(model,model_data,sim_data,model);
%  model.plot_sim_data_state(model,model_data,sim_data,model);
  %compute sum for plausibility check of time range/dimension 
%  l1norm = sum(sim_data.X);
%  figure, plot(l1norm); xlabel('time'); title('l1norm of X over time');
  fprintf('\n')
  disp(['time for detailed_simulation :',num2str(t)]);
  if model.debug
    disp(['finished, please interact with data']);  
    keyboard;
  end;
  res = sim_data;
 
 case 1.1 % call step 1 with changed parameters
  params.range = 20;
  params.output_range = 20;
%  params.output_plot_indices = [1];
  model = follicle_model(params);
  model_data = gen_model_data(model);
  follicle_experiments(1,model,model_data);

 case 1.2 % call step 1 with changed parameters
  params.range = 50;
  params.output_range = 20;
%  params.output_plot_indices = [1,2];
  model = follicle_model(params);
  model_data = gen_model_data(model);
  follicle_experiments(1,model,model_data);

 case 1.3 % call step 1 with changed parameters
  params.range = 50;
  params.output_range = 50;
%  params.output_plot_indices = [1,2];
  model = follicle_model(params);
  model_data = gen_model_data(model);
  follicle_experiments(1,model,model_data);
  
 case 1.4 % call step 1 with changed parameters
  params.range = 70;
  params.output_range = 50;
  %  params.output_plot_indices = [1,2];
  model = follicle_model(params);
  model_data = gen_model_data(model);
  follicle_experiments(1,model,model_data);

 case 1.5 % call step 1 with changed parameters
  tic;
  params.range = 200;
  params.output_range =  50;
  %  params.output_plot_indices = [1,2];
  model = follicle_model(params);
  model.clim = [0,0.00005];
  % model_data = gen_model_data(model);
  % mu = (u1, u2), mu_ranges : {[0.005,0.02],[0.005,0.02]};
%  mu = [0.005,0.005]; % => blob to 3e-6, yout to 0.3
%  mu = [0.01,0.01]; %   => blob to 1e-6, yout to 0.1
%  mu = [0.02,0.02]; %   => blob to 0, yout to almost 0, left corner visible
%  mu = [0.005,0.02]; %  => blob to 0, yout to almost 0, left corner invisible
  mu = [0.02,0.005]; %  => blob non existing, yout 0.9, linear decay
                    
  model = model.set_mu(model,mu);
  %  model.save_time_indices = 0:(model.nt/100):model.nt;  
  %model.save_time_indices = 1:100:1001;
  %model.nt = 1000;
  %model.T = 1000*model.nt;
%  sim_data = follicle_experiments(1,model,model_data);
  sim_data = follicle_experiments(1,model);
  save(fullfile(rbmatlabtemp,'follicle_n200_example.mat'),...
       'sim_data','model','model_data');
  t = toc;
  disp(['computation time t = ',num2str(t),' sec.']);
  keyboard;
  
 case 1.6 
  % test nt behaviour

  params.range = 200;
  params.output_range =  50;
  %  params.output_plot_indices = [1,2];
  model = follicle_model(params);
  if isempty(model_data)
    model_data = gen_model_data(model);
  end;
  model.clim = [0,0.00005];

  % nt = 2000, return every 4th slice
  % nt = 1000, return every second slice
  % nt = 500

  disp('nt = 500');
  mu = model.get_mu(model);
  sim_data500 = filecache_function(@detailed_simulation,model, ...
                                   model_data);
%  model.plot_title = ['state probabilities for nt = 500'];
%  plot_state_probabilities(model,model_data,sim_data500,model);
  
  disp('nt = 1000');
  model.nt = 1000;
  model.save_time_indices = 0:2:1000;
  sim_data1000 = filecache_function(@detailed_simulation,model, ...
                                model_data);
%  model.plot_title = ['state probabilities for nt = 1000'];
%  plot_state_probabilities(model,model_data,sim_data1000,model);

  sim_data_err = sim_data500;
  sim_data_err.X = abs(sim_data_err.X- sim_data1000.X);
  model.plot_title = ['absolute value of error'];
%  plot_state_probabilities(model,model_data,sim_data_err,model);

  nz = find(abs(sim_data1000.X)>0);
  z = find(sim_data1000.X==0);
  sim_data_err_rel = sim_data500;
  sim_data_err_rel.X(nz)= sim_data_err.X(nz)./abs(sim_data1000.X(nz)); 
  sim_data_err_rel.X(z) = sim_data_err.X(z) * 1e10; 
  model.plot_title = ['relative error'];
  model.clim =[0,0.1];
  
  model.plot_sim_data_state(model,model_data,sim_data_err_rel,model);

  keyboard;
  
 case 1.7
   %generate movie from 1.6

   % currently problem in movie generation: mpgwrite produces
   % non-readable file, explorer crash, matlab-dll warning, etc.
   % perhaps use Markus method, see basisgen_animation
   
   load(fullfile(rbmatlabtemp,'follicle_n200_example.mat'));
   model.clim = [0,0.00005];
   outputfn = fullfile(rbmatlabtemp,'follicle_trajectory_small2.mpg');
   if exist(outputfn)
     error(['file ',outputfn,' existing!']);
   end;
   fps = 15;
   %params.show_colorbar = 1
   params.xnumintervals = model.range+1;
   params.ynumintervals = model.range+1;
   params.xrange = [-0.5,model.range+0.5];
   params.yrange = [-0.5,model.range+0.5];
   params.verbose = 0;
   grid = rectgrid(params);
   model.no_lines = 1;
   %grid = Grids.rectgrid(params);
   plot_state_probability_slice(model,model_data,grid,sim_data.X(:,1),model);
%   set(gcf,'Position',[10,10,300,200]);
   set(gcf,'Position',[10,10,800,800]);
   for j = 1:10
%   for j = 1:size(sim_data.X,2)
     cla;
     plot_state_probability_slice(model,model_data,grid,sim_data.X(:,j),model);
%     text(0.1,0.9,text_sequence{i},'FontSize',30,'Color',[0,0,1])
     F(j) = getframe(gca);
   end;
   
   C = colormap;
   keyboard;
   % play movie
   movie(F);
   %mov = close(mov);
   mpgwrite(F,C,outputfn,...
            [1, 0, 1, 0, 10, 1, 1, 1]);

   keyboard;
   
 case 1.8    
  tic;
  params.range = 200;
  params.output_range =  50;
  %  params.output_plot_indices = [1,2];
  params.init_range = 20;
%  params.output_type = 'separatrix';
  params.output_type = 'expectation_N';
%  params.reg_epsilon = 0.00009;
  params.reg_epsilon = 0.00000;
  model = follicle_model(params);
  model.clim = [0,0.00005];
  % model_data = gen_model_data(model);
  % mu = (u1, u2), mu_ranges : {[0.005,0.02],[0.005,0.02]};
%  mu = [0.005,0.005]; % => blob to 5e-5, yout to 3.7
%  mu = [0.01,0.01]; %   => blob to 5e-5, yout to 3.8
%  mu = [0.02,0.02]; %   => blob to 0, yout to almost 3.7,
%  mu = [0.005,0.02]; %  => blob to 0, yout to almost 9.2
  mu = [0.02,0.005]; %  => blob non existing, yout 0.3
                    
  model = model.set_mu(model,mu);
  %  model.save_time_indices = 0:(model.nt/100):model.nt;  
  %model.save_time_indices = 1:100:1001;
  %model.nt = 1000;
  %model.T = 1000*model.nt;
%  sim_data = follicle_experiments(1,model,model_data);
  sim_data = follicle_experiments(1,model);
  save(fullfile(rbmatlabtemp,'follicle_n200_example.mat'),...
       'sim_data','model','model_data');
  t = toc;
  disp(['computation time t = ',num2str(t),' sec.']);
  keyboard;
   
 case 2 % reduced simulation
  %define reduced basis of first two coordinates (i.e. error only
  %by u!)

  model.debug = 1;

  %    produce high-dimensional data, e.g. reduced basis, etc.
  %    mu is not known
  detailed_data     = gen_detailed_data(model,model_data); 
  disp('gen_detailed_data successful');
  
  %    produce low-dimensional data, e.g. gram matrices between
  %    reduced basis vectors, subgrid extraction, etc.
  %    mu is not known (should comprise previous offline and online steps)

  reduced_data = gen_reduced_data(model,detailed_data);
  disp('gen_reduced_data successful');
  model.N = reduced_data.N;
  
  % rb_sim_data  = rb_simulation(reduced_data,model)
  %    perform reduced simulation, i.e. mu known
  %    produces some structure with suitable fields, e.g DOF vector
  
  disp('entering reduced simulation');
  rb_sim_data = rb_simulation(model,reduced_data);
  disp('reduced simulation successful');
  
  rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);

  model.plot_title = 'reduced trajectory';
  %model.plot_indices = 1:3;   => default
  p = plot_sim_data(model,model_data,rb_sim_data,[]);

  disp(['finished, please interact with data']);  
  keyboard;

 case 2.1
  par.range = 20;
  par.output_range = 20;  
  model = follicle_model(par);
  model_data = gen_model_data(model);
  model_data.RB = eye(model.dim_x,231); % i.e. first 231 states 
  model.RB_generation_mode = 'from_detailed_data';
  follicle_experiments(2,model,model_data);  
   
 case 3 % detailed and reduced simulation
  
  %model.x0_shift=0;
  %model.u_function_ptr = @(model) 0; % define u constant 0
  % => leads to error 0 :-) if x0_shift = 0
  
  %model.u_function_ptr = @(model) 0.1; % define u constant 
  %=> error linearly increasing
  %model.u_function_ptr = @(model) 0.3*sin(2*model.t)+0.002*model.t;     

  tic;
  
  % remove some fields such that caching works
  tmp_model_data = model_data;
  tmp_model = model;
  if isfield(tmp_model_data,'RB')
    tmp_model_data = rmfield(tmp_model_data,'RB');
  end;
  if isfield(tmp_model_data,'RB_generation_mode')
    tmp_model = rmfield(tmp_model,'RB_generation_mode');
  end;
  sim_data = filecache_function(@detailed_simulation,tmp_model,tmp_model_data);
  t1 = toc;
  fprintf('\n');
  disp(['time for detailed simulation: ',num2str(t1)]);

%  model.plot_title = 'full simulation';
%  model.plot_sim_data_state(model,model_data,sim_data,model);

  if isempty(detailed_data)
    detailed_data     = gen_detailed_data(model,model_data); 
  end;
  
  reduced_data = gen_reduced_data(model,detailed_data);
  model.N = reduced_data.N;
  tic;
  rb_sim_data = rb_simulation(model,reduced_data);
  t2 = toc;
  disp(['time for reduced simulation: ',num2str(t2)]);
  
  rb_sim_data = rb_reconstruction(model,detailed_data, ...
                                         rb_sim_data);

%  model.plot_title = 'reduced simulation';
%  model.plot_sim_data_state(model,model_data,rb_sim_data,model);

%  figure, plot3(sim_data.X(1,:),...
%               sim_data.X(2,:),...
%               sim_data.X(3,:), ...
%                'color','b' ...
%               );
%  hold on
%  plot3(rb_sim_data.X(1,:),...
%               rb_sim_data.X(2,:),...
%               rb_sim_data.X(3,:),...
%               'color','r');
%  axis equal;
%  title('exact and reduced trajectories');
%  legend({'exact','reduced'});

  er = sim_data;
  er.X = sim_data.X-rb_sim_data.X;
  er.Y = sim_data.Y-rb_sim_data.Y;

  model.plot_title = 'error';
  model.plot_sim_data_state(model,model_data,er,model);

  err = sqrt(sum((detailed_data.G * er.X).*er.X,1)); 
  max_t_err = err;
  
  % somehow replace by matrix version!!!
  for i = 2:length(max_t_err)
    if max_t_err(i)<max_t_err(i-1)
      max_t_err(i) = max_t_err(i-1);
    end;
  end;
  
  %  figure, plot(sim_data.time,[rb_sim_data.DeltaX(:),max_t_err(:),err(:)]), 
  if model.enable_error_estimator
    figure, plot([rb_sim_data.DeltaX(:),max_t_err(:),err(:)]), 
    Ylim = get(gca,'Ylim');
    set(gca,'Ylim',[0,max(Ylim)]);
    title('rb state error and estimator');
    legend({'estimator \Delta_x(t)','||e||_{L\infty([0,t],L2)}',...
            '||e(t)||_{L2}'});
  else
    figure, plot([max_t_err(:),err(:)]), 
    Ylim = get(gca,'Ylim');
    set(gca,'Ylim',[0,max(Ylim)]);
    title('rb state error');
    legend({'||e||_{L\infty([0,t],L2)}',...
            '||e(t)||_{L2}'});
  end;
  
  if model.dim_y > 1
    er.Ynorm = sum(er.Y.^2)';
  else
    er.Ynorm = abs(er.Y(:));
  end;
  
  if model.enable_error_estimator
    figure, plot([rb_sim_data.DeltaY(:),er.Ynorm]); 
    title('rb output error and estimator');
    Ylim = get(gca,'Ylim');
    set(gca,'Ylim',[0,max(Ylim)]);
    legend({'estimator','error'});
    figure, plot(rb_sim_data.time_indices,rb_sim_data.Rnorms);
    title('Residual norms')
  else
    figure, plot([er.Ynorm]); 
    title('rb output error');
    Ylim = get(gca,'Ylim');
    set(gca,'Ylim',[0,max(Ylim)]);
    legend({'error'});
  end;

%  plot_sim_data(model,model_data,sim_data,model);
  if model.debug
    disp(['finished, please interact with data']);  
    keyboard;
  end;

 case 3.1
  
  model = follicle_model([]);
  model_data.RB = eye(model.dim_x,231); % i.e. first 231 states 
  model.RB_generation_mode = 'from_detailed_data';
  follicle_experiments(3,model,model_data);  

 case 3.2
  
  model = follicle_model([]);
  % length i = 136 (range <=15)
  E = eye(model.dim_x,136); 
  model_data.RB = E;
  model.RB_generation_mode = 'from_detailed_data';
  follicle_experiments(3,model,model_data);  
 
 case 3.3

  par = [];
  par.output_type = 'expectation_N';
  model = follicle_model(par);
  model = model.set_mu(model,[0.01,0.01]);
  model_data = gen_model_data(model);
  % length i = 136 (range <=15)
  i = find((model_data.Ns+model_data.Ms)<=15); % i.e. lower states
                                               
  E = eye(model.dim_x); 
  model_data.RB = E(:,i);
  model.RB_generation_mode = 'from_detailed_data';
  follicle_experiments(3,model,model_data);  
  keyboard;

 case 3.31
  fn = fullfile(rbmatlabtemp,'follicle_PODGreedy_detailed_data_proj_err_9x9.mat');
  if ~exist(fn)
    error('please run step 3.6 (or any similar) first!!');
  end;
  load(fn);
  % length i = 136 (range <=15)
  %i = find((model_data.Ns+model_data.Ms)<=15); % i.e. lower states
                                               
  model_data.RB = detailed_data.V;
  model.RB_generation_mode = 'from_detailed_data';
  follicle_experiments(3,model,model_data);  

 case 3.4
  % generate POD basis of single trajectory
  params.range = 200;
  params.output_range =  50;
  %  params.output_plot_indices = [1,2];
  params.init_range = 20;
%  params.output_type = 'separatrix';
  params.output_type = 'expectation_N';
%  params.reg_epsilon = 0.00009;
  params.reg_epsilon = 0.00000;
  model = follicle_model(params);
  %  model.RB_generation_mode = 'PCA_trajectory';
  %  model.RB_stop_Nmax = 10;
  model_data = gen_model_data(model);
  model.RB_generation_mode = 'from_detailed_data';
  tmpmodel = model;
  tmpmodel.save_time_indices = 0:model.nt; 
  disp('computing trajectory');
  sim_data = filecache_function(@detailed_simulation,tmpmodel,model_data);
  disp('computing  POD');
  RBinit = orthonormalize_gram_schmidt(sim_data.X(:,1),model_data.G,1e-20);
  Xorth = sim_data.X(:,2:end)- ...
          RBinit * ((RBinit' * model_data.G) * sim_data.X(:,2:end));
  %  keyboard;
  clear('sim_data');
  [U,S,V,flag] = svds(Xorth,20); % single svds(..20) takes about
                                 % 5 min!!
  i = find(diag(S)>=eps);
  U = U(:,i);
  %[U,S] = svds_batch(sim_data.X,5,100);
  % warning: basis not fully orthogonal!!!
  %  U = orth(U);
  disp('computing  gram-schmidt orthogonalization');
  model_data.RB = orthonormalize_gram_schmidt([RBinit,U],...
                                              model_data.G,eps);
  %  keyboard;
  detailed_data = gen_detailed_data(model,model_data);
  save(fullfile(rbmatlabtemp,'follicle_POD_detailed_data.mat'),...
       'model','model_data','detailed_data');

  model.clim = [-0.00005,0.00005];
  model.plot_detailed_data(model,detailed_data,[]);
%  follicle_experiments(3,model,model_data);  

 case 3.41
   % generate POD basis of multiple trajectories
   
   %  model.RB_generation_mode = 'PCA_trajectory';
  %  model.RB_stop_Nmax = 10;
  par.range = 200;
  par.output_range = 50;  
  params.init_range = 20;
%  params.output_type = 'separatrix';
  params.output_type = 'expectation_N';
%  params.reg_epsilon = 0.00009;
  params.reg_epsilon = 0.00000;
  model = follicle_model(par);
  model_data = gen_model_data(model);
  model.RB_generation_mode = 'from_detailed_data';

  u1s = [0.005,0.007,0.01,0.014,0.02];
  u2s = [0.005,0.007,0.01,0.014,0.02];
%  u2s = [0.005,0.01,0.02];

  Us = [];
  Ss = [];
  
  tic;
  for u1i = 1:length(u1s);
    for u2i = 1:length(u2s);
      tmpmodel = model;
      tmpmodel.save_time_indices = 0:model.nt; 
      mu = [u1s(u1i),u2s(u2i)];
      disp(['processing trajectory for mu=',num2str(mu)])
      tmpmodel = tmpmodel.set_mu(tmpmodel,mu);
      disp('computing trajectory');
      sim_data = filecache_function(@detailed_simulation,tmpmodel,model_data);
      disp('computing  POD');
      RBinit = orthonormalize_gram_schmidt(sim_data.X(:,1),model_data.G,1e-20);
      Xorth = sim_data.X(:,2:end)- ...
              RBinit * ((RBinit' * model_data.G) * sim_data.X(:,2:end));
      %  keyboard;
      clear('sim_data');
      [U,S,V,flag] = svds(Xorth,15);
      clear('Xorth');
      i = find(diag(S)>=eps);
      U = U(:,i);
      disp('computing  gram-schmidt orthogonalization');
      U = orthonormalize_gram_schmidt(U,model_data.G,eps);
      S = diag(S);
      Ss = [Ss;S(i)];
      Us = [Us,U];
    end;
  end;
  
  % save temporary result
  save(fullfile(rbmatlabtemp,'partial_PODs'));
  
  % one final SVD merging the partial bases:
  Us = Us*diag(Ss);
  disp('computing  final POD');
  [U,S,V,flag] = svds(Us,min(100,size(Us,2)));
  i = find(diag(S)>=eps);
  Utot = Us(:,i);
  disp('computing final gram-schmidt orthogonalization');
  model_data.RB = orthonormalize_gram_schmidt([RBinit,Utot],model_data.G,eps);
  
  detailed_data = gen_detailed_data(model,model_data);
  save(fullfile(rbmatlabtemp,'follicle_multi_POD_detailed_data.mat'),...
       'model','model_data','detailed_data');

  model.clim = [-0.00005,0.00005];
  model.plot_detailed_data(model,detailed_data,[]);
%  follicle_experiments(3,model,model_data);  
  t = toc;
  disp(['runtime ',num2str(t)]);
  
 case 3.5

  % test POD basis, only part of basis vectors
%  fn = fullfile(rbmatlabtemp,'follicle_multi_POD_detailed_data.mat');
  fn = fullfile(rbmatlabtemp, ...
                'follicle_PODgreedy_detailed_data_non_orth.mat');
  if ~exist(fn)
    error('please run step 3.4 first!!');
  end;
  load(fn);
  model.clim = [-0.00005,0.00005];
  plot_detailed_data(model,detailed_data,[]);
  return;
  %  model_data.RB = model_data.RB(:,1:2); => Acceleration 1000 !!
%  RB = orthonormalize_qr(model_data.RB(:,1:20));
%  RB = orth(model_data.RB(:,1:2));
  RB = model_data.RB(:,1:4);
  model_data.RB = RB;
  model.debug = 1;
%  model.enable_error_estimator = 1;
  model.enable_error_estimator = 1;
  follicle_experiments(3,model,model_data);  
  keyboard;
  
 case 3.6
  
  % POD-Greedy basis with true error as indicator (default settings
  % are in model)
  tic;
  
  par.range = 200;
  par.output_range = 50;
  par.init_range = 20;
  model = follicle_model(par);
  model.RB_generation_mode = 'greedy_log_uniform_fixed';
  model.RB_detailed_train_savepath = 'follicle_model_detailed_train_eps0';
  
  model_data = gen_model_data(model);
  detailed_data = gen_detailed_data(model,model_data);
  save(fullfile(rbmatlabtemp,'follicle_PODGreedy_detailed_data_eps0.mat'),...
       'model','model_data','detailed_data');

  t = toc;
  disp('computation time: ')
  t
%  model.N = 20;
%  follicle_experiments(3,model,model_data);  

 case 3.61
  
  % POD-Greedy basis with true error as indicator (default settings
  % are in model) with regularization 
  par.range = 200;
  par.output_range = 50;  
  par.init_range = 20;
  model = follicle_model(par);
  model.RB_generation_mode = 'greedy_log_uniform_fixed';
  model.RB_detailed_train_savepath = 'follicle_model_detailed_train_eps0.0186714';
%  model.RB_numintervals = [1,1];
%  model.RB_stop_Nmax = 20;
%  model.reg_epsilon = 0.00009;
  model.reg_epsilon = 0.0186714; % maximum eigenvalue over 25 points
%  model.debug = 1;
  model_data = gen_model_data(model);
  detailed_data = gen_detailed_data(model,model_data);
  save(fullfile(rbmatlabtemp,...
                'follicle_PODGreedy_detailed_data_eps0.0186714.mat'),...
       'model','model_data','detailed_data');

%  model.N = 20;
%  follicle_experiments(3,model,model_data);  

 case 3.62 % search for several mu largest required regularization parameter
  par.range = 200;
  par.output_range = 50;  
  par.init_range = 20;
  model = follicle_model(par);
  
  par.numintervals = model.RB_numintervals;
  par.range = model.mu_ranges;
  for i = 1:length(par.range);
    par.range{i} = log(par.range{i});
  end;
  %  enable restart of basis generation!!!  
  %  detailed_data.RB_info = [];
  MMesh0 = cubegrid(par);
  Mtrain_log = get(MMesh0,'vertex')';
  mus = exp(Mtrain_log);
% mus = rand_uniform(nmus,model.mu_ranges);
  model.reg_epsilon = 0.0186714; % maximum eigenvalue over 25 points
  model.decomp_mode = 0;
  model_data = gen_model_data(model);
  lmax = -10000;
%  for i = 1:size(mus,2);
  for i = 25;
    mu = mus(:,i)
    model = model.set_mu(model,mu);
    A = model.A_function_ptr(model,model_data);
    As = 0.5*(A+A');
    l_current = eigs(As-10*speye(size(As)),1,'SM');
    l_current = 10+l_current;
    if l_current > lmax
      lmax = l_current;
      disp(['detected eigenvalue ',num2str(lmax)]);
    end; 
  end;

  
 case 3.63
  
  % POD-Greedy basis with true projection error as indicator (default settings
  % are in model) 9x9 parameter grid
  tic;
  
  par.range = 200;
  par.output_range = 50;
  par.init_range = 20;
  model = follicle_model(par);
  model.RB_generation_mode = 'greedy_log_uniform_fixed';
  model.RB_detailed_train_savepath = 'follicle_model_detailed_train_eps0_9x9';
  model.RB_error_indicator = 'projection-error';
  model.RB_stop_epsilon = 1e-16;
  model.RB_numintervals = [8,8];
  model_data = gen_model_data(model);
  detailed_data = gen_detailed_data(model,model_data);
  save(fullfile(rbmatlabtemp,'follicle_PODGreedy_detailed_data_proj_err_9x9.mat'),...
       'model','model_data','detailed_data');

  t = toc;
  disp('computation time: ')
  t
%  model.N = 20;
%  follicle_experiments(3,model,model_data);  
  
 case 3.7
    
  % test POD-Greedy basis
%  fn = fullfile(rbmatlabtemp,'follicle_PODGreedy_detailed_data.mat');
  fn = fullfile(rbmatlabtemp,'follicle_PODGreedy_detailed_data_proj_err_9x9.mat');
  if ~exist(fn)
    error('please run step 3.6 (or any similar) first!!');
  end;
  load(fn);
%  model_data.RB = model_data.RB(:,1:2); => Acceleration 1000 !!
%  model_data.RB = model_data.RB(:,1:20);
%  model.debug = 1;
  model.enable_error_estimator = 0;
  follicle_experiments(3,model,model_data,detailed_data);  
  plot_detailed_data(model,detailed_data,[]);
  keyboard;

 case 3.8 % generate BT export matrices

  params.range = 200;
%  params.range = 70;
  params.output_range =  50;
  %  params.output_plot_indices = [1,2];
  model = follicle_model(params);
  model.clim = [0,0.00005];
  model_data = gen_model_data(model);
  
  fn = 'follicle_model_system_matrices.mat';
  fn2 = 'follicle_model_BT.m';
  fn3 = 'follicle_model_BT_basis.mat';
  fullfn = fullfile(rbmatlabtemp, ...
                fn);
  fullfn2 = fullfile(rbmatlabtemp, ...
                fn2);
  if exist(fullfn)
    delete(fullfn);
%    error(['file ',fullfn,' existing, please remove.']);
  end;
  if exist(fullfn2)
    delete(fullfn2);
%    error(['file ',fullfn2,' existing, please remove.']);
  end;
  
  % matrix export for remote BT basis computation
  model.decomp_mode = 0;
  A = filecache_function(model.A_function_ptr,model,model_data);
  B = model.B_function_ptr(model,model_data);
  C = model.C_function_ptr(model,model_data);
  D = model.D_function_ptr(model,model_data);
  save(fullfn,'model','A','B','C','D');
  disp(['system matrices exported to file ',fullfn]);

  diary(fullfile(rbmatlabtemp,fn2));
  disp( [' load(''',fn,''');  '...
         ' sys = ss(A,B,C,D); ',...
         ' tic; ', ...
         ' [sysb,G,T,Tinv] = balreal(sys); ', ...
         ' time = toc; ', ...
         ' W = transpose(T(1:500,:)); ',...
         ' V = Tinv(:,1:500); ',...
         ' H = G(1:500); ',...
         ' save([''~/',fn3,'''],''V'',''W'',''H'',''time'');']);
  diary off;  
  disp(['BT script exported to file ',fullfn2]);
  
  disp(['please copy files to remote host into your matlab startup' ...
        ' directory and execute'])
  disp(['  matlab -nodesktop -nosplash -r ',fn2]);
  
  disp(['then copy remote file ~/',fn3,' to local ',rbmatlabtemp,' and press enter'])

  keyboard;

  tmp = load(fullfile(rbmatlabtemp,fn3));
  if isfield(tmp,'V') && isfield(tmp,'W')
    disp('BT basis OK')
  else
    error('error in file generation');
  end;
   
 case 3.9
   % test BT basis
   
 case 4.0 % reduced basis in demo_rb_gui;

  % test POD basis only part of basis vectors
%   fn = fullfile(rbmatlabtemp,'follicle_POD_detailed_data.mat');
%  fn = fullfile(rbmatlabtemp,'follicle_PODgreedy_detailed_data_eps0.0186714.mat');
%fn = fullfile(rbmatlabtemp, ...
%             'follicle_PODGreedy_detailed_data_eps0.0186714.mat');
%fn = fullfile(rbmatlabtemp,'follicle_PODgreedy_detailed_data_eps0.mat');
%fn = fullfile(rbmatlabtemp,'follicle_PODgreedy_detailed_data_proj_err.mat');
fn = fullfile(rbmatlabtemp,'follicle_PODgreedy_detailed_data_proj_err_9x9.mat');
%  fn = fullfile(rbmatlabtemp,'follicle_multi_POD_detailed_data.mat');
%  fn = fullfile(rbmatlabtemp,'follicle_PODGreedy_detailed_data_non_orth.mat');
  if ~exist(fn)
    error('please run step 3.4 (or similar) first!!');
  end;
  load(fn);
  model.N = size(detailed_data.V,2);
  model.clim = [0.00000,0.00005];
  demo_rb_gui(model,detailed_data);
  keyboard;
  
 case 5.0
  % solve optimization problem: target output curve for mu = (0.01,0.01)
  
  % load target output curve
  params.range = 200;
  params.output_range =  50;
  %  params.output_plot_indices = [1,2];
  params.init_range = 20;
  params.output_type='expectation_N';
  model = follicle_model(params);
  model.clim = [0,0.00005];
  model_data = gen_model_data(model);
  mu = [0.01,0.01];
  model = model.set_mu(model,mu);
%  mu = model.get_mu(model);
  sim_data = filecache_function(@detailed_simulation,model,model_data);
  Ytarget = sim_data.Y;

  % load reduced model
%  fn = fullfile(rbmatlabtemp,'follicle_multi_POD_detailed_data.mat');
%  fn = fullfile(rbmatlabtemp,'follicle_POD_detailed_data.mat');
%  fn = fullfile(rbmatlabtemp,'follicle_PODGreedy_detailed_data_hand_orth.mat');
fn = fullfile(rbmatlabtemp,'follicle_PODgreedy_detailed_data_proj_err_9x9.mat');
%  fn = fullfile(rbmatlabtemp,'follicle_PODgreedy_detailed_data_eps0.mat');
%  fn = fullfile(rbmatlabtemp,'follicle_PODGreedy_detailed_data.mat');
  if ~exist(fn)
    error('please run basis generation step 3.* first!!');
  end;
%  keyboard;
  model_old = model;
  model_data_old = model_data;
  load(fn);
  model = model_old;
  model_data.C_components = model_data_old.C_components;
  detailed_data.C_components = model_data_old.C_components;

  sim_data_check = filecache_function(@detailed_simulation,model,model_data);
  Ytarget_check = sim_data_check.Y;

  if ~isequal(Ytarget,Ytarget_check)
    error('model and stored data inconsistent!');
  end;
  %  model.C_function_ptr = model_old.C_function_ptr;
  %  keyboard; 
%  model.N = floor(size(detailed_data.V,2)/2);
%  model.N = 12;
  model.N = size(detailed_data.V,2);
  reduced_data = gen_reduced_data(model,detailed_data);
  
  % set up optimization function
  fun = @(mu) my_opt_target(mu,Ytarget,model,reduced_data)
  
  if 1
    disp('generating surface');
    % plot target surface
    logr = log([0.005,0.02]);
    logrange = logr(1):(logr(2)-logr(1))/20:logr(2);
    range = exp(logrange)
    [X,Y] = meshgrid(range,range);
    Z = zeros(size(X));
    for i = 1:length(X(:))
      Z(i) = fun([X(i),Y(i)]);
%      if (Z(i))>1e10
%       keyboard;
%      end;
    end;
    surf(X,Y,Z);
    set(gca,'Zscale','log')
    keyboard;
  end;
  
%  mu0 = [0.005,0.005];
%  mu0 = [0.02,0.005];
  mu0 = [0.02,0.02];
  lb = [0.005,0.005];
  ub = [0.02,0.02];
  [mumin,fval,exitflag,output,lambda,grad] = fmincon(fun,mu0,[],[],[],[],lb,ub);    

  disp('The original value is')
  
  mu;
  
  disp('And target value J(mu)')
  fun(mu);
  
  disp('optimization finished, please inspect workspace');
  keyboard;
  
 otherwise
  error('step unknown');
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function J = my_opt_target(mu,Ytarget,model,reduced_data);
% function J = my_opt_target(mu,Ytarget,model,reduced_data);
% my_opt_target

tmpmodel = model.set_mu(model,mu);
rb_sim_data = rb_simulation(tmpmodel,reduced_data);
J = sum((rb_sim_data.Y-Ytarget).^2);
disp(['J(',num2str(mu(:)'),')=',num2str(J)]);