Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
300 lines (260 sloc) 11.5 KB
function varargout=maternprc(theta,prc,meth,xver)
% [Kprc,krange,Sk,kSkgues]=MATERNPRC(theta,prc,meth,xver)
%
% Given a two-dimensional Matern parameter vector theta [s2 nu rho]
% returns the wavenumber K for which $int_\0^K S(k;\theta) k dk$
% equals a target percentage of $int_\0^\inf S(k;\theta) k dk$
%
% INPUT:
%
% theta vector of Matern parameters [s2 nu rho]
% prc target percentile of overall area under Sk [default: 50]
% meth 0 symbolically build up to numerical/analytic methods (slow!)
% 1 numerical (using an initial guess)
% 2 analytical [default, fast, and preferred]
% xver 0 no extra verification [default]
% 1 extra verification via the symbolic method, and after the
% fact for the other methods
%
% OUTPUT:
%
% Kprc wavenumber corresponding to target percentage
% krange some trial wavenumber range, used to create kSkgues
% Sk the actual Matern function
% kSkgues the half-Sk wavenumber, used as a guess for meth==1
%
% NOTE: The FZERO function used requires an initial guess, picked from a
% logarithmically spaced set of wavenumbers that should be good for
% most Matern Sk curves, but this is not guaranteed.
%
% EXAMPLE:
%
% maternprc('demo1') % Just a bunch of default evaluations
% maternprc('demo2') % A bunch of random evaluations with some plots
% maternprc('demo3') % A bunch of sequential evaluations
% maternprc('demo4') % Gabe no-like values
%
% TO DO: shunt the 0th and 100th
%
% SEE ALSO: MATERNOS, MATERNOSP, MATERNOSY, FZERO, SYM, SYMS, TRAPEZE
%
% Last modified by gleggers-at-princeton.edu, 03/19/2014
% Last modified by fjsimons-at-princeton.edu, 06/23/2018
% Default values
defval('prc',50)
defval('meth',2)
defval('xver',0)
defval('theta',[0.022 2.5 35000])
% Still need one for the recursive part where xver is set to 0
defval('krange',NaN);
defval('Sk',NaN);
defval('kSkgues',[NaN NaN]);
% Not a demo!
if ~isstr(theta)
% Regardless of the method, shunt the 0 and 100 percentiles but note
% that in those cases you only ever get one output
if prc==0
Kprc=0;
varns={Kprc,krange,Sk,kSkgues};
varargout=varns(1:nargout);
return
elseif prc==100
Kprc=theta(1);
varns={Kprc,krange,Sk,kSkgues};
varargout=varns(1:nargout);
return
end
% If you'll be needing a krange output or input, make it here
if nargout>1 || meth==1 || xver==1
% Start with a trial wavenumber range
krange=logspace(-9,-1,1e4);
end
% And this is the theoretical y-axis intercept
S0=[theta(1)*theta(3)^2/4*pi];
% Symbolic math method, and extra verification
if meth==0 || xver==1
% Matern isotropic 2D spectral density, Simons & Olhede (2013), eq. (72)
syms s2 nu rho positive
syms k K pee prcnt positive
Sk=sym(...
's2*nu^(nu+1)*4^nu/pee/(pee*rho)^(2*nu)*(4*nu/(pee*rho)^2+k^2)^(-nu-1)');
% The area under the curve is the TWO-dimensional integral that
% in polar coordinates is a ONE-dimensional integral
% A = \int_0^K \int_0^K S(k) d\bk
% = 2\pi \int_0^K S(k) k dk
% The integrand of the equivalent ONE-dimensional integral
Skk=2*pee*Sk*k;
% The total area under the curve from 0 to inf simply IS s2...
Sintinf=simplify(int(Skk,'k',0,Inf));
% ... which we should explicitly check by subbing anything
difer(subs(Sintinf-s2),[],[],NaN)
% The area under the curve from 0 to K ranges between 0 and s2
SintK=simplify(int(Skk,'k',0,K));
% This is the evaluated Sintk spelled out...
spl1=2*4^nu*nu^(nu + 1)*s2/(pee*rho)^(2*nu)*...
[((pee^2*rho^2)/(4*nu))^nu/(2*nu)-1/(2*nu*(K^2+(4*nu)/(pee^2*rho^2))^nu)];
% ... which we should explicitly check by subbing anything (even
% nothing!) into what should be exactly, symbolically zero
difer(subs(SintK-spl1,symvar(spl1)),[],[],NaN)
% And this is how we simplify it even further, by hand
spl2=s2*[1-2^(2*nu)*nu^nu/(K^2*pee^2*rho^2 + 4*nu)^nu];
% Divide by s2 since we'll have a right hand side containing s2
spl3= [1-2^(2*nu)*nu^nu/(K^2*pee^2*rho^2 + 4*nu)^nu];
% you can verify that these two representations are identical by
% substitution, e.g.
% ... this should be exactly, symbolically zero
difer(subs(symvar(spl2)-symvar(spl1)),[],[],NaN)
% ... this will be exactly, numerically zero (I need to evaluate)
% since Matlab isn't smart enough to realize they are actually,
% symbolically zero, except perhaps after simplification
difer(subs(spl1-spl2,[K s2 nu rho pee],[rand theta pi]),[],[],NaN)
% ... turns out that simplify brings out the symbolic zero
difer(subs(simplify(spl1-spl2)),[],[],NaN)
% or you can turn them into symbolic functions, e.g.
spl1f=symfun(spl1,[K s2 nu rho pee]);
spl2f=symfun(spl2,[K s2 nu rho pee]);
% Watch the order of the input variables - symvar lists
% alphabetically
rando=rand;
difer(eval(spl1f(rando,theta(1),theta(2),theta(3),pi))-...
eval(spl2f(rando,theta(1),theta(2),theta(3),pi)),...
[],[],NaN)
% Solve for a certain percentile of s2; divide by s2 then solve this
K2solv=symfun(spl3-prcnt,[K prcnt nu rho pee]);
% We can write the solution out by hand - RB X page 25
Ksolved=symfun(2/(pee*rho)*sqrt(nu*[(1-prcnt)^(-1/nu)-1]),[prcnt nu rho pee]);
end
% Now cycle through the methods
switch meth
case 0
% Symbolic derivation of the numerical and the analytical method
% If you made it all the way here, you can just evaluate the
% analytical solution...
Kprc=eval(Ksolved(prc/100,theta(2),theta(3),pi));
% Verfiy that you have just found the solution by plugging in
difer(eval(K2solv(Kprc,prc/100,theta(2),theta(3),pi)),[],[],NaN)
% ... or you can also solve the symbolic expression numerically
% after converting it to an anonymous function handle
K2solvemf=matlabFunction(K2solv);
% Then solve numerically using FZERO to be sure
% in this case we simply take Kprc as the initial guess...
Kprc2=fzero(@(K) K2solvemf(K,theta(2),pi,prc/100,theta(3)),Kprc);
% Check that they're both consistent
difer(Kprc-Kprc2,[],[],NaN)
case 1
% Direct route to the numerical method...
K2solvemf=@(K,nu,pee,prcnt,rho) ...
-prcnt-2.0.^(nu.*2.0).*nu.^nu.*(nu.*4.0+K.^2.*pee.^2.*rho.^2).^(-nu)+1.0;
% ... but now need proper guess
% Here is the actual Matern function
Sk=maternos(krange,theta);
% Check that the intercept is correct to 1e-7 or thereabouts
difer(maternos(0,theta)-S0,7,[],NaN);
% Find the halfway point of the curve to use as a starting guess K
% kSkgues=krange(find(diff(sign(Sk/S0-0.5))));
% Find the largest wavenumber where the normalized kSk curve
% reaches a low value (to being the interesting nook) 0.01
kSkgues=kguess(krange,Sk,prc);
% What is the solution in the end? One could force it all to be
% positive - in earlier experiments we found the right number, only
% the wrong sign (we doing other things differently...)
Kprc=abs(fzero(@(K) K2solvemf(K,theta(2),pi,prc/100,theta(3)),kSkgues(1)));
case 2
% Direct route to the analytical method
Kprc=2/(pi*theta(3))*sqrt(theta(2)*[(1-prc/100)^(-1/theta(2))-1]);
end
% More output requested, but you hadn't needed them: make them
if nargout>2 && [exist('Sk','var')==0 || strcmp(class(Sk),'sym') || any(isnan(Sk))]
Sk=maternos(krange,theta);
end
if nargout>3 && [exist('kSkgues','var')==0 || any(isnan(kSkgues))]
% kSkgues=krange(find(diff(sign(Sk/S0-0.5))));
kSkgues=kguess(krange,Sk,prc);
end
% If extra verification is requested, compare the calculated results
if xver == 1
% Check "Sintinf" via explicit numerical integration on the krange
cSintinf=trapeze(krange(:),2*pi*maternos(krange(:),theta).*krange(:));
disp(sprintf('\nS2 %g TRAPEZE %g Difference %g',...
theta(1),cSintinf,theta(1)-cSintinf))
if meth~=0
% Check "Kprc"via the "other" method (of course meth==0 is OK)
Kprc2=maternprc(theta,prc,meth+(-1)^(meth+1),0);
disp(sprintf('METHOD %g OTHER METHOD %g Difference %g',...
Kprc,Kprc2,Kprc-Kprc2))
end
% Check "Kprc": Compare the calculated area under the curve from 0 to K
% to the area that K corresponds to as a percentage of the total area
ktoK=krange(krange(:)<=Kprc); ktoK=ktoK(:);
cSintK=trapeze(ktoK,2*pi*maternos(ktoK,theta).*ktoK);
disp(sprintf('%s0K %g TRAPEZE %g Difference %g\n',...
'%',prc/100*theta(1),cSintK,prc/100*theta(1)-cSintK))
end
% Provide as much output as requested
varns={Kprc,krange,Sk,kSkgues};
varargout=varns(1:nargout);
elseif strcmp(theta,'demo1')
prc=rand*100;
[Kprc,krange,Sk,kguess]=maternprc([],prc,0,1); disp(sprintf('%g',Kprc))
[Kprc,krange,Sk,kguess]=maternprc([],prc,1,1); disp(sprintf('%g',Kprc))
[Kprc,krange,Sk,kguess]=maternprc([],prc,2,1); disp(sprintf('%g',Kprc))
[Kprc,krange,Sk,kguess]=maternprc([],prc,0,0); disp(sprintf('%g',Kprc))
[Kprc,krange,Sk,kguess]=maternprc([],prc,1,0); disp(sprintf('%g',Kprc))
[Kprc,krange,Sk,kguess]=maternprc([],prc,2,0); disp(sprintf('%g',Kprc))
elseif strcmp(theta,'demo2')
prc=rand*100; meth=randi(3)-1; xver=round(rand);
thet=[max(round(rand*1000)/1000,0.0001) max(round(rand*10),0.5) max(round(rand*10000),0.00001)];
disp(sprintf('s2 %g nu %g rho %g',thet))
disp(sprintf('Trying the %gth percentile using method %i',prc,meth))
[a,b,c,d]=maternprc(thet,prc,meth,xver);
plotsk(a,b,c,d,prc,meth,xver)
elseif strcmp(theta,'demo3')
thet=[max(round(rand*1000)/1000,0.0001) max(round(rand*10),0.5) max(round(rand*10000),0.00001)];
disp(sprintf('s2 %g nu %g rho %g',thet))
meth=randi(3)-1; xver=round(rand);
for prc=1:3:99
[a,b,c,d]=maternprc(thet,prc,meth,xver);
plotsk(a,b,c,d,prc,meth,xver)
pause
end
elseif strcmp(theta,'demo4')
prc=rand*100; meth=randi(3)-1; xver=round(rand);
prc=50;
disp(sprintf('2-D Matern %ith percentile method %i xver %i',...
round(prc),meth,xver))
[Kprc,b,c,d]=maternprc([0.1636 0.0049 8.3618e+03],prc,meth,xver);
disp(sprintf('%g',Kprc)); plotsk(Kprc,b,c,d,prc,meth,xver); pause
[Kprc,b,c,d]=maternprc([0.0160 0.1146 5.9485e+04],prc,meth,xver);
disp(sprintf('%g',Kprc)); plotsk(Kprc,b,c,d,prc,meth,xver); pause
[Kprc,b,c,d]=maternprc([0.0046 0.2629 3.6451e+04],prc,meth,xver);
disp(sprintf('%g',Kprc)); plotsk(Kprc,b,c,d,prc,meth,xver); pause
[Kprc,b,c,d]=maternprc([7.5641e-04 0.3199 1.5281e+05],prc,meth,xver);
disp(sprintf('%g',Kprc)); plotsk(Kprc,b,c,d,prc,meth,xver); pause
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plotsk(a,b,c,d,prc,meth,xver)
clf
Skp=semilogx(b,c/max(c),'LineW',2); hold on
Skkp=semilogx(b,c(:).*b(:)/max(c(:).*b(:)),'LineW',2); hold on
set(Skkp,'Color','r')
plot([a a],ylim,'r-','LineW',1); grid on
if meth~=0 && [meth~=2 && xver~=1]
plot(d(1),d(2),'o','MarkerF','r','MarkerE','r'); hold off
end
title(sprintf('2-D Matern %5.2fth percentile method %i xver %i',prc,meth,xver))
ylim([-0.1 1.1])
xlabel('wavenumber k')
ylabel('S(k)/S(0) and normalized S(k)k')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function kg=kguess(kr,Sk,prc)
% There is art in picking this level - if FZERO complains, change this
thepoint=0.01;
% Here's a trick
cands=kr(find(diff(...
sign(kr(:).*Sk/max(kr(:).*Sk)-thepoint))));
if prc<50
kg(1)=min(cands);
else
kg(1)=max(cands);
end
kg(2)=thepoint;