0001 function [v,x,t,m,ze]=v_quadpeak(z)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035 persistent wz a
0036
0037 sz=size(z);
0038 psz=prod(sz);
0039 dz=numel(sz);
0040 mz=find(sz>1);
0041 nm=numel(mz);
0042 vz=sz(mz);
0043 dx=max(mz);
0044 if ~nm
0045 error('Cannot find peak of a scalar');
0046 end
0047 nc=(nm+1)*(nm+2)/2;
0048 if min(vz)<3
0049 error('Need at least 3 points in each non-singleton dimension');
0050 end
0051 if isempty(wz) || numel(wz)~=numel(vz) || ~all(wz==vz)
0052 wz=vz;
0053 a=ones(psz,nc);
0054 ix=(0:psz-1)';
0055 for i=1:nm
0056 jx=floor(ix/sz(mz(i)));
0057 a(:,i+nc-nm-1)=1+ix-jx*sz(mz(i));
0058 ix=jx;
0059 a(:,(i^2-i+2)/2:i*(i+1)/2)=a(:,nc-nm:i+nc-nm-1).*repmat(a(:,i+nc-nm-1),1,i);
0060 end
0061 a=(a'*a)\a';
0062 end
0063
0064
0065
0066 c=a*z(:);
0067 w=zeros(nm+1,nm+1);
0068 i=1:(nm+1)*(nm+2)/2;
0069 j=floor((sqrt(8*i-7)-1)/2);
0070 w(i+j.*(2*nm+1-j)/2)=c;
0071 w=(w+w.')/2;
0072 mr=w(1:nm,1:nm);
0073 we=w(1:nm,nm+1);
0074 y=-(mr\we);
0075 v=y'*we+w(nm+1,nm+1);
0076
0077
0078
0079 x=zeros(dx,1);
0080 x(mz)=y;
0081 m=zeros(dx+1,dx+1);
0082 mz(nm+1)=dx+1;
0083 m(mz,mz)=w;
0084 if nargout>2
0085 ev=eig(mr);
0086 t=all(ev>0)-all(ev<0);
0087 end
0088 if nargout>4
0089 ze=zeros(sz);
0090 scp=cumprod([1 sz(1:end-1)]);
0091 ivec=fix(repmat((0:psz-1)',1,dz)./repmat(scp,psz,1));
0092 xe=[1+ivec-repmat(sz,psz,1).*fix(ivec./repmat(sz,psz,1)) ones(psz,1)];
0093 ze=reshape(sum((xe*m).*xe,2),sz);
0094 end
0095
0096
0097 if ~nargout && nm<=2
0098
0099 desc={'Maximum','Saddle Point','Minimum'};
0100 if nargout<=2
0101 ev=eig(mr);
0102 t=all(ev>0)-all(ev<0);
0103 end
0104 if nm==1
0105 xax=linspace(1,psz,100);
0106 plot(xax,c(1)*xax.^2+c(2)*xax+c(3),'-r',1:psz,z(:),'ob',x,v,'^k');
0107 set(gca,'xlim',[0.9 psz+0.1]);
0108 ylabel('z');
0109 xlabel(sprintf('x%d',mz(1)));
0110 title(sprintf('\\Delta = %s: z(%.2g) = %.2g',desc{t+2},y(1),v));
0111 else
0112 ngr=17;
0113 xax=repmat(linspace(1,vz(1),ngr)',1,ngr);
0114 yax=repmat(linspace(1,vz(2),ngr),ngr,1);
0115 zq=(c(1)*xax+c(2)*yax+c(4)).*xax+(c(3)*yax+c(5)).*yax+c(6);
0116 hold off
0117 mesh(xax,yax,zq,'EdgeColor','r');
0118 hold on
0119 plot3(repmat((1:vz(1))',1,vz(2)),repmat(1:vz(2),vz(1),1),reshape(z,vz),'ob',y(1),y(2),v,'^k');
0120 hold off
0121 set(gca,'xlim',[0.9 vz(1)+0.1],'ylim',[0.9 vz(2)+0.1]);
0122 xlabel(sprintf('x%d',mz(1)));
0123 ylabel(sprintf('x%d',mz(2)));
0124 zlabel('z');
0125 title(sprintf('\\Delta = %s: z(%.2g,%.2g) = %.2g',desc{t+2},y(1),y(2),v));
0126 end
0127 end
0128
0129