gusucode.com > 基于互信息的图像配准matlab源码程序 > 基于互信息的图像配准matlab源码程序/PSO.m
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];