clear all
format long g

% Opt1.m is a file used to find the optimum solutions 
% for a minimizaton problem in Test Examples, EG.111.
%
% 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=75; 
rho=0.1; alpha=0.5; 

c=[-6.089 -17.164 -34.054 -5.914 -24.721 -14.986 -24.100 -10.708 -26.662 -22.179];
n=size(c,2);

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

% Initial mean
mu=-4*ones(1,10);

% Best known solution by genetic algorithm
xstar1=[-3.201212 -1.912060 -0.2444413 -6.537489 -0.7231524 -7.267738 -3.596711 -4.017769 -3.287462 -2.335582];
fx1=S(xstar1);
mu = xstar1; sigma=0.00001*ones(1,10);
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],11);
gamhat=sortSX(order,11); gam=-47.7;
Sstar = 1000;
% main loop
%while gamhat<gam
for t=0:100
    t=t+1;
    x=gen(N,mu,sigma.^2);
    SX=S(x)';
    sortSX=sortrows([x SX],11);
    mu1=mean(sortSX(1:order,1:10)); % update the mean vector
    %mu1 = sortSX(1,1:10);
    mu=alpha.*mu1+(1-alpha).*mu; % smoothed updating of mu
    sigma1=std(sortSX(1:order,1:10)); % update sigma
    sigma=alpha.*sigma1+(1-alpha).*sigma; %smoothed updating of sigma
    if (sortSX(1,11) < Sstar)
      Sstar = sortSX(1,11);
      xbest = sortSX(1,1:10);
    end
    [t, Sstar, xbest]
end
time=toc

fx=S(mu)
fx1

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