%% computing RGB input values for a display to obtain desired XYZ outputs
%clear all;
close all;
addpath(genpath('fmid'));


%Get data ready for training and testing
% Load Data 
%load('DataSets_MSIQled.mat')

% TS Modeling via Clustering
%Set U inputs and Y outputs
Dat.U = TrainingMeasurementXYZn; %we will need this struct for training later
Dat.Y = TrainingMeasurementRGB;
%Set testing/validation data
u = ValidationMeasurementXYZn; %no need for struct for validation/testing
y = ValidationMeasurementRGB;

%Cluster parameters: optimization needed!
Par.ante = 1; %Antedecent type: means no dynamic (other values are used for dynamic systems that use as inputs previous outputs of the system itself). 
Par.tol = 0.005; % tolerance for cluster divisions
Par.m = 1.25; %overlapping parameter for the fuzzy clusters (higher more fuzziness and overlapping). 2 implies 50% overlapping approx.
%are outputs independent or correlated?? Depending on this choose next
Par.unifycenter = 0; % 1-Only one clustering for all outputs, 0-Different clustering for each output
%choose type of local approximation per output
%Lineal
Par.cons = 1*ones(1,3); % 1-Linear Model + Offset, 6-Parabolic Model + offset, 7-Logaritmic Model + offset, 8-Sqrt Model + offset, ... 

h = waitbar(0,'Initializing waitbar...');
% Init BestModel var
% Try different number of clusters/rules i
%Set default error
for i = 2:1:15
    BestModel{i}.MSE = 10^6;
    BestModel{i}.MSE_Learn = 10^6;
end

% Monte Carlo Simulation from 2 clusters/rules to 10
for i = 2:1:15
    Par.c = i*ones(1,3); %i clusters per output whose initial values the method will randomize
    for j=1:10000 %10000 monte carlo simulations 
        clear FM; %just in case
        % Learning
        [FM,Mu] = colorfmclust(Dat,Par); %function to cluster data using the provided parameters
        [ym_learning,VAF_Learn] = colorfmsim(Dat.U,Dat.Y,FM); %compute training error VAF in [0-100]        
        MSE_Learn = mean(sum((ym_learning'-Dat.Y').^2)./3); %MSE error is other figure of merit that we use in a NN for same purpose
        % Validation
        [ym,VAF] = colorfmsim(u,y,FM); %compare this VAF or other figure of merit to detect potential overfitting
        MSE = mean(sum((ym'-y').^2)./3); %MSE error is other figure of merit that we use in a NN for same purpose
        
        if MSE < BestModel{i}.MSE %is best so far? we check mean because VAF has one value per output
            BestModel{i}.MSE = MSE;
            BestModel{i}.MSE_Learn = MSE_Learn;
            BestModel{i}.FM = FM;
            BestModel{i}.Mu = Mu;
        end
        if mod(j,2) == 0 %progress bar
            waitbar(((i-2)*100000+j)/(9*10000),h,strcat('Progress:','Clusters = ',num2str(i),'- Iteration = ',num2str(j)));   
        end
    end
end
save('FMLineal','BestModel')

%%
%Pot2
Par.cons = 6*ones(1,3); % 1-Linear Model + Offset, 6-Parabolic Model + offset, 7-Logaritmic Model + offset, 8-Sqrt Model + offset, ... 

h = waitbar(0,'Initializing waitbar...');
% Init BestModel var
% Try different number of clusters/rules i
%Set default error
for i=2:1:15
    BestModel{i}.MSE = 10^6;
    BestModel{i}.MSE_Learn = 10^6;
end

% Monte Carlo Simulation from 2 clusters/rules to 10
for i=2:1:15
    Par.c = i*ones(1,3); %i clusters per output whose initial values the method will randomize
    for j=1:10 %10000 monte carlo simulations 
        clear FM; %just in case
        % Learning
        [FM,Mu] = colorfmclust(Dat,Par); %function to cluster data using the provided parameters
        [ym_learning,VAF_Learn] = colorfmsim(Dat.U,Dat.Y,FM); %compute training error VAF in [0-100]        
        MSE_Learn = mean(sum((ym_learning'-Dat.Y').^2)./3); %MSE error is other figure of merit that we use in a NN for same purpose
        % Validation
        [ym,VAF] = colorfmsim(u,y,FM); %compare this VAF or other figure of merit to detect potential overfitting
        MSE = mean(sum((ym'-y').^2)./3); %MSE error is other figure of merit that we use in a NN for same purpose
        
        if MSE < BestModel{i}.MSE %is best so far? we check mean because VAF has one value per output
            BestModel{i}.MSE = MSE;
            BestModel{i}.MSE_Learn = MSE_Learn;
            BestModel{i}.FM = FM;
            BestModel{i}.Mu = Mu;
        end
        if mod(j,2) == 0 %progress bar
            waitbar(((i-2)*10+j)/(9*15),h,strcat('Progress:','Clusters = ',num2str(i),'- Iteration = ',num2str(j)));   
        end
    end
end
save('FMPot2','BestModel')

%%
%Sqrt
Par.cons = 8*ones(1,3); % 1-Linear Model + Offset, 6-Parabolic Model + offset, 7-Logaritmic Model + offset, 8-Sqrt Model + offset, ... 

h = waitbar(0,'Initializing waitbar...');
% Init BestModel var
% Try different number of clusters/rules i
%Set default error
for i=2:1:15
    BestModel{i}.MSE = 10^6;
    BestModel{i}.MSE_Learn = 10^6;
end

% Monte Carlo Simulation from 2 clusters/rules to 10
for i=2:1:15
    Par.c = i*ones(1,3); %i clusters per output whose initial values the method will randomize
    for j=1:10 %10000 monte carlo simulations 
        clear FM; %just in case
        % Learning
        [FM,Mu] = colorfmclust(Dat,Par); %function to cluster data using the provided parameters
        [ym_learning,VAF_Learn] = colorfmsim(Dat.U,Dat.Y,FM); %compute training error VAF in [0-100]        
        MSE_Learn = mean(sum((ym_learning'-Dat.Y').^2)./3); %MSE error is other figure of merit that we use in a NN for same purpose
        % Validation
        [ym,VAF] = colorfmsim(u,y,FM); %compare this VAF or other figure of merit to detect potential overfitting
        MSE = mean(sum((ym'-y').^2)./3); %MSE error is other figure of merit that we use in a NN for same purpose
        
        if MSE < BestModel{i}.MSE %is best so far? we check mean because VAF has one value per output
            BestModel{i}.MSE = MSE;
            BestModel{i}.MSE_Learn = MSE_Learn;
            BestModel{i}.FM = FM;
            BestModel{i}.Mu = Mu;
        end
        if mod(j,2) == 0 %progress bar
            waitbar(((i-2)*10+j)/(9*15),h,strcat('Progress:','Clusters = ',num2str(i),'- Iteration = ',num2str(j)));   
        end
    end
end
save('FMSqrt','BestModel')

%% Plot Results 

load FMLineal
%Extract MSEs per cluster number
FMLinealMSEs = [];
FMLinealMSEs_Learn = [];
for i=2:1:15
    FMLinealMSEs = [FMLinealMSEs BestModel{i}.MSE]; 
    FMLinealMSEs_Learn = [FMLinealMSEs_Learn BestModel{i}.MSE_Learn]; 
end

load FMPot2
FMPot2MSEs = [];
FMPot2MSEs_Learn = [];
for i=2:1:15
    FMPot2MSEs = [FMPot2MSEs BestModel{i}.MSE]; 
    FMPot2MSEs_Learn = [FMPot2MSEs_Learn BestModel{i}.MSE_Learn]; 
end

load FMSqrt
FMSqrtMSEs = [];
FMSqrtMSEs_Learn = [];
for i=2:1:15
    FMSqrtMSEs = [FMSqrtMSEs BestModel{i}.MSE]; 
    FMSqrtMSEs_Learn = [FMSqrtMSEs_Learn BestModel{i}.MSE_Learn]; 
end


%close all
 

% Plot validation and training errors
figure, hold on
%title('Best Fuzzy Models - 10 simulations for each cluster')
axis([2 15 -25 100]);
xlabel('Numer of Clusters - Fuzzy Rules', 'FontSize',12,'FontWeight','bold');
ylabel('MSE','FontSize',12,'FontWeight','bold');

plot(FMLinealMSEs,'r--','LineWidth',2)
plot(FMLinealMSEs_Learn,'r-','LineWidth',2)
plot(FMPot2MSEs,'g--','LineWidth',2)
plot(FMPot2MSEs_Learn,'g-','LineWidth',2)
plot(FMSqrtMSEs,'b--','LineWidth',2)
plot(FMSqrtMSEs_Learn,'b-','LineWidth',2)

legend('Linear model MSE','Linear model MSE Learn', 'Parabolic model MSE','Parabolic model MSE Learn', 'Squared model MSE','Squared model MSE Learn', 'FontSize',12,'FontWeight','bold')
legend boxoff
exportgraphics(gca, 'Plot1.png', 'Resolution',1200)

%% Extract information from the trained model
% Plot Clustering
%Cluster centers
load FMlineal
close all
KRules = 5; %choose model with KRules
output = 3
centers = BestModel{KRules}.FM.V{output}; %output is relevant if clusterization is different per output

%plot
figure
strCenter=['r*','ro','r+','rx','rd','rs'];
plot3(u(:,1),u(:,2),u(:,3),'.')
for i=1:1
   centers = BestModel{KRules}.FM.V{i};
   for j=1:2
       hold on
       plot3(centers(j,1),centers(j,2),centers(j,3),strCenter(j));
       text(centers(j,1),centers(j,2),centers(j,3),strcat(num2str(j),'-',num2str(i)));
   end
end
grid on
%%
%Membership functions
plotmfs(BestModel{8}.FM) %one figure per output (will be identical if output clusterization is common)
% To plot the MFs for a publication setting 
plotmfs1(BestModel{8}.FM)
%Show ouputs for centroids is helpful
[centersOutput,VAF] = colorfmsim(centers,centers,BestModel{8}.FM); %second input is irrelevant as it is desired output, unknown in real cases 
%%
%After choosing best clusterization and output type, obtain predictions
% load FMLineal
% load Training_Data
% load Validation_Data

% u = ValidationXYZ; %no need for struct for validation/testing
% y = ValidationRGB;
[ym,VAF] = colorfmsim(u,y,BestModel{5}.FM);
