function pareto = moead(mop)
%MOEAD run moea/d algorithms for the given mop.
    %global variable definition.
    global subproblems  itrCounter;
    evolve(mop);
    itrCounter=itrCounter+1;
    pareto=[subproblems.curpoint];
end

% The evoluation setp in MOEA/D
function evolve(mop)
  global params subproblems idealpoint ;
   selindex = 1:length(subproblems);  
   selindex =random_shuffle(selindex);
  for i=1:length(subproblems)
    index = selindex(i);
    r = rand;
    useneighbour = r < params.neighbormating;
    ind = genetic_op(index,mop.domain,useneighbour);
    if params.genevolve
       inds(1,i)=ind;
    else
    [obj,ind] = evaluate(mop, ind);
    idealpoint = min(idealpoint, obj);
      if params.globupdate
         globupdate(index, ind);
      else
         origupdate(index, ind, useneighbour);
      end
    end     
    clear ind obj useneighbour;
  end
  if params.genevolve
  [V, INDS] = arrayfun(@evaluate, repmat(mop, size(inds)), inds, 'UniformOutput', 0);
   v = cell2mat(V); 
  idealpoint = min(idealpoint, min(v,[],2));
  oldinds=[subproblems.curpoint];
  oldinds=num2cell(oldinds, 1);
  indiv=get_structure('individual');
  allind = repmat(indiv, 1,length(oldinds)+length(INDS));
  allind(1:length(oldinds))=cell2mat(oldinds);
  allind(length(oldinds)+1:length(oldinds)+length(INDS))=cell2mat(INDS);
  genupdate(allind);
  clear inds INDS oldinds V v indiv allind;
  end
end
% update the population.
% index is the subproblem's index in the main population.
% ind is the individual structure.
% useneighbour is a bool determine whether the neighbourhood of index, or the whole population should be updated.
% this procedure is also governed by a parameter from params: params.updatesize, which determine how many subproblem
% should be updated at most by this new individual in MOEAD-DE.
function globupdate(index, ind)
 global subproblems idealpoint params itrCounter neighborrelation;
 allobj=subobjective( [subproblems.weight], ind.objective,  idealpoint, params.dmethod);
 [~,pos]=min(allobj);
 repl=ind;
 if ~isempty(params.adaptmethod)&&~isempty(params.adaptparams)
 subproblems(index).Tr = agr(params.popsize,params.adaptmethod,params.adaptparams,params.iteration,itrCounter);
 end
 neighborindex= neighborrelation(1:subproblems(index).Tr,pos);
 neighborpoints=[subproblems(neighborindex).curpoint];   
 neighbor_sub=subobjective( [subproblems(neighborindex).weight], [neighborpoints.objective],  idealpoint, params.dmethod);
 repl_sub=subobjective( [subproblems(neighborindex).weight],  repl.objective,  idealpoint, params.dmethod);
 C       =repl_sub <  neighbor_sub;  
 if (length(C(C==1)) <= params.updatesize)
    toupdate = neighborindex(C);
 else
    toupdate = randsample(neighborindex(C), params.updatesize);
 end
 [subproblems(toupdate).curpoint]=deal(repl);   
end
function origupdate(index, ind, updateneighbour)
  global subproblems idealpoint params neighborrelation;
  
  % collect the updation index
  if (updateneighbour)
    updateindex = neighborrelation(1:subproblems(index).Tr,index);
  else
    updateindex = 1:length(subproblems);
  end
  
  updateindex = random_shuffle(updateindex);
  neighborpoints=[subproblems(updateindex).curpoint];  
  newobj  = subobjective([subproblems(updateindex).weight], ind.objective, idealpoint, params.dmethod);    %objective values of current solution
  oldobj  = subobjective([subproblems(updateindex).weight], [neighborpoints.objective], idealpoint, params.dmethod);        %previous objective values
  C       = newobj < oldobj;    %new solution is better?
    if (length(C(C==1)) <= params.updatesize)
        toupdate = updateindex(C);
    else
        toupdate = randsample(updateindex(C), params.updatesize);
    end    
  [subproblems(toupdate).curpoint]=deal(ind);      %repace with the new one
end
function genupdate(allind)
global subproblems idealpoint dirvectors itrCounter params;
  objs=[allind.objective];
  y1=objs-idealpoint(:,ones(1, length(allind)));
  y1long=sqrt(sum(y1.*y1));
  Tr = agr(params.popsize,params.adaptmethod,params.adaptparams,params.iteration,itrCounter);
 for i=1:length(subproblems)
    y2=dirvectors(:,i);
    y2long=sqrt(sum(y2.*y2));
    cos=y2'*y1./(y1long*y2long);
    [~,pos]=sort(cos,'descend');
    inds=allind(pos(1:Tr));
    allobj=subobjective(subproblems(i).weight,[inds.objective], idealpoint, 'te');
    [new,pos_in]=min(allobj);
    old=subobjective( subproblems(i).weight,  subproblems(i).curpoint.objective,  idealpoint, 'te');
    if new<old
    subproblems(i).curpoint=inds(pos_in);
    end
    clear y2 y2long cos pos inds allobj pos_in;
 end
 clear objs y1 y1long;
end
