function [err,sim_output,output]=simulation_error(m,varargin)
%SIMULATION_ERROR computes the mean square error of the
%    simulated model
%    
%    This function performs a simulation of the model using
%    the SIMULATE method and then calculates the mean square
%    error by comparing the simulated output with the value of
%    the output contained in the output data set. 
%    
%    M=SIMULATION_ERROR(M,J,'PLOT') get the simulation error
%    of the model M respect	the output J. If the instead of
%    the number of the output the label 'all' is given to the
%    function the error respect all the outputs is returned in
%    an array. The optional label PLOT may be specified if a
%    plot of the original data and the simulation is desired. 
%     
%    
%    M=SIMULATION_ERROR(M,DATA,J) get the error of the
%    simulated model M respect   the output J and the data
%    DATA input-output set. DATA must be a structure defining
%    the 'dat_in' and 'dat_out'   fields which contain the
%    input output data. If the instead of the number of the
%    output J the label 'all' is given to the function the
%    error respect all the outputs is returned in an array.
%    The optional label <code>'plot'</code> may be specified
%    if a plot of the original data and the simulation is
%    desired.   
%    
%    See also: FLAT_ERROR
%    
%    

%MAN_PAGE_BEGIN
%
%   @purpose  computes the mean square error of the simulated model. This function
%   performs a simulation of the model using the <code>simulate</code> method and
%   then calculates the mean square error by comparing the simulated output
%   with the value of the output contained in the output data set. 
%
%   @synopsis	m=simulation_error(m,j,'plot')
%   @description    get the simulation error of the model <code>m</code> respect 
%   the output <code>j</code>. If the instead of the number of the output
%   the label 'all' is given to the function the error respect all the outputs is
%   returned in an array. The optional label <code>plot</code> may be
%   specified if a plot of the original data and the simulation is
%   desired.   
%
%   @synopsis	m=simulation_error(m,data,j)
%   @description    get the error of the simulated model <code>m</code> respect  
%   the output <code>j</code> and the data <code>data</code> input-output set.
%   <code>data</code> must be a structure defining the 'dat_in' and 'dat_out'  
%   fields which contain the input output data. If the instead of the
%   number of the output <code>j</code>
%   the label 'all' is given to the function the error respect all the outputs is
%   returned in an array. The optional label <code>'plot'</code> may be
%   specified if a plot of the original data and the simulation is
%   desired.   
%
%   @see flat_error
%
%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
  
if nargin > 4,
  error('bad number of inputs');
elseif ~ isa(m,'external'),
  error('the first argument is not a external object');
end
  
if isa(varargin{1},'dataset'),
  old_data=get(m,'data');
  m=add_data(m,varargin{1},get(m,'T'));
  j=varargin{2};
else 
  j=varargin{1};
end

error(checkval(j,'intpos','j',{1 1:get(m,'n_out') 0},[1 get(m, ...
						  'n_out')]));
do_plot=0;
if isa(varargin{end},'char')
  do_plot=1;
end  

%get the output of the system
dd = get(get(m,'data'),'contdata');
out=dd(:,get(m,'n_in')+1:end);
 

%check if the simulation data has already been cached 
opt=get(m,'opt');
res=[];
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};
    end
  end
end
%if we did not find anything 
if isempty(res), 
  res=simulate(m);
end

err=sum((out(:,j) - res(:,j)).^2)./size(out,1);
sim_output=res;
output=out;

if (do_plot),
  hold on;
  max_y=max(max([out(:,j) res(:,j)]));
  max_y=max_y+max_y*0.05;
  min_y=min(min([out(:,j) res(:,j)]));
  min_y=min_y-max_y*0.05;
  axis([0 size(out,1) min_y max_y]);
  for i=1:size(j,2);
    plot([out(:,j(i)),res(:,j(i))]);
  end
  hold off;
end


%update the caller
if isa(varargin{1},'dataset'),
  m=add_data(m,old_data,get(m,'T'));
end

assignin('caller',inputname(1),m);


