function [Output] = ModelCalculation(data,NbInput,NbRuleMax,ErrorMax,FileName,OutputRelation,m,Thresold,ErrorCriteria,Debug)
%*******************************************************************************%
% File : ModelCalculation.m
% Last update: 15/10/98
% Author : SIEMENS\Dominique PASSAQUAY
% Subject: Build fuzzy models
%*******************************************************************************%
% ModelCalculation(data,NbInput,NbRuleMax,ErrorMax,FileName,OutputRelation,m,Thresold,ErrorCriteria,Debug)
%   data is an array where:
%    each rows represents a data
%    each columns represents an input or an output 
%   NbInput is the number of input of the fuzzy model
%   NbRuleMax is the maximal number of rules
%   ErrorMax is the maximal modelisation error accepted for all the outputs
%   FileName is an optional parameter to specify the name of the file where the fuzzy models are saved (default='model')
%   OutputRelation is an optional matrix to specify the relations between inputs and outputs (matrix[NbOutput,NbInput])
%    OutputRelation(i,j)=0 - The ith output and the jth input are independents 
%    OutputRelation(i,j)=1 - The ith output and the jth input are dependents 
%   m is an optional parameter to specify the fuzziness exponent (default=1.5)
%   Thresold is an optional parameter (default=1.0)
%    0 <= Thresold <= 1
%   ErrorCriteria is an optional parameter to define the maximal error criteria used (default=1)
%    1 - The maximal error is equal to the biggest error on the outputs
%    2 - The maximal error is the biggest euclidian distance between the output datas and the modelised datas
%   Debug is an optional boolean to plot the output errors after each modelisation step (default=0)
%
%   The fuzzy model structure is:
%    NbInput                  ... number of inputs
%    NbOutput                 ... number of outputs
%    NbRule                   ... number of fuzzy rules
%    m                        ... fuzziness exponent
%    ScaleFactor              ... scale factors of inputs and outputs (vector[NbInput+NbOutput])
%    Range                    ... Inputs and outputs ranges (vector[NbInput+NbOutput])
%    MembershipFunctionCenter ... memberships function centers (matrix[NbRule,NbInput+NbOutput])
%    LocalLinearModel         ... local linear model (matrix[NbRule,NbInput+1,NbOutput])
%    NbData                   ... number of datas used to build the fuzzy model
%    Thresold                 ... Thresold used to build the fuzzy model
%    ErrorCriteria            ... Error criteria used to build the fuzzy model

 if nargin<4, error('Incorrect number of inputs.'); end
 if not(NbRuleMax>0 & (NbInput>0 & NbInput<size(data,2)) & size(data,1)>NbInput+1), error('Incorrect parameters.'); end

 global Models;
 Models=[];
 
 NbData=size(data,1);
 NbOutput=size(data,2)-NbInput;

 if nargin<10, Debug=0; end
 if nargin<9, ErrorCriteria=1; end
 if nargin<8, Thresold=1; end
 if nargin<7, m=1.5; end
 if (nargin<6 | isempty(OutputRelation)), OutputRelation=ones(NbOutput,NbInput); end
 if nargin<5, FileName='model'; end
 
 if not((Debug==0 | Debug==1) & (ErrorCriteria==1 | ErrorCriteria==2) & ...
        (0<=Thresold & Thresold<=1) & ...
		(size(OutputRelation,1)==NbOutput & size(OutputRelation,2)==NbInput & OutputRelation==not(not(OutputRelation))))
  error('Incorrect parameters.'); end

 %General information display
 clc;
 fprintf(1,'Data number used    : %6d\n',NbData);
 fprintf(1,'Input number        : %6d\n',NbInput);
 fprintf(1,'Output number       : %6d\n',NbOutput);
 fprintf(1,'Maximum rule number : %6d\n',NbRuleMax);
 fprintf(1,'Maximal error       : %2.4f\n',ErrorMax);
 fprintf(1,'Fuzzy parameter (m) : %2.4f\n',m);
 fprintf(1,'Thresold            : %2.4f\n',Thresold);
 fprintf(1,'Error criteria      : %6d\n',ErrorCriteria);
 fprintf('************************************************************************\n');
 drawnow;

 %Normalisation
 Range=[min(data);max(data)];
 ScaleFactor=repmat(1,1,NbInput+NbOutput)./(Range(2,:)-Range(1,:));
 NormalisedData=(data-repmat(Range(1,:),NbData,1)).*repmat(ScaleFactor,NbData,1);
 
 %Initialisation with a global linear model
 Model=ModelInitialisation(NormalisedData,NbData,NbInput,NbOutput,ScaleFactor,Range,m,OutputRelation,Thresold,ErrorCriteria); 
 [Error,Indice]=ModelSave(NormalisedData,Model,FileName,ErrorCriteria,Debug); 

 %Growth of the model
 while Error>ErrorMax & Model.NbRule<NbRuleMax
  Model=ModelGrowth(NormalisedData,Model,Indice,Thresold,OutputRelation); 
  [Error,Indice]=ModelSave(NormalisedData,Model,FileName,ErrorCriteria,Debug); 
 end;
 Output=Models;
 Models=[];
%*******************************************************************************%
function [Model] = ModelInitialisation(NormalisedData,NbData,NbInput,NbOutput,ScaleFactor,Range,m,OutputRelation,Thresold,ErrorCriteria)
%*******************************************************************************%
% Subject: Initialize the fuzzy model whit one rule
%*******************************************************************************%
 for Output=NbInput+1:NbInput+NbOutput
  LinearModelInput=find(OutputRelation(Output-NbInput,:));
  if size(LinearModelInput,2)>0
   LinearModel=[NormalisedData(:,LinearModelInput) repmat(1,NbData,1)]\NormalisedData(:,Output);
   GlobalLinearModel(:,Output-NbInput)=zeros(NbInput+1,1);
   GlobalLinearModel(LinearModelInput,Output-NbInput)=LinearModel(1:size(LinearModelInput));
   GlobalLinearModel(NbInput+1,Output-NbInput)=LinearModel(size(LinearModelInput,2)+1);
  else
   GlobalLinearModel(:,Output-NbInput)=zeros(NbInput+1,1);
   GlobalLinearModel(NbInput+1,Output-NbInput)=mean(NormalisedData(:,Output));
  end
 end
 ErrorOutput=...
  abs(NormalisedData(:,NbInput+1:NbInput+NbOutput)-[NormalisedData(:,1:NbInput) repmat(1,NbData,1)]*GlobalLinearModel);

 if ErrorCriteria==1
  Error=max(ErrorOutput,[],2);
 else
  Error=(sum(ErrorOutput.^2,2)).^0.5;
 end
 [ErrorMin,IndiceError]=min(Error);

 DeltaData=NormalisedData-repmat(NormalisedData(IndiceError,:),NbData,1);

 for Output=NbInput+1:NbInput+NbOutput
  LinearModelInput=find(OutputRelation(Output-NbInput,:));
  if size(LinearModelInput,2)>0
   LinearModel=zeros(1,NbInput+1);
   LinearModel(LinearModelInput)=DeltaData(:,LinearModelInput)\DeltaData(:,Output);
   LinearModel(NbInput+1)=...
    NormalisedData(IndiceError,Output)-...
	NormalisedData(IndiceError,1:NbInput)*LinearModel(1:NbInput)';
   else
    LinearModel=zeros(1,NbInput+1);
    LinearModel(NbInput+1)=NormalisedData(IndiceError,Output);
   end
   LocalLinearModel(1,:,Output-NbInput)=LinearModel;
  end

 Model=struct(...
  'NbInput',NbInput,...
  'NbOutput',NbOutput,...
  'NbRule',1,...
  'm',m,...
  'ScaleFactor',ScaleFactor,...
  'Range',Range,...
  'MembershipFunctionCenter',[[NormalisedData(IndiceError,:)]],...
  'LocalLinearModel',LocalLinearModel,...
  'NbData',NbData,...
  'Thresold',Thresold,...
  'ErrorCriteria',ErrorCriteria...
              );
%*******************************************************************************%
function [MembershipDegree] = ModelMembership(NormalisedData,Model);
%*******************************************************************************%
% Subject: Calculate the membership degrees
%*******************************************************************************%
 MembershipDegree=zeros(Model.NbData,Model.NbRule);
 
 %Membership degree calculation
 x1=MembershipDegree;
 x2=0;
 for Rule=1:Model.NbRule
  MembershipFunctionCenter=...
   repmat(Model.MembershipFunctionCenter(Rule,:),Model.NbData,1);
  SquareDistance=...
   sum(((NormalisedData(:,1:Model.NbInput)-MembershipFunctionCenter(:,1:Model.NbInput))).^2,2);
  x1(:,Rule)=...
    SquareDistance.^(1/(Model.m-1));
  if not(all(SquareDistance))
   %No division by zero
   position=find(SquareDistance==0);
   SquareDistance(position)=ones(size(position));
  end   
  x2=...
   x2+ones(Model.NbData,1)./(SquareDistance.^(1/(Model.m-1)));
 end
 for Data=1:Model.NbData
  if all(x1(Data,:))
   MembershipDegree(Data,:)=ones(1,Model.NbRule)./((x1(Data,:).*repmat(x2(Data),1,Model.NbRule)));
  else %NormalisedData(Data,:) is equal to a membership function center
   MembershipDegree(Data,:)=not(x1(Data,:));
   MembershipDegree(Data,:)=MembershipDegree(Data,:)/sum(MembershipDegree(Data,:));
  end
 end
%*******************************************************************************%
function [ErrorMax,IndiceMax] = ModelSave(NormalisedData,Model,FileName,ErrorCriteria,Debug);
%*******************************************************************************%
% Subject: Save the current and previous fuzzy models to the disk and calculate the model errors
%*******************************************************************************%
 global Models;

 ErrorOutput=abs(NormalisedData(:,Model.NbInput+1:Model.NbInput+Model.NbOutput)-ModelOutput(NormalisedData,Model));
 ErrorMeanOutput=mean(ErrorOutput);
 ErrorMaxOutput=max(ErrorOutput,[],1);

 if ErrorCriteria==1
  Error=max(ErrorOutput,[],2);
 else
  Error=(sum(ErrorOutput.^2,2)).^0.5;
 end

 [ErrorMax,IndiceMax]=max(Error);
 ErrorMean=mean(ErrorMeanOutput);
 
 fprintf('Number of rules: %4d\tMaximal error: %4.4f \tMean error: %4.4f\n',Model.NbRule,ErrorMax,ErrorMean);

 %Save the current fuzzy model and the previous ones 
 CurrentModel=struct(...
   'MeanErrorOutput',ErrorMeanOutput,...
   'MeanError',ErrorMean,...
   'MaxErrorOutput',ErrorMaxOutput,...
   'MaxError',ErrorMax,...
   'Model',ModelDenormalize(Model)...
               );
 if isempty(Models)
  Models=[CurrentModel];
 else
  Models=[Models CurrentModel];
 end
 save(FileName,'Models');

 if Debug
  %Plot the mean and maximal error evolution
  clf;
  for Output=1:Model.NbOutput
   subplot(Model.NbOutput,1,Output)
   hold on;  
   ylabel(sprintf('Output %1d',Output));
   axis([0,length(Error),0,1]);
   plot(ErrorOutput(:,Output));
   plot([1 length(ErrorOutput)],[ErrorMeanOutput(1,Output) ErrorMeanOutput(1,Output)],'r');
   plot(IndiceMax,abs(ErrorOutput(IndiceMax,Output)),'ro');
  end
 end
 drawnow;
%*******************************************************************************%
function [NewModel] = ModelGrowth(NormalisedData,Model,Center,Thresold,OutputRelation)
%*******************************************************************************%
% Subject: Add a new rule to the current fuzzy model
%*******************************************************************************%
 NewModel=Model;
 %New rule
 NewModel.NbRule=Model.NbRule+1;
 NewModel.MembershipFunctionCenter=[Model.MembershipFunctionCenter;[NormalisedData(Center,:)]];
 MembershipDegree=ModelMembership(NormalisedData,NewModel); 

 DataUsedMap=MembershipDegree>=repmat(max(MembershipDegree,[],2)*Thresold,1,NewModel.NbRule);

 for Rule=1:NewModel.NbRule
  %Local linear model calculation
  DataUsed=find(DataUsedMap(:,Rule));
  if size(DataUsed,1)<max(sum(OutputRelation,2))+1
   % Not enough datas to calculate the local linear model
   warning('Not enough datas to calculate the local linear model');
   [Degree,Membership]=sort(MembershipDegree(:,Rule));
   DataUsed=Membership(NewModel.NbData-max(sum(OutputRelation,2)):NewModel.NbData);
  end
  DeltaData=NormalisedData(DataUsed,:)-repmat(NewModel.MembershipFunctionCenter(Rule,:),size(DataUsed,1),1);
  DeltaData=DeltaData.*repmat(MembershipDegree(DataUsed,Rule),1,NewModel.NbInput+NewModel.NbOutput);

  for Output=NewModel.NbInput+1:NewModel.NbInput+NewModel.NbOutput
   LinearModelInput=find(OutputRelation(Output-NewModel.NbInput,:));
   if size(LinearModelInput,2)>0
    LinearModel=zeros(1,NewModel.NbInput+1);
    LinearModel(LinearModelInput)=DeltaData(:,LinearModelInput)\DeltaData(:,Output);
    LinearModel(NewModel.NbInput+1)=...
     NewModel.MembershipFunctionCenter(Rule,Output)-...
	  NewModel.MembershipFunctionCenter(Rule,1:NewModel.NbInput)*LinearModel(1:NewModel.NbInput)';
   else
    LinearModel=zeros(1,NewModel.NbInput+1);
    LinearModel(NewModel.NbInput+1)=NewModel.MembershipFunctionCenter(Rule,Output);
   end
   NewModel.LocalLinearModel(Rule,:,Output-NewModel.NbInput)=LinearModel;
  end
 end
%*******************************************************************************%
function [ModelOutput] = ModelOutput(NormalisedData,Model);
%*******************************************************************************%
% Subject: Calculate the estimated outputs with the fuzzy model
%*******************************************************************************%
 MembershipDegree=zeros(Model.NbData,Model.NbRule);
 
 %Membership degree calculation
 x1=MembershipDegree;
 x2=0;
 for Rule=1:Model.NbRule
  MembershipFunctionCenter=...
   repmat(Model.MembershipFunctionCenter(Rule,:),Model.NbData,1);
  SquareDistance=...
   sum(((NormalisedData(:,1:Model.NbInput)-MembershipFunctionCenter(:,1:Model.NbInput))).^2,2);
  x1(:,Rule)=...
    SquareDistance.^(1/(Model.m-1));
  if not(all(SquareDistance))
   %No division by zero
   position=find(SquareDistance==0);
   SquareDistance(position)=ones(size(position));
  end   
  x2=...
   x2+ones(Model.NbData,1)./(SquareDistance.^(1/(Model.m-1)));
 end
 for Data=1:Model.NbData
  if all(x1(Data,:))
   MembershipDegree(Data,:)=ones(1,Model.NbRule)./((x1(Data,:).*repmat(x2(Data),1,Model.NbRule)));
  else %NormalisedData(Data,:) is equal to a membership function center
   MembershipDegree(Data,:)=not(x1(Data,:));
   MembershipDegree(Data,:)=MembershipDegree(Data,:)/sum(MembershipDegree(Data,:));
  end
 end
 
 %Output calculation
  ModelOutput=zeros(Model.NbData,Model.NbOutput);
  for Output=1:Model.NbOutput
   ModelOutput(:,Output)=...
    sum((MembershipDegree*Model.LocalLinearModel(:,:,Output)).*[NormalisedData(:,1:Model.NbInput) ones(Model.NbData,1)],2);
  end
%*******************************************************************************%
function [DenormalizedModel] = ModelDenormalize(Model)
%*******************************************************************************%
% Subject: Denormalize the fuzzy model
%*******************************************************************************%
 DenormalizedModel=Model;
 
 % Denormalize membership function centers
 for Rule=1:Model.NbRule
  DenormalizedModel.MembershipFunctionCenter(Rule,:)= ...
   (Model.MembershipFunctionCenter(Rule,:)./Model.ScaleFactor)+Model.Range(1,:);
 end

 % Denormalize loacl linear models
 ModelOutput=zeros(Model.NbInput+1,Model.NbOutput);
 for Rule=1:Model.NbRule
  for Output=1:Model.NbOutput
   ModelOutput(1:Model.NbInput+1,Output)=Model.LocalLinearModel(Rule,1:Model.NbInput+1,Output)';
  end
  %Denormalized local linear model
  for Input=1:Model.NbInput
   ModelOutput(Input,:)=ModelOutput(Input,:).*...
    repmat(Model.ScaleFactor(Input),1,Model.NbOutput)./Model.ScaleFactor(Model.NbInput+1:Model.NbInput+Model.NbOutput);
  end
  ModelOutput(Model.NbInput+1,:)=...
   ModelOutput(Model.NbInput+1,:)./Model.ScaleFactor(Model.NbInput+1:Model.NbInput+Model.NbOutput)+...
   Model.Range(1,Model.NbInput+1:Model.NbInput+Model.NbOutput);
  for Input=1:Model.NbInput
   ModelOutput(Model.NbInput+1,:)=...
    ModelOutput(Model.NbInput+1,:)-ModelOutput(Input,:)*Model.Range(1,Input);
  end
  DenormalizedModel.LocalLinearModel(Rule,:,:)=ModelOutput;
 end
%*******************************************************************************%
