-
Notifications
You must be signed in to change notification settings - Fork 3
/
maternosy.m
328 lines (317 loc) · 14 KB
/
maternosy.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
function varargout=maternosy(y,th,dth,meth)
% CyordCydth=MATERNOSY(y,th,dth,meth)
%
% Calculates d-dimensional isotropic Matern correlation, which is
% independent of d, also the counterpart to the spectral covariance.
% See Olhede & Simons (2013), doi: 10.1093/gji/ggt056.x, eq. (72)
% Additionally can calculate the partial derivatives of the isotropic Matern
% correlation with respect to the Matern parameters.
%
% INPUT:
%
% y Lag parameter, the distance between spatial positions,
% e.g. from XXPDIST or SSPDIST [m], can be any dimension
% th The unscaled parameter vector, with, in the last three slots:
% s2 The first Matern parameter [variance in units^2]
% nu The second Matern parameter [differentiability]
% rh The third Matern parameter [range in units]
% dth When empty, the spatial Matern covariance is calculated; otherwise,
% 1, 2, or 3 specifies which element of th gets differentiated
% meth 1 straight up calculation
% 2 analytical simplification for special values of th(end-1)
%
% OUTPUT:
%
% CyordCydth The spatial Matern covariance at all the requested lags, possibly
% calculated using the Abramowitz & Stegun (1965) simplifications to
% the Bessel function for half-integer values of nu, the infinite limit
% of nu, OR returns the DERIVATIVE in the dth element of th, to feed
% into BLUROSY
%
% SEE ALSO:
%
% MATERNOS, MATERNOSP, MATERNPRC
%
% EXAMPLE:
%
% Discrete approximation of the integral compared to exact result:
%
% [Hx,th0,p]=simulosl; bigN=1000;
% y=linspace(0,sqrt(prod(p.dydx))*sqrt(prod(p.NyNx)),bigN);
% % S(0)=1/(2*pi)\int C(r)rdr
% [sum(y.*maternosy(y,th0))*(y(2)-y(1))/(2*pi) maternos(0,th0)]
%
% Confirmation that the special case simplifications are equivalent to the
% straight calculation:
%
% nu=[1/3 1/2 1 3/2 5/2 31/2]; th0(2)=nu(randi(length(nu)));
% Cy=maternosy(y,th0,[],1); Cy2=maternosy(y,th0,[],2); diferm(Cy/(Cy(1)),Cy2/Cy2(1))
%
% Calculate the infinite smoothness, Gaussian, squared exponential case:
%
% th0(2)=Inf; Cy=maternosy(y,th0);
%
% Calculate and visualize the partial derivatives:
%
% maternosy('demo2')
%
% Demo1 provides a visual comparison of the general output of MATERNOS and
% MATERNOSY as Fourier pairs:
%
% maternosy('demo1')
%
% Last modified by fjsimons-at-alum.mit.edu, 06/17/2024
% Last modified by olwalbert-at-princeton.edu, 06/17/2024
if ~isstr(y)
% Defaults
defval('meth',1)
defval('dth',[])
% The Matern parameters are always the last three elements of the TH input
s2=th(end-2);
nu=th(end-1);
rh=th(end );
% The argument, make sure it is a distance
argu=2*sqrt(nu)/pi/rh*abs(y);
% Calculate the spatial covariance?
if isempty(dth)
% Check whether we are asking for the case that nu approaches
% infinity; if so, by the regular method Cy would have Inf/NaN entries
if isinf(nu)
% Squared exponential
% This analytic form of the isotropic Matern covariance
% was solved from the inverse Fourier transform of the limit as
% nu approaches infinity of the spectral density (see MATERNOS
% for details) with use of Eq. 3.323.2 of Gradshteyn & Ryzhik (1980).
%%%OLW: (Consistency with BLUROSY demo?)
Cy=s2/(pi*rh)*exp(-abs(y).^2/(pi^2*rh^2));
else
% Switch the calculation method
switch meth
case 1
% The evaluation
Cy=2^(1-nu)*s2/gamma(nu)*argu.^nu.*besselk(nu,argu);
% Supply the smallest arguments
Cy(y==0)=s2;
case 2
% By selecting the second method of calculation, we are seeking to
% evaluate Cy from the simplified analytic expression of the
% isotropic Matern covariance for a special value of nu. The
% following analytic forms of half-integer nu are calculated from
% substitution of Eqs 10.1.9 and 10.2.15 of Abramowitz & Stegun
% (1965) for modified Bessel functions of half-integer orders.
% See also 10.1093/biomet/93.4.989 for comparison.
if nu==1/3
% von Karman
Cy=s2*2^(2/3)/gamma(nu).*(2*sqrt(3)/(3*pi*rh)*abs(y)).^...
(nu).*besselk(nu,2*sqrt(3)/(3*pi*rh)*abs(y));
% Compute the value at zero lag
Cy(y==0)=s2;
elseif nu==1/2
% Exponential
Cy=s2*exp(-sqrt(2)/(pi*rh)*abs(y));
elseif nu==1
% Whittle
Cy=s2*2/(pi*rh)*abs(y).*besselk(nu,2/(pi*rh)*abs(y));
% Compute the value at zero lag
Cy(y==0)=s2;
elseif nu==3/2
% Second-order autoregressive
Cy=s2*exp(-sqrt(6)/(pi*rh)*abs(y)).*(1+sqrt(6)/(pi*rh)*abs(y));
elseif nu==5/2
% Third-order autoregressive
Cy=s2*exp(-sqrt(10)/(pi*rh)*abs(y)).*(1+sqrt(10)/(pi*rh)*...
abs(y)+10/(3*pi^2*rh^2)*abs(y).^2);
elseif mod(nu+0.5,1)==0
% A general half-integer case beyond 1/2, 3/2, or 5/2,
% nu=n+1/2 for n = 1, 2, 3, ...
k=0:nu-0.5;k=k(:);
n=nu-0.5;
Cy=s2*exp(-2/(pi*rh)*abs(y)*sqrt(nu))*factorial(n)/factorial(2*n).*...
sum(factorial(n+k)./(factorial(k).*factorial(n-k)).*...
(4*sqrt(nu)*abs(y)/(pi*rh)).^(n-k),1);
else
% However, if the nu provided to MATERNOSY is not one of the
% special values of nu that we have considered so far, we should
% throw an error
error('This is not a special case of nu. Ask for meth=1.')
end
end
end
% Or, calculate the derivatives?
else
% Calculate the partial derivative wrt the requested parameter index
% disp(sprintf('Calculating %i%s parameter derivative',dth,ith(dth)))
% Abuse of Cy nomenclature, now you're getting a specific derivative
if dth==1
% The partial derivative of Cy with respect to the variance, dCyds2
Cy=2^(1-nu)/gamma(nu)*argu.^nu.*besselk(nu,argu);
% Supply the smallest arguments
Cy(y==0)=1;
elseif dth==2
% The partial derivative of Cy with respect to the smoothness,
% dCydnu; the derivative of the gamma function is provided by Eq.
% 6.3.1 of Abramowitz & Stegun (1965);
if nu==0
% the partial derivative of the modified Bessel function of the
% second kind with respect to order when the order is evaluated
% at 0 is given by Eq. 9.6.46 of A&S, which we do not consider for
% Matern
dKdnuo=0;
elseif mod(nu,1)==0
% for integer order, the partial derivative of the modified
% Bessel function of the second kind with respect to order is
% provided by Eq. 9.6.45 of A&S
syms k;
dKdnuo=(factorial(nu).*(argu/2).^(-nu)/2).*...
vpasum((argu/2).^k.*besselk(k,argu)./...
((nu-k).*factorial(k)),0,nu-1);
elseif nu==1/2
% for order of 1/2, we use Eq. 10.2.34 of A&S; note that
% dCydnu(y=0,nu=1/2)=NaN because Ei(0)=-Inf; there also appear to
% be numerical instability issues -- many NaNs and Infs for small
% parameter values, e.g., th0=[1 0.5 1]; Olivia will study this
% further. Default s2,rh for th0=[1e6 0.5 2e4] seems to behave.
warning('derivative wrt nu is still in development for nu given')
dKdnuo=-sqrt(pi/2*argu).*ei(-2*argu).*exp(argu);
%elseif mod(nu-1/2,1)==0
% for half-integer order, Brychkov & Geddes (2005) presented a
% solution for n-1/2 and n+1/2 order. We will use the n-1/2
% order (Eq. 26)
%warning('derivative wrt nu is still in development for nu given')
% I need to look into this expression further before we include
% it.
%n=nu+0.5;
%argun=2*sqrt(n)/pi/rh*abs(y);
%syms k p;
%dKdnuo=(-1)^n*(pi/2)*(besseli(n-0.5,argun)+...
% besseli(0.5-n,argun)).*(coshint(2*argun)-sinhint(2*argun))+0.5*...
% vpasum((nchoosek(n,k)).*factorial(n-k-1).*(argun/2).^(k-n).*...
% besselk(k-0.5,argun),0,n-1)-...
% (-1)^n*sqrt(pi)*sum((-1).^k.*(nchoosek(n,k)).*2.^(k-1).*...
% (besseli(n-k-0.5,argun)+besseli(k-n+0.5,argun)).*...
% vpasum((nchoosek(k-1,p)).*factorial(k-p-1).*argun.^(p-k+0.5).*...
% besselk(p-0.5,2*argun),0,k-1),...
% 1,n);
else
% for non-integer order, substitute DLMF Eq. 10.38.1 into
% Eq. 9.6.43 of A&S; we need to make an approximation to the
% infinite sum terms by choosing a large sumlim
warning('derivative wrt nu is still in development for nu given')
sumlim=1;
syms k;
dIdnup=besseli(nu,argu).*log(argu/2)-(argu/2).^nu.*...
vpasum((psi(nu+k+1)./gamma(nu+k+1).*...
(argu/2).^(k*2)./factorial(k)),0,sumlim);
dIdnun=-besseli(-nu,argu).*log(argu/2)+(argu/2).^(-nu).*...
vpasum((psi(-nu+k+1)./gamma(-nu+k+1).*...
(argu/2).^(k*2)./factorial(k)),0,sumlim);
dKdnuo=pi/2*csc(nu*pi)*(dIdnun-dIdnup)-...
pi*cot(nu*pi)*besselk(nu,argu);
keyboard % solution blowing up for large y when th0=[1e6 1.5 1e6],
% NaNs and Inf when [1 1.5 1]?
figure();plot(y(1:end-50),dKdnuo(1:end-50))
end
% Now put things together
Cy=s2*2^(1-nu)/gamma(nu)*argu.^nu.*...
(besselk(nu,argu).*(0.5+log(argu/2)-psi(nu))-...
argu.^nu/(2*nu).*(argu.*besselk(nu-1,argu)+...
nu*besselk(nu,argu)).*double(dKdnuo));
% Supply the smallest arguments
Cy(y==0)=0;
% From the asymptotic expansion of Kn for small z, differentiated in Mathematica
%Cy=2^(1-nu)*s2/gamma(nu)*argu.^nu.*gamma(nu)*(argu/2)^(-nu)/2;
%bla=-1/4*pi^nu*(sqrt(nu/rh*y)).^(-nu)*gamma(nu).*(1+log(nu)-2*log(pi*rh)+2*log(y)-2*psi(nu));
elseif dth==3
% The partial derivative of Cy with respect to the range, dCydrho;
% simplification of the derivative of the Bessel term with respect
% to argument is made through Eq. 3.71.3 of Watson (1962)
Cy=(s2/rh)*2^(1-nu)/gamma(nu)*argu.^(nu+1).*besselk(nu-1,argu);
% Supply the smallest arguments
Cy(y==0)=0;
else
error('Not a valid partial derivative index. Ask for dth=1,2,or 3.')
end
end
% Serve output
varns={Cy};
varargout=varns(1:nargout);
elseif strcmp(y,'demo1')
% The Fourier relation between the correlation and the covariance
% is hard to verify... since we can never observe all lags on a
% fine and large enough grid... and that is the point of estimating
% parameters differently. See the comparison in SIMULOSL1('demo1')
% Matern parameters... very sensitive to the last one
th=[2000 0.5 2/3];
% Grid size, also in physical units, keep it even for this example
p.NyNx=[4300 5500]+randi(1000,[1 2])*(-1)^round(rand);
p.NyNx=p.NyNx+mod(p.NyNx,2);
lY=13; lX=24;
p.dydx=[lY/(p.NyNx(1)-1) lX/(p.NyNx(2)-1)];
% Wavenumber grid
[k,kx,ky,dci]=knum2(p.NyNx,[lY lX]);
% Space grid
x=[-floor(p.NyNx(2)/2):1:+floor(p.NyNx(2)/2)-1]*p.dydx(2);
y=[-floor(p.NyNx(1)/2):1:+floor(p.NyNx(1)/2)-1]*p.dydx(1);
[X,Y]=meshgrid(x,y); yy=sqrt(X.^2+Y.^2);
% Evaluate the Matern spectral covariance
Sbb=v2s(maternos(k,th,[],2),p);
% Evaluate the Matern correlation
Cy=maternosy(yy,th,[],2); difer(Cy(dci(1),dci(2))-th(1))
% Fourier transform the correlation to check its relation to the covariance
Skk=fftshift(v2s(tospace(Cy,p),p));
% Compare two profiles somehow
m=median(Sbb(dci(1),:)./abs(Skk(dci(1),:)));
m=max(max(Sbb))./max(max((Skk)));
plot(log10(abs(Skk(dci(1),:))),'Color','r');
hold on
plot(log10(Sbb(dci(1),:))-log10(m),'LineWidth',2,'Color','m');
hold off
legend('FFT(MATERNOSY)','MATERNOS')
axis tight
grid on
% Annotate
title(sprintf('p.NyNx = [%i %i]',p.NyNx))
xlabel('wavenumber')
ylabel('MATERNOS 2s FFT(MATERNOSY)')
elseif strcmp(y,'demo2')
clf
th0 = [1 2 27]; p = []; p.NyNx = [256 256]; p.dydx = [1 1];
y = linspace(0,sqrt(prod(p.dydx)*prod(p.NyNx)),1000);
xels=[0 min(6*pi*th0(3),max(y))];
xels=xels+[-1 1]*range(xels)/20;
% Calculate Matern Covariance
labs={'\sigma^2','\nu','\rho'};
% Calculate Matern Covariance derivatives
for index=1:3
ah1(index)=subplot(2,3,index);
pc(index)=plot(y,maternosy(y,th0,1),'k','LineWidth',1); hold on
% Let's add a few neighbors
% Percentage perturbation in each of the three parameters
pc1(index)=plot(y,maternosy(y,th0.*(1+rindeks(eye(3),index)/20)),'r');
pc2(index)=plot(y,maternosy(y,th0.*(1-rindeks(eye(3),index)/20)),'b'); hold off
ylim([0 th0(1)]+[-1 1]*th0(1)/20)
% ylim([-0.1 0.1])
xlim(xels)
if index==2
title(sprintf('%s = %g | %s = %g | %s = %g',...
'\sigma^2',th0(1),'\nu',th0(2),'\rho',th0(3)))
end
grid on
ylabel('C(y)')
xlabel('y')
ah2(index)=subplot(2,3,index+3);
zdiff=maternosy(y,th0,index);
pd(index)=plot(y,zdiff,'k');
ylim(halverange(zdiff,105,NaN))
% if index==2
% hold on
% plot([1 1]*pi*th0(3)/2/sqrt(th0(2)),ylim)
% hold off
% end
grid on
ylabel(sprintf('dC(y)/d%s',labs{index}))
xlabel('y')
xlim(xels)
end
longticks([ah1 ah2],2)
end