clear all
format long g

% Opt1.m is a file used to find the optimum solutions 
% for a minimizaton problem in Test Examples, EG.63.
%
% Assume that x_1, x_2, ...,x_n comes from a multi-variate normal distribution.
%
% Inputs: 1. mu     - initial mean
%         2. sigma  - initial standard deviation
%         3. rho    - rarity parameter
%         4. alpha  - smoothing parameter
%         5. N      - sample size
%         6. S      - performance function
%
% Outputs: 1. xhat    - optimal solution for x_j, j=1,2,...,n.
%          2. fxhat - value of function with optimal x_j, j=1,2,...,n.
%          3. time    - CPU time

% Initial CE parameters
N=50; n=3;
rho=0.2; alpha=0.3; 

%c=[3.512118414 0.2169881741 3.552174034];

% Initial covariance matrix
sigma=ones(1,3);

% Initial mean
mu=zeros(1,3);

% Best known solution by genetic algorithm
xstar1=[3.512118414 0.2169881741 3.552174034];
fx1=S(xstar1);

order=ceil(rho*N);

tic
t=0; %iteration counter
x=gen(N,mu,sigma.^2); % Generate multi-variate truncated normal distribution X.
SX = S(x)';
sortSX=sortrows([x SX],n+1);
gamhat=sortSX(order,n+1);
Sstar = 1000;
% main loop
%while gamhat<gam
for t=0:200
    t=t+1;
    x=gen(N,mu,sigma.^2);
    SX=S(x)';
    sortSX=sortrows([x SX],n+1);
    %mu1=mean(sortSX(1:order,1:10)); % update the mean vector
    mu1 = sortSX(1,1:n);
    mu=alpha.*mu1+(1-alpha).*mu; % smoothed updating of mu
    sigma1=std(sortSX(1:order,1:n)); % update sigma
    sigma=alpha.*sigma1+(1-alpha).*sigma; %smoothed updating of sigma
    if (sortSX(1,n+1) < Sstar)
      Sstar = sortSX(1,n+1);
      xbest = sortSX(1,1:n);
    end
    [t, Sstar, xbest]
end
time=toc

fx=S(mu)
fx1

format short g
% Compute relative experimental error
RE=(fx-fx1)/fx1