clc;
clear;
close all;

%% Load Data

data=CreateData();

%% Generate Basic FIS

fis=CreateInitialFIS(data,10);

%% Tarin Using PSO

fis=TrainUsingPSO(fis,data);

%% Results

% Train Data
TrainOutputs=evalfis(data.TrainInputs,fis);
PlotResults(data.TrainTargets,TrainOutputs,'Train Data');

% Test Data
TestOutputs=evalfis(data.TestInputs,fis);
PlotResults(data.TestTargets,TestOutputs,'Test Data');

function data=CreateData()
    load('trainingdata.mat'); 
    entrance.data.total =trainingdata(: , [1 2  13 9]); 
    [Inputs, Targets]=abalone_dataset();
    Inputs= entrance.data.total(:,1:3);
    Targets=entrance.data.total(:,4);
    
    nSample=size(Inputs,1);
    
    S=randperm(nSample);
    Inputs=Inputs(S,:);
    Targets=Targets(S,:);
    
    % Train Data
    pTrain=0.7;
    nTrain=round(pTrain*nSample);
    TrainInputs=Inputs(1:nTrain,:);
    TrainTargets=Targets(1:nTrain,:);
    
    % Test Data
    TestInputs=Inputs(nTrain+1:end,:);
    TestTargets=Targets(nTrain+1:end,:);
    
    % Export
    data.TrainInputs=TrainInputs;
    data.TrainTargets=TrainTargets;
    data.TestInputs=TestInputs;
    data.TestTargets=TestTargets;

end
function fis=CreateInitialFIS(data,nCluster)

    if ~exist('nCluster','var')
        nCluster='auto';
    end

    x=data.TrainInputs;
    t=data.TrainTargets;
    
    fcm_U=2;
    fcm_MaxIter=1000;
    fcm_MinImp=1e-5;
    fcm_Display=false;
    fcm_options=[fcm_U fcm_MaxIter fcm_MinImp fcm_Display];
    fis=genfis3(x,t,'sugeno',nCluster,fcm_options);

end
function bestfis=TrainUsingPSO(fis,data)
       
    p0=GetFISParams(fis);
    Problem.CostFunction=@(x) TrainFISCost(x,fis,data);
    Problem.nVar=numel(p0);    
    alpha=1;
    Problem.VarMin=-(10^alpha);
    Problem.VarMax=10^alpha;  
    Params.MaxIt=1000;
    Params.nPop=40;    
    results=RunPSO(Problem,Params);  
    p=results.BestSol.Position.*p0;
    bestfis=SetFISParams(fis,p);
    
end
function results=RunPSO(Problem,Params)

    disp('Starting PSO ...');
    CostFunction=Problem.CostFunction;       
    nVar=Problem.nVar;          
    VarSize=[1 nVar];           
    VarMin=Problem.VarMin;      
    VarMax=Problem.VarMax;      
    MaxIt=Params.MaxIt;     
    nPop=Params.nPop;        
    w=1;            
    wdamp=0.99;     
    c1=1;           
    c2=2;           

    VelMax=0.1*(VarMax-VarMin);
    VelMin=-VelMax;

    empty_particle.Position=[];
    empty_particle.Cost=[];
    empty_particle.Velocity=[];
    empty_particle.Best.Position=[];
    empty_particle.Best.Cost=[];

    particle=repmat(empty_particle,nPop,1);

    BestSol.Cost=inf;

    for i=1:nPop
     
        if i>1
            particle(i).Position=unifrnd(VarMin,VarMax,VarSize);
        else
            particle(i).Position=ones(VarSize);
        end
        
        particle(i).Velocity=zeros(VarSize);
       
        particle(i).Cost=CostFunction(particle(i).Position);
        
        particle(i).Best.Position=particle(i).Position;
        particle(i).Best.Cost=particle(i).Cost;

        if particle(i).Best.Cost<BestSol.Cost

            BestSol=particle(i).Best;
        end
    end

    BestCost=zeros(MaxIt,1);

    for it=1:MaxIt

        for i=1:nPop

           
            particle(i).Velocity = w*particle(i).Velocity ...
                +c1*rand(VarSize).*(particle(i).Best.Position-particle(i).Position) ...
                +c2*rand(VarSize).*(BestSol.Position-particle(i).Position);

            particle(i).Velocity = max(particle(i).Velocity,VelMin);
            particle(i).Velocity = min(particle(i).Velocity,VelMax);          
            particle(i).Position = particle(i).Position + particle(i).Velocity;          
            IsOutside=(particle(i).Position<VarMin | particle(i).Position>VarMax);
            particle(i).Velocity(IsOutside)=-particle(i).Velocity(IsOutside);            
            particle(i).Position = max(particle(i).Position,VarMin);
            particle(i).Position = min(particle(i).Position,VarMax);            
            particle(i).Cost = CostFunction(particle(i).Position);           
            if particle(i).Cost<particle(i).Best.Cost

                particle(i).Best.Position=particle(i).Position;
                particle(i).Best.Cost=particle(i).Cost;

               
                if particle(i).Best.Cost<BestSol.Cost

                    BestSol=particle(i).Best;

                end
            end

        end

        BestCost(it)=BestSol.Cost;

        disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);

        w=w*wdamp;

    end

    disp('End of PSO.');
    disp(' ');
    results.BestSol=BestSol;
    results.BestCost=BestCost;
    
end
function p=GetFISParams(fis)

    p=[];

    nInput=numel(fis.input);
    for i=1:nInput
        nMF=numel(fis.input(i).mf);
        for j=1:nMF
            p=[p fis.input(i).mf(j).params];
        end
    end

    nOutput=numel(fis.output);
    for i=1:nOutput
        nMF=numel(fis.output(i).mf);
        for j=1:nMF
            p=[p fis.output(i).mf(j).params];
        end
    end
    
end
function PlotResults(Targets,Outputs,Name)

    Errors=Targets-Outputs;

    figure;

    subplot(2,2,[1 2]);
    plot(Targets,'r:');
    hold on;
    plot(Outputs,'b');
    legend('Targets','Outputs');
    title([Name ' - Targets and Outputs']);
    ylabel('Targets and Outputs');
    
    ErrorMean=mean(Errors);
    ErrorStD=std(Errors);
    MSE=mean(Errors.^2);
    RMSE=sqrt(MSE);
    
    subplot(2,2,3);
    plot(Errors);
    ylabel('Errors');
    title(['MSE = ' num2str(MSE) ', RMSE = ' num2str(RMSE)]);
    
    nBin=max(round(numel(Errors)/20),10);
    
    subplot(2,2,4);
    histfit(Errors,nBin);
    xlabel('Errors');
    title(['Mean = ' num2str(ErrorMean) ', StD = ' num2str(ErrorStD)]);
    
end
function [z, out]=TrainFISCost(x,fis,data)

    MinAbs=1e-5;
    if any(abs(x)<MinAbs)
        S=(abs(x)<MinAbs);
        x(S)=MinAbs.*sign(x(S));
    end

    p0=GetFISParams(fis);

    p=x.*p0;
    
    fis=SetFISParams(fis,p);
    
    x=data.TrainInputs;
    t=data.TrainTargets;
    y=evalfis(x,fis);
    
    e=t-y;
    
    MSE=mean(e(:).^2);
    RMSE=sqrt(MSE);
    
    z=RMSE;
    
    out.fis=fis;
   
    out.MSE=MSE;
    out.RMSE=RMSE;
    
end
function fis=SetFISParams(fis,p)

    nInput=numel(fis.input);
    for i=1:nInput
        nMF=numel(fis.input(i).mf);
        for j=1:nMF
            k=numel(fis.input(i).mf(j).params);
            fis.input(i).mf(j).params=p(1:k);
            p(1:k)=[];
        end
    end

    nOutput=numel(fis.output);
    for i=1:nOutput
        nMF=numel(fis.output(i).mf);
        for j=1:nMF
            k=numel(fis.output(i).mf(j).params);
            fis.output(i).mf(j).params=p(1:k);
            p(1:k)=[];
        end
    end

end