follicle_model.m

Go to the documentation of this file.

function model = follicle_model(params)
%function model = follicle_model(params)
% model of the human follicle growth
%
% two chemical components are assumed to be characteristic for the
% growth of human follicles. 
%
% A deterministic model for the two species is given by a positive
% feedback loop, a bistable switch, i.e. 2-d nonlinear ODE.
% 
% An approximation to the Chemical Master Equation gives a linear
% dynamical system.
%
% here, the state vector is the probability of states (pair of number of two
% chemical components) and the matrix A describes the change in
% these probabilities over time. The number of states is
% p = (range+1)*(range+2)/2
%
% The output is given by params.output_type:
% 'separatrix': an average over a subset of states, the ones below
%   the "separatrix", which are assumed to be still "in rest",
%   i.e. waiting for their growth specified by model.output_range
% 'expectation_M': the expectation of M state is computed
% 'expectation_N': the expectation of N state is computed

% Q: D.h. die Zellen oberhalb der separatrix sind einfach abgestorben, oder zur
%   Ovulation gekommen?
%
% Explanations from S. Waldherr "A stochastic switch for the
% primordial to primary transition in follicle growth", IST
% internal report (2009-04-09) v1

% Bernard Haasdonk 22.9.2009
if (nargin < 1) 
  params = [];
end;

model = lin_ds_model_default;
model.data_const_in_time = 1;
model.verbose = 0;
model.debug = 0;

%%%%%%%%%%%%%%%%%%%%%%% parameters from steffen
model.M1 = 25;
model.h = 3;
model.M2 = 25;
%inomi=6;
if isfield(params,'range')
  model.range = params.range;
else
  model.range = 50; % connected with end time, increase e.g. to 50 later
end;
if isfield(params,'output_range')
  model.output_range = params.output_range;
else
  model.output_range = model.range; % connected with end time, increase e.g. to 50 later
end;
if (model.output_range > model.range)
  error('please choose output range <= range');
end;
% determines the initial distribution in a lower left triangle
if isfield(params,'init_range')
  model.init_range = params.init_range;
else
  model.init_range = 0; 
end;
% 
if isfield(params,'reg_epsilon')
  model.reg_epsilon = params.reg_epsilon;
else
  model.reg_epsilon = 0; 
end;
if ~isfield(params,'output_type')
  params.output_type = 'separatrix'; 
end;
%for i = 1:10 
%     u1 = 0.01+(i-6)*0.001;
%     u2 = 0.01+(i-6)*0.0005;
%     k1 = 4*u1+(i-6)*0.0001;
%     V1 = 75*u1;
%     V2 = 75*0.01;
%     parameterlist(i,1:5)= [k1 u1 u2 V1 V2];
%end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%model_data = my_ds_gen_model_data(model);

model.u1 = 0.01;
model.u2 = 0.01;
model.k1fact = 4;
model.V1fact = 75;
model.V2fact = 75;

% used as components later:
%model.k1 = model.k1fact * model.u1;
%model.V1 = model.V1fact * model.u1;
%model.V2 = model.V2fact * model.u2;

model.mu_names = {'u1','u2'};
%model.mu_names = {'u1','u2','k1fact','V1fact','V2fact'};
model.mu_ranges = {[0.005,0.02],[0.005,0.02]};

%model.T = 1000000;  % Minutes, i.e. 1.9 Jahre 
%model.nt = 500;
%model.T = 1000000;  % Minutes, i.e. 1.9 Jahre 
%model.nt = 10000;
model.T = 10000000;  % Minutes, i.e. 19 Jahre 
model.nt = 500;
%model.T = 5000000;  % Minutes, i.e. 9.5 Jahre 
%model.nt = 500;
%model.T = 50000000;  % Minutes, i.e. 95 Jahre 
%model.nt = 10000;
%disp('nt to be adjusted later!!!') => 500 OK!!
model.affinely_decomposed = 1;
model.A_function_ptr = @A_func;
% ratio of first to second ellipses axes
%model.A_axes_ratio = 2;
%model.A_axes_ratio = 2;
model.B_function_ptr = @B_func;
switch params.output_type
 case 'separatrix'
  model.C_function_ptr = @C_func_separatrix;
 case 'expectation_M'
  model.C_function_ptr = @C_func_expectation_M;
 case 'expectation_N'
  model.C_function_ptr = @C_func_expectation_N;
 otherwise
  error('params.output_type not known');
end;
model.D_function_ptr = @D_func;
model.x0_function_ptr = @x0_func;
% shift of x0 in third coordinate
%model.x0_shift=1;
model.u_function_ptr = @(model) 0; % define u constant 
%model.u_function_ptr = @(model) sin(4*model.t);  
%model.u_function_ptr = @(model) 0; % define u constant 0
%model.u_function_ptr = @(model) model.t; % define u = 1 * t
model.G_matrix_function_ptr = ...
    @(model,model_data) speye(model.dim_x,model.dim_x);
% the following lead to both larger C1 and C2!!!
%model.G_matrix_function_ptr = @(model) diag([2,0.5,1]);
%model.G_matrix_function_ptr = @(model) diag([0.5,2,1]);
 
% set function pointers to main model methods
model.gen_model_data = @my_ds_gen_model_data;
%model.detailed_simulation = @lin_ds_detailed_simulation;
model.gen_detailed_data = @my_ds_gen_detailed_data;
%model.gen_reduced_data = @lin_ds_gen_reduced_data;
%model.rb_simulation = @lin_ds_rb_simulation;
%model.plot_sim_data = @my_ds_plot_sim_data;
model.plot_sim_data_state = @plot_state_probabilities;
model.plot_slice = @my_plot_slice;
model.axis_tight = 1;

%model.rb_reconstruction = @lin_ds_rb_reconstruction;
% determin model data, i.e. matrix components, matrix G

model.orthonormalize = @model_orthonormalize_qr_no_perm;

model_data = gen_model_data(model); % matrix components
model.dim_x = size(model_data.X,2);
model.dim_y = size(model_data.C_components{1},1);

if isfield(params,'output_plot_indices')
  model.output_plot_indices = params.output_plot_indices;
else
  model.output_plot_indices = 1:model.dim_y;
end;

model.theta = 1; % implicit discretization

% be sure to determine correct constants once, if anything in the
% problem setting has changed!!

%model.error_estimation = 1; % model is capable of error_estimation
model.enable_error_estimator = 1;
if ~isfield(params,'estimate_bounds')
  params.estimate_bounds = 0;
end;
% the following only needs to be performed, if the model is changed.
if (params.estimate_bounds == 0)
  % for G = Id
  model.state_bound_constant_C1 = 1 ;
  if ~isequal(params.output_type,'separatrix')
    warning('please adjust the computation of the error estimator constants!')
  end;
  model.output_bound_constant_C2 = sqrt(model.dim_x);
elseif params.estimate_bounds == 1
  model.estimate_lin_ds_nmu = 3;
  model.estimate_lin_ds_nX = 10;
  model.estimate_lin_ds_nt = 10;
  format long;
  [C1,C2] = lin_ds_estimate_bound_constants(model,model_data)
  model.state_bound_constant_C1 = C1;
  model.output_bound_constant_C2 = C2;
else
  model.estimate_lin_ds_nmu = 3;
  model.estimate_lin_ds_nX = 10;
  model.estimate_lin_ds_nt = 10;
  format long;
  [C1,C2] = lin_ds_estimate_bound_constants(model,model_data)
  model.state_bound_constant_C1 = C1;
  model.output_bound_constant_C2 = C2;  
end;


model.inner_product_matrix_algorithm =@(model,model_data) ...
     speye(model.dim_x,model.dim_x);

model.error_algorithm = @(u1,u2,dummy1,dummy2) ...
    sqrt(max(sum((u1-u2).^2),0));
%    [4,3]* [3,2];

%model.error_norm = 'l2';
%%                  'RB_extension_max_error_snapshot'}; 
model.RB_stop_timeout = 4*60*60; % 4 hours              
model.RB_stop_epsilon = 1e-10; 
model.RB_error_indicator = 'error'; % true error
                                    
%model.RB_error_indicator = 'estimator'; % Delta from rb_simulation
                                    
model.RB_stop_Nmax = 200;
model.RB_generation_mode = 'greedy_uniform_fixed';
model.RB_numintervals = [4,4];
model.RB_detailed_train_savepath = 'follicle_model_detailed_train';
model.orthonormalize  = @model_orthonormalize_qr;
model.RB_extension_algorithm = @RB_extension_PCA_fixspace;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% functions defining the dynamical system and input/initial data
function A = A_func(model,model_data);
  % A_func
A = eval_affine_decomp(@A_components,@A_coefficients,...
                       model,model_data);

function Acomp = A_components(model,model_data)
  % A_components
p = model.dim_x;
X = model_data.X;
h = model.h;
M1 = model.M1;
M2 = model.M2;
A1 = sparse([],[],[],p,p,7*p); % for k1
A2 = sparse([],[],[],p,p,7*p); % for u1              
A3 = sparse([],[],[],p,p,7*p); % for u2
A4 = sparse([],[],[],p,p,7*p); % for V1
A5 = sparse([],[],[],p,p,7*p); % for V2
A6 = -speye(p,p); % for regularization
for m=1:p
  A1(m,m)= -1;
  j = find ((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
  A1(j,m) = 1;
  A2(m,m) = -X(1,m);
  j = find((X(1,:)==X(1,m)-1) & (X(2,:)==X(2,m)));
  A2(j,m) = X(1,m);
  A3(m,m) = -X(2,m);
  j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)-1));
  A3(j,m) = X(2,m);
  A4(m,m) = -X(2,m)^h/(M1^h+X(2,m)^h);
  j = find((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
  A4(j,m) = X(2,m)^h/(M1^h+X(2,m)^h);
  A5(m,m)= -X(1,m)^h/(M2^h+X(1,m)^h);
  j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)+1));
  A5(j,m) = X(1,m)^h/(M2^h+X(1,m)^h);
end;
%keyboard;
Acomp = {A1, A2, A3, A4, A5, A6};

function Acoeff = A_coefficients(model)
  % A_coefficients

% Acoeff = [model.k1, model.u1, model.u2, model.V1, model.V2];
k1 = model.k1fact * model.u1;
V1 = model.V1fact * model.u1;
V2 = model.V2fact * model.u2;
Acoeff = [k1, model.u1, model.u2, V1, V2, model.reg_epsilon];

function B = B_func(model,model_data)
  % B_func
B = eval_affine_decomp(@B_components,@B_coefficients,model,model_data);

function Bcomp = B_components(model,model_data)
  % B_components
Bcomp = {zeros(model.dim_x,1)};

function Bcoeff = B_coefficients(model)
  % B_coefficients
Bcoeff = 1;

function C = C_func_separatrix(model,model_data)
  % C_func_separatrix
C = eval_affine_decomp(@C_components_separatrix,...
                       @C_coefficients,model,model_data);

function C = C_func_expectation_M(model,model_data)
  % C_func_expectation_M
C = eval_affine_decomp(@C_components_expectation_M,...
                       @C_coefficients,model,model_data);

function C = C_func_expectation_N(model,model_data)
  % C_func_expectation_N
C = eval_affine_decomp(@C_components_expectation_N,...
                       @C_coefficients,model,model_data);

function Ccomp = C_components_separatrix(model,model_data)
  % C_components_separatrix
i = find((model_data.Ms+model_data.Ns)<=model.output_range);
C = zeros(1,model.dim_x);
C(i) = 1;
Ccomp = {C};
%  Ccomp = {ones(1,model.dim_x)};
%Ccomp = {[2,0,1],[0,2,1]};

function Ccomp = C_components_expectation_M(model,model_data)
  % C_components_expectation_M
C = model_data.Ms(:)';
Ccomp = {C};

function Ccomp = C_components_expectation_N(model,model_data)
  % C_components_expectation_N
C = model_data.Ns(:)';
Ccomp = {C};

function Ccoeff = C_coefficients(model)
  % C_coefficients
Ccoeff = [1];
%Ccoeff = [model.C_weight,1-model.C_weight];

function D = D_func(model,model_data)
  % D_func

%D = [0;0]; % no decomposition required for scheme
D = [0]; % no decomposition required for scheme

function x0 = x0_func(model,model_data)
  % x0_func
x0 = eval_affine_decomp(@x0_components,@x0_coefficients,model,model_data);

function x0comp = x0_components(model,model_data)
  % x0_components

%if model.init_range == 0
%  x0comp = {eye(model.dim_x,1)};
%else
i = ((model_data.Ms+model_data.Ns)<=model.init_range);
%i = (model_data.Ms<=1);
x0comp = {i(:)/sum(i)}; % sum should be 1!
%end;

function x0coeff = x0_coefficients(model)
  % x0_coefficients
x0coeff = 1;

function model_data = my_ds_gen_model_data(model)
  % model data generation
s = 1; % state counter
% generate list of states X and output weight vector
Ms = zeros(1,(model.range+1)*model.range/2);
Ns = zeros(1,(model.range+1)*model.range/2);
for m=0:model.range
  for n=0:(model.range-m)
    X(:,s) = [m;n];
    Ns(s) = n;
    Ms(s) = m;
    s = s+1;
  end;
end;
p = s-1; % number of states
model.dim_x = p;
%model_data = lin_ds_gen_model_data(model);
model_data.X = X;
model_data.Ns = Ns;
model_data.Ms = Ms;

if model.affinely_decomposed
  model.decomp_mode = 1;
  % the following call is extremely expensive: 30 seconds!!!!
  % therefore caching is activated
  %  model_data.A_components = model.A_function_ptr(model,model_data);
  model_data.A_components = filecache_function(...
      model.A_function_ptr,model,model_data);
  model_data.B_components = model.B_function_ptr(model,model_data);
  model_data.C_components = model.C_function_ptr(model,model_data);
  % no D components
  %  model_data.D_components = model.D_function_ptr(model);
  model_data.x0_components = model.x0_function_ptr(model,model_data);
end;
model_data.G = model.G_matrix_function_ptr(model,[]);

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;

% only for plotting:
model_data.plot_grid = rectgrid(params);
model_data.plot_nm_index_vector = ...
    model_data.Ms + model_data.Ns*(model.range+1)+1;

function detailed_data = my_ds_gen_detailed_data(model,model_data)
  % detailed data generation
detailed_data = lin_ds_gen_detailed_data(model,model_data);
% only for plotting:
detailed_data.plot_grid = model_data.plot_grid;
detailed_data.plot_nm_index_vector = model_data.plot_nm_index_vector;

function p = plot_state_probabilities(model,model_data,sim_data,params)
% plot of trajectory and output

%p = lin_ds_plot_sim_data(model,model_data,sim_data,params);
%close(gcf-1);
p = []; % no handles
%grid = Grids.rectgrid(params);
U = nan((model.range+1)*(model.range+1),size(sim_data.X,2));
U(model_data.plot_nm_index_vector,1:end)= sim_data.X;
% for debug:
%U(nm_index_vector,1) = model_data.Ns;
%U(nm_index_vector,2) = model_data.Ms;
%params.plot_function = @plot_element_data;
merged_plot_params = merge_model_plot_params(model, params);
merged_plot_params.plot_function = 'plot_element_data';
plot_data_sequence(model,model_data.plot_grid,U,merged_plot_params)
axis tight;
xlabel('species m');
ylabel('species n');
if ~isfield(model,'plot_title')
  title('state probabilities');
else
  title(model.plot_title);
end;

function p = my_plot_slice(model,model_data,sim_data,params)
  % plots only a slice
U = nan((model.range+1)*(model.range+1),1);
U(model_data.plot_nm_index_vector)= sim_data.X(:,params.timestep+1);
% for debug:
%U(nm_index_vector,1) = model_data.Ns;
%U(nm_index_vector,2) = model_data.Ms;
%params.plot_function = @plot_element_data;
params.plot_function = 'plot_element_data';
if ~isfield(params,'clim') && isfield(model,'clim')
  params.clim = model.clim;
end;
p = plot_element_data(U,model_data.plot_grid,params);
%axis tight;
xlabel('species m');
ylabel('species n');
title('state probabilities');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CODE FROM STEFFEN, RESP HANWEI LI:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% $$$ 
% $$$ 
% $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ % define parameter 
% $$$ clear;
% $$$ M1 = 25;
% $$$ h = 3;
% $$$ M2 = 25;
% $$$ 
% $$$ for i = 1:10 
% $$$      u1 = 0.01+(i-6)*0.001;
% $$$      u2 = 0.01+(i-6)*0.0005;
% $$$      k1 = 4*u1+(i-6)*0.0001;
% $$$      V1 = 75*u1;
% $$$      V2 = 75*0.01;
% $$$      parameterlist(i,1:5)= [k1 u1 u2 V1 V2];
% $$$ end;
% $$$ inomi=6;
% $$$ 
% $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ % build the nominal model
% $$$ 
% $$$ range = 20;
% $$$ s = 1; % state counter
% $$$ 
% $$$ % generate list of states X
% $$$ clear X;
% $$$ for m=0:range
% $$$   for n=0:(range-m)
% $$$           X(:,s) = [m;n];
% $$$           s = s+1;
% $$$   end;
% $$$ end;
% $$$ p = s-1; % number of states
% $$$ 
% $$$ % defining propensity functions 
% $$$ a1 = @(x) (parameterlist(inomi,1)+parameterlist(inomi,4)*x(2)^h/(M1^h+x(2)^h));
% $$$ a2 = @(x) parameterlist(inomi,2)*x(1);
% $$$ a3 = @(x) parameterlist(inomi,5)*x(1)^h/(M2^h+x(1)^h);
% $$$ a4 = @(x) parameterlist(inomi,3)*x(2);
% $$$      
% $$$ A = sparse([],[],[],p,p,7*p);
% $$$   for m=1:p
% $$$       A(m,m) = -(a1(X(:,m))+a2(X(:,m))+a3(X(:,m))+a4(X(:,m)));
% $$$       j = find((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
% $$$     A(j,m) = a1(X(:,m));
% $$$       j = find((X(1,:)==X(1,m)-1) & (X(2,:)==X(2,m)));
% $$$     A(j,m) = a2(X(:,m));
% $$$     j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)+1));
% $$$     A(j,m) = a3(X(:,m));
% $$$     j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)-1));
% $$$     A(j,m) = a4(X(:,m));
% $$$   end;
% $$$     
% $$$ 
% $$$ %  KoeffizientenZerlegung
% $$$ clear A1;      % for k1
% $$$ A1 = sparse([],[],[],p,p,7*p);
% $$$   for m=1:p
% $$$       A1(m,m)= -1;
% $$$       j = find ((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
% $$$       A1(j,m) = 1;
% $$$   end;
% $$$ 
% $$$ clear A2;       % for u1              
% $$$ A2 = sparse([],[],[],p,p,7*p);
% $$$   for m=1:p;
% $$$       A2(m,m) = -X(1,m);
% $$$       j = find((X(1,:)==X(1,m)-1) & (X(2,:)==X(2,m)));
% $$$       A2(j,m) = X(1,m);
% $$$   end;
% $$$ 
% $$$ clear A3;       % for u2
% $$$ A3 = sparse([],[],[],p,p,7*p);
% $$$   for m=1:p;
% $$$       A3(m,m) = -X(2,m);
% $$$       j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)-1));
% $$$       A3(j,m) = X(2,m);
% $$$   end;
% $$$ 
% $$$ clear A4;       % for V1
% $$$ A4 = sparse([],[],[],p,p,7*p);
% $$$   for m=1:p;
% $$$       A4(m,m) = -X(2,m)^h/(M1^h+X(2,m)^h);
% $$$       j = find((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
% $$$       A4(j,m) = X(2,m)^h/(M1^h+X(2,m)^h);
% $$$   end;
% $$$ 
% $$$ clear A5;       % for V2
% $$$ A5 = sparse([],[],[],p,p,7*p);
% $$$   for m=1:p;
% $$$       A5(m,m)= -X(1,m)^h/(M2^h+X(1,m)^h);
% $$$       j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)+1));
% $$$       A5(j,m) = X(1,m)^h/(M2^h+X(1,m)^h);
% $$$   end;
% $$$    
% $$$ % state space of the nominal system
% $$$ A = full(A);
% $$$ B = eye(p);
% $$$ C = ones(1,p);
% $$$ D = zeros(1,p);
% $$$ sys=ss(A,B,C,D);
% $$$ 
% $$$ %Model reduction, sysb is a balanced realization for sys.sysb=T*sys*T^-1
% $$$ [sysb,g,T,Ti]= balreal(sys);
% $$$ n = 3;
% $$$ sysr = modred(sysb,[n:p],'del');
% $$$ a1 = eye(n-1);
% $$$ a2 = zeros(p-n+1,n-1);
% $$$ V = T^-1*[a1;a2];
% $$$ W = T'*[a1;a2];
% $$$ 
% $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ %  compare the reduced system with different parameters with the original system
% $$$ 
% $$$ for i = 1:10
% $$$     % Ar = k1*A1r+u1*A2r+u2*A3r+V1*A4r+V2*A5r
% $$$     % Air = W'*Ai*V
% $$$     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;
% $$$     Cr = C*V;
% $$$    
% $$$     % original system
% $$$     a1 = @(x) (parameterlist(i,1)+parameterlist(i,4)*x(2)^h/(M1^h+x(2)^h));
% $$$     a2 = @(x) parameterlist(i,2)*x(1);
% $$$     a3 = @(x) parameterlist(i,5)*x(1)^h/(M2^h+x(1)^h);
% $$$     a4 = @(x) parameterlist(i,3)*x(2);
% $$$      
% $$$     A = sparse([],[],[],p,p,7*p);
% $$$     for m=1:p
% $$$        A(m,m) = -(a1(X(:,m))+a2(X(:,m))+a3(X(:,m))+a4(X(:,m)));
% $$$        j = find((X(1,:)==X(1,m)+1) & (X(2,:)==X(2,m)));
% $$$        A(j,m) = a1(X(:,m));
% $$$            j = find((X(1,:)==X(1,m)-1) & (X(2,:)==X(2,m)));
% $$$        A(j,m) = a2(X(:,m));
% $$$        j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)+1));
% $$$        A(j,m) = a3(X(:,m));
% $$$        j = find((X(1,:)==X(1,m)) & (X(2,:)==X(2,m)-1));
% $$$        A(j,m) = a4(X(:,m));
% $$$     end;
% $$$    
% $$$     
% $$$     subplot(4,3,i);
% $$$     x0=zeros(p,1);
% $$$     x0(1)=1;
% $$$     [t,x]=ode15s(@(t,x) A*x, [0 1000000], x0);
% $$$     plot(t,sum(x'),'r');
% $$$     grid on;
% $$$     hold on;
% $$$    
% $$$     % reduced system
% $$$     x1=T*x0;
% $$$     x2=x1(1:(n-1),:);
% $$$     [tr,xr]=ode15s(@(tr,xr) Ar*xr, [0 1000000], x2);
% $$$     
% $$$     clear y;
% $$$     for m = 1:size(xr,1)
% $$$         y(m)=sysr.c*xr(m,:)';
% $$$     end;
% $$$     plot(tr,y);
% $$$     hold on;
% $$$     xlabel('zeit t');
% $$$     ylabel('y');
% $$$ 
% $$$ end;
% $$$ 
% $$$ 
% $$$ 
% $$$ 
% $$$ 
%| \docupdate