function [res, m]=simulate(m,time)
%SIMULATE simulates the model
%    
%    This function performs a simulation of the  model M using
%    the input data and feeding back the results of the
%    simulation to calculate at each step the regressors. The
%    result of  the simulation RES is a SIMULATION_TIME *
%    NUM_OUTPUTS array. 
%    
%    [RES,M]=SIMULATE(M) simulate the model M using the input
%    data.
%    
%    [RES,M]=SIMULATE(M,TIME) simulate the model  M for the
%    time horizon TIME
%    
%    

%MAN_PAGE_BEGIN
%
%   @purpose    simulates the model. This function performs a simulation of the 
%   model <code>m</code> using the input data and feeding back the results
%   of the simulation to calculate at each step the regressors. The result of 
%   the simulation <code>res</code> is a <code>simulation_time * num_outputs</code>
%   array. 
%
%   @synopsis   [res,m]=simulate(m)
%   @description  simulate the model <code>m</code> using the input data.
%
%   @synopsis    [res,m]=simulate(m,time)
%   @description  simulate the model  <code>m</code> for the time horizon <code>time</code>
%
%MAN_PAGE_END 

check(m); 

if get(m,'n_out') <=0,
  error ('the number of the outputs has not yet been initialized');
end
if get(m,'n_in') <=0,
  error ('the number of the outputs has not yet been initialized');
end

error(checkargs({m},1,'external')); 

%check if the simulation data has already been cached 
opt=get(m,'opt');
if isfield(opt,'simulation_database')  & ~isempty(opt.simulation_database.links),
  database=opt.simulation_database;
  links=get(m,'links');
  for i=1:size(database.links,2),
    %if all the links filename are the same
    if all(links==database.links{i}) & ...
	  strcmp(get(get(m,'data'),'filename'),database.file{i}),
      res=database.simulated_outputs{i};
      return
    end
  end
end

%since no cached simulation is available perform it    
%get the input and output data sets
dd = get(get(m,'data'),'contdata'); 
if isempty(dd)
  error('no data attached to the object')
end

in=dd(:,1:get(m,'n_in'));
out=dd(:,get(m,'n_in')+1:end);

%if no time horizon is specified simulate for all the inputs 
if nargin<2, time=size(in,1); end

%gets the dynamic information from the system
%how many inputs are used in the dynamics
dynIn =size(m.dynamics.nu,2);

%the time horizon of the model
maxt = max(max(max(m.dynamics.ny+(m.dynamics.ny>0))),...
	   max(max(m.dynamics.nu+m.dynamics.nd)));

%initialises the system. the outputs precedent to the 
%start of the simulation are supposed to be equal to
%the first output, while the inputs are all 0
x=[out(ones(maxt,1),:) zeros(maxt,dynIn)];

%checks if the mex-file exists
if exist('simulate') == 3,
	res = simulate(struct(m),in',x(1:end-1,[size(out,2)+(1:dynIn) 1:size(out,2)])')';
else
	limits = get(m,'limits');
	inlim = limits(:,1:get(m,'n_in'));
	outlim = limits(:,get(m,'n_in')+(1:get(m,'n_out')));

	%saturate the inputs
	in = min(max(in,inlim(ones(size(in,1),1),:)),inlim(1+ones(size(in,1),1),:));

	%creates the 'work in progress bar' window
	h = my_waitbar(-1,0,'Simulating...');

	%inizializes the results of the simulation
	res=zeros(time,get(m,'n_out'));

	%simulates the system       
	for i=1:time,
	  %put the newest input data in the last row of the matrix,
	  %leaving the output part to zero 
	  xx=[x(2:end,:);[zeros(1,get(m,'n_out')) in(i,:)]];
	  %build the needed regressor vectors
	  [r,yy] = build_regressors(xx,m.dynamics.ny, m.dynamics.nu,m.dynamics.nd,1);
	  %evaluates the output of the model
	  for output=1:get(m,'n_out'),
	    y(output)=min(max(eval(m,r{output},output),outlim(1,output)),outlim(2,output));
	  end
	  %updates the matrix with the past data
	  x=[x(2:end,:);[y in(i,:)]];
	  %saves the result
	  res(i,:)=y;
	  %update the 'work in progress bar' window 
	  h= my_waitbar(h,i/time);
	end
	%closes the 'work in progress bar' window
	close(h);

end 

%update the simulation database
opt=get(m,'opt');
if isfield(opt,'simulation_database')  & ~isempty(opt.simulation_database.links),
  database=opt.simulation_database;
else,
  database=struct('links',cell(1,1),'file',cell(1,1),...
		  'simulation_errors',cell(1,1),...
		  'outputs',cell(1,1),...
		  'simulated_outputs',cell(1,1));
end 
data=get(m,'data');

%the new points have the following structure:
% {links, data_file, simulation_error outputs simulated_outputs}
database.links{end+1}=get(m,'links');
database.file{end+1}=get(get(m,'data'),'filename');
database.simulation_errors{end+1}=sum((out- res).^2)./size(out,1);
database.outputs{end+1}=out;
database.simulated_outputs{end+1}=res;

%write the results in the cache
opt.simulation_database=database; 
set(m,'opt',opt);

%update the database on the calling object
assignin('caller',inputname(1),m);

