rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
follicle_rect_model.m
Go to the documentation of this file.
1 function model = follicle_rect_model(params)
2 %function model = follicle_model(params)
3 % model of the human follicle growth
4 %
5 % two chemical components are assumed to be characteristic for the
6 % growth of human follicles.
7 %
8 % A deterministic model for the two species is given by a positive
9 % feedback loop, a bistable switch, i.e. 2-d nonlinear ODE.
10 %
11 % An approximation to the Chemical Master Equation gives a linear
12 % dynamical system.
13 %
14 % here, the state vector is the probability of states (pair of number of two
15 % chemical components) and the matrix A describes the change in
16 % these probabilities over time. The number of states is
17 % p = (range+1)*(range+1) hence rectangular domain
18 %
19 % The output is given by params.output_type:
20 % 'separatrix': an average over a subset of states, the ones below
21 % the diagonal "separatrix", which are assumed to be still "in rest",
22 % i.e. waiting for their growth specified by model.output_range
23 % 'expectation_M': the expectation of M state is computed
24 % 'expectation_N': the expectation of N state is computed
25 
26 % Q: D.h. die Zellen oberhalb der separatrix sind einfach abgestorben, oder zur
27 % Ovulation gekommen?
28 %
29 % Explanations from S. Waldherr "A stochastic switch for the
30 % primordial to primary transition in follicle growth", IST
31 % internal report (2009-04-09) v1
32 
33 % Bernard Haasdonk 22.9.2009
34 if (nargin < 1)
35  params = [];
36 end;
37 
38 model = lin_ds_model_default;
39 model.data_const_in_time = 1;
40 model.verbose = 0;
41 model.debug = 0;
42 
43 %%%%%%%%%%%%%%%%%%%%%%% parameters from steffen
44 model.M1 = 25;
45 model.h = 3;
46 model.M2 = 25;
47 %inomi=6;
48 if isfield(params,'range')
49  model.range = params.range;
50 else
51  model.range = 150; % connected with end time, increase e.g. to 50 later
52 end;
53 if isfield(params,'output_range')
54  model.output_range = params.output_range;
55 else
56  model.output_range = 50; % roughly between two bulbs
57 end;
58 if (model.output_range > model.range)
59  error('please choose output range <= range');
60 end;
61 % determines the initial distribution in a lower left triangle
62 if isfield(params,'init_range')
63  model.init_range = params.init_range;
64 else
65  model.init_range = 20;
66 end;
67 %
68 if isfield(params,'reg_epsilon')
69  model.reg_epsilon = params.reg_epsilon;
70 else
71  model.reg_epsilon = 0;
72 end;
73 if ~isfield(params,'output_type')
74  params.output_type = 'separatrix';
75 end;
76 %for i = 1:10
77 % u1 = 0.01+(i-6)*0.001;
78 % u2 = 0.01+(i-6)*0.0005;
79 % k1 = 4*u1+(i-6)*0.0001;
80 % V1 = 75*u1;
81 % V2 = 75*0.01;
82 % parameterlist(i,1:5)= [k1 u1 u2 V1 V2];
83 %end;
84 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
85 
86 %model_data = my_ds_gen_model_data(model);
87 
88 model.u1 = 0.01;
89 model.u2 = 0.01;
90 model.k1fact = 4;
91 model.V1fact = 75;
92 model.V2fact = 75;
93 
94 % used as components later:
95 %model.k1 = model.k1fact * model.u1;
96 %model.V1 = model.V1fact * model.u1;
97 %model.V2 = model.V2fact * model.u2;
98 
99 model.mu_names = {'u1','u2'};
100 %model.mu_names = {'u1','u2','k1fact','V1fact','V2fact'};
101 model.mu_ranges = {[0.005,0.02],[0.005,0.02]};
102 
103 %model.T = 1000000; % Minutes, i.e. 1.9 Jahre
104 %model.nt = 500;
105 %model.T = 1000000; % Minutes, i.e. 1.9 Jahre
106 %model.nt = 10000;
107 model.T = 10000000; % Minutes, i.e. 19 Jahre
108 model.nt = 500;
109 %model.T = 5000000; % Minutes, i.e. 9.5 Jahre
110 %model.nt = 500;
111 %model.T = 50000000; % Minutes, i.e. 95 Jahre
112 %model.nt = 10000;
113 %disp('nt to be adjusted later!!!') => 500 OK!!
114 model.affinely_decomposed = 1;
115 model.A_function_ptr = @A_func;
116 % ratio of first to second ellipses axes
117 %model.A_axes_ratio = 2;
118 %model.A_axes_ratio = 2;
119 model.B_function_ptr = @B_func;
120 switch params.output_type
121  case 'separatrix'
122  model.C_function_ptr = @C_func_separatrix;
123  case 'expectation_M'
124  model.C_function_ptr = @C_func_expectation_M;
125  case 'expectation_N'
126  model.C_function_ptr = @C_func_expectation_N;
127  otherwise
128  error('params.output_type not known');
129 end;
130 model.D_function_ptr = @D_func;
131 model.x0_function_ptr = @x0_func;
132 % shift of x0 in third coordinate
133 %model.x0_shift=1;
134 model.u_function_ptr = @(model) 0; % define u constant
135 %model.u_function_ptr = @(model) sin(4*model.t);
136 %model.u_function_ptr = @(model) 0; % define u constant 0
137 %model.u_function_ptr = @(model) model.t; % define u = 1 * t
138 model.G_matrix_function_ptr = ...
139  @(model,model_data) speye(model.dim_x,model.dim_x);
140 % the following lead to both larger C1 and C2!!!
141 %model.G_matrix_function_ptr = @(model) diag([2,0.5,1]);
142 %model.G_matrix_function_ptr = @(model) diag([0.5,2,1]);
143 
144 % set function pointers to main model methods
145 model.gen_model_data = @my_ds_gen_model_data;
146 %model.detailed_simulation = @lin_ds_detailed_simulation;
147 model.gen_detailed_data = @my_ds_gen_detailed_data;
148 %model.gen_reduced_data = @lin_ds_gen_reduced_data;
149 %model.rb_simulation = @lin_ds_rb_simulation;
150 %model.plot_sim_data = @my_ds_plot_sim_data;
151 model.plot_sim_data_state = @plot_state_probabilities;
152 model.plot_slice = @my_plot_slice;
153 model.axis_tight = 1;
154 
155 %model.rb_reconstruction = @lin_ds_rb_reconstruction;
156 % determin model data, i.e. matrix components, matrix G
157 
158 model.orthonormalize = @model_orthonormalize_qr;
159 
160 %model_data = gen_model_data(model); % matrix components
161 %model.dim_x = size(model_data.X,2);
162 model.dim_x = (model.range+1)*(model.range+1);
163 
164 %model.decomp_mode = 1;
165 %tmp = model.C_function_ptr(model,[]);
166 model.dim_y = 1;
167 %model.dim_y = size(tmp{1},1);
168 %model.decomp_mode = 0;
169 
170 %model.dim_y = size(model_data.C_components{1},1);
171 
172 if isfield(params,'output_plot_indices')
173  model.output_plot_indices = params.output_plot_indices;
174 else
175  model.output_plot_indices = 1:model.dim_y;
176 end;
177 
178 model.theta = 1; % implicit discretization
179 
180 % be sure to determine correct constants once, if anything in the
181 % problem setting has changed!!
182 
183 %model.error_estimation = 1; % model is capable of error_estimation
184 model.enable_error_estimator = 1;
185 if ~isfield(params,'estimate_bounds')
186  params.estimate_bounds = 0;
187 end;
188 % the following only needs to be performed, if the model is changed.
189 if (params.estimate_bounds == 0)
190  % for G = Id
191  switch params.output_type
192  case 'separatrix'
193  model.state_bound_constant_C1 = 1 ;
194  model.output_bound_constant_C2 = sqrt(model.dim_x);
195  case 'expectation_N'
196  model.state_bound_constant_C1 = 1 ;
197  N = model.range;
198  model.output_bound_constant_C2 = 100000;
199  otherwise
200  warning('please adjust the computation of the error estimator constants!')
201  model.state_bound_constant_C1 = 100000 ;
202  model.output_bound_constant_C2 = 100000;
203  end;
204 %elseif params.estimate_bounds == 1
205 % model.estimate_lin_ds_nmu = 3;
206 % model.estimate_lin_ds_nX = 10;
207 % model.estimate_lin_ds_nt = 10;
208 % format long;
209 % [C1,C2] = lin_ds_estimate_bound_constants(model,model_data)
210 % model.state_bound_constant_C1 = C1;
211 % model.output_bound_constant_C2 = C2;
212 else
213  model.estimate_lin_ds_nmu = 3;
214  model.estimate_lin_ds_nX = 10;
215  model.estimate_lin_ds_nt = 10;
216  format long;
217  [C1,C2] = lin_ds_estimate_bound_constants(model,model_data)
218  model.state_bound_constant_C1 = C1;
219  model.output_bound_constant_C2 = C2;
220 end;
221 
222 model.inner_product_matrix_algorithm =@(model,model_data) ...
223  speye(model.dim_x,model.dim_x);
224 
225 model.error_algorithm = @(u1,u2,dummy1,dummy2) ...
226  sqrt(max(sum((u1-u2).^2),0));
227 % [4,3]* [3,2];
228 
229 %model.error_norm = 'l2';
230 %% 'RB_extension_max_error_snapshot'};
231 model.RB_stop_timeout = 12*60*60; % 4 hours
232 model.RB_stop_epsilon = 1e-16;
233 model.RB_error_indicator = 'projection-error'; % true projection error
234 %model.RB_error_indicator = 'estimator'; % Delta from rb_simulation
235 
236 model.RB_stop_Nmax = 200;
237 model.RB_generation_mode = 'greedy_uniform_fixed';
238 model.RB_numintervals = [20,20];
239 model.RB_detailed_train_savepath = ...
240  'follicle_rect_model_detailed_train';
241 model.RB_save_temp_detailed_data = 1;
242 model.orthonormalize = @model_orthonormalize_qr;
243 model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
245 model.transformed_colormap = @my_transformed_colormap;
246 
247 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
248 % functions defining the dynamical system and input/initial data
249 function A = A_func(model,model_data);
250 A = eval_affine_decomp(@A_components,@A_coefficients,...
251  model,model_data);
252 
253 function Acomp = A_components(model,model_data)
254 p = model.dim_x;
255 X = model_data.X;
256 h = model.h;
257 M1 = model.M1;
258 M2 = model.M2;
259 A1 = sparse([],[],[],p,p,7*p); % for k1
260 A2 = sparse([],[],[],p,p,7*p); % for u1
261 A3 = sparse([],[],[],p,p,7*p); % for u2
262 A4 = sparse([],[],[],p,p,7*p); % for V1
263 A5 = sparse([],[],[],p,p,7*p); % for V2
264 A6 = -speye(p,p); % for regularization
265 for m=1:p
266  A1(m,m)= -1;
267  j = find ((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
268  A1(j,m) = 1;
269  A2(m,m) = -X(1,m);
270  j = find((X(1,:)==X(1,m)-1) & (X(2,:)==X(2,m)));
271  A2(j,m) = X(1,m);
272  A3(m,m) = -X(2,m);
273  j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)-1));
274  A3(j,m) = X(2,m);
275  A4(m,m) = -X(2,m)^h/(M1^h+X(2,m)^h);
276  j = find((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
277  A4(j,m) = X(2,m)^h/(M1^h+X(2,m)^h);
278  A5(m,m)= -X(1,m)^h/(M2^h+X(1,m)^h);
279  j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)+1));
280  A5(j,m) = X(1,m)^h/(M2^h+X(1,m)^h);
281 end;
282 %keyboard;
283 Acomp = {A1, A2, A3, A4, A5, A6};
284 
285 function Acoeff = A_coefficients(model)
286 %Acoeff = [model.k1, model.u1, model.u2, model.V1, model.V2];
287 k1 = model.k1fact * model.u1;
288 V1 = model.V1fact * model.u1;
289 V2 = model.V2fact * model.u2;
290 Acoeff = [k1, model.u1, model.u2, V1, V2, model.reg_epsilon];
291 
292 function B = B_func(model,model_data)
293 B = eval_affine_decomp(@B_components,@B_coefficients,model,model_data);
294 
295 function Bcomp = B_components(model,model_data)
296 Bcomp = {zeros(model.dim_x,1)};
297 
298 function Bcoeff = B_coefficients(model)
299 Bcoeff = 1;
300 
301 function C = C_func_separatrix(model,model_data)
302 C = eval_affine_decomp(@C_components_separatrix,...
303  @C_coefficients,model,model_data);
304 
305 function C = C_func_expectation_M(model,model_data)
306 C = eval_affine_decomp(@C_components_expectation_M,...
307  @C_coefficients,model,model_data);
308 
309 function C = C_func_expectation_N(model,model_data)
310 C = eval_affine_decomp(@C_components_expectation_N,...
311  @C_coefficients,model,model_data);
312 
313 function Ccomp = C_components_separatrix(model,model_data)
314 i = find((model_data.Ms+model_data.Ns)<=model.output_range);
315 C = zeros(1,model.dim_x);
316 C(i) = 1;
317 Ccomp = {C};
318 % Ccomp = {ones(1,model.dim_x)};
319 %Ccomp = {[2,0,1],[0,2,1]};
320 
321 function Ccomp = C_components_expectation_M(model,model_data)
322 C = model_data.Ms(:)';
323 Ccomp = {C};
324 
325 function Ccomp = C_components_expectation_N(model,model_data)
326 C = model_data.Ns(:)';
327 Ccomp = {C};
328 
329 function Ccoeff = C_coefficients(model)
330 Ccoeff = [1];
331 %Ccoeff = [model.C_weight,1-model.C_weight];
332 
333 function D = D_func(model,model_data)
334 %D = [0;0]; % no decomposition required for scheme
335 D = [0]; % no decomposition required for scheme
336 
337 function x0 = x0_func(model,model_data)
338 x0 = eval_affine_decomp(@x0_components,@x0_coefficients,model,model_data);
339 
340 function x0comp = x0_components(model,model_data)
341 %if model.init_range == 0
342 % x0comp = {eye(model.dim_x,1)};
343 %else
344 i = ((model_data.Ms+model_data.Ns)<=model.init_range);
345 %i = (model_data.Ms<=1);
346 x0comp = {i(:)/sum(i)}; % sum should be 1!
347 %end;
348 
349 function x0coeff = x0_coefficients(model)
350 x0coeff = 1;
351 
352 function model_data = my_ds_gen_model_data(model)
353 s = 1; % state counter
354 % generate list of states X and output weight vector
355 Ms = zeros(1,(model.range+1)*(model.range+1));
356 Ns = zeros(1,(model.range+1)*(model.range+1));
357 for m=0:model.range
358  for n=0:model.range
359  X(:,s) = [m;n];
360  Ns(s) = n;
361  Ms(s) = m;
362  s = s+1;
363  end;
364 end;
365 p = s-1; % number of states
366 model.dim_x = p;
367 %model_data = lin_ds_gen_model_data(model);
368 model_data.X = X;
369 model_data.Ns = Ns;
370 model_data.Ms = Ms;
371 
372 if model.affinely_decomposed
373  model.decomp_mode = 1;
374  % the following call is extremely expensive: 30 seconds!!!!
375  % therefore caching is activated
376  % model_data.A_components = model.A_function_ptr(model,model_data);
377  model_data.A_components = filecache_function(...
378  model.A_function_ptr,model,model_data);
379  model_data.B_components = model.B_function_ptr(model,model_data);
380  model_data.C_components = model.C_function_ptr(model,model_data);
381  % no D components
382  % model_data.D_components = model.D_function_ptr(model);
383  model_data.x0_components = model.x0_function_ptr(model,model_data);
384 end;
385 model_data.G = model.G_matrix_function_ptr(model,[]);
386 
387 params.xnumintervals = model.range+1;
388 params.ynumintervals = model.range+1;
389 params.xrange = [-0.5,model.range+0.5];
390 params.yrange = [-0.5,model.range+0.5];
391 params.verbose = 0;
392 
393 % only for plotting:
394 model_data.plot_grid = rectgrid(params);
395 model_data.plot_nm_index_vector = ...
396  model_data.Ms + model_data.Ns*(model.range+1)+1;
397 
398 function detailed_data = my_ds_gen_detailed_data(model,model_data)
399 detailed_data = lin_ds_gen_detailed_data(model,model_data);
400 % only for plotting:
401 detailed_data.plot_grid = model_data.plot_grid;
402 detailed_data.plot_nm_index_vector = model_data.plot_nm_index_vector;
403 
404 function p = plot_state_probabilities(model,model_data,sim_data,params)
405 % plot of trajectory and output
406 %p = lin_ds_plot_sim_data(model,model_data,sim_data,params);
407 %close(gcf-1);
408 p = []; % no handles
409 %grid = Grids.rectgrid(params);
410 U = nan((model.range+1)*(model.range+1),size(sim_data.X,2));
411 U(model_data.plot_nm_index_vector,1:end)= sim_data.X;
412 % for debug:
413 %U(nm_index_vector,1) = model_data.Ns;
414 %U(nm_index_vector,2) = model_data.Ms;
415 %params.plot_function = @plot_element_data;
416 merged_plot_params = merge_model_plot_params(model, params);
417 merged_plot_params.plot_function = 'plot_element_data';
418 plot_data_sequence(model,model_data.plot_grid,U,merged_plot_params)
419 axis tight;
420 xlabel('species m');
421 ylabel('species n');
422 if ~isfield(model,'plot_title')
423  title('state probabilities');
424 else
425  title(model.plot_title);
426 end;
427 
428 function p = my_plot_slice(model,model_data,sim_data,params)
429 U = nan((model.range+1)*(model.range+1),1);
430 U(model_data.plot_nm_index_vector)= sim_data.X(:,params.timestep+1);
431 % for debug:
432 %U(nm_index_vector,1) = model_data.Ns;
433 %U(nm_index_vector,2) = model_data.Ms;
434 %params.plot_function = @plot_element_data;
435 params.plot_function = 'plot_element_data';
436 %colormap = jet(2^20);
437 params.colormap = my_transformed_colormap;
438 %keyboard;
439 if ~isfield(params,'clim') && isfield(model,'clim')
440  params.clim = model.clim;
441 end;
442 p = plot_element_data(model_data.plot_grid,U,params);
443 
444 %axis tight;
445 xlabel('species m');
446 ylabel('species n');
447 title('state probabilities');
448 
449 function cm = my_transformed_colormap;
450 colormap = jet(65536);
451 gamma = 0.4;
452 %ind = round((double(0:65535).^gamma)/(65535^gamma)*(2^20-1)+1);
453 ind = round((double(0:255).^gamma)/(255^gamma)*65535+1);
454 cm = colormap(ind,:);
455 
456 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
457 % CODE FROM STEFFEN, RESP HANWEI LI:
458 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
459 
460 
461 % $$$
462 % $$$
463 % $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
464 % $$$ % define parameter
465 % $$$ clear;
466 % $$$ M1 = 25;
467 % $$$ h = 3;
468 % $$$ M2 = 25;
469 % $$$
470 % $$$ for i = 1:10
471 % $$$ u1 = 0.01+(i-6)*0.001;
472 % $$$ u2 = 0.01+(i-6)*0.0005;
473 % $$$ k1 = 4*u1+(i-6)*0.0001;
474 % $$$ V1 = 75*u1;
475 % $$$ V2 = 75*0.01;
476 % $$$ parameterlist(i,1:5)= [k1 u1 u2 V1 V2];
477 % $$$ end;
478 % $$$ inomi=6;
479 % $$$
480 % $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
481 % $$$ % build the nominal model
482 % $$$
483 % $$$ range = 20;
484 % $$$ s = 1; % state counter
485 % $$$
486 % $$$ % generate list of states X
487 % $$$ clear X;
488 % $$$ for m=0:range
489 % $$$ for n=0:(range-m)
490 % $$$ X(:,s) = [m;n];
491 % $$$ s = s+1;
492 % $$$ end;
493 % $$$ end;
494 % $$$ p = s-1; % number of states
495 % $$$
496 % $$$ % defining propensity functions
497 % $$$ a1 = @(x) (parameterlist(inomi,1)+parameterlist(inomi,4)*x(2)^h/(M1^h+x(2)^h));
498 % $$$ a2 = @(x) parameterlist(inomi,2)*x(1);
499 % $$$ a3 = @(x) parameterlist(inomi,5)*x(1)^h/(M2^h+x(1)^h);
500 % $$$ a4 = @(x) parameterlist(inomi,3)*x(2);
501 % $$$
502 % $$$ A = sparse([],[],[],p,p,7*p);
503 % $$$ for m=1:p
504 % $$$ A(m,m) = -(a1(X(:,m))+a2(X(:,m))+a3(X(:,m))+a4(X(:,m)));
505 % $$$ j = find((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
506 % $$$ A(j,m) = a1(X(:,m));
507 % $$$ j = find((X(1,:)==X(1,m)-1) & (X(2,:)==X(2,m)));
508 % $$$ A(j,m) = a2(X(:,m));
509 % $$$ j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)+1));
510 % $$$ A(j,m) = a3(X(:,m));
511 % $$$ j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)-1));
512 % $$$ A(j,m) = a4(X(:,m));
513 % $$$ end;
514 % $$$
515 % $$$
516 % $$$ % KoeffizientenZerlegung
517 % $$$ clear A1; % for k1
518 % $$$ A1 = sparse([],[],[],p,p,7*p);
519 % $$$ for m=1:p
520 % $$$ A1(m,m)= -1;
521 % $$$ j = find ((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
522 % $$$ A1(j,m) = 1;
523 % $$$ end;
524 % $$$
525 % $$$ clear A2; % for u1
526 % $$$ A2 = sparse([],[],[],p,p,7*p);
527 % $$$ for m=1:p;
528 % $$$ A2(m,m) = -X(1,m);
529 % $$$ j = find((X(1,:)==X(1,m)-1) & (X(2,:)==X(2,m)));
530 % $$$ A2(j,m) = X(1,m);
531 % $$$ end;
532 % $$$
533 % $$$ clear A3; % for u2
534 % $$$ A3 = sparse([],[],[],p,p,7*p);
535 % $$$ for m=1:p;
536 % $$$ A3(m,m) = -X(2,m);
537 % $$$ j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)-1));
538 % $$$ A3(j,m) = X(2,m);
539 % $$$ end;
540 % $$$
541 % $$$ clear A4; % for V1
542 % $$$ A4 = sparse([],[],[],p,p,7*p);
543 % $$$ for m=1:p;
544 % $$$ A4(m,m) = -X(2,m)^h/(M1^h+X(2,m)^h);
545 % $$$ j = find((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
546 % $$$ A4(j,m) = X(2,m)^h/(M1^h+X(2,m)^h);
547 % $$$ end;
548 % $$$
549 % $$$ clear A5; % for V2
550 % $$$ A5 = sparse([],[],[],p,p,7*p);
551 % $$$ for m=1:p;
552 % $$$ A5(m,m)= -X(1,m)^h/(M2^h+X(1,m)^h);
553 % $$$ j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)+1));
554 % $$$ A5(j,m) = X(1,m)^h/(M2^h+X(1,m)^h);
555 % $$$ end;
556 % $$$
557 % $$$ % state space of the nominal system
558 % $$$ A = full(A);
559 % $$$ B = eye(p);
560 % $$$ C = ones(1,p);
561 % $$$ D = zeros(1,p);
562 % $$$ sys=ss(A,B,C,D);
563 % $$$
564 % $$$ %Model reduction, sysb is a balanced realization for sys.sysb=T*sys*T^-1
565 % $$$ [sysb,g,T,Ti]= balreal(sys);
566 % $$$ n = 3;
567 % $$$ sysr = modred(sysb,[n:p],'del');
568 % $$$ a1 = eye(n-1);
569 % $$$ a2 = zeros(p-n+1,n-1);
570 % $$$ V = T^-1*[a1;a2];
571 % $$$ W = T'*[a1;a2];
572 % $$$
573 % $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
574 % $$$ % compare the reduced system with different parameters with the original system
575 % $$$
576 % $$$ for i = 1:10
577 % $$$ % Ar = k1*A1r+u1*A2r+u2*A3r+V1*A4r+V2*A5r
578 % $$$ % Air = W'*Ai*V
579 % $$$ Ar = parameterlist(i,1)*W'*A1*V + parameterlist(i,2)*W'*A2*V + parameterlist(i,3)*W'*A3*V + parameterlist(i,4)*W'*A4*V + parameterlist(i,5)*W'*A5*V;
580 % $$$ Cr = C*V;
581 % $$$
582 % $$$ % original system
583 % $$$ a1 = @(x) (parameterlist(i,1)+parameterlist(i,4)*x(2)^h/(M1^h+x(2)^h));
584 % $$$ a2 = @(x) parameterlist(i,2)*x(1);
585 % $$$ a3 = @(x) parameterlist(i,5)*x(1)^h/(M2^h+x(1)^h);
586 % $$$ a4 = @(x) parameterlist(i,3)*x(2);
587 % $$$
588 % $$$ A = sparse([],[],[],p,p,7*p);
589 % $$$ for m=1:p
590 % $$$ A(m,m) = -(a1(X(:,m))+a2(X(:,m))+a3(X(:,m))+a4(X(:,m)));
591 % $$$ j = find((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
592 % $$$ A(j,m) = a1(X(:,m));
593 % $$$ j = find((X(1,:)==X(1,m)-1) & (X(2,:)==X(2,m)));
594 % $$$ A(j,m) = a2(X(:,m));
595 % $$$ j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)+1));
596 % $$$ A(j,m) = a3(X(:,m));
597 % $$$ j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)-1));
598 % $$$ A(j,m) = a4(X(:,m));
599 % $$$ end;
600 % $$$
601 % $$$
602 % $$$ subplot(4,3,i);
603 % $$$ x0=zeros(p,1);
604 % $$$ x0(1)=1;
605 % $$$ [t,x]=ode15s(@(t,x) A*x, [0 1000000], x0);
606 % $$$ plot(t,sum(x'),'r');
607 % $$$ grid on;
608 % $$$ hold on;
609 % $$$
610 % $$$ % reduced system
611 % $$$ x1=T*x0;
612 % $$$ x2=x1(1:(n-1),:);
613 % $$$ [tr,xr]=ode15s(@(tr,xr) Ar*xr, [0 1000000], x2);
614 % $$$
615 % $$$ clear y;
616 % $$$ for m = 1:size(xr,1)
617 % $$$ y(m)=sysr.c*xr(m,:)';
618 % $$$ end;
619 % $$$ plot(tr,y);
620 % $$$ hold on;
621 % $$$ xlabel('zeit t');
622 % $$$ ylabel('y');
623 % $$$
624 % $$$ end;
625 % $$$
626 % $$$
627 % $$$
628 % $$$
629 % $$$
630 %| \docupdate
function model = follicle_rect_model(params)
model of the human follicle growth
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...
Definition: verbose.m:17
function [ RBext , dummy ] = RB_extension_PCA_fixspace(model, detailed_data)
function computing a RB basis extension for given parameters by the POD-Greedy algorithm.
function varargout = filecache_function(funcptr, varargin)
function used for file-caching other function calls.
function model = follicle_model(params)
model of the human follicle growth
A cartesian rectangular grid in two dimensions with axis parallel elements.
Definition: rectgrid.m:17
function save_detailed_simulations(model, model_data, M, savepath)
perform loop over detailed simulations and save results or check consistency with existing saved resu...