% evaluate fitnesses
        try
            fitnesses = feval(func, pop);
        catch 
            error('swarm:function_error', ...
                 ['Swarm cannot continue because the supplied cost function ',...
                  'gave the following error:\n %s'], lasterr);
        end
        if (numel(fitnesses) ~= swarmsize)
            error('swarm:function_not_vectorized_or_incorrect_dimensions', ...
                 ['The user-supplied cost function does not appear to get enough arguments,\n',...
                  'or is not vectorized properly. Make sure the dimensions of the limits [Lb]\n',...
                  'and [Ub] are the same as the required input vector for the function, or that\n',...
                  'the function is vectorized properly.'])            
        end
        
        % update evaluations
        evals = evals + numel(fitnesses);
        
        % global best 
        [globalbestfit, index] = min(fitnesses);
        globalbest = pop(index, :);        
       
        % improving solutions
        if (globalbestfit < previousbestfit)
            
            % new best individual
            improvement     = maxiterinpop;
            previousbestfit = globalbestfit;
            previousbest    = globalbest;
            
            % assign also the globals
            SWARM_bestfval = globalbestfit;
            SWARM_bestind  = globalbest;
            
            % display improving solution
            if display                
                if converging1
                   fprintf(1, converge1bck);
                   pause(0.05)
                end
                if converging2
                   fprintf(1, converge2bck);
                   pause(0.05)
                end 
                if (globalbestfit < 0)
                    fprintf('\b')
                end
                fprintf(1, fitbackspace);
                fprintf(1, '%1.8e', globalbestfit);
                converging1 = false;
                firsttime   = true;
                pause(0.05)
            end 
        end      
        
        % Check for convergence
        if (evals > maxfevals)

            % maximum allowed function evaluations has been superceded
            if display
                fprintf(1, overfevals1);
                fprintf(1, overfevals2);
            end
            break
            
        end
        switch convergence
            
            % exhaustive 
            case 1
                if (globalbestfit == previousbestfit) && ...
                   (previousbestfit ~= inf)
                   if (improvement <= 100) && display
                        if ~converging1
                           fprintf(1, convergestr);
                           fprintf(1, '%3.0f', improvement-1);                   
                           converging1 = true;
                           pause(0.05)
                        else
                            fprintf(1, '\b\b\b%3.0f', improvement-1);
                            pause(0.05)
                        end    
                    end           
                    improvement = improvement - 1;
                end
                
            % max iterations 
            case 2                   
                if display && (maxiters-iters < 100)            
                    if firsttime                
                        fprintf(1, itersstr);
                        fprintf(1, '%3.0f', maxiters-iters);                
                        pause(0.05)
                    else
                        fprintf(1, '\b\b\b%3.0f', maxiters-iters);
                        pause(0.05)
                    end
                    firsttime   = false;
                    converging2 = true;
                end  
                iters = iters + 1;
                
            % goal-attain
            case 3
                if (globalbestfit <= convvalue)
                    break;
                end
        end 
        
        % update individuals' local fitnesses
        localbetterfit = (fitnesses < localbestfit);
        localbestfit(localbetterfit) = fitnesses(localbetterfit);
        localbest(localbetterfit, :) = pop(localbetterfit, :);
        
        % update positions
        prevpop = pop;
        pop     = pop + velocity;

        % enforce boundaries (bounce against bounds)
        outsiders1 = (pop < mins);
        outsiders2 = (pop > maxs);
        if any(outsiders1(:)) || any(outsiders2(:))
            pop(outsiders1) = pop(outsiders1) - velocity(outsiders1);
            pop(outsiders2) = pop(outsiders2) - velocity(outsiders2);
            velocity(outsiders1) = -velocity(outsiders1);
            velocity(outsiders2) = -velocity(outsiders2);
        end
        
        % random vectors
        rs = rand(swarmsize, 3);      
        r1 = rs(:, 1*dimsrep);
        r2 = rs(:, 2*dimsrep);
        r3 = rs(:, 3*dimsrep);
        
        % get best fitness from neighbor           
        [bestneighborfit, neighbor] = max(localbestfit(neighbors), [], 2);                
        
        % find actual neighbor
        neighbor     = neighbor + neighinds*numneighbors;
        flipneighbor = neighbors';
        neighbor     = flipneighbor(neighbor);        
        
        % catch best neighbors 
        neighborbest = localbest(neighbor, :);
        
        % update swarms's velocities
        globalbest = globalbest(swarmrep, :);
        velocity   = omega*velocity ...
                     + eta1 * r1.*(neighborbest - pop)...   
                     + eta2 * r2.*(globalbest   - pop)...
                     + eta3 * r3.*(localbest    - pop);                     
                 
         % limit velocity
         Velsign    = sign(velocity);
         outsiders1 = abs(velocity) > abs(velUb);
         outsiders2 = abs(velocity) < abs(velLb);
         if any(outsiders1(:)) || any(outsiders2(:))
             velocity(outsiders1) = Velsign(outsiders1).*velUb(outsiders1);
             velocity(outsiders2) = Velsign(outsiders2).*velLb(outsiders2);
         end         
    end
    
%%  (pre-) end values

    % if solution has been found
    if isfinite(previousbestfit)
        
        % when called normally
        if (~skippop)
            fval = previousbestfit;
            sol  = previousbest;
        else
            % when called from GODLIKE
            fval = fitnesses;
            sol  = prevpop;
        end
        if display
            fprintf(1, '\nSwarm optimization algorithm has converged.\n');
        end
        
    % all trials might be infeasible 
    else
        fprintf(1,'\n');
        warning('swarm:no_solution',...
              'SWARM was unable to find any solution that gave finite values.')
        fval = previousbestfit;
        sol  = NaN;
    end
    
%%  Grace function evaluations

    if (grace > 0)
        
        % display progress
        if display
            fprintf(1, 'Performing direct-search...');
            pause(0.05)
        end
        
        % perform direct search
        options = optimset('TolX', eps, 'MaxFunEvals', grace, 'TolFun', ...
                           eps, 'MaxIter', inf, 'Display', 'off');
        [soltry, fvaltry] = fminsearch(func, sol, options);
        
        % enforce boundaries
        if ~any(soltry >= ub | soltry <= lb)
           sol  = soltry;
           fval = fvaltry;
        end    
        evals = evals + grace;        
    end
    
%%  finalize
    
    % display progress
    if display
        fprintf(1, 'All done.\n\n');
        pause(0.05)
    end 
    
    % clear the temp globals
    clear global SWARM_bestfval SWARM_bestind
    
end

% parser function to easily parse the input arguments
function [pop, func, swarmsize, lb, ub, grace, display, maxfevals, convmethod, convvalue, ...
          eta1, eta2, eta3, omega, numneighbors] = parseprob(problem)
           
        func       = problem.costfun;
        swarmsize  = problem.popsize; 
        lb         = problem.lb;
        ub         = problem.ub; 
        grace      = problem.grace;
        display    = problem.display;   
        maxfevals  = problem.maxfevals;
        convmethod = problem.conv.method;
        convvalue  = problem.conv.value;
        pop           = problem.PS.pop;
        eta1          = problem.PS.eta1; 
        eta2          = problem.PS.eta2; 
        eta3          = problem.PS.eta3; 
        omega         = problem.PS.omega;  
        numneighbors  = problem.PS.numneigh; 
        
end