function [OUT,varargout]=PSO(structure)
D=3;
rand('state',sum(100*clock));
if nargin < 1
   error('Not enough arguments.');
end

    
% PSO PARAMETERS
   VRmin=ones(D,1)*-20; 
   VRmax=ones(D,1)*20;    
   VR=[VRmin,VRmax];
   minmax = 1;

P =[1 2000 20 4 2 2 0.9 0.2 1500 2 1e-5 20 1];
df  = P(1);
me  = P(2);
ps  = P(3);
mv  = P(4);
ac1 = P(5);
ac2 = P(6);
iw1 = P(7);
iw2 = P(8);
iwe = P(9);
flagg=P(10);
ergrd=P(11);
ergrdep=P(12);
plotflg=P(13);

% PLOTTING
 message = sprintf('PSO: %%g/%g iterations, GBest = %%g.\n',me);

pos=40*rand(ps,D)-20;
vel=8*rand(ps,D)-4;

 
% initial pbest positions vals
 pbest=pos;
 
 for j=1:ps  % start particle loop
    numin='0';
    for i=1:D
        numin=strcat(numin,',',num2str(pos(j,i)));
    end
%     evstrg=strcat('feval(''',functname,'''',numin(2:end),',structure',')');  
evstrg=strcat('feval(''myMI''',numin(2:end),',structure',')');
    out(j)=eval(evstrg);     % evaluate desired function with particle j      
 end
 
 pbestval=out;   % initially, pbest is same as pos

% assign initial gbest here also (gbest and gbestval)
 if minmax==1
    [gbestval,idx1]=max(pbestval);  % this picks gbestval when we want to maximize the function
 elseif minmax==0
    [gbestval,idx1]=min(pbestval);  % this works for straight minimization
 end
 gbest=pbest(idx1,:);  % this is gbest position
 tr(1)=gbestval;       % save for output

% start PSO iterative procedures
cnt=0; % counter used for updating display according to df in the options
cnt2=0; % counter used for the stopping subroutine based on error convergence


for i=1:me  % start epoch loop (iterations)
   if flagg==0   % randimization control, one random set for each epoch
       rannum1=rand(1);  
       rannum2=rand(2);
   end
   
   for j=1:ps  % start particle loop
       
     if flagg==1   % randomization control, one random set for each particle at each epoch
         rannum1=rand(1);
         rannum2=rand(1);
     end

     numin='0';
     for dimcnt=1:D
         numin=strcat(numin,',',num2str(pos(j,dimcnt)));
     end
%      evstrg=strcat('feval(''',functname,'''',numin(2:end),',structure',')'); 
evstrg=strcat('feval(''myMI''',numin(2:end),',structure',')');
     out(j)=eval(evstrg);     % evaluate desired function with particle j  
     e(j) = out(j);              % use to minimize or maximize function to unknown values

     %SSEhist(j) = sumsqr(e);    % sum squared 'error' for jth particle (averages if there is more than one output)
     
    % update pbest to reflect whether searching for max or min of function
     if minmax==0
       if pbestval(j)>=e(j);
          pbestval(j)=e(j);
          pbest(j,:)=pos(j,:);
       end
     elseif minmax==1
       if pbestval(j)<=e(j);
           pbestval(j)=e(j);
           pbest(j,:)=pos(j,:);
       end
     end
 
          
    % assign gbest by finding minimum of all particle pbests 
     if minmax==1
       [iterbestval,idx1]=max(pbestval);  % this picks gbestval when we want to maximize the function
       if gbestval<=iterbestval
          gbestval=iterbestval;
          gbest=pbest(idx1,:);
       end       
     elseif minmax==0 
       [iterbestval,idx1]=min(pbestval);  % this works for straight minimization and for minimizing error to target
       if gbestval>=iterbestval
          gbestval=iterbestval;
          gbest=pbest(idx1,:);
       end       
     end    
    
     tr(i+1)=gbestval; % keep track of global best val
     te=i;           % this will return the epoch number to calling program when done

    % get new velocities, positions (this is the heart of the PSO algorithm)
     if i<=iwe
        iwt(i)=((iw2-iw1)/(iwe-1))*(i-1)+iw1; % get inertia weight, just a linear funct w.r.t. epoch parameter iwe
     else
        iwt(i)=iw2;
     end
    
     if flagg==2              % each component of each particle gets a different random number set
        for dimcnt=1:D
           rannum1=rand(1);
           rannum2=rand(1);
           vel(j,dimcnt)=iwt(i)*vel(j,dimcnt)...
                    +ac1*rannum1*(pbest(j,dimcnt)-pos(j,dimcnt))...
                    +ac2*rannum2*(gbest(1,dimcnt)-pos(j,dimcnt));
        end
     else                     % this velocity update is for flagg= 0 or 1
        vel(j,:)=iwt(i)*vel(j,:)...
                 +ac1*rannum1*(pbest(j,:)-pos(j,:))...
                 +ac2*rannum2*(gbest(1,:)-pos(j,:));
     end

    % update new position
     pos(j,:)=pos(j,:)+vel(j,:);    

    % limit velocity/position components to maximums
     for dimcnt=1:D
        if vel(j,dimcnt)>mv
           vel(j,dimcnt)=mv;
        end
       
        if vel(j,dimcnt)<-mv
          vel(j,dimcnt)=-mv;
        end    
        
        if pos(j,dimcnt)>=VR(dimcnt,2)
           pos(j,dimcnt)=VR(dimcnt,2);
        end
       
        if pos(j,dimcnt)<=VR(dimcnt,1)
           pos(j,dimcnt)=VR(dimcnt,1);
        end       
        
        
     end
       
   end         % end particle loop
   
  %---------------------------------------------------------------------------------------------------------------------- 
   
   
   % check for stopping criterion based on speed of convergence to desired error   
    tmp1=abs(tr(i)-gbestval);
    if tmp1>ergrd
       cnt2=0;
    elseif tmp1<=ergrd
       cnt2=cnt2+1;
       if cnt2>=ergrdep
         if plotflg==1
          fprintf(message,i,gbestval);           
%           disp(' ');
%           disp('***** global error gradient too small for too long');
%           disp('***** this means you''ve got the solution or it got stuck');
         end
         break
       end       
    end     
end  % end epoch loop
      OUT=[gbest';gbestval];
 OUT(1:3)=round(OUT(1:3));
 varargout{1}=[0:te];
 varargout{2}=[tr];