oscillator_model.m

function model = oscillator_model(params);
%function model = oscillator_model;
%
% Steffen's oscillator model

% B.Haasdonk and S. Waldherr 5.8.2010
% Time-stamp: <Last change 2010-09-16 15:21:48 by Steffen Waldherr>

if (nargin < 1) 
  params = [];
end;

model = lin_ds_model_default;
model.data_const_in_time = 1;
model.verbose = 0;
model.debug = 0;
model.clim = [0,0.00005];

% useful range and number of parameters:
model.mu_names = {'k3'};
model.mu_ranges = {[0.001,1]};
%model.T = 20;  
%model.nt = 200; % number of timesteps
model.T = 10;  
model.nt = 100; % number of timesteps

% model parameters
model.epsilon = 0.01;
model.range = [300, 300];
%model.range = [51, 51];
model.kp = 1;
model.k1 = 15;
model.kq = 1;
model.k2 = 0.2;
model.k3 = 0.1;
model.k4 = 6.5;
model.s = 10;

model.affinely_decomposed = 1;
model.u_function_ptr = @(model) 0; % define u constant 
model.A_function_ptr = @A_func;
model.x0_function_ptr = @x0_func;
model.B_function_ptr = @B_func;
model.C_function_ptr = @C_func;
model.D_function_ptr = @D_func;

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;

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 time 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 = 0;
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
  warning('please adjust the computation of the error estimator constants!')
  model.state_bound_constant_C1 = 1 ;
  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];
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 = eval_affine_decomp(@A_components,@A_coefficients,...
                       model,model_data);

function Acomp = A_components(model,model_data)
Acomp = filecache_function(...
    @Acomp_assembly,model.s,model.k2, model.k4,model_data.X, model_data.numstates);

function Acomp = Acomp_assembly(s,k2,k4,X,numstates);
% define reaction rate functions without linear parameters
% $$$ reactions
% $$$ degradep, P -> , kp;
% $$$ producep, -> P, [k1*s^2/(Q+s*k2)];
% $$$ degradeq, Q -> , kq/epsilon;
% $$$ produceq, -> Q, [(Q^2/(Q^2+k4*s^2)*P + k3*s)/epsilon];

a{1} = @(x) x(1); % degrade x1
a{2} = @(x) s^2/(x(2)+s*k2); % produce x1 regulated
a{3} = @(x) x(2); % degrade x2
a{4} = @(x) x(1)*x(2)^2/(x(2)^2+k4*s^2); % produce x2 regulated
a{5} = @(x) s; % produce x2 constitutively
% define stoichiometric matrix
S = [-1 1 0 0 0; 0 0 -1 1 1];
% define the initial condition
% define range of discrete network configurations

% Acoeff contains the coefficient matrices (without parameters)
%Acoeff = cell(5);
for j=1:size(S,2)
  Acoeff{j} = sparse([],[],[],numstates,numstates,...
                     size(S,2)*numstates);
end;
for i=1:numstates
  for j=1:size(S,2)
    Acoeff{j}(i,i) = - a{j}(X(:,i));
    % find the index of the state that the system goes to with reaction j
    ind = find((X(1,:) == X(1,i)+S(1,j)) & ...
               (X(2,:) == X(2,i)+S(2,j)));
%    assert(isscalar(ind) || isempty(ind), 'Found multiple occurences of target state in reaction, network structure may be incorrect.');
    Acoeff{j}(ind,i) = a{j}(X(:,i));
  end;
end;
Acomp = {Acoeff{1},Acoeff{2},Acoeff{3},Acoeff{4},Acoeff{5}};



function Acoeff = A_coefficients(model)
Acoeff = [model.kp,model.k1,model.kq/model.epsilon,1/model.epsilon,model.k3/model.epsilon];

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

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

function Bcoeff = B_coefficients(model)
Bcoeff = 1;

function C = C_func(model,model_data)
C = eval_affine_decomp(@C_components,...
                       @C_coefficients,model,model_data);

function Ccomp = C_components(model,model_data)
%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_data.numstates)};
%Ccomp = {[2,0,1],[0,2,1]};

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

function D = D_func(model,model_data)
%D = [0;0]; % no decomposition required for scheme
D = [0]; % no decomposition required for scheme

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

function x0comp = x0_components(model,model_data)
% each row in ic is of the form [x1, x2, weight of (x1,x2) in ic]
ic = [0 0 1];
Pic = zeros(model_data.numstates,1);
i = find((model_data.X(1,:)<=50) & (model_data.X(2,:)<=50));
Pic(i) = 1/length(i);
x0comp = {Pic};
%weights = sum(ic(:,3));
%for i=1:size(ic,1)
%  ind = find((model_data.X(1,:) == ic(i,1)) & (X(2,:) == ic(i,2)));
%  Pic(ind) = ic(i,3)/weights;
%end;

function x0coeff = x0_coefficients(model)
x0coeff = 1;

function model_data = my_ds_gen_model_data(model)
% compute state vectors using a rectangular shape
numstates = (model.range(1)+1)*(model.range(2)+1);
X = zeros(2,numstates);
s = 1;
for x1=0:model.range(1)
  for x2=0:model.range(2)
    X(:,s) = [x1; x2];
    s = s+1;
  end;
end;
model_data.X = X;
model_data.numstates = numstates;

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;

warning('please use  G_matrix_function_ptr')
%model_data.G = model.G_matrix_function_ptr(model,[]);
model_data.G = speye(model_data.numstates);

% only required for plotting:
params.xnumintervals = model.range(1)+1;
params.ynumintervals = model.range(2)+1;
params.xrange = [-0.5,model.range+0.5];
params.yrange = [-0.5,model.range+0.5];
params.verbose = 0;
model_data.plot_grid = rectgrid(params);
%model_data.plot_nm_index_vector = ...
%    model_data.Ms + model_data.Ns*(model.range+1)+1;
plot_permind = reshape(1:(model.range(1)+1)*(model.range(2)+1),...
                  model.range(2)+1,model.range(1)+1)';
model_data.plot_permind = plot_permind(:);

function detailed_data = my_ds_gen_detailed_data(model,model_data)
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
% 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,sim_data.X(model_data.plot_permind,:),merged_plot_params)
axis tight;
xlabel('species p');
ylabel('species q');
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)
params.plot_function = 'plot_element_data';
if ~isfield(params,'clim') && isfield(model,'clim')
  params.clim = model.clim;
end;
p = plot_element_data(sim_data.X(model_data.plot_permind,params.timestep+1),model_data.plot_grid,params);
%axis tight;
xlabel('species p');
ylabel('species q');
title('state probabilities');

return;


%%%%% code from Steffen:

% generate and reduce CME model for stochastic resonance oscillator
% El-Samad and Khammash, CDC 2006
%
% File initiated 2010-07-29
% Time-stamp: <Last change 2010-07-29 17:44:34 by Steffen Waldherr>

% define linear reaction parameters
clear p
epsilon = 0.01;
p{1} = 1; % kp
p{2} = 15; % k1
p{3} = 1/epsilon; % kq
p{4} = 1/epsilon; % implicit parameter in original model
p{5} = 0.1/epsilon; % k3

% define non-linear (constant) reaction parameters
k2 = 0.2;
k4 = 6.5;

% define reaction rate functions without linear parameters
clear a;
a{1} = @(x) x(1);
a{2} = @(x) 1/(x(2)+k2);
a{3} = @(x) x(2);
a{4} = @(x) x(1)*x(2)^2/(x(2)^2+k4);
a{5} = @(x) 1;

% define stoichiometric matrix
S = [-1 1 0 0 0; 0 0 -1 1 1];

% define the initial condition
% each row in ic is of the form [x1, x2, weight of (x1,x2) in ic]
ic = [0 0 1];

% define range of discrete network configurations
range = [50, 50];
% range = [3, 3];

% compute state vectors using a rectangular shape
numstates = (range(1)+1)*(range(2)+1);
X = zeros(2,numstates);
s = 1;
for x1=0:range(1)
  for x2=0:range(2)
    X(:,s) = [x1; x2];
    s = s+1;
  end;
end;

% TODO: define relevant outputs: average, probability of being near deterministic equilibrium?

% Acoeff contains the coefficient matrices (without parameters)
for j=1:size(S,2)
  Acoeff{j} = sparse([],[],[],numstates,numstates,size(S,2)*numstates);
end;
for i=1:numstates
  for j=1:size(S,2)
    Acoeff{j}(i,i) = - a{j}(X(:,i));
    % find the index of the state that the system goes to with reaction j
    ind = find((X(1,:) == X(1,i)+S(1,j)) & (X(2,:) == X(2,i)+S(2,j)));
    assert(isscalar(ind) || isempty(ind), 'Found multiple occurences of target state in reaction, network structure may be incorrect.');
    Acoeff{j}(ind,i) = a{j}(X(:,i));
  end;
end;

% compute transition matrix A
A = sparse([],[],[],numstates,numstates,size(S,2)*numstates);
for j=1:size(S,2)
  A = A + p{j}*Acoeff{j};
end;

Pic = zeros(numstates,1);
weights = sum(ic(:,3));
for i=1:size(ic,1)
  ind = find((X(1,:) == ic(i,1)) & (X(2,:) == ic(i,2)));
  Pic(ind) = ic(i,3)/weights;
end;

[t,P] = ode15s(@(t,x) A*x,[0 1], Pic);