Home > atmlab > gridcreation > uniformsphere > ParticleSampleSphere.m

ParticleSampleSphere

PURPOSE ^

Create an approximately uniform triangular tessellation of the unit

SYNOPSIS ^

function [V,Tri,Ue_i,Ue]=ParticleSampleSphere(varargin)

DESCRIPTION ^

 Create an approximately uniform triangular tessellation of the unit 
 sphere by minimizing generalized electrostatic potential energy 
 (aka Reisz s-energy) of the system of charged particles. Effectively, 
 this function produces a locally optimal solution to the problem that 
 involves finding a minimum Reisz s-energy configuration of N equal 
 charges confined to the surface of the unit sphere (s=1 corresponds to 
 the problem originally posed by J. J. Thomson). 

 SYNTAX:
 [V,Tri,Ue_i,Ue]=ParticleSampleSphere(option_i,value_i,option_j,value_j,...)

 OPTIONS:
   - 'N'    : desired number of particles. Corresponding value of N must 
              be a positive interger greater than 9 and less than 1001.
              N=200 particles is the default setting. 
   - 'Vo'   : particle positions used to initialize the search.
              Corresponding value of Vo must be a N-by-3 array, where n is
              the number of particles. N=10 is the lowest permissible 
              number of particles. Initializations consisting of more than 
              1E3 particles are admissible but may lead to unreasonably 
              long optimization times.
   - 's'    : Reisz s-energy parameter used to control the strength of
              particle interactions. Corresponding value of s must be
              a real number greater than zero. s=1 is the default setting.
   - 'Etol' : covergence tolerance. Coresponding value of Etol must be
              a real, positive number. Etol=1E-5 is the default setting.
   - 'Nitr' : Maximum number of iterations. Corresponding value of Nitr
              must be a non-negative interger. Nitr=1E3 is the default 
              setting.

 OUTPUT:
   - V     : N-by-3 array of vertex positions.
   - Tri   : M-by-3 list of face-vertex connectivities. 
   - Ue_i  : N-by-1 array of particle energies.
   - Ue    : K-by-1 array of energy scores, where K-1 was the total number 
             of iterations. Ue(1) corresponds to the energy of the initial
             configuration.

 EXAMPLE: 
 -------------------------------------------------------------------------
 % Uniformly distribute 100 particles across the surface of the unit sphere

  [V,Tri,~,Ue]=ParticleSampleSphere('N',100);
  figure, plot(0:numel(Ue)-1,Ue,'.-')
  set(get(gca,'Title'),'String','Optimization Progress','FontSize',20)
  xlabel('Iteration #','FontSize',15)
  ylabel('Reisz s-enrgy','FontSize',15)
  TR=TriRep(Tri,V);
  figure, h=trimesh(TR); set(h,'EdgeColor','b'), axis equal
 -------------------------------------------------------------------------

 AUTHOR: Anton Semechko (a.semechko@gmail.com)
 DATE: June.2012

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

ParticleSampleSphere.m

SOURCE CODE ^

0001 function [V,Tri,Ue_i,Ue]=ParticleSampleSphere(varargin)
0002 % Create an approximately uniform triangular tessellation of the unit
0003 % sphere by minimizing generalized electrostatic potential energy
0004 % (aka Reisz s-energy) of the system of charged particles. Effectively,
0005 % this function produces a locally optimal solution to the problem that
0006 % involves finding a minimum Reisz s-energy configuration of N equal
0007 % charges confined to the surface of the unit sphere (s=1 corresponds to
0008 % the problem originally posed by J. J. Thomson).
0009 %
0010 % SYNTAX:
0011 % [V,Tri,Ue_i,Ue]=ParticleSampleSphere(option_i,value_i,option_j,value_j,...)
0012 %
0013 % OPTIONS:
0014 %   - 'N'    : desired number of particles. Corresponding value of N must
0015 %              be a positive interger greater than 9 and less than 1001.
0016 %              N=200 particles is the default setting.
0017 %   - 'Vo'   : particle positions used to initialize the search.
0018 %              Corresponding value of Vo must be a N-by-3 array, where n is
0019 %              the number of particles. N=10 is the lowest permissible
0020 %              number of particles. Initializations consisting of more than
0021 %              1E3 particles are admissible but may lead to unreasonably
0022 %              long optimization times.
0023 %   - 's'    : Reisz s-energy parameter used to control the strength of
0024 %              particle interactions. Corresponding value of s must be
0025 %              a real number greater than zero. s=1 is the default setting.
0026 %   - 'Etol' : covergence tolerance. Coresponding value of Etol must be
0027 %              a real, positive number. Etol=1E-5 is the default setting.
0028 %   - 'Nitr' : Maximum number of iterations. Corresponding value of Nitr
0029 %              must be a non-negative interger. Nitr=1E3 is the default
0030 %              setting.
0031 %
0032 % OUTPUT:
0033 %   - V     : N-by-3 array of vertex positions.
0034 %   - Tri   : M-by-3 list of face-vertex connectivities.
0035 %   - Ue_i  : N-by-1 array of particle energies.
0036 %   - Ue    : K-by-1 array of energy scores, where K-1 was the total number
0037 %             of iterations. Ue(1) corresponds to the energy of the initial
0038 %             configuration.
0039 %
0040 % EXAMPLE:
0041 % -------------------------------------------------------------------------
0042 % % Uniformly distribute 100 particles across the surface of the unit sphere
0043 %
0044 %  [V,Tri,~,Ue]=ParticleSampleSphere('N',100);
0045 %  figure, plot(0:numel(Ue)-1,Ue,'.-')
0046 %  set(get(gca,'Title'),'String','Optimization Progress','FontSize',20)
0047 %  xlabel('Iteration #','FontSize',15)
0048 %  ylabel('Reisz s-enrgy','FontSize',15)
0049 %  TR=TriRep(Tri,V);
0050 %  figure, h=trimesh(TR); set(h,'EdgeColor','b'), axis equal
0051 % -------------------------------------------------------------------------
0052 %
0053 % AUTHOR: Anton Semechko (a.semechko@gmail.com)
0054 % DATE: June.2012
0055 %
0056 
0057 % Check the inputs
0058 [V,s,Etol,Nitr]=VerifyInputArgs(varargin);
0059 N=size(V,1);    % number of particles
0060 clear varargin
0061     
0062 % Compute geodesic distances between the points
0063 DOT=V*V';       % dot product
0064 DOT(DOT<-1)=-1; DOT(DOT>1)=1;
0065 GD=acos(DOT);   % geodesic distance
0066 
0067 % Evaluate the energy functional
0068 GD(1:(N+1):end)=Inf; % set diagonal entries of GD to Inf
0069 Ue_ij=1./((GD.^s)+eps);
0070 Ue_i=sum(Ue_ij,2);
0071 Ue=sum(Ue_i);
0072 
0073 % Iteratively optimize the position of the particles along the negative
0074 % gradient of the energy functional using an adaptive Gauss-Seidel
0075 % update scheme -----------------------------------------------------------
0076 
0077 fprintf('\nWait while particle positions are being optimized ...\n')
0078 t0=clock;
0079 i=0; 
0080 
0081 a=ones(N,1);    % step sizes used during position updates
0082 a_min=1E-14;    % minimum step size
0083 a_max=0.1;      % maximum step size
0084 dE=Inf;
0085 while i<Nitr && dE>Etol
0086     
0087     i=i+1;
0088     
0089     % Sort the particles according to their energy contribution
0090     [~,idx_sort]=sort(Ue_i,'descend');
0091     
0092     % Update the position of individual particles
0093     for k=1:N
0094         
0095         j=idx_sort(k);
0096         
0097         idx_j=true(N,1); % particle indices, except the current one
0098         idx_j(j)=false;
0099         
0100         DOTj=DOT(idx_j,j); 
0101         GDj=GD(idx_j,j); 
0102         
0103         % Compute the gradient for the j-th particle
0104         dVj=bsxfun(@times,s./(sqrt(1-DOTj.^2)+eps),V(idx_j,:));                
0105         dVj=bsxfun(@rdivide,dVj,(GDj.^(s+1))+eps);
0106         if i<5 % remove contributions from particles which are too close
0107             idx_fail=sum(dVj.^2,2)>(1E4)^2;
0108             dVj(idx_fail,:)=[];
0109             if isempty(dVj) % perturb initial positions
0110                 V=V+randn(size(V))/1E5;
0111                 V_L2=sqrt(sum(V.^2,2));
0112                 V=bsxfun(@rdivide,V,V_L2);
0113                 DOT=V*V';
0114                 GD=acos(DOT);
0115                 break
0116             end 
0117         end
0118         dVj=sum(dVj,1);
0119         
0120         % Only retain the tangential component of the gradient
0121         dVj_n=(dVj*V(j,:)')*V(j,:);
0122         dVj_t=dVj-dVj_n;        
0123         
0124         % Adaptively update position of the j-th particle
0125         m=0; 
0126         Uj_old=sum(Ue_ij(j,:));
0127         while true
0128         
0129             % Update the position
0130             Vj_new=V(j,:)-a(j)*dVj_t;
0131             
0132             % Constrain the point to surface of the sphere
0133             Vj_new=Vj_new/norm(Vj_new);
0134                      
0135             % Recompute the dot products and the geodesics
0136             DOTj=sum(bsxfun(@times,V,Vj_new),2);
0137             DOTj(DOTj<-1)=-1; DOTj(DOTj>1)=1;
0138             GDj=acos(DOTj);
0139             GDj(j)=Inf;
0140             
0141             Ue_ij_j=1./((GDj.^s)+eps);
0142             
0143             % Check if the system potential decreased
0144             if sum(Ue_ij_j)<=Uj_old
0145                 
0146                 V(j,:)=Vj_new;
0147                 
0148                 DOT(j,:)=DOTj';
0149                 DOT(:,j)=DOTj;
0150                 
0151                 GD(j,:)=GDj';
0152                 GD(:,j)=GDj;
0153                 
0154                 Ue_ij(j,:)=Ue_ij_j';
0155                 Ue_ij(:,j)=Ue_ij_j;
0156                 
0157                 if m==1, a(j)=a(j)*2.5; end
0158                 if a(j)>a_max, a(j)=a_max; end
0159                 
0160                 break
0161                 
0162             else
0163                 
0164                 if a(j)>a_min
0165                     a(j)=a(j)/2.5;
0166                     if a(j)<a_min, a(j)=a_min; end
0167                 else
0168                     break
0169                 end
0170                 
0171             end
0172             
0173         end
0174         
0175     end
0176     
0177     % Evaluate the total energy of the system
0178     Ue_i=sum(Ue_ij,2);
0179     Ue(i+1)=sum(Ue_i);
0180 
0181     % Progress update
0182     if i==1
0183         fprintf('Iteration#\tEnergy Score\n')
0184         fprintf('%u\t        %.3f\t\n',0,Ue(1))
0185     else
0186         if mod(i,50)==0
0187         fprintf('%u\t        %.3f\t\n',i,Ue(end))
0188         end
0189     end
0190     
0191     % Change in energy
0192     if i>=10, dE=(Ue(i-2)-Ue(i+1))/10; end
0193  
0194     % Reset the step sizes
0195     if mod(i,20)==0, a=a_max*ones(N,1); end
0196     
0197 end
0198 clear DOT GD Ue_ij
0199 if Nitr==0, fprintf('%u\t        %.3f\t\n',i,Ue(1)); end
0200 if mod(i,50)~=0, fprintf('%u\t        %.3f\t\n',i,Ue(end)); end
0201 fprintf('Optimization completed after %u iterations. Elapsed time: %5.1f sec\n',i,etime(clock,t0))
0202 
0203 try
0204     Tri=fliplr(convhulln(V));
0205 catch err
0206     msg=sprintf('Unable to triangulate the points. %s',err.message);
0207     disp(msg)
0208     Tri=[];
0209 end
0210 
0211 
0212 %==========================================================================
0213 function [V,s,Etol,Nitr]=VerifyInputArgs(VarsIn)
0214 % Make sure all supplied input argumets have valid format
0215 
0216 % Default settings
0217 V=RandSampleSphere(200,'stratified');
0218 s=1; Etol=1E-5; Nitr=1E3;
0219 if isempty(VarsIn), return; end
0220 
0221 % First check that there is an even number of inputs
0222 Narg=numel(VarsIn);
0223 if mod(Narg,2)~=0
0224     error('Incorrect number of input arguments')
0225 end
0226 
0227 % Get the properties
0228 FNo={'N','Vo','s','Etol','Nitr'};
0229 flag=false(1,5); exit_flag=false;
0230 for i=1:Narg/2
0231     
0232     % Make sure the input is a string
0233     str=VarsIn{2*(i-1)+1};
0234     if ~ischar(str)
0235         error('Input argument #%u is not valid',2*(i-1)+1)
0236     end
0237     
0238     % Get the value
0239     Val=VarsIn{2*i};
0240     
0241     % Match the string against the list of avaliable options
0242     chk=strcmpi(str,FNo);
0243     id=find(chk,1);
0244     if isempty(id), id=0; end
0245     
0246     switch id
0247         case 1 % number of particles
0248             
0249             % Check if 'initialization' option has also been specified
0250             if flag(2) 
0251                 error('Ambigious combination of options. Specify option ''%s'' or option ''%s'', but not both.',FNo{2},FNo{1})
0252             end
0253             
0254             % Check the format
0255             if sum(Val(:)<10)>0 || sum(Val(:)>1E3)>0 || ~isreal(Val) || ~isnumeric(Val) || numel(Val)~=1 
0256                 error('Incorrect entry for the ''%s'' option. N must be a positive interger greater than 9 and less than 1001.',FNo{1})
0257             end
0258             V=RandSampleSphere(round(Val),'stratified'); %#ok<*NASGU>
0259             
0260         case 2 % initialization
0261         
0262             % Check if 'number' option has also been specified
0263             if flag(1) 
0264                 error('Ambigious combination of options. Specify option ''%s'' or option ''%s'', but not both.',FNo{1},FNo{2})
0265             end
0266             
0267             % Check the format
0268             if ~isreal(Val) || ~isnumeric(Val) || size(Val,2)~=3 || size(Val,1)<10 
0269                 error('Incorrect entry for the ''%s'' option. Vo must be a N-by-3 array, where N is the number of particles. N=10 is the lowest number allowed.',FNo{2})
0270             end
0271                         
0272             % Make sure the particles are constrained to the surface of the unit sphere
0273             V=Val;
0274             V_L2=sqrt(sum(V.^2,2));
0275             V=bsxfun(@rdivide,V,V_L2);
0276             clear Val V_L2 
0277             
0278             % Check if there are more than 1E3 particles
0279             if size(V,1)>1E3
0280                 
0281                 % Construct a 'yes'/'no' questdlg
0282                 choice = questdlg('Default particle limit exceeded. Would you like to continue?', ...
0283                     'Particle Limit Exceeded','   YES   ','   NO   ','   NO   ');
0284                 
0285                 % Handle response
0286                 if strcmpi(choice,'   NO   '), exit_flag=true; end
0287                 
0288             end
0289             
0290         case 3 % s parameter
0291             
0292             % Check the format
0293             if Val<0 || ~isreal(Val) || ~isnumeric(Val) || numel(Val)~=1 
0294                 error('Incorrect entry for the ''%s'' option. s must be a positive real number.',FNo{3})
0295             end
0296             s=Val;
0297             
0298         case 4 % energy tolerance parameter
0299         
0300             % Check the format
0301             if Val<0 || ~isreal(Val) || ~isnumeric(Val) || numel(Val)~=1 
0302                 error('Incorrect entry for the ''%s'' option. Etol must be a positive real number.',FNo{4})
0303             end
0304             Etol=Val;
0305             
0306         case 5 % maximum number of iterations
0307             
0308             % Check the format
0309             if Val<0 || ~isreal(Val) || ~isnumeric(Val) || numel(Val)~=1 
0310                 error('Incorrect entry for the ''%s'' option. Nitr must be a non-negative integer.',FNo{5})
0311             end
0312             Nitr=Val;
0313             
0314         otherwise
0315             error('''%s'' is not a valid option',str)
0316     end
0317     flag(id)=true;
0318     
0319 end
0320 
0321 if exit_flag, Nitr=0; end %#ok<*UNRCH>
0322

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005