%% Algorithmic Parameter Setting
AlgorithmParams.NumOfCountries = 80;               % Number of initial countries.
AlgorithmParams.NumOfInitialImperialists = 8;      % Number of Initial Imperialists.
AlgorithmParams.NumOfAllColonies = AlgorithmParams.NumOfCountries - AlgorithmParams.NumOfInitialImperialists;
AlgorithmParams.NumOfDecades = 30;
AlgorithmParams.RevolutionRate = 0.3;               % Revolution is the process in which the socio-political characteristics of a country change suddenly.
AlgorithmParams.AssimilationCoefficient = 2;        % In the original paper assimilation coefficient is shown by "beta".
AlgorithmParams.AssimilationAngleCoefficient = .5;  % In the original paper assimilation angle coefficient is shown by "gama".
AlgorithmParams.Zeta = 0.02;                        % Total Cost of Empire = Cost of Imperialist + Zeta * mean(Cost of All Colonies);
AlgorithmParams.DampRatio = 0.99;
AlgorithmParams.StopIfJustOneEmpire = false;         % Use "true" to stop the algorithm when just one empire is remaining. Use "false" to continue the algorithm.
AlgorithmParams.UnitingThreshold = 0.02;            % The percent of Search Space Size, which enables the uniting process of two Empires.

%% Start of Program
clc
clear
close all

%% Input Output Data
Data = xlsread('Data.xls');

X = Data(:,1:end-1);
Y = Data(:,end);

NumOfInputs = size(X,2);
NumOfOutputs = size(Y,2);
NumOfData = size(X,1);

%% Normalization
XN = X;
YN = Y;

MinX = min(X); MaxX = max(X);
MinY = min(Y); MaxY = max(Y);

for ii = 1:size(X,2)
    XN(:,ii) = Normalize_Fcn(X(:,ii),MinX(ii) , MaxX(ii));
end

for ii = 1:size(Y,2)
    YN(:,ii) = Normalize_Fcn(Y(:,ii),MinY(ii) , MaxY(ii));
end

%% Test and Train Data
TrPercent = 80;
DataNum = size(X,1);
TrNum = round(TrPercent * DataNum/100);
R = randperm(DataNum);
TrInx = R(1:TrNum);
TsInx = R(TrNum+1:end);

Xtr = XN(TrInx,:);
Ytr = YN(TrInx,:);

Xts = XN(TsInx,:);
Yts = YN(TsInx,:);

%% Network Structure
pr = [-1 1];
PR = repmat(pr,NumOfInputs,1);
Network = newff(PR,[5 NumOfOutputs],{'tansig' 'tansig'});
Network.trainparam.goal = .0001;

%% Network Training
[Network BestCost] = TrainUsingICA_Fcn(Network,Xtr,Ytr);

%% Network Assessment
YtrNet = sim(Network,Xtr')';
YtsNet = sim(Network,Xts')';

MSEtr = mse(YtrNet-Ytr)
MSEts = mse(YtsNet-Yts)

CorrelatioTr = corrcoef(YtrNet,Ytr);
CorrelatioTs = corrcoef(YtsNet,Yts);

CorrelatioTr = CorrelatioTr(2)
CorrelatioTs = CorrelatioTs(2)

%% Display
figure
plot(Ytr,'-or'); hold on;
plot(YtrNet,'-sk'); hold off;
xlabel('Index'); ylabel('Output')
title('Real and Network output for Train Data');
legend('Real Output','Network Output')

figure
plot(Yts,'-or'); hold on;
plot(YtsNet,'-sk'); hold off;
xlabel('Index'); ylabel('Output')
title('Real and Network output for Test Data');
legend('Real Output','Network Output')

figure
t = -1:.1:1;
plot(t,t,'linewidth',2)
hold on
plot(Ytr,YtrNet,'ok')
xlabel('Real Output'); ylabel('Network Output')
title('Network output vs. Real Output for Train Data');

figure
t = -1:.1:1;
plot(t,t,'linewidth',2)
hold on
plot(Yts,YtsNet,'ok')
xlabel('Real Output'); ylabel('Network Output')
title('Network output vs. Real Output for Test Data');

MSEtr =
    0.0055

MSEts =
    0.0057

CorrelatioTr =
    0.9895

CorrelatioTs =
    0.9863