function [fb,xb,con] = ESGHS(alg,prob)
% HS: Harmoney Search for optimization problem with continuous variables

% the 'alg' struct may include the following fields:
% 1.hms:  size of harmony memory
% 2.HM:   Harmony memory
% 3.NI:   number of improvisations

% the struct of problem includes the following fields:
% 1. dim:  dimension of problem to be solved
% 2. lb:  lower bound for each decision variable
% 3. ub:  uper bound for each decision variable
% 4. fun:   functional handle

% fb: the best value of the objective function
% xb: the best  solution
% con:	 archived convergence

% Edited by: Kaiping Luo, Beihang University, lkpyonghengyiran@163.com
% Updated time 2018-5-30

n = prob.dim;
lb = prob.lb';
ub = prob.ub';

hms = alg.hms;
HM = rand(hms,n).*repmat(ub-lb,hms,1)+repmat(lb,hms,1);
FHM = zeros(hms,1);
con = zeros(alg.NFE,1);
FHM(1) = prob.fun(HM(1,:)');
con(1) = FHM(1);
for i = 2:hms
    FHM(i) = prob.fun(HM(i,:)');
    con(i) = min(FHM(i),con(i-1));
end
t = 1+hms;
[fw,w] = max(FHM);
[fb,b] = min(FHM);
xb = HM(b,:);

NI = alg.NI;
xnew = zeros(1,n);
for k = 1:NI % NI+hms=NFE
    % ------------improvise a new harmoney, begin
    hmcr = (n+randn)/(n+1);
    par = 1-(k-1)/NI;
    delta = (1-((k-1)/NI)^2);
    for j = 1:n
        if rand < hmcr %  select some values from memory 
            xnew(j) = HM(ceil(rand*hms), j);
            if rand<par
                bw = abs(HM(ceil(rand*hms),j)-xnew(j));
                if bw==0
                    bw = (ub(j)-lb(j))*exp((lb(j)-ub(j))*k/NI);
                end
                xnew(j) = xnew(j)+(2*rand-1)*bw;
            end
        else 
            xnew(j) = xb(j)+delta*randn;
        end
        if xnew(j) < lb(j) || xnew(j) > ub(j)
            xnew(j) = lb(j)+rand*(ub(j)-lb(j));
        end
    end
    % ------------improvise a new harmoney, end
    
    % -----update memory, begin
    fnew = prob.fun(xnew');
    if fnew < fw 
        HM(w,:) = xnew;
        FHM(w) = fnew;
        if fnew < fb
            xb = xnew;
            fb = fnew;
        end
        [fw,w] = max(FHM);
    end
    con(t) = fb;
    t = t+1;
    if t>alg.NFE
        return
    end
    % -----update memory, end
end
